text:
function writeParaviewMooring(moorDyn, model,t,simdate,nline,nnode,paraviewPath,TimeBodyParav,NewTimeParaview)
% Method to write ``vtp`` Paraview visualization files for the
% mooringClass. Executed by paraviewVisualization.m when simu.paraview.option=1
% in the wecSimInputFile.m
%
% Parameters
% ------------
%   moorDyn : obj
%       The moorDyn object
%   model : string
%       The simMechanics ``.slx`` filename
%   t : float vector
%       Moordyn time vector
%   simdate : string
%       Date and time of the simulation
%   nline : integer
%       Number of mooring lines
%   nnode : integer
%       Number of nodes on mooring lines
%   paraviewPath : directory
%       Directory the Paraview files were saved
%   TimeBodyParav : float vector
%       Paraview time vector
%   NewTimeParaview : float vector
%       Total simulation time interpolated by the Paraview time step
%

nsegment = nnode -1;
for iline = 1:nline
    for inode = 0:nnode(iline)-1
        moorDyn.(['Line' num2str(iline)]).(['Node' num2str(inode) 'px']) = interp1(t,moorDyn.(['Line' num2str(iline)]).(['Node' num2str(inode) 'px'])(:),NewTimeParaview);
        moorDyn.(['Line' num2str(iline)]).(['Node' num2str(inode) 'py']) = interp1(t,moorDyn.(['Line' num2str(iline)]).(['Node' num2str(inode) 'py'])(:),NewTimeParaview);
        moorDyn.(['Line' num2str(iline)]).(['Node' num2str(inode) 'pz']) = interp1(t,moorDyn.(['Line' num2str(iline)]).(['Node' num2str(inode) 'pz'])(:),NewTimeParaview);
    end
end
for iline = 1:nline
    for isegment = 0:nsegment(iline)-1
        moorDyn.(['Line' num2str(iline)]).(['Seg' num2str(isegment) 'Te']) = interp1(t,moorDyn.(['Line' num2str(iline)]).(['Seg' num2str(isegment) 'Te'])(:),NewTimeParaview);
    end
end
for it = 1:length(TimeBodyParav)
    % open file
    filename = [paraviewPath, filesep 'mooring' filesep 'mooring_' num2str(it) '.vtp'];
    fid = fopen(filename, 'w');
    % write header
    fprintf(fid, '<?xml version="1.0"?>\n');
    fprintf(fid, ['<!-- WEC-Sim Visualization using ParaView -->\n']);
    fprintf(fid, ['<!--   model: ' model ' - ran on ' simdate ' -->\n']);
    fprintf(fid, ['<!--   mooring:  MoorDyn -->\n']);
    fprintf(fid, ['<!--   time:  ' num2str(TimeBodyParav(it)) ' -->\n']);
    fprintf(fid, '<VTKFile type="PolyData" version="0.1">\n');
    fprintf(fid, '  <PolyData>\n');
    % write line info
    for iline = 1:nline
        fprintf(fid,['    <Piece NumberOfPoints="' num2str(nnode(iline)) '" NumberOfLines="' num2str(nsegment(iline)) '">\n']);
        % write points
        fprintf(fid,'      <Points>\n');
        fprintf(fid,'        <DataArray type="Float32" NumberOfComponents="3" format="ascii">\n');
        for inode = 0:nnode(iline)-1
            pt = [moorDyn.(['Line' num2str(iline)]).(['Node' num2str(inode) 'px'])(it), moorDyn.(['Line' num2str(iline)]).(['Node' num2str(inode) 'py'])(it), moorDyn.(['Line' num2str(iline)]).(['Node' num2str(inode) 'pz'])(it)];
            fprintf(fid, '          %5.5f %5.5f %5.5f\n', pt);
        end; clear pt inode
        fprintf(fid,'        </DataArray>\n');
        fprintf(fid,'      </Points>\n');
        % write lines connectivity
        fprintf(fid,'      <Lines>\n');
        fprintf(fid,'        <DataArray type="Int32" Name="connectivity" format="ascii">\n');
        count = 0;
        for isegment = 1:nsegment(iline)
            fprintf(fid, ['          ' num2str(count) ' ' num2str(count+1) '\n']);
            count = count +1;
        end; clear count isegment
        fprintf(fid,'        </DataArray>\n');
        fprintf(fid,'        <DataArray type="Int32" Name="offsets" format="ascii">\n');
        fprintf(fid, '         ');
        for isegment = 1:nsegment(iline)
            n = 2*isegment;
            fprintf(fid, ' %i', n);
        end; clear n isegment
        fprintf(fid, '\n');
        fprintf(fid,'        </DataArray>\n');
        fprintf(fid, '      </Lines>\n');
        % write cell data
        fprintf(fid,'      <CellData>\n');
        % Segment Tension
        fprintf(fid,'        <DataArray type="Float32" Name="Segment Tension" NumberOfComponents="1" format="ascii">\n');
        for isegment = 0:nsegment(iline)-1
            fprintf(fid, '          %i', moorDyn.(['Line' num2str(iline)]).(['Seg' num2str(isegment) 'Te'])(it));
        end;
        fprintf(fid, '\n');
        fprintf(fid,'        </DataArray>\n');
        fprintf(fid,'      </CellData>\n');
        % end file
        fprintf(fid, '    </Piece>\n');
    end;
    % close file
    fprintf(fid, '  </PolyData>\n');
    fprintf(fid, '</VTKFile>');
    fclose(fid);
end;
clear it iline nline nnode nsegment model t
end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function [Theta0,Re0] = lme_mass_fit_EMinit(X,Zcols,Y,ni,maskvtx,nit,prs)
% [Theta0,Re0] = lme_mass_fit_EMinit(X,Zcols,Y,ni,maskvtx,nit,prs)
%
% Starting values at each location for the linear mixed-effects iterative
% estimation. These starting values are based on the ordinary least squares
% estimators plus some fast expectation maximization iterations. This
% function should give better estimators than the lme_mass_fit_init
% function but at the cost of extra computational time.
%
% Input
% X: Ordered (according to time for each subject) design matrix (nmxp, nm
% total # of maps, p # of fixed effects).
% Zcols: Vector with the indices of the colums of X that will be considered
% as random effects.
% Y: Data matrix (nmxnv, nv #vertices) whos colums represent the ordered
% data vector (according to X) for each voxel/vertex.
% ni: Vector whose entries are the number of repeated measures for each
% subject (ordered according to X).
% maskvtx: Mask's vertices (1-based). Default [] (all vertices included).
% nit: Number of Expectation Maximization iterations. Default 5.
% prs: Number of workers for parallel computing. Default numcores;
%
% Output
% Theta0: Matrix whose colums are estimators of the covariance components at.
% each location.
% Re0: Matrix of residual errors at each location.
%
% $Revision: 1.3 $  $Date: 2016/04/08 19:39:24 $
% Original Author: Jorge Luis Bernal Rusiel
% CVS Revision Info:
%    $Author: mreuter $
%    $Date: 2016/04/08 19:39:24 $
%    $Revision: 1.3 $
% References: Bernal-Rusiel J.L., Greve D.N., Reuter M., Fischl B., Sabuncu
% M.R., 2012. Statistical Analysis of Longitudinal Neuroimage Data with Linear
% Mixed Effects Models, NeuroImage, doi:10.1016/j.neuroimage.2012.10.065.
%
tic;
if nargin < 4
    error('Too few inputs');
elseif nargin < 7
    prs = feature('numcores');
    if nargin < 6
        nit = 5;
        if nargin < 5
            maskvtx = [];
        end;
    end;
end;
nv0 = size(Y,2);
if isempty(maskvtx)
    maskvtx = 1:nv0;
end;
Y = Y(:,maskvtx);
m = length(ni);
p = size(X,2);
q = length(Zcols);
n = sum(ni);
nimax = max(ni);
nv = size(Y,2);
nth = q*(q+1)/2+1;
nDth = q*q;
ind = false(nDth,1);
for k=1:q
    for j=1:k
        ind((k-1)*q+j) = true;
    end;
end;
[Theta1,Re1] = lme_mass_fit_init(X,Zcols,Y,ni,[],prs);
display(' ');
display('Starting Expectation Maximization iterations ...');
fname = parfor_progress('init',nv);
parfor j=1:nv
    % Iterations
    W = zeros(n,nimax);
    ths = Theta1(:,j);
    D = zeros(nDth,1);
    D(ind) = ths(1:end-1);
    D = reshape(D,q,q);
    D = D + triu(D,1)';
    phisq = ths(end);
    for k=1:nit
        %Computation of W = SIGMA^-1 and H.
        posi = 1; H = 0; Term = 0;
        y = Y(:,j);
        scInvD = D\eye(q)*phisq;
        for i=1:m
            posf = posi+ni(i)-1;
            Zi = X(posi:posf,Zcols);
            Wi = (eye(ni(i))-Zi/(Zi'*Zi+scInvD)*Zi')/phisq;
            W(posi:posf,1:ni(i)) = Wi;
            Xi = X(posi:posf,:);
            Ti = Xi'*Wi;
            H = H + Ti*Xi;
            Term = Term + Ti*y(posi:posf);
            posi = posf+1;
        end;
        invH = H\eye(p);
        %Estimation
        Bhat = invH*Term;
        r = y-X*Bhat;
        Re1(:,j) = r;
        Dhat = zeros(q,q);
        phisqhat = 0;
        posi = 1;
        for i=1:m
            posf = posi+ni(i)-1;
            Zi = X(posi:posf,Zcols);
            Wi = W(posi:posf,1:ni(i));
            Xi = X(posi:posf,:);
            ri = r(posi:posf);
            bihat = D*Zi'*Wi*ri;
            Pi = Wi - Wi*Xi*invH*Xi'*Wi;
            phisqhat = phisqhat + (ri-Zi*bihat)'*(ri-Zi*bihat)+phisq*...
                trace(eye(ni(i))-phisq*Pi);
            Dhat = Dhat+bihat*bihat'+D-D*Zi'*Pi*Zi*D;
            posi = posf+1;
        end;
        D = Dhat/m;
        phisq = phisqhat/n;
    end
    vecD = vec(D);
    Theta1(:,j) = [vecD(ind);phisq];
    if mod(j,2000) == 0
        parfor_progress(fname,2000);
    end;
end
parfor_progress(fname,0);
Theta0 = zeros(nth,nv0);
Theta0(:,maskvtx) = Theta1;
Re0 = zeros(n,nv0);
Re0(:,maskvtx) = Re1;
if license('test','distrib_computing_toolbox')
    if verLessThan('matlab','8.2.0.29')
        if (matlabpool('size') > 0)
            matlabpool close;
        end;
    else
        if ~isempty(gcp('nocreate'))
            delete(gcp('nocreate'))
        end;
    end;
end;
et = toc;
display(['Total elapsed time is ' num2str(et/60) ' minutes.']);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  battistn@tcnj.edu
% Date Simulation Created: August 28, 2019
% Institution: TCNJ
%
% Date IB2d Was Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%   .
%   .
%   .
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific fiber model or example, please let Nick (battistn@tcnj.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the PENDULUM POINT MASS EXAMPLE geometry and prints associated input files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Give_Pendulum_Geometry_PointMass(Mass)

%
% Grid Parameters (MAKE SURE MATCHES IN input2d !!!)
%
Nx =  1024;      % # of Eulerian Grid Pts. in x-Direction (MUST BE EVEN!!!)
Ny =  1024;      % # of Eulerian Grid Pts. in y-Direction (MUST BE EVEN!!!)
Lx = 1.0;        % Length of Eulerian Grid in x-Direction
Ly = 1.0;        % Length of Eulerian Grid in y-Direction

% Immersed Structure Geometric / Dynamic Parameters %
ds= min(Lx/(2*Nx),Ly/(2*Ny));  % Lagrangian spacing
Lp = 0.2;                      % Length of Pendulum
struct_name = 'pendulum';      % Name for .vertex, .spring, etc files.

% x and y vectors represent pairs of coordinates describing pendulum and its point mass bob.
x = -Lp/2:Lp:Lp/2;
y = zeros(length(x));

% Angle of rotation from x-axis
t = pi/5;

% Carries out rotation
for j=1:length(x)
    xN(j) = x(j)*cos(t) - y(j)*sin(t);
    yN(j) = x(j)*sin(t) + y(j)*cos(t);
end

%
% Shifts to appropriate place in domain
%
xN = xN - xN(end) + Lx/2;
yN = yN - yN(end) + 2*Ly/3;

%
% Plots the Geometry
%
plot(x, y, 'b.-'); hold on;
plot(xN, yN, 'r.-'); hold on;
axis([0 1 0 1]);

%
% Defines xLag and yLag for printing input files below
%
xLag = xN;
yLag = yN;

% Prints .vertex file!
print_Lagrangian_Vertices(xLag,yLag,struct_name);

% Prints .mass file!
k_Mass = 2.5e6;            % 'spring' stiffness parameter for mass tethering
%Mass = 1e2;               % "MASS" value for point mass (input to this script)
print_Lagrangian_Mass_Pts(xLag,k_Mass,Mass,struct_name);

% Prints .spring file!
k_Spring = 5e7; ds_Rest = sqrt( ( xLag(1)-xLag(2)  )^2 + ( yLag(1)-yLag(2)  )^2    );
print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Rest,struct_name);

% Prints .target file!
k_Target = 5e7;
print_Lagrangian_Target_Pts(xLag,k_Target,struct_name);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints VERTEX points to a file called rubberband.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Vertices(xLag,yLag,struct_name)

    N = length(xLag);

    vertex_fid = fopen([struct_name '.vertex'], 'w');

    fprintf(vertex_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        X_v = xLag(s);
        Y_v = yLag(s);
        fprintf(vertex_fid, '%1.16e %1.16e\n', X_v, Y_v);
    end

    fclose(vertex_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints TARGET points to a file called rubberband.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name)

    N = 1; %for 3 channels

    target_fid = fopen([struct_name '.target'], 'w');

    fprintf(target_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.

    fprintf(target_fid, '%d %1.16e\n', 2, k_Target); %to choose correct lag. pt.

    fclose(target_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints TARGET points to a file called struct_name.mass
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Mass_Pts(xLag,kMass,Mass,struct_name)

    %LOOP OVER LAG PTS FOR MASS PTS IN LAGRANGIAN INDEXING

    N = 1;

    mass_fid = fopen([struct_name '.mass'], 'w');

    fprintf(mass_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.

    fprintf(mass_fid, '%d %1.16e %1.16e\n', 1, kMass, Mass);

    fclose(mass_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints BEAM (Torsional Spring) points to a file called rubberband.beam
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Beams(xLag,yLag,k_Beam,C,struct_name)

    % k_Beam: beam stiffness
    % C: beam curvature

    N = length(xLag); % NOTE: Total number of beams = Number of Total Lag Pts. - 2

    beam_fid = fopen([struct_name '.beam'], 'w');

    fprintf(beam_fid, '%d\n', 3*N );

    %spring_force = kappa_spring*ds/(ds^2);

    %BEAMS BETWEEN VERTICES
    for s = 1:N
            if s==1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',N, s, s+1,   k_Beam, C);
            elseif  s <= N-1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C);
            elseif s==N
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, 1,   k_Beam, C);
            end
    end

    for s = N+1:2*N
            if s==N+1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',2*N, s, s+1,   k_Beam, C);
            elseif  s <= 2*N-1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C);
            elseif s==2*N
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, N+1,   k_Beam, C);
            end
    end

    for s = 2*N+1:3*N
            if s==2*N+1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',3*N, s, s+1,   k_Beam, C);
            elseif  s <= 3*N-1
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C);
            elseif s==3*N
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, 2*N+1,   k_Beam, C);
            end
    end
    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints SPRING points to a file called rubberband.spring
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Springs(xLag,yLag,k_Spring,ds_Rest,struct_name)

    N = length(xLag);

    spring_fid = fopen([struct_name '.spring'], 'w');

    fprintf(spring_fid, '%d\n', 1 ); %N MUST BE EVEN!

    %spring_force = kappa_spring*ds/(ds^2);

    %SPRINGS BETWEEN VERTICES

    fprintf(spring_fid, '%d %d %1.16e %1.16e\n', 1, 2, 2*k_Spring, ds_Rest);

    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry for cylinder
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Circular_Immersed_Boundary_Geometry(ds,r,x0,y0)

% The immsersed structure is a circle %

dtheta = ds/ (2*r);
theta = 0; i=1;
while theta < 2*pi
   xLag(i) = x0 - r*cos(theta);
   yLag(i) = y0 - r*sin(theta);
   theta = theta + dtheta;
   i=i+1;
end

author: nickabattista


text:
function test_suite=test_cosmo_notify_test_skipped
% tests for cosmo_align
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function s=randstr()
    s=char(ceil(rand(1,10)*26+64));

function test_cosmo_notify_test_skipped_basics()
    notify_state=cosmo_notify_test_skipped();
    warning_state=warning();

    notify_resetter=onCleanup(@()cosmo_notify_test_skipped(notify_state));
    warning_resetter=onCleanup(@()warning(warning_state));

    % empty notified tests
    cosmo_notify_test_skipped('on');
    s=cosmo_notify_test_skipped();
    assert(isempty(s));

    % switch of warnings
    warning('off');

    db=dbstack();

    if is_moxunit(db);
        try
            moxunit_throw_test_skipped_exception('foo');
        catch
            e=lasterror();
        end
        identifier=e.identifier;
        caller=@(varargin)assertExceptionThrown(@()...
                        cosmo_notify_test_skipped(varargin{:}),identifier);
    else
        caller=@(varargin)cosmo_notify_test_skipped(varargin{:});
    end

    % one notificiation
    reason1=['a' randstr()];
    caller(reason1);
    s=cosmo_notify_test_skipped();
    assert_has_substrings(s,{reason1});

    reason2=['b' randstr()];
    caller(reason2);
    s=cosmo_notify_test_skipped();
    assert_has_substrings(s,{reason1,reason2});

function assert_has_substrings(haystack,needle)
    n=numel(haystack);
    assertEqual(n,numel(needle));
    for k=1:n
        assert(~isempty(findstr(haystack{k},needle{k})));
    end

function tf=is_moxunit(db)
    tf=false;

    nm='initTestSuite';
    main='moxunit_runtests';

    if isempty(strmatch(nm,{db.name},'exact'))
        return;
    end

    rootdir=fileparts(which(nm));
    pth=fullfile(rootdir,[main '.m']);

    tf=isequal(which(main),pth);

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
%% polyContourThick
% Below is a demonstration of the features of the |polyContourThick| function

%%
clear; close all; clc;

%% Syntax
% |[C]=polyContourThick(varargin);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function M = perform_varying_blurring(M, Sigma, options)

% perform_varying_blurring - perform a spacially varying gaussian blurring
%
%   M = perform_varying_blurring(M, Sigma, options);
%
%   Sigma(i,j) is the width of blurring (in pixels) around pixel (i,j).
%
%   Copyright (c) 2007 Gabriel Peyre

options.null = 0;
if size(M,3)>1
    for i=1:size(M,3)
        M(:,:,i) = perform_varying_blurring(M(:,:,i), Sigma, options);
    end
    return;
end

n = size(M,1);
if isfield(options, 'nb_filters')
    m = nb_filters;
else
    m = 15;
end

if std(Sigma(:))<1e-8
    % special case: in fact constant blurring
    m = 1;
end

sigma = linspace( min(Sigma(:)), max(Sigma(:)), m );
Mh = zeros(n,n,m);
for i=1:m
    Mh(:,:,i) = perform_blurring(M, sigma(i), options);
end

% combine together the filters
[tmp,I] = min( abs( repmat(Sigma,[1 1 m]) - repmat( reshape(sigma, [1 1 m]), [n n 1]) ), [], 3 );
J = (1:n^2) + (I(:)'-1)*n^2; J = reshape(J, n,n);
M = Mh(J);

author: gpeyre

repo: matlab-toolboxes


text:
function combo_test12 ( )

%*****************************************************************************80
%
%% COMBO_TEST12 tests KNAPSACK_REORDER and KNAPSACK_RATIONAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 January 2011
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  mass_limit = 26.0;
  p = [ 24.0, 13.0, 23.0, 15.0, 16.0 ]';
  w = [ 12.0,  7.0, 11.0,  8.0,  9.0 ]';

  fprintf ( 1, ' \n' );
  fprintf ( 1, 'COMBO_TEST12\n' );
  fprintf ( 1, '  KNAPSACK_REORDER reorders the knapsack data.\n' );
  fprintf ( 1, '  KNAPSACK_RATIONAL solves the rational knapsack problem.\n' );

  fprintf ( 1, ' \n' );
  fprintf ( 1, '  Object, Profit, Mass, "Profit Density"\n' );
  fprintf ( 1, ' \n' );
  for i = 1 : n
    fprintf ( 1, '  %6d  %7.3f  %7.3f  %7.3f\n', i, p(i), w(i), p(i) / w(i) );
  end

  [ p, w ] = knapsack_reorder ( n, p, w );

  fprintf ( 1, ' \n' );
  fprintf ( 1, '  After reordering by Profit Density:\n' );
  fprintf ( 1, ' \n' );
  fprintf ( 1, '  Object, Profit, Mass, "Profit Density"\n' );
  fprintf ( 1, ' \n' );
  for i = 1 : n
    fprintf ( 1, '  %6d  %7.3f  %7.3f  %7.3f\n', i, p(i), w(i), p(i) / w(i) );
  end

  fprintf ( 1, ' \n' );
  fprintf ( 1, '  Total mass restriction is %f\n', mass_limit );

  [ x, mass, profit ] = knapsack_rational ( n, mass_limit, p, w  );

  fprintf ( 1, ' \n' );
  fprintf ( 1, '  Object, Density, Choice, Profit, Mass\n' );
  fprintf ( 1, ' \n' );
  for i = 1 : n
    fprintf ( 1, '%6d  %7.3f  %7.3f  %7.3f  %7.3f\n', ...
      i, p(i) / w(i), x(i),  x(i) * p(i), x(i) * w(i) );
  end

  fprintf ( 1, ' \n' );
  fprintf ( 1, '  Total:                  %7.3f  %7.3f\n', profit, mass );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [x,y]=gpos(h_axes)
%GPOS Get current position of cusor and return its coordinates in axes with handle h_axes
% h_axes - handle of specified axes
% [x,y]  - cursor coordinates in axes h_aexs
%
% -------------------------------------------------------------------------
% Note:
%  1. This function should be called in the figure callback WindowButtonMotionFcn.
%  2. It works like GINPUT provided by Matlab,but it traces the position
%       of cursor without click and is designed for 2-D axes.
%  3. It can also work even the units of figure and axes are inconsistent,
%       or the direction of axes is reversed.
% -------------------------------------------------------------------------

% Written by Kang Zhao,DLUT,Dalian,CHINA. 2003-11-19
% E-mail:kangzhao@student.dlut.edu.cn

h_figure=gcf;

units_figure = get(h_figure,'units');
units_axes   = get(h_axes,'units');

if_units_consistent = 1;

if ~strcmp(units_figure,units_axes)
    if_units_consistent=0;
    set(h_axes,'units',units_figure); % To be sure that units of figure and axes are consistent
end

% Position of origin in figure [left bottom]
pos_axes_unitfig    = get(h_axes,'position');
width_axes_unitfig  = pos_axes_unitfig(3);
height_axes_unitfig = pos_axes_unitfig(4);

xDir_axes=get(h_axes,'XDir');
yDir_axes=get(h_axes,'YDir');

% Cursor position in figure
pos_cursor_unitfig = get( h_figure, 'currentpoint'); % [left bottom]

if strcmp(xDir_axes,'normal')
    left_origin_unitfig = pos_axes_unitfig(1);
    x_cursor2origin_unitfig = pos_cursor_unitfig(1) - left_origin_unitfig;
else
    left_origin_unitfig = pos_axes_unitfig(1) + width_axes_unitfig;
    x_cursor2origin_unitfig = -( pos_cursor_unitfig(1) - left_origin_unitfig );
end

if strcmp(yDir_axes,'normal')
    bottom_origin_unitfig     = pos_axes_unitfig(2);
    y_cursor2origin_unitfig = pos_cursor_unitfig(2) - bottom_origin_unitfig;
else
    bottom_origin_unitfig = pos_axes_unitfig(2) + height_axes_unitfig;
    y_cursor2origin_unitfig = -( pos_cursor_unitfig(2) - bottom_origin_unitfig );
end

xlim_axes=get(h_axes,'XLim');
width_axes_unitaxes=xlim_axes(2)-xlim_axes(1);

ylim_axes=get(h_axes,'YLim');
height_axes_unitaxes=ylim_axes(2)-ylim_axes(1);

x = xlim_axes(1) + x_cursor2origin_unitfig / width_axes_unitfig * width_axes_unitaxes;
y = ylim_axes(1) + y_cursor2origin_unitfig / height_axes_unitfig * height_axes_unitaxes;

% Recover units of axes,if original units of figure and axes are not consistent.
if ~if_units_consistent
    set(h_axes,'units',units_axes);
end

author: Sable

repo: mcbench-benchmarks


text:
classdef xt_gizmo < mp.extension

%   MATPOWER
%   Copyright (c) 2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end     %% properties

    methods
        function nm_class = network_model_class(obj, nm_class, task_tag, mpopt)
            nm_class = @mp.net_model_acp_test;
        end

        function dmc_elements = dmc_element_classes(obj, dmc_class, fmt, mpopt)
            dmc_elements = { @mp.dmce_gizmo_mpc2 };
        end

        function dm_elements = dm_element_classes(obj, dm_class, task_tag, mpopt)
            dm_elements = { @mp.dme_gizmo };
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [tc,relres,iter,frec,cd] = franalasso(F,f,lambda,varargin)
%FRANALASSO  Frame LASSO regression
%   Usage: tc = franalasso(F,f,lambda)
%          tc = franalasso(F,f,lambda,C,tol,maxit)
%          [tc,relres,iter,frec,cd] = franalasso(...)
%
%   Input parameters:
%       F        : Frame definition
%       f        : Input signal
%       lambda   : Regularisation parameter, controls sparsity of the solution
%       C        : Step size of the algorithm.
%       tol      : Reative error tolerance.
%       maxit    : Maximum number of iterations.
%   Output parameters:
%       tc       : Thresholded coefficients
%       relres   : Vector of residuals.
%       iter     : Number of iterations done.
%       frec     : Reconstructed signal
%       cd       : The min ||c||_2 solution using the canonical dual frame.
%
%   `franalasso(F,f,lambda)` solves the LASSO (or basis pursuit denoising)
%   regression problem for a general frame: minimize a functional of the
%   synthesis coefficients defined as the sum of half the $l^2$ norm of the
%   approximation error and the $l^1$ norm of the coefficient sequence, with
%   a penalization coefficient `lambda` such that
%
%   .. argmin lambda||c||_1 + 1/2||Fc - f||_2^2
%
%   .. math:: \arg\!\min_c \ \lambda ||c||_1 +  \frac{1}{2}||Fc - f||_2^2
%
%   The solution is obtained via an iterative procedure, called Landweber
%   iteration, involving iterative soft thresholdings.
%
%   The following flags determining an algorithm to be used are recognized
%   at the end of the argument list:
%
%   'ista'
%     The basic (Iterative Soft Thresholding) algorithm given *F* and *f*
%     and parameters *C*$>0$, `lambda` $>0$ acts as follows::
%
%        Initialize c
%        repeat until stopping criterion is met
%          c <- soft(c + F*(f - Fc)/C,lambda/C)
%        end
%
%   'fista'
%     The fast version of the previous. This is the default option.
%     ::
%
%        Initialize c(0),z,tau(0)=1,n=1
%        repeat until stopping criterion is met
%           c(n) <- soft(z + F*(f - Fz)/C,lambda/C)
%           tau(n) <- (1+sqrt(1+4*tau(n-1)^2))/2
%           z   <- c(n) + (c(n)-c(n-1))(tau(n-1)-1)/tau(n)
%           n   <- n + 1
%        end
%
%   `[tc,relres,iter] = franalasso(...)` returns the residuals *relres* in
%   a vector and the number of iteration steps done *iter*.
%
%   `[tc,relres,iter,frec,cd] = franalasso(...)` returns the reconstructed
%   signal from the coefficients, *frec* and coefficients *cd* obtained by
%   analysing using the canonical dual system.
%   Note that this requires additional computations.
%
%   The relationship between the output coefficients and *frec* is
%   given by ::
%
%     frec = frsyn(F,tc);
%
%   The function takes the following optional parameters at the end of
%   the line of input arguments:
%
%   'C',cval
%      Landweber iteration parameter: must be larger than square of upper
%      frame bound. Default value is the upper frame bound.
%
%   'tol',tol
%      Stopping criterion: minimum relative difference between norms in
%      two consecutive iterations. Default value is 1e-2.
%
%   'maxit',maxit
%      Stopping criterion: maximal number of iterations to do. Default
%      value is 100.
%
%   'print'
%      Display the progress.
%
%   'quiet'
%      Don't print anything, this is the default.
%
%   'printstep',p
%      If 'print' is specified, then print every p'th iteration. Default
%      value is 10;
%
%   'debias'
%      Performs debiasing (improves approximation accuracy) of the solution
%      using the conjugate gradient algorithm. This amounts to doing
%      best-fit with respect to a subdictionary selected by ISTA.
%
%   'pcgmaxit',pcgmaxit
%      Maximum allowed number of iterations for pcg. Only used together
%      with the 'debias' flag. The default value is 100.
%
%   'pcgtol',pcgtol
%      Tolerance of pcg. Only used together with the 'debias' flag.
%      The default value is 1e-6.
%
%   The parameters *C*, *itermax* and *tol* may also be specified on the
%   command line in that order: `franalasso(F,x,lambda,C,tol,maxit)`.
%
%   **Note**: If you do not specify *C*, it will be obtained as the upper
%   framebound. Depending on the structure of the frame, this can be an
%   expensive operation.
%
%   Examples:
%   ---------
%
%   The following example shows how |franalasso| produces a sparse
%   representation of a test signal *greasy*:::
%
%      f = greasy;
%      % Gabor frame with redundancy 8
%      F = frame('dgtreal','gauss',64,512);
%      % Choosing lambda (weight of the sparse regularization param.)
%      lambda = 0.1;
%      % Solve the basis pursuit problem
%      [c1,~,~,frec1,cd] = franalasso(F,f,lambda);
%      % Solve again with debiasing enabled
%      [c2,~,~,frec2,cd] = franalasso(F,f,lambda,'debias');
%      % Plot sparse coefficients
%      figure(1); plotframe(F,c1,'dynrange',50);
%      % Plot debiased coefficients
%      figure(2); plotframe(F,c2,'dynrange',50);
%      % Plot coefficients obtained by applying an analysis operator of a
%      % dual Gabor system to f
%      figure(3); plotframe(F,cd,'dynrange',50);
%
%      % Normalized MSE of the approximation
%      % frec is obtained by applying the synthesis operator of frame F
%      % to sparse coefficients c.
%      fprintf('Normalized MSE:                  %.2f dB\n',20*log10(norm(f-frec1)/norm(f)));
%      fprintf('Normalized MSE (with debiasing): %.2f dB\n',20*log10(norm(f-frec2)/norm(f)));
%
%      % Compare decay of coefficients sorted by absolute values
%      % (compressibility of coefficients)
%      figure(3);
%      semilogx([sort(abs(c1),'descend')/max(abs(c1)),...
%      sort(abs(c2),'descend')/max(abs(c2)),...
%      sort(abs(cd),'descend')/max(abs(cd))]);
%      legend({'sparsified coefficients','with debiasing','dual system coefficients'});
%
%   See also: frame, frsyn, framebounds, franabp, franagrouplasso
%
%   References: dademo04 beck09 finowr07

%   AUTHOR : Bruno Torresani.
%   TESTING: OK

%   XXX Removed Remark: When the frame is an orthonormal basis, the solution
%   is obtained by soft thresholding of the basis coefficients, with
%   threshold lambda.  When the frame is a union of orthonormal bases, the
%   solution is obtained by applying soft thresholding cyclically on the
%   basis coefficients (BCR algorithm)
%

complainif_notenoughargs(nargin,2,'FRANALASSO');
complainif_notvalidframeobj(F,'FRANALASSO');

if sum(size(f)>1)>1
  error('%s: Too many input channels.',upper(mfilename));
end

% Define initial value for flags and key/value pairs.
definput.keyvals.C=[];
definput.keyvals.tol=1e-2;
definput.keyvals.maxit=100;
definput.keyvals.printstep=10;
definput.flags.print={'print','quiet'};
definput.flags.algorithm={'fista','ista'};
definput.flags.debias = {'nodebias','debias'};
definput.keyvals.pcgmaxit=100;
definput.keyvals.pcgtol=1e-6;
[flags,kv]=ltfatarghelper({'C','tol','maxit'},definput,varargin);

% Accelerate frame, we will need it repeatedly
Ls = size(f,1);
F=frameaccel(F,Ls);
L=F.L;
f = postpad(f,L);

% Use the upper framebound as C
if isempty(kv.C)
  [~,kv.C] = framebounds(F,L);
end;

% Initialization of thresholded coefficients
% frana is used instead of F.frana to get the correct zero padding of f
c0 = frana(F,f);

% Various parameter initializations
threshold = lambda/kv.C;

tc0 = c0;
relres = 1e16;
iter = 0;

if flags.do_ista
   % Main loop
   while ((iter < kv.maxit)&&(relres >= kv.tol))
       tc = c0 - F.frana(F.frsyn(tc0));
       tc = tc0 + tc/kv.C;
       tc = thresh(tc,threshold,'soft');
       relres = norm(tc(:)-tc0(:))/norm(tc0(:));
       tc0 = tc;
       iter = iter + 1;
       if flags.do_print
         if mod(iter,kv.printstep)==0
           fprintf('Iteration %d: relative error = %f\n',iter,relres);
         end;
       end;
   end
elseif flags.do_fista
   tz0 = c0;
   tau0 = 1;
   % Main loop
   while ((iter < kv.maxit)&&(relres >= kv.tol))
       tc = c0 - F.frana(F.frsyn(tz0));
       tc = tz0 + tc/kv.C;
       tc = thresh(tc,threshold,'soft');

       tau = 1/2*(1+sqrt(1+4*tau0^2));
       tz0 = tc + (tau0-1)/tau*(tc-tc0);
       relres = norm(tc(:)-tc0(:))/norm(tc0(:));
       tc0 = tc;
       tau0 = tau;
       iter = iter + 1;
       if flags.do_print
         if mod(iter,kv.printstep)==0
           fprintf('Iteration %d: relative error = %f\n',iter,relres);
         end;
       end;
   end
end

if flags.do_debias
   mask = double(abs(tc)>0);

   % Simulate the pseudoinverse of the synthesis operator of the
   % reduced system
   if numel(find(mask>0)) >= L
       % Invert the frame matrix
       A=@(x) F.frsyn(mask.*F.frana(x));

       [fout,flag,~,iter,relres]=pcg(A,f,kv.pcgtol,kv.pcgmaxit);
       tc=mask.*F.frana(fout);
   else
       % Invert the Gram matrix
       A=@(x) mask.*F.frana(F.frsyn(mask.*x));
       [tc,flag,~,iter,relres]=pcg(A,mask.*F.frana(f),kv.pcgtol,kv.pcgmaxit);
       tc=mask.*tc;
   end
end

% Optional reconstruction
if nargout>3
   frec = postpad(F.frsyn(tc),Ls);
end;

% Calculate coefficients using the canonical dual system
% May be conviniently used for comparison
if nargout>4
  try
     Fd = framedual(F);
     cd = frana(Fd,f);
  catch
     warning(sprintf(['%s: Dual frame is not available. Using franaiter'],...
                     upper(mfilename)));
     cd = franaiter(F,f);
  end
end;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% Function to estimate glottal closure instants (GCIs) using an extension
% of Thomas Drugman's SEDREAMS algorithm, optimised to better handle
% non-modal phonation types
%
% Octave compatible
%
% Description
%  Function to extract GCIs using an adapted version of the SEDREAMS
%  algorithm which is optimised for non-modal voice qualities. Ncand maximum
%  peaks are selected from the LP-residual signal in the interval defined by
%  the mean-based signal. A dynamic programming algorithm is then used to
%  select the optimal path of GCI locations. Then a post-processing method,
%  using the output of a resonator applied to the residual signal, is
%  carried out to remove false positives occurring in creaky speech regions
%
% Inputs
%  x               : [samples] [Nx1] Speech signal
%  fs              : [Hz]      [1x1] sampling frequency
%  F0mean          : [Hz]      [1x1] Mean fundamental frequency
%  creak           : [binary]  [Nx1] Creak decision (OPTIONAL)
%
% Outputs
%  GCI             : [s] [Mx1] Glottal closure instants
%
% Example
%  GCI = se_vq(x,fs);
%
% References
%  [1] Kane, J., Gobl, C., (2013) `Evaluation of glottal closure instant
%       detection in a range of voice qualities', Speech Communication
%       55(2), pp. 295-314.
%
% Copyright (c) 2013 Trinity College Dublin
%
% License
%  This code is a part of the Voice Analysis Toolkit with the following
%  licence:
%  The software product (and any modifications thereof) is distributed under
%  a dual licence: an open source license for individual, non-commercial
%  purposes, and a commercial license. The opensource licence under which
%  the product is distributed is GNU GPL v2. For individual users, this
%  licence suits their use as these are not distributing proprietary
%  modifications, additions to, or derivatives of the product and don't
%  require legal protection of a commercial licence. For commercial users,
%  where open source does not meet their requirements, we offer commercial
%  licensing of the product. A commercial license permits customers to
%  modify, add or produce derivative products without the obligation of
%  making the subsequent code open source. For more information regarding
%  our commercial licence, please contact john.whelan@tcd.ie
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  John Kane kanejo@tcd.ie
%

function GCI = se_vq(x,fs,F0mean,creak)

%% Settings
F0min=20;
F0max=500;
T0mean = fs/F0mean; % Rough period length for mean-based signal

winLen = 25; % window length in ms
winShift = 5; % window shift in ms
LPC_ord = round(fs/1000)+2; % LPC order
Ncand=5; % Number of candidate GCI residual peaks to be considered in the dynamic programming

trans_wgt=1; % Transition cost weight
relAmp_wgt=0.3; % Local cost weight

repNum=2;
removeThresh=0.4; % Threshold for removing false GCIs
search_reg=1.3/1000*fs;

%% Calculate LP-residual and extract N maxima per mean-based signal determined intervals
res = lpcresidual(x,winLen/1000*fs,winShift/1000*fs,LPC_ord); % Get LP residual
rep = RCVD_reson_GCI(res,fs,F0mean); % Get resonator output
MBS = get_MBS(x,fs,T0mean); % Extract mean based signal
interval = get_MBS_GCI_intervals(MBS,fs,T0mean,F0max); % Define search intervals
[GCI_N,GCI_relAmp] = search_res_interval_peaks(res,interval,Ncand); % Find residual peaks
GCI = RESON_dyProg_mat(GCI_relAmp',GCI_N',F0mean,x,fs,trans_wgt,relAmp_wgt); % Do dynamic programming

GCI(GCI>length(x))=[];

%% Remove false alarms as weak peaks in resonator output
if nargin > 3
    GCI = GCI_creak_postproc(GCI,creak,search_reg,rep,removeThresh,repNum);
end

GCI = (GCI-1)/fs;

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function prior = uniformPriorSetBounds(prior, bounds)

% UNIFORMPRIORSETBOUNDS Set uniform prior bounds.

% PRIOR

prior.a = bounds(1);
prior.b = bounds(2);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
filename='Cantilever_tetrahedra_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'IPOPT';
incrementFactor = 1;
designVariable = 'Density';
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.05;
Perimeter_target = 1;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-2;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%Pad The Original Image
function [r,c,padIMG]=padO(img)
[iLength,iWidth]=size(img);
iDiag=sqrt(iLength^2 + iWidth^2);
LengthPad=ceil(iDiag - iLength) + 2;
WidthPad=ceil(iDiag - iWidth) + 2;
padIMG = zeros(iLength+LengthPad, iWidth+WidthPad);
%figure,imshow(padIMG);
padIMG(ceil(LengthPad/2):(ceil(LengthPad/2)+iLength-1), ...
       ceil(WidthPad/2):(ceil(WidthPad/2)+iWidth-1)) = img;
%Size of Original padIMG
   [r,c]=size(padIMG);
end


text:
function rank = kel2rank(kel)
%KEL2RANK Convert temperature from kelvin to rankine
%
%  rank = KEL2RANK(kel) converts temperature from kelvin to rankine.
%
%  See also RANK2KEL, KEL2CEL, KEL2FAR.

rank = kel * 1.8;


text:
% [INPUT]
% data = A float t-by-n matrix (-Inf,Inf) representing the model input.
% a = A float [0.01,0.10] representing the target quantile.
% k = An integer [1,60] representing the target lag.
% cis = A float {0.005;0.010;0.025;0.050;0.100} representing the significance level of confidence intervals (optional, default=0.050).
% cif = A float {0.00;0.01;0.03;0.05;0.10;0.15;0.20;0.30} representing the minimum subsample size, as a fraction, for confidence intervals (optional, default=0.10).
%
% [OUTPUT]
% cq = A float (-Inf,Inf) representing the cross-quantilogram.
% ci = A row vector of floats (-Inf,Inf) of length 2 representing the lower and upper confidence intervals.
%
% [NOTES]
% The model computes partial cross-quantilograms when n is greater than 2 using exogenous variables from 3 to n.

function [cq,ci] = cross_quantilograms_sn(varargin)

    persistent ip;

    if (isempty(ip))
        ip = inputParser();
        ip.addRequired('data',@(x)validateattributes(x,{'double'},{'real' 'finite' '2d' 'nonempty'}));
        ip.addRequired('a',@(x)validateattributes(x,{'double'},{'real' 'finite' '>=' 0.01 '<=' 0.10 'scalar'}));
        ip.addRequired('k',@(x)validateattributes(x,{'double'},{'real' 'finite' 'integer' '>=' 1 '<=' 60 'scalar'}));
        ip.addOptional('cis',0.050,@(x)validateattributes(x,{'double'},{'real' 'finite' 'scalar'}));
        ip.addOptional('cif',0.10,@(x)validateattributes(x,{'double'},{'real' 'finite' 'scalar'}));
    end

    ip.parse(varargin{:});

    ipr = ip.Results;
    data = validate_input_data(ipr.data);
    a = ipr.a;
    k = ipr.k;
    [civ,cif] = validate_input_ci(ipr.cis,ipr.cif);

    nargoutchk(2,2);

    [cq,ci] = cross_quantilograms_sn_internal(data,a,k,civ,cif);

end

function [cq,ci] = cross_quantilograms_sn_internal(x,a,k,civ,cif)

    [t,n] = size(x);
    len = max(round(cif * t,0),1);
    partial = n > 2;

    cq_sn = zeros(t,1);

    if (partial)
        for i = len:t
            x_i = x(1:i,:);
            x_t = size(x_i,1);
            x_t_mk = x_t - k;

            q_sn = (x_i <= repmat(gumbel_quantile(x_i,a),x_t,1)) - (ones(x_t,n) .* a);

            d_sn = zeros(x_t_mk,n);
            d_sn(:,1) = q_sn(k+1:x_t,1);
            d_sn(:,2:n) = q_sn(1:x_t_mk,2:n);

            h_sn = d_sn.' * d_sn;

            if (det(h_sn) <= 1e-08)
                hi_sn = pinv(h_sn);
            else
                hi_sn = inv(h_sn);
            end

            cq_sn(i) = -hi_sn(1,2) / sqrt(hi_sn(1,1) * hi_sn(2,2));
        end
    else
        for i = len:t
            x_i = x(1:i,:);
            x_t = size(x_i,1);
            x_t_mk = x_t - k;

            q_sn = (x_i <= repmat(gumbel_quantile(x_i,a),x_t,1)) - (ones(x_t,n) .* a);

            d_sn = zeros(x_t_mk,n);
            d_sn(:,1) = q_sn(k+1:x_t,1);
            d_sn(:,2:n) = q_sn(1:x_t_mk,2:n);

            h_sn = d_sn.' * d_sn;

            cq_sn(i) = h_sn(1,2) / sqrt(h_sn(1,1) * h_sn(2,2));
        end
    end

    cq = cq_sn(end);

    cqc = (cq_sn - repmat(cq,t,1)) .* (1:t).';
    cv0_bb = cqc(len:t);
    cv0_sn = (t * cq^2) / ((1 / t^2) * (cv0_bb.' * cv0_bb));
    cv0 = sqrt(civ * (cq^2 / cv0_sn));
    ci = [-cv0 cv0];

end

function q = gumbel_quantile(x,p)

    index = 1 + ((size(x,1) - 1) * p);
    low = floor(index);
    high = ceil(index);

    x = sort(x);
    x_low = x(low,:);
    x_high = x(high,:);

    h = max(index - low,0);
    q = (h .* x_high) + ((1 - h) .* x_low);

end

function [civ,cif] = validate_input_ci(cis,cif)

    persistent cif_allowed;
    persistent cis_allowed;
    persistent v;

    if (isempty(cif_allowed))
        cif_allowed = [0.00 0.01 0.03 0.05 0.10 0.15 0.20 0.30];
    end

    if (isempty(cis_allowed))
        cis_allowed = [0.005 0.010 0.025 0.050 0.100];
    end

    if (isempty(v))
        v = [
            129.15490  99.44085  66.00439 45.43917 28.06313;
            131.82880 101.21590  66.49058 45.48538 28.31850;
            131.83560 101.31700  67.55465 46.02739 28.51908;
            135.24000 103.32090  68.20319 46.48723 28.82658;
            139.34750 106.80970  71.07829 48.55132 30.01695;
            150.60740 115.90190  76.13708 51.38898 31.66900;
            166.53440 127.71000  83.45021 55.85094 34.03793;
            206.01210 155.00120 101.10960 67.53906 40.84930
        ];
    end

    [cis_ok,j] = ismember(cis,cis_allowed);

    if (~cis_ok)
        cis_allowed_text = [sprintf('%.3f',cis_allowed(1)) sprintf(', %.3f',cis_allowed(2:end))];
        error(['The value of ''cis'' is invalid. Expected input to have one of the following values: ' cis_allowed_text '.']);
    end

    [cif_ok,i] = ismember(cif,cif_allowed);

    if (~cif_ok)
        cif_allowed_text = [sprintf('%.2f',cif_allowed(1)) sprintf(', %.2f',cif_allowed(2:end))];
        error(['The value of ''cif'' is invalid. Expected input to have one of the following values: ' cif_allowed_text '.']);
    end

    civ = v(i,j);

end

function data = validate_input_data(data)

    n = size(data,2);

    if (n < 2)
        error('The value of ''data'' is invalid. Expected input to be a matrix with at least 2 columns.');
    end

end

author: TommasoBelluzzo

repo: SystemicRisk


text:
function [ fea, out ] = ex_eddycurrents2( varargin )
%EX_EDDYCURRENTS2 3D Eddy currents test example.
%
%   [ FEA, OUT ] = EX_EDDYCURRENTS2( VARARGIN ) 3D Eddy currents test
%   example for vector elements (Nedelec). Accepts the following
%   property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       icase       integer 1/{1}          Predefined test case
%       hmax        scalar {0.05}          Grid cell size
%       sfun        string {sf_simp_N1}    Vector shape function (Nedelec)
%       iplot       scalar 0/{1}           Plot solution (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = { 'icase',    1;
            'hmax',     0.05;
            'sfun',     'sf_simp_N1';
            'iplot',    1;
            'tol',      0.01;
            'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});

% Geometry and grid generation.
fea.sdim = {'x','y','z'};
fea.geom.objects = { gobj_rectangle() };
fea.grid = hex2tet( blockgrid(ceil(1/opt.hmax)), 2 );

% Problem definition.
fea = addphys( fea, @customeqn, {'E'} );
fea.phys.ce.eqn.seqn = {'Ec_c + E_t = 0'};
fea.phys.ce.sfun = {opt.sfun};
if( opt.icase==2 )   % Set homogenous Neumann BCs.
  [fea.phys.ce.bdr.coef{5}{:}] = deal(0);
end

% Parse problem.
fea = parsephys( fea );
fea = parseprob( fea );

% Define manual source term with rank 3 to each vector compnent.
switch( opt.icase )
  case 1
    f1 = '0';
    f2 = '0';
    f3 = '(2*pi^2 + 1)*sin(pi*x)*sin(pi*y)';
end
fea.eqn.f.form{1} = 1;
fea.eqn.f.coef{1} = {cat(3,{f1},{f2},{f3})};

% Solve problem.
fea.sol.u = solvestat( fea, 'fid', opt.fid );

% Postprocessing.
if( opt.iplot )
  postplot( fea, 'sliceexpr', 'Ec#1' )
  title( 'Curl solution #1')
end

% Error checking.
switch( opt.icase )
  case 1
    refsol = { '0', ...
               '0', ...
               'sin(pi*x)*sin(pi*y)', ...
               'pi*cos(pi*y)*sin(pi*x)', ...
               '-pi*cos(pi*x)*sin(pi*y)', ...
               '0' };
end
out.err = [ intsubd(['(',refsol{1},'-E#1)^2'],fea), ...
            intsubd(['(',refsol{2},'-E#2)^2'],fea), ...
            intsubd(['(',refsol{3},'-E#3)^2'],fea), ...
            intsubd(['(',refsol{4},'-Ec#1)^2'],fea), ...
            intsubd(['(',refsol{5},'-Ec#2)^2'],fea), ...
            intsubd(['(',refsol{6},'-Ec#3)^2'],fea) ];
out.pass = all( out.err < opt.tol );
if( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
function page3dEBSD( api )

getEBSD = api.Export.getOptEBSD;
setEBSD = api.Export.setOptEBSD;

gui = localCreatePage();

api.setWizardTitle('Import Wizard EBSD');
api.setWizardDescription('3d Data Z-Values');

api.setLeavePageCallback(@leavePage);
api.setGotoPageCallback(@gotoPage);

% set(gui.h3dRadio,'Callback',@localUpdate3d)
set(gui.h3dGroup  ,'SelectionChangeFcn',@localUpdate3d)
set(gui.hUpdate   ,'Callback'          ,@localUpdate3d)
% set(gui.hTable,'DataChangedCallback',@localUpdate3d)

  function nextPage = leavePage

    % store z-values, while DataChangedCallback fails somehow
    localUpdate3d([],[]);

    nextPage = @pageCS;

  end

  function gotoPage

    api.Progress.enableNext(true);

    localUpdateGUI();

  end

  function localUpdate3d(source,event)

    setEBSD('is3d',get(gui.h3dRadio(2),'value'));

    tdata = cell(get(gui.hTable,'Data'));
    Z = tdata(:,2);
    x = cellfun('isclass',Z,'char');
    Z(x) = cellfun(@str2num,Z(x),'UniformOutput',false);
    Z(cellfun('isempty',Z)) = {0};

    setEBSD('Z',Z);

    localUpdateGUI();

  end

  function localUpdateGUI()

    set(gui.h3dRadio(1),'Value',~getEBSD('is3d'));
    set(gui.h3dRadio(2),'Value',getEBSD('is3d'));

    files = api.getFiles();
    [pathes,fnames,ext] = cellfun(@fileparts,files,'Uniformoutput',false);
    fnames = strcat(fnames,ext);

    Z = getEBSD('Z');
    if isempty(Z) || numel(Z) ~= numel(files)
      Z = num2cell(1:numel(files));
    end

    state = {'off','on'};

    set(gui.hTable,  'data',[fnames(:),Z(:)]);
    set(gui.hTable,  'Visible',getEBSD('is3d'));
    set(gui.hUpdate, 'Visible',state{1+getEBSD('is3d')});

  end

  function gui = localCreatePage()

    page = api.hPanel;

    h    = api.Spacings.PageHeight;
    w    = api.Spacings.PageWidth;
    m    = api.Spacings.Margin;
    bH   = api.Spacings.ButtonHeight;
    bW   = api.Spacings.ButtonWidth;
    fs   = api.Spacings.FontSize;

    group3d = uibuttongroup(...
      'parent',page,...
      'Units','pixels',...
      'FontSize',fs,...
      'position',[1,0,w,h],...
      'Title','3d EBSD Data');

    is3d(1) = uicontrol(...
      'Parent',group3d,...
      'Style','radio',...
      'FontSize',fs,...
      'String',' All data files lie on the same surface',...
      'Value',1,...
      'position',[m h-m-1.5*bH w-2*m bH]);

    is3d(2) = uicontrol(...
      'Parent',group3d,...
      'Style','radio',...
      'FontSize',fs,...
      'String',' Each data file is a serial section',...
      'Value',0,...
      'position',[m h-m-2.5*bH w-2*m bH]);

    update = uicontrol(...
      'Parent',group3d,...
      'style','pushbutton',...
      'FontSize',fs,...
      'String','Update Table',...
      'Position',[w-bW-m h-m-2.5*bH bW bH]);

    if verLessThan('matlab','7.6'),  v0 = {}; else  v0 = {'v0'}; end

    cdata =  {'Z-Layer Data Source','Z-Value'};

    table = uitable(v0{:},'parent',page,...
      'Data',{'filename',0},...
      'Position',[m,m,w-2*m,h-3*bH-2*m],...
      'ColumnNames',cdata,...
      'ColumnWidth',[w/2-10-2*m]);

    gui.h3dGroup = group3d;
    gui.h3dRadio = is3d;
    gui.hTable   = table;
    gui.hUpdate  = update;

  end

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [  ] = rpn_detect(  )

    if p.rpn_detect
        objpos = annolist(i).annorect.objpos;
        objpos = single([objpos.x objpos.y]);

        pos = int32((objpos*scale_factor - half_stride)/stride) + 1;
        row = pos(2);
        col = pos(1);

        anchors = squeeze(rpn_prob(row, col, :));
        [~,anch] = max(anchors);

        t = squeeze(rpn_bbox(row, col, anch, :));

        anchor_types = single([ 1, 130; 1, 211; 2, 153; 3, 125; 4, 97]);
        w_a = anchor_types(anch, 2);
        h_a = w_a * anchor_types(anch, 1);
        x_a = objpos(1);
        y_a = objpos(2);

        x = x_a + w_a*t(1);
        y = y_a + h_a*t(2);
        w = w_a * exp(t(3));
        h = h_a * exp(t(4));

        extra_margin = 0;
        w = w + extra_margin;
        h = h + extra_margin;

        x_r = x - w/2;
        y_r = y - h/2;

        rect_orig = single([x_r y_r x_r+w y_r+h]);
        rect = int32((rect_orig*scale_factor - half_stride)/stride) + 1;
        sm_height = size(rpn_prob, 1);
        sm_width = size(rpn_prob, 2);
        rect(1) = max(rect(1), 1);
        rect(2) = max(rect(2), 1);
        rect(3) = min(rect(3), sm_width);
        rect(4) = min(rect(4), sm_height);

        if bVis
            scmap = visualise_scoremap( rpn_prob );
            figure(3);
            imshow(scmap);

            figure(4);
            imagesc(im_orig);
            rectangle('Position', [x_r y_r w h]);
        end
    end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
% std_prepare_neighbors() - prepare Fieldtrip channel neighbor structure.
%                           Only prepare the structure if necessary based
%                           on statistical options in STUDY.etc.statistics.
%                           Use the 'force' option to force preparing the
%                           matrix.
%
% Usage:
%   >> [STUDY neighbors] = std_prepare_neighbors( STUDY, ALLEEG, 'key', val)
%
% Inputs:
%   STUDY        - an EEGLAB STUDY set of loaded EEG structures
%   ALLEEG       - ALLEEG vector of one or more loaded EEG dataset structures
%
% Optional inputs:
%   'force'        - ['on'|'off'] force generating the structure irrespective
%                    of the statistics options selected. Default is 'off'.
%   'channels'     - [cell] list of channels to include in the matrix
%   'method'       - [string] method for preparing. See ft_prepare_neighbors
%   'neighbordist' - [float] max distance. See ft_prepare_neighbors
%
% Note: other ft_prepare_neighbors fields such as template, layout may
%       also be used as optional keywords.
%
% Outputs:
%   STUDY       - an EEGLAB STUDY set of loaded EEG structures
%   neighbors   - Fieldtrip channel neighbour structure
%
% Author: Arnaud Delorme, SCCN, UCSD, 2012-
%
% See also: statcondfieldtrip()

% Copyright (C) Arnaud Delorme
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [STUDY neighbors] = std_prepare_neighbors(STUDY, ALLEEG, varargin);

neighbors = [];
if nargin < 2
    return;
end;

[opt addopts] = finputcheck( varargin, {  'force'    'string'  { 'on','off' }   'off';
                                          'channels' 'cell'    {}               {} }, 'std_stat', 'ignore');

if strcmpi(opt.force, 'on') || (strcmpi(STUDY.etc.statistics.fieldtrip.mcorrect, 'cluster') && ...
        strcmpi(STUDY.etc.statistics.mode, 'fieldtrip') && (strcmpi(STUDY.etc.statistics.groupstats, 'on') || strcmpi(STUDY.etc.statistics.condstats, 'on')))

    EEG = eeg_emptyset;
    EEG.chanlocs = eeg_mergelocs(ALLEEG.chanlocs);
    if isempty(EEG.chanlocs)
        disp('std_prepare_neighbors: cannot prepare channel neighbour structure because of empty channel structures');
        return;
    end;

    if ~isempty(STUDY.etc.statistics.fieldtrip.channelneighbor) && isempty(addopts) && ...
        length(STUDY.etc.statistics.fieldtrip.channelneighbor) == length(EEG.chanlocs)

        disp('Using stored channel neighbour structure');
        neighbors = STUDY.etc.statistics.fieldtrip.channelneighbor;

    else

        if ~isempty(opt.channels)
            indChans = eeg_chaninds(EEG, opt.channels);
            EEG.chanlocs = EEG.chanlocs(indChans);
        end;

        EEG.nbchan   = length(EEG.chanlocs);
        EEG.data     = zeros(EEG.nbchan,100,1);
        EEG.trials   = 1;
        EEG.pnts     = 100;
        EEG.xmin     = 0;
        EEG.srate    = 1;
        EEG.xmax     = 99;
        EEG = eeg_checkset(EEG);
        tmpcfg = eeglab2fieldtrip(EEG, 'preprocessing', 'none');

        % call the function that find channel neighbors
        % ---------------------------------------------
        addparams = eval( [ '{' STUDY.etc.statistics.fieldtrip.channelneighborparam '}' ]);
        for index = 1:2:length(addparams)
            tmpcfg = setfield(tmpcfg, addparams{index}, addparams{index+1});
        end;
        for index = 1:2:length(addopts)
            tmpcfg = setfield(tmpcfg, addopts{index}, addopts{index+1});
        end;
        warning off;
        cfg.neighbors = ft_prepare_neighbours(tmpcfg, tmpcfg);
        warning on;
        neighbors = cfg.neighbors;

    end;

    STUDY.etc.statistics.fieldtrip.channelneighbor = neighbors;
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
clc;
close all;
clearvars;

rng default;
A = mat_selector();
k = 30;
propack_lanbdpro(A, k);

author: indigits

repo: sparse-plex


text:
function Ntot = PQloud (Ehs, Ver, Mod)
% Calculate the loudness

% P. Kabal $Revision: 1.1 $  $Date: 2003/12/07 13:35:09 $

e = 0.23;

persistent Nc s Et Ets Version Model

if (~strcmp (Ver, Version) | ~strcmp (Mod, Model))
    Version = Ver;
    Model = Mod;
    if (strcmp (Model, 'FFT'))
        [Nc, fc] = PQCB (Version);
        c = 1.07664;
    else
        [Nc, fc] = PQFB;
        c = 1.26539;
    end
    E0 = 1e4;
    Et = PQ_enThresh (fc);
    s = PQ_exIndex (fc);
    for (m = 0:Nc-1)
        Ets(m+1) = c * (Et(m+1) / (s(m+1) * E0))^e;
    end
end

sN = 0;
for (m = 0:Nc-1)
    Nm = Ets(m+1) * ((1 - s(m+1) + s(m+1) * Ehs(m+1) / Et(m+1))^e - 1);
    sN = sN + max(Nm, 0);
end
Ntot = (24 / Nc) * sN;

%====================
function s = PQ_exIndex (f)
% Excitation index

N = length (f);
for (m = 0:N-1)
    sdB = -2 - 2.05 * atan(f(m+1) / 4000) - 0.75 * atan((f(m+1) / 1600)^2);
    s(m+1) = 10^(sdB / 10);
end

%--------------------
function Et = PQ_enThresh (f)
% Excitation threshold

N = length (f);
for (m = 0:N-1)
    EtdB = 3.64 * (f(m+1) / 1000)^(-0.8);
    Et(m+1) = 10^(EtdB / 10);
end

author: stephencwelch

repo: Perceptual-Coding-In-Python


text:
function numgrad = computeNumericalGradient(J, theta)
%COMPUTENUMERICALGRADIENT Computes the gradient using "finite differences"
%and gives us a numerical estimate of the gradient.
%   numgrad = COMPUTENUMERICALGRADIENT(J, theta) computes the numerical
%   gradient of the function J around theta. Calling y = J(theta) should
%   return the function value at theta.

% Notes: The following code implements numerical gradient checking, and
%        returns the numerical gradient.It sets numgrad(i) to (a numerical
%        approximation of) the partial derivative of J with respect to the
%        i-th input argument, evaluated at theta. (i.e., numgrad(i) should
%        be the (approximately) the partial derivative of J with respect
%        to theta(i).)
%

numgrad = zeros(size(theta));
perturb = zeros(size(theta));
e = 1e-4;
for p = 1:numel(theta)
    % Set perturbation vector
    perturb(p) = e;
    loss1 = J(theta - perturb);
    loss2 = J(theta + perturb);
    % Compute Numerical Gradient
    numgrad(p) = (loss2 - loss1) / (2*e);
    perturb(p) = 0;
end

end

author: schneems

repo: Octave

save_path: github-repos/MATLAB/schneems-Octave


text:
function writeDMAT(filename,W,ascii)
  % WRITEDMAT  writes a matrix from a dmat file.  first line is <# columns> <#
  % rows>, then values with columns running faster
  %
  % writeDMAT(filename,W)
  % writeDMAT(filename,W,asci)
  %
  % Input:
  %   filename  name of .dmat file
  %   W  matrix read from file
  %   ascii  write ascii file {true}
  %
  % Example:
  %   % weights for bones (C,BE), size #V by #BE
  %   WB;
  %   % rearrange weights to match sorted bone edge order used by lbs.app
  %   [~,I] = sortrows(BE,[1 2]);
  %   WB = WB(:,I);
  %   % Append zeros for dummy point handle weights
  %   W = [zeros(size(V,1),size(C,1)) WB];
  %   % finally, write weights to file, usable by lbs.app
  %   writeDMAT(filename,W);
  %
  % See also: readDMAT
  %

  %disp(['writing: ',filename]);
    % open file
    fp = fopen(filename,'w');
  if ~exist('ascii','var') || ascii
    % write header
    fprintf(fp,'%d %d\n',fliplr(size(W)));
    % write data
    fprintf(fp,'%.16g\n',W);
  else
    % write header for ascii
    fprintf(fp,'0 0\n');
    % write header for binary
    fprintf(fp,'%d %d\n',fliplr(size(W)));
    if ~isa(W,'double')
      warning('Input will be cast from %s to double',class(W));
    end
    fwrite(fp,double(W),'double');
  end
  % close file
  fclose(fp);

end

author: alecjacobson

repo: gptoolbox


text:
function ZYahtzee(np)
% Play text-based Yahtzee!
%
% ZYAHTZEE
% ZYAHTZEE(n), where n is the number of players.
%
% All standard Yahtzee rules apply.
%
% On your turn, you may roll the dice up to 3 times, although you may stop
% and score after your first or second roll. The 1st roll will be done
% automatically and displayed next to the "your roll" field. The 1st roll
% of each turn is of 5 6-sided dice. For your 2nd or 3rd rolls you may keep
% any or all of the dice of your previous roll. Select which rolls to keep
% by entering their die number(s) in the prompt. (You are now no longer
% required to enter spaces when deciding which dice to hold, i.e.,
% entering 125 will hold dice numbers 1, 2, and 5.)
%
% In addition to entering numbers in the prompt, you may enter "all" to
% hold all dice, enter to reroll all dice or "quit" to exit the game.
%
% Once you have held all dice or taken 3 rolls on this turn, you select the
% category to apply the roll to, and earn points accordingly. Your score
% sheet will be updated according to standard Yahtzee rules. You will be
% allowed 13 total turns per game.
%
% Good luck, and good dice.
%
% For more information see the official rules at:
% http://www.hasbro.com/common/instruct/Yahtzee.pdf
%
% Technical Note: You may wish to randomize the seed of MATLAB's random
% number generators prior to playing to avoid games that are the same, or
% even insert a randomization line at the beginning of this code.
%
%
% %%% ZCD June 2010 %%%
%

% Edits (ZCD) July 2010:
%   * fixed spelling error in game board
%   * do not accept bad die number inputs (instead of ignoring them)
%   * added statistics "easter egg"
%   * added multi-player versus capability
%   * automatically score final roll
% Edits (ZCD) August 2010:
%   * fixed improper Yahtzee points allocations in versus mode
%   * adjusted figure location on stats display
% Edits (ZCD) September 2010:
%   * added features to stats bonus display
%   * will save scores to appropriate directory (even when not in that dir)
% Edits (ZCD) February 2011:
%   * there is now no need to include spaces when entering held die numbers

% check inputs
if nargin == 0
    np=1;
elseif ~(np>0 && real(np) && numel(np)==1)
    error('Invalid input for number of players.')
end
% display variable for multiple players
pz = [np 1:np-1];

% set up initial points array
points = ones(13,np)*NaN;
% roll variable
roll = zeros(1,5);

% main game loop runs 13 rounds, with 3 throws per round
for r = 1:13
    % run one round loop per player
    for p = 1:np
    % when playing with multiple players, briefly show the previous
    % player's updated board
    if np>1 && ~(p==1 && r==1)
        clc
        dispPoints(points(:,pz(p)),zeros(1,5),0,pz(p));
        pause(1.5)
    end

    % initialize our rerolls variable
    rerolls = 1:5;
    % set up for a bonus yatzee
    bflag = 0; strbonus = '';

    for t = 1:3
        % throw the dice
        roll(rerolls) = randi(6,[1 length(rerolls)]);
        % show the user the game state
        dispPoints(points(:,p),roll,t,p);

        % --- do special yahtzee things ---
        % check for yahtzee
        fq = hist(roll,1:6);
        if any(fq==5)
            if isnan(points(13,p))
                points(13,p) = 50;
                dispPoints(points(:,p),roll,t,p);
                % be excited about the yahtzee
                disp(' ');disp('YAHTZEE!');disp(' ');pause(3);
                break;
            elseif points(13,p)>0
                points(13,p) = points(13,p)+100;
                dispPoints(points(:,p),roll,t,p);
                disp(' ');
                disp(['BONUS YAHTZEE' repmat('!!',[1 (points(13,p)-50)/100])]);
                disp(' ');pause(3);
                bflag = 1;
                strbonus = 'BONUS ';
                action = 'all';
            elseif points(13,p)==0
                bflag = 1;
                action = 'all';
                disp(' ')
                disp('Seems like you shouldn''t have taken a zero for Yahtzees.')
                disp(' ');pause(3);
            end
        end
        % there are special restrictions when the joker scoring
        % for bonus yahtzees are in effect
        if bflag
            % you must fill the top if its available
            if isnan(points(roll(1)))
                points(roll(1),p) = scoreLogic(roll,roll(1));
                break;
            end
        end
        % --- end special yahtzee things ---

        % obtain the user action
        if t<3 && ~bflag
            keepflag = 1;
            while keepflag
                action = input('Which die numbers do you keep?   ','s');
                % separate each character by a space
                sAction = blanks(2*length(action)); % long row of blanks
                sAction(1:2:end) = action;          % interleave with choies
                % check for bad inputs
                [actstr cf] = str2num(sAction);
                if ( isempty(setdiff(actstr,1:5)) && cf ) || any(strcmp(action,{'','all','quit'}))
                    % accept this input if it is a (legit die # AND we can
                    % do a succesful number conversion on it) OR its a
                    % special word ('' all or quit).
                    keepflag = 0;
                else
                    disp('invalid die number')
                end
            end
        end
        % take user action
        switch action
            case 'all'
                rerolls = [];
            case 'quit'
                return;
            otherwise
                % obtain which die indicies we are to reroll
                % rerolls = setxor(str2num(action),1:5);
                rerolls = setxor(str2num(sAction),1:5);
        end

        % get the score placement
        if t==3 || isempty(rerolls) || bflag
            catflag = 1;        % viable category choice
            % make sure this is an appropriate category
            while catflag
                % if this is the final turn, we know the category
                if sum(isnan(points(:,p))) == 1
                    catflag = 0;
                    catg = find(isnan(points(:,p)));
                    points(catg,p) = scoreLogic(roll,catg,bflag);
                    pause(1.5);
                else
                    % ask user for category
                    response = input(['Mark ' strbonus 'score in which category?   '],'s');
                    % the only valid string here is 'quit' which exits the
                    % game now
                    if strcmpi(response,'quit')
                        return;
                    else
                        catg = str2double(response);
                    end
                    % if we are in the right category range and do not yet have
                    % a score for that category we allow the assignment
                    if catg>0 && catg<14 && isnan(points(catg,p))
                        catflag = 0;
                        % get the category score
                        points(catg,p) = scoreLogic(roll,catg,bflag);
                    else
                        disp('invalid category')
                    end
                end
            end
            break;
        end

    end % throw loop
    end % player loop
end     % round loop

% --- End game operations ---
clc
gtotal = zeros(1,np);
for p=1:np
    % give final accounting to the player
    gtotal(p) = dispPoints(points(:,p),zeros(1,5),0,p);
end
for p=1:np
    % run score records
    highScores(gtotal(p),points(13,p),[np p])
end
if np>1
    % interleave player no and score
    c = [1:np;gtotal];
    s = sprintf('Player%d: %d,  ',c(:));
    disp(s(1:end-3)); disp(' ')
    disp(['PLAYER ' num2str(find(gtotal==max(gtotal))) ' WINS' repmat('!',[1 np])])
    disp(' ')
end
% ---------------------------

    function score = scoreLogic(roll,cat,joker)
        score = 0;
        % how many points is this roll worth to this category?
        switch cat
            case {1 2 3 4 5 6}
                score = sum(roll==cat)*cat;
            case 7
                fq = hist(roll,1:6);
                if any(fq>2) || joker, score = sum(roll); end
            case 8
                fq = hist(roll,1:6);
                if any(fq>3) || joker, score = sum(roll); end
            case 9
                fq = hist(roll,1:6);
                if all(diff(sort(fq))==[0 0 0 2 1]) || joker
                    score = 25;
                end
            case 10
                if ~isempty( strfind(diff(unique(roll)),[1 1 1]) ) || joker
                    score = 30;
                end
            case 11
                if ~isempty( strfind(diff(unique(roll)),[1 1 1 1]) ) || joker
                    score = 40;
                end
            case 12
                score = sum(roll);
            case 13
                % case dealt with earlier in code
        end
    end

    function gtotal = dispPoints(points,roll,t,p)
    % display the current game state
        % clear the command window
        if t~=0 clc; end

        % displays current points to the user
        slots = {'1: ones';'2: twos';'3: threes';'4: fours';'5: fives';...
            '6: sixes';'-->bonus (35/63)';'-->Top Subtotal';'7: set';'8: quads';...
            '9: boat';'10: 4-straight';'11: 5-straight';'12: Chance';...
            '13: YAHTZEE';'-->Bottom Subtotal';'-->Grand Total'};
        % top subtotal
        top = nansum(points(1:6));
        % do we have a bonus?
        if top<63, b=0; else b=35; end
        % bottom subtotal
        bot = nansum(points(7:13));

        gtotal = top+b+bot;

        % concatenate all our points information according to the
        % description in the slots variable
        disp(' ')
        disp(['PLAYER: ' num2str(p)])
        disp([slots,num2cell( [points(1:6);b;top+b;points(7:13);bot;gtotal] )])
        disp(' ')
        if t~=0
            disp(['Roll  ' num2str(t) '/3'])
            fprintf('Your Roll:    %d   %d   %d   %d   %d\n',roll)
            fprintf('Die Number:  (1) (2) (3) (4) (5)\n')
        end
    end

    function highScores(total,yah,np)
        dirFiles = dir(strrep(mfilename('fullpath'),mfilename,''));
        if ~any(strcmp({dirFiles.name},'YahtzeeRecords.mat'))
            scores = {'Player' 'Score' 'Date' 'Yahtzees' 'Rank'};
        else
            load('YahtzeeRecords')
        end

        name = inputdlg('Enter Your Name For Record Keeping',['Player ' num2str(np(2))]);
        if isempty(name), return; end
        name = name{1};
        if length(name)>15
            name = name(1:15);
        end

        if yah~=0, yah = mod(yah,99)-49; end

        scores(end+1,:) = {name total datestr(now,1) yah ''};

        sc = scores(2:end,2);
        [B IX] = sort([sc{:}]','descend');

        scores(2:end,:) = scores(IX+1,:);
        scores(2:end,end) = num2cell(1:size(scores,1)-1);
        save([strrep(mfilename('fullpath'),mfilename,'') 'YahtzeeRecords'],'scores');

        if size(scores,1)>10
            dispLengthTop = 10;
        else
            dispLengthTop = size(scores,1);
        end
        curRank = find(B==total,1,'last');
        if curRank>=dispLengthTop
            dispLengthNow = [curRank-3:curRank+3];
            if ~isempty(find(dispLengthNow==size(scores,1)))
                dispLengthNow = dispLengthNow(1:find(dispLengthNow==size(scores,1)));
            end
        else
            dispLengthNow = [];
        end

        disp(' ');disp(' ');
        disp(scores(1:dispLengthTop,:));
        if ~isempty(dispLengthNow)
            disp('     ...............................');
            disp(scores([1 dispLengthNow],:));
        end
        disp(' ');

        % occasionally do some interesting stats
        if rand>0.85 && size(scores,1)>5 && np(1)==1
            % print stats to screen
            disp(' '); disp('Your Statistics:'); disp(' ')
            tg = size(scores,1)-1;
            fprintf('Total games played: %g, Average games per day: %g\n\n',...
                [tg, tg/(etime(datevec(max(datenum(scores(2:end,3)))),...
                datevec(min(datenum(scores(2:end,3)))))/(60*60*24))]);
            disp('Correlation between scores and Yahtzees')
            rho = corrcoef([scores{2:end,2}],[scores{2:end,4}]);
            disp(rho(2))
            disp('Correlation between scores and the date')
            rho = corrcoef([scores{2:end,2}],datenum(scores(2:end,3)));
            disp(rho(2))
            disp('Correlation between scores and rank')
            rho = corrcoef([scores{2:end,2}],[scores{2:end,end}]);
            disp(rho(2))
            disp(' ')
            disp(['Mean of scores: ' num2str(mean([scores{2:end,2}]))])
            disp(['Standard deviation of scores: ' num2str(std([scores{2:end,2}]))])
            disp(' ')
            ny = repmat(unique([scores{2:end,4}]),[2 1]);
            eval(sprintf('disp([''Mean for %d Yahtzees: '' num2str( mean( [scores{[false %d==[scores{2:end,4}]],2}] ) )]), ',ny(:)'))
            disp(' ')
            eval(sprintf('disp( [''Best score for %d Yahtzees: '' num2str(max([scores{[false scores{2:end,4}]==%d,2}]))] ), ',ny(:)'))
            disp(' ')
            disp(['Average Yahtzee''s per game: ' num2str(sum([scores{2:end,4}])/scores{end,end})]);
            disp(' ')

            % print stats to figure
            hh = figure;
            subplot(3,1,1)
            hist([scores{2:end,2}])
            xlabel('Scores'); ylabel('Frequency')
            subplot(3,1,2)
            d = datenum(scores(2:end,3));
            du = unique(d);
            x = zeros(1,length(du)); xs=x;
            for i=1:length(du)
                x(i) = mean([scores{[false; d==du(i)],2}]);
                xs(i) = std([scores{[false; d==du(i)],2}]);
            end
            if any(xs)
                [bds ibds] = max(x(xs~=0)); bdd = du(xs~=0);
                fprintf(' Best day (with more than 1 game played), %s, has mean score of %g\n',datestr(bdd(ibds)),bds)
                [bds ibds] = min(x(xs~=0)); bdd = du(xs~=0);
                fprintf('Worst day (with more than 1 game played), %s, has mean score of %g\n',datestr(bdd(ibds)),bds)
            end
            ln = pinv([ones(length(x),1) (1:length(x))'])*x';
            plot(x,'-.sb','linewidth',2), hold on
            line(1:length(x),ln(1)+(1:length(x)).*ln(2),'color','r')
            set(gca,'Xtick',[])
            xlabel('Increasing Date (days) \rightarrow'); ylabel('Mean Score')
            legend('day score','linear fit','location','best')
            subplot(3,1,3)
            plist = unique(scores(2:end,1));
            x = zeros(1,length(plist));
            for i=1:length(plist)
                tmp = [scores{boolean([0; strcmpi(plist{i},scores(2:end,1))]),2}];
                x(i) = mean(tmp);
                plist{i} = [plist{i} ': ' num2str(length(tmp))];
            end
            [xb IX] = sort(x(1,:),'descend');
            if length(plist)<=6, lp=length(plist); else lp = 6; end
            bar(xb(1:lp),'r'); xlim([0.5 lp+0.5])
            set(gca,{'xtick','fontsize','xticklabel'},{1:lp,8,plist(IX(1:lp))})
            ylabel('Mean Score'); xlabel('Player: Games Played')
            disp(' ');
        end
    end

end


text:
function vw = makeGrayROI(vw,name,select,color)
%
% vw = makeGrayROI(vw,[name],[select],[color])
%
% Makes an ROI consisting of all of the gray.coords
%
% name: name (string) for the ROI (default = 'gray')
% select: if non-zero, chooses the new ROI as the selectedROI
%         (default=1).
% color: sets color for drawing the ROI (default 'b').
%
% djh, 2/15/2001
% 09/2005 ras, if non-gray view, makes it in a hidden gray and xforms over

% needs image size and nSlices
mrGlobals;

if notDefined('name'),      name='gray'; end
if notDefined('select'),    select=1;    end
if notDefined('color'),     color=[.6 .6 .6]; end

if ~strcmp(vw.viewType,'Gray')
    % make a gray ROI in a hidden gray view,
    % and xform it over
    % (but, first check if a segmentation is installed)
    if ~exist(fullfile('Gray','coords.mat'), 'file')
        error('No Segmentation currently installed.');
    end

    dt = viewGet(vw,'curdt');
    scan = viewGet(vw,'curscan');
    hG = initHiddenGray(dt,scan);
    hG = makeGrayROI(hG,name,select,color);
    switch vw.viewType
        case 'Inplane', vw = vol2ipAllROIs(hG,vw);
        case 'Flat', vw = vol2flatAllROIs(hG,vw);
        case 'Volume', vw = addROI(vw,hG.ROIs);
        otherwise, error('Unkown view type.');
    end

    return
end

ROI.name=name;
ROI.viewType=vw.viewType;
ROI.coords=vw.coords;
ROI.color=color;

vw = addROI(vw,ROI,select);

return;

author: vistalab

repo: vistasoft


text:
function result = analyzeLocalizationQuality(debug, ground_truth, PerformanceSpec, visualize)

% Author: Alexander Wischnewski     Date: 24-02-2019
%
% Description:
%   scans the debug data and ground truth and evalutes the localization
%   quality
%
% Input:
%   debug:                  Debug structure of the vehicle
%   ground_truth:           Physical ground truth data from simulation
%   PerformanceSpec:        Structure with fields:
%       sDevMax_m               Maximum deviation position in longitudinal direction
%       dDevMax_m               Maximum deviation position in lateral direction
%       psiDevMax_rad           Maximum deviation heading
%       vxDevMax_mps            Maximum deviation longitudinal velocity
%       vyDevMax_mps            Maximum deviation lateral velocity
%       dPsiDevMax_radps        Maximum deviation yaw rate
%   visualize:              Set to true if debug plots shall be done
%
% Outputs:
%   result:       1 if the localization quality is sufficient

result = 1;

disp('------------------------------');
disp('Check Localization Quality: ');
disp('------------------------------');
% check whether sensor fusion is accurate enough during the whole dataset
% after it was valid for the first time
idxStart = find(uint16(debug.debug_mvdc_state_estimation_debug_StateEstimate_SE_Status.Data) == 2, 1);
if(isempty(idxStart))
    disp('No valid state estimate was found in the dataset at all');
    result = 0;
    return
end
tStart = debug.debug_mvdc_state_estimation_debug_StateEstimate_SE_Status.Time(idxStart);
disp(['Localization valid after ' num2str(tStart) 's']);

% localization quality has to be assessed dependent on the current vehicle
% orientation, as the quality requirements in lateral are way higher than
% in longitudinal

t_vec = tStart:0.04:ground_truth.SimRealState_x_m.Time(end);
x_m_gt = resample(ground_truth.SimRealState_x_m, t_vec);
y_m_gt = resample(ground_truth.SimRealState_y_m, t_vec);
x_m_se = resample(debug.debug_mvdc_state_estimation_debug_StateEstimate_Pos_x_m, t_vec);
y_m_se = resample(debug.debug_mvdc_state_estimation_debug_StateEstimate_Pos_y_m, t_vec);

if(visualize)
    figure; hold on; grid on;
    plot(x_m_gt.Data, y_m_gt.Data);
    plot(x_m_se.Data, y_m_se.Data);
    axis equal;
    legend('Ground Truth', 'State Estimation');
end

psi_rad_gt = -1*(resample(ground_truth.SimRealState_psi_rad, t_vec) + pi/2);
dx = x_m_se.Data - x_m_gt.Data;
dy = y_m_se.Data - y_m_gt.Data;
p_s = zeros(length(dx), 1);
p_d = zeros(length(dx), 1);
% transform localization deviation to vehicle coordinate system
for i = 1:1:length(dx)
    p_s(i) = cos(psi_rad_gt.Data(i))*dx(i) - sin(psi_rad_gt.Data(i))*dy(i);
    p_d(i) = sin(psi_rad_gt.Data(i))*dx(i) + cos(psi_rad_gt.Data(i))*dy(i);
end
s_m = timeseries(p_s, t_vec);
d_m = timeseries(p_d, t_vec);
zero_ts = timeseries(zeros(length(dx), 1), t_vec);

if(compareSignals(zero_ts, s_m,...
        tStart, 'Longitudinal localization error', PerformanceSpec.sDevMax_m, false, visualize) == 0)
    result = 0;
end

if(compareSignals(zero_ts, d_m,...
        tStart, 'Lateral localization error', PerformanceSpec.dDevMax_m, false, visualize) == 0)
    result = 0;
end

if(compareSignals(ground_truth.SimRealState_psi_rad, debug.debug_mvdc_state_estimation_debug_StateEstimate_Pos_psi_rad,...
        tStart, 'Heading', PerformanceSpec.psiDevMax_rad, true, visualize) == 0)
    result = 0;
end

if(compareSignals(ground_truth.SimRealState_vx_mps, debug.debug_mvdc_state_estimation_debug_StateEstimate_vx_mps,...
        tStart, 'Longitudinal velocity', PerformanceSpec.vxDevMax_mps, false, visualize) == 0)
    result = 0;
end

if(compareSignals(ground_truth.SimRealState_vy_mps, debug.debug_mvdc_state_estimation_debug_StateEstimate_vy_mps,...
        tStart, 'Lateral velocity', PerformanceSpec.vyDevMax_mps, false, visualize) == 0)
    result = 0;
end

if(result == 1)
   disp(' ');
   disp('All localization checks passed.');
end

end

function dev_ok = compareSignals(signal1, signal2, tStart, name, maxDev, normalize, visualize)
    s1 = resample(signal1, tStart:0.004:signal1.Time(end));
    s2 = resample(signal2, tStart:0.004:signal1.Time(end));
    s1Dev = s1.Data-s2.Data;
    if(normalize)
        s1Dev = normalizeAngle(s1Dev);
    end
    s1Dev = abs(s1Dev);
    disp(' ');
    disp(name);
    display(['Maximum deviation: ' num2str(max(s1Dev))]);
    if(~(max(s1Dev) < maxDev))
        disp('Deviation too high');
        dev_ok = 0;
    else
        disp('Deviation within spec');
        dev_ok = 1;
    end
    if(visualize)
        figure;
        subplot(2, 1, 1);
        if(normalize)
            plot(signal1.Time, normalizeAngle(signal1.Data));
        else
            plot(signal1);
        end
        hold on; grid on;
        if(normalize)
            plot(signal2.Time, normalizeAngle(signal2.Data));
        else
            plot(signal2);
        end
        title(name);
        subplot(2, 1, 2);
        plot(s1.Time, s1Dev);
        grid on; hold on;
        plot([s1.Time(1), s1.Time(end)], [maxDev, maxDev], '--')
        ylabel('Deviation');
        xlabel('Time in s');
    end
end

author: TUMFTM


text:
function [C,alpha_C] = imdarken(A,alpha_A,B,alpha_B)
  % IMDARKEN Blend A over B, replacing pixels in B if the cooresponding pixel in
  % A is darker.
  %
  % Inputs:
  %   A  w by h by c image A
  %   alpha_A  alpha mask for image A, either:
  %     w by h by c
  %     w by h  (same across channels)
  %     s  (same for every pixel)
  %   B  w by h by c image B
  %   alpha_B  alpha mask for image B, either:
  %     w by h by c
  %     w by h  (same across channels)
  %     s  (same for every pixel)
  % Outputs:
  %   C  w by h by c result image, *ALWAYS of type double*
  %   alpha_C  alpha mask for resulting image

  % https://developer.android.com/reference/android/graphics/BlendMode#DARKEN
  alpha_C = alpha_A+alpha_B-alpha_B.*alpha_A;
  C = (1-alpha_B).*A + (1-alpha_A).*B+min(B,A);
end


text:
function b = r8plu_mul ( n, pivot, lu, x )

%*****************************************************************************80
%
%% R8PLU_MUL computes A * x using the PLU factors of A.
%
%  Discussion:
%
%    It is assumed that R8MAT_TO_R8PLU has computed the PLU factors of
%    the matrix A.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer PIVOT(N), the pivot vector computed by R8MAT_TO_R8PLU.
%
%    Input, real LU(N,N), the matrix factors computed by
%    R8MAT_TO_R8PLU.
%
%    Input, real X(N), the vector to be multiplied.
%
%    Output, real B(N), the result of the multiplication.
%
  b(1:n) = x(1:n);
%
%  Y = U * X.
%
  for j = 1 : n
    b(1:j-1) = b(1:j-1) + lu(1:j-1,j)' * b(j);
    b(j) = lu(j,j) * b(j);
  end
%
%  B = PL * Y = PL * U * X = A * x.
%
  for j = n-1 : -1 : 1

    b(j+1:n) = b(j+1:n) - lu(j+1:n,j)' * b(j);

    k = pivot(j);

    if ( k ~= j )
      temp = b(k);
      b(k) = b(j);
      b(j) = temp;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
addpath(genpath('../mutils/My/'));
addpath(genpath('../ptv'));
%%
data = load('data/SA_ppu.mat');
%%
imgs2 = abs(double(squeeze(data.data)));
imgs = imgs2(60:140,80:170, 40:end);
imgs = reshape(imgs, [size(imgs, 1), size(imgs,2), 1,1, size(imgs, 3)]);
imgs(imgs> 1000) = 1000;

opts = [];
opts.pix_resolution = [1, 1];
opts.interp_type = 0;

maxpf= max(max(imgs, [], 1), [],  2);
imgs = bsxfun(@rdivide, imgs, maxpf);

opts.border_mask = 5;

opts.singular_coefs = ones(1, size(imgs, 5));
opts.singular_coefs(1) = 0;

opts.T_D2L1 = 1e-2;
opts.D2L1 = 1e-2;
opts.isoTV = 0.000;

opts.max_iters = 150;
opts.grid_spacing = [4, 4];
opts.cp_refinements = 0;

opts.display = 'iter';
opts.metric = 'nuclear';
opts.k_down = 0.75;
opts.check_gradients = 0;
opts.spline_order = 3;
opts.display = 'off';

tic
[voldef, Tmin_out] = ptv_register(imgs, [], opts);
toc
szs = size(imgs);
%% we can mad displacements back to the first reference frame
refid = 1;
[T_1, imgs_1] = remap_displacements(Tmin_out, refid, imgs, opts.pix_resolution);

%%
imgs_exp = [];
for i = 1 : size(imgs, 5)
    hold off;
    imagesc(squeeze(imgs(:,:, 1, :, i)));
    colormap gray;
    hold on;
    plot_wiremesh(squeeze(T_1(:,:, 1, :, i)), 4, 'r-');
    daspect([1,1,1]);
    F = getframe;
    imgs_exp = cat(4, imgs_exp, double(F.cdata));
%     pause(0.1);
end
%%
savegif('results/heart_SA_register.gif', squeeze([imgs, voldef]), 1/20);
savegif('results/heart_SA_track.gif', squeeze([imgs_exp]/250), 1/20);

repo: pTVreg


text:
function blas1_s_test02 ( )

%*****************************************************************************80
%
%% TEST02 tests ISAMAX, SAXPY, SDOT, SSCAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2006
%
%  Author:
%
%    John Burkardt
%
  n = 10;
  lda = n;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  Use ISAMAX, SAXPY, SDOT and SSCAL\n' );
  fprintf ( 1, '  in a Gauss elimination routine.\n' );
  fprintf ( 1, '\n' );
%
%  Set the matrix.
%
  for i = 1 : n
    for j = 1 : n

      if ( i == j )
        a(i,j) = 2.0;
      elseif ( i == j + 1 )
        a(i,j) = - 1.0;
      elseif ( i == j - 1 )
        a(i,j) = - 1.0;
      else
        a(i,j) = 0.0;
      end

    end
  end
%
%  Set the right hand side.
%
  b(1:n-1) = 0.0;
  b(n) = n + 1.0;

  info = 0;

  for k = 1 : n - 1

    l = isamax ( n-k+1, a(k:n,k), 1 ) + k - 1;
    ipvt(k) = l;

    if ( a(l,k) == 0.0 )

      info = k;

    else

      if ( l ~= k )
        t = a(l,k);
        a(l,k) = a(k,k);
        a(k,k) = t;
      end

      t = -1.0 / a(k,k);
      a(k+1:n,k) = sscal ( n-k, t, a(k+1:n,k), 1 );

      for j = k+1 : n
        t = a(l,j);
        if ( l ~= k )
          a(l,j) = a(k,j);
          a(k,j) = t;
        end
        a(k+1:n,j) = saxpy ( n-k, t, a(k+1:n,k), 1, a(k+1:n,j), 1 );
      end

    end

  end

  ipvt(n) = n;
  if ( a(n,n) == 0.0 )
    info = n;
  end

  if ( info ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The matrix is singular.\n' );
    return
  end

  for k = 1 : n-1
    l = ipvt(k);
    t = b(l);
    if ( l ~= k )
      b(l) = b(k);
      b(k) = t;
    end
    b(k+1:n) = saxpy ( n-k, t, a(k+1:n,k)', 1, b(k+1:n), 1 );
  end

  for k = n : -1 : 1
    b(k) = b(k) / a(k,k);
    t = - b(k);
    b(1:k-1) = saxpy ( k-1, t, a(1:k-1,k)', 1, b(1:k-1), 1 );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Solution (should be (1,2,3,...,n)):\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %12f\n', b(i) );
  end

  return
end


text:
function pts =densify(p1,spacing)
np=length(p1);
p2=zeros(np,2); padd=zeros(floor(np/3),2);
k=0;
for i = 1:np-1
    if isnan(p1(i,:))
        k=k+1;
        p2(k,1:2)=NaN;
        continue;
    end
    if isnan(p1(i+1,:))
        k=k+1;
        p2(k,1:2)=NaN;
        continue
    end;
    d = sqrt((p1(i,1)-p1(i+1,1)).^2 + (p1(i,2)-p1(i+1,2)).^2);
    nsplit = ceil(d/spacing);
    k=k+1;
    p2(k,:)=p1(i,:);
    padd(1:nsplit,1)=(linspace(p1(i,1),p1(i+1,1),nsplit))';
    padd(1:nsplit,2)=(linspace(p1(i,2),p1(i+1,2),nsplit))';
    k=k+nsplit;
    p2(k-nsplit+1:k,:)=padd(1:nsplit,:);
    k=k+1;
    p2(k,:)=p1(i+1,:);
    padd=padd*0;
end
x=p2(1:k,1); y=p2(1:k,2);
pts = [x,y];
end

author: CHLNDDEV


text:
%% Policy evaluation
%
% V(s) = \sum_a \pi(s,a) \sum_s' P_ss'^a [R_ss'^a + \gamma V(s')]
%

clear all;

gw_options = [];
gw_options.dims = [15,15];
grid_world = Grid_world(gw_options);
actions    = [1 0; -1 0;  0 1;  0 -1];

%% Create reward

r_options               = [];
r_options.x_goal        = [10,10];
r_options.x_bad         = rl_reward_firepit( 'central_pit' );

r_options.type           = 'discrete';
r_options.r_goal         =   10; % reward for reaching the goal.
r_options.r_pit          =   -10; % reward for falling into a pit.
r_options.r_step         =     -1; % reward for each time step.

r_options.grid_world    = grid_world;

reward                  = Reward_gridworld(r_options);

% Discrete reward function

reward_f                = @(s,x)reward_grid_world_goal(s,x,reward.s_goal,reward.s_bad,reward.r_goal,reward.r_pit,reward.r_step);

%% Create Transition and Reward Model

p_noise   = 0.0;

[P,N,R]   = rl_create_PR_models(grid_world,actions,reward_f,p_noise);

%% Value Iteration
bRecord    = true;
V          = zeros(grid_world.num_states,1);
gamma      = 0.9;
stop_c     = 0.001;

Q          = rand(grid_world.num_states,4);
tau        = 10;
policy     = @(s)rl_softmax_policy(s,Q,tau);

num_iterations = 1;
dtau1          = @(i) -(tau)/(num_iterations+1) * i + tau;
s              = tau;
dtau2          = @(i) 10 * exp(- 0.1 * i*i) + 0.05;
dtau2(1)
dtau2(2)

% Policy Evaluation -> Policy Improvement -> Policy Evaluation -> .... -> until convergence

data_save = [];
data_save.V   = cell(num_iterations,1);
data_save.Q   = cell(num_iterations,1);
data_save.tau = [];
%%
close all;
hf = figure('Position',[195         359        1143         419]);

subplot(1,2,1);
[~,idx] = max(Q,[],2);
U       = actions(idx,:);hold on;
h_P     = rl_2D_plot_policy_uw(0.5 .* U,grid_world.X(:),grid_world.Y(:),[]);
reward.plot_firpit();
reward.plot_goal_state();
xlim([-1 11]);
ylim([-1 11]);
box on;

subplot(1,2,2);
h_Q = plot_qtable(gca,V,grid_world.X,grid_world.Y);

title(['Policy Evaluation # ' num2str(1) ]);
pause(0.1);

for i=1:num_iterations

    disp(['Policy Evaluation (' num2str(i) ')']);
    [V,Q,data] = rl_policy_evaluation(policy,V,actions,P,reward_f,N,gamma,stop_c,bRecord);

    disp('Policy Improvement');
    tau        = dtau2(i);
    policy     = @(s)rl_softmax_policy(s,Q,tau);

    data_save.V{i} = V;
    data_save.Q{i} = Q;
    data_save.tau  = [data_save.tau,tau];

    % Plot Update
    [~,idx] = max(Q,[],2);
    U       = actions(idx,:);hold on;
    rl_2D_plot_policy_uw(0.5 .* U,grid_world.X(:),grid_world.Y(:),h_P);
    plot_qtable(gca,V,grid_world.X,grid_world.Y,h_Q);
    title(['Policy Evaluation # ' num2str(i) ]);
    pause(0.1);

end

%% Play Sequce of Evaulation

bSaveVideo = true;

if bSaveVideo
   writer = VideoWriter('PE_multiple.avi');
   open(writer);
end

V =  data_save.V{1};
Q =  data_save.Q{1};
close all;
hf = figure('Position',[195         359        1143         419]);

subplot(1,2,1);
[~,idx] = max(Q,[],2);
U       = actions(idx,:);hold on;
h_P     = rl_2D_plot_policy_uw(0.5 .* U,grid_world.X(:),grid_world.Y(:),[]);
reward.plot_firpit();
reward.plot_goal_state();
xlim([-1 11]);
ylim([-1 11]);
box on;

subplot(1,2,2);
h_Q = plot_qtable(gca,V,grid_world.X,grid_world.Y);
title(['Policy Evaluation # ' num2str(1) ]);

bRun     = true;
i        = 1;
count    = 0;
num_iter = 20;% size(data_save.V,1);

while(bRun)

    V =  data_save.V{i};
    Q =  data_save.Q{i};

    % Plot Policy
    [~,idx] = max(Q,[],2);
    U       = actions(idx,:);hold on;
    rl_2D_plot_policy_uw(0.5 .* U,grid_world.X(:),grid_world.Y(:),h_P);

    % Plot Value function
    plot_qtable(gca,V,grid_world.X,grid_world.Y,h_Q);
    title(['Policy Evaluation # ' num2str(i) ]);

    if bSaveVideo == true
        F = getframe(hf);
        writeVideo(writer,F);
        pause(0.01);
        count = count + 1;
        if(count == 10)
            count = 0;
            i = i + 1;
        end
    else
        i = i + 1;
        pause(1);
    end

    if i == num_iter
        bRun = false;
    end

end
disp('Finished');

if bSaveVideo
 close(writer);
end

%% Record video for Policy Evaluation Iterations (single)

bSaveVideo = true;

if bSaveVideo
   writer = VideoWriter('PE_single.avi');
   open(writer);
end

close all;
hf = figure;
set(gcf,'color','w');
h_Q = plot_qtable(gca,data.Vs(:,1),grid_world.X,grid_world.Y);
title(['Policy Evaluation #1 iter(' num2str(1) ')']);

bRun     = true;
i        = 1;
count    = 0;
num_iter = size(data.Vs,1);

while(bRun)

    % Plot Policy

    % Plot Value function
    plot_qtable(gca,data.Vs(:,i),grid_world.X,grid_world.Y,h_Q);
    title(['Policy Evaluation #1 iter(' num2str(i) ')']);

    if bSaveVideo == true
        F = getframe(hf);
        writeVideo(writer,F);
        pause(0.01);
        count = count + 1;
        if(count == 10)
            count = 0;
            i = i + 1;
        end
    else
        i = i + 1;
        pause(1);
    end

    if i == num_iter
        bRun = false;
    end

end
disp('Finished');

if bSaveVideo
 close(writer);
end

%%

author: epfl-lasa


text:
function [cum_ret, cumprod_ret, daily_ret, daily_portfolio] = ...
    cwmr_var_run(fid, data, phi, epsilon, tc, opts)
% This program simulates the CWMR-Var algorithm
%
% function [cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ...
%    = cwmr_var_run(fid, data, phi, epsilon, tc, opts)
%
% cum_ret: a number representing the final cumulative wealth.
% cumprod_ret: cumulative return until each trading period
% daily_ret: individual returns for each trading period
% daily_portfolio: individual portfolio for each trading period
%
% data: market sequence vectors
% fid: handle for write log file
% phi: confidence parameter
% epsilon: mean reversion threshold
% tc: transaction cost rate parameter
% opts: option parameter for behvaioral control
%
% Example: [cum_ret, cumprod_ret, daily_ret, daily_portfolio, exp_ret] ...
%            = cwmr_var_run(fid, data, phi, epsilon, tc, opts)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n, m]=size(data);

% Variables for return, start with uniform weight
cum_ret = 1;
cumprod_ret = ones(n, 1);
daily_ret = ones(n, 1);
day_weight = ones(m, 1)/m;  %#ok<*NASGU>
day_weight_o = zeros(m, 1);
daily_portfolio = zeros(n, m);

% print file head
fprintf(fid, '-------------------------------------\n');
fprintf(fid, 'Parameters [phi:%f, epsilon:%f, tc:%f]\n', phi, epsilon, tc);
fprintf(fid, 'day\t Daily Return\t Total return\n');

fprintf(1, '-------------------------------------\n');
if(~opts.quiet_mode)
    fprintf(1, 'Parameters [tc:%f]\n', tc);
    fprintf(1, 'day\t Daily Return\t Total return\n');
end

%% Initialization
alpha = 1; lambda = 0;
mu = alpha*ones(m, 1)/m;
sigma = alpha* eye(m)/(m^2);

%% Online
if (opts.progress)
	progress = waitbar(0,'Executing Algorithm...');
end
for t = 1:1:n,
    % Calculate t's portfolio at the beginning of t-th trading day
    if (t >= 2)
         [day_weight, mu, sigma] ...
             = cwmr_var_kernel(data(1:t-1, :), mu, sigma, phi, epsilon);
    end

    day_weight = day_weight./sum(day_weight);
    daily_portfolio(t, :) = day_weight';

    % Step 3: Cal t's return and total return
    daily_ret(t, 1) = (data(t, :)*day_weight)*(1-tc/2*sum(abs(day_weight-day_weight_o)));
    cum_ret = cum_ret * daily_ret(t, 1);
    cumprod_ret(t, 1) = cum_ret;

    % Adjust weight(t, :) for the transaction cost issue
    day_weight_o = day_weight.*data(t, :)'/daily_ret(t, 1);

    % Time consuming part, other way?
    fprintf(fid, '%d\t%f\t%f\n', t, daily_ret(t, 1), cumprod_ret(t, 1));
    if (~opts.quiet_mode)
        if (~mod(t, opts.display_interval)),
            fprintf(1, '%d\t%f\t%f\n', t, daily_ret(t, 1), cumprod_ret(t, 1));
        end
    end
    if (opts.progress)
		if mod(t, 50) == 0
			waitbar((t/n));
		end
	end
end

% Debug Information
if(opts.log_mode)
    fprintf(fid, 'CWMR-Var(%.2f, %.2f, %.4f), Final return: %.2f\n', ...
        phi, epsilon, tc, cum_ret);
    fprintf(fid, '-------------------------------------\n');
end

fprintf(1, 'CWMR(%.2f, %.2f, %.4f), Final return: %.2f\n', ...
    phi, epsilon, tc, cum_ret);
fprintf(1, '-------------------------------------\n');

	if (opts.progress)
		close(progress);
	end

end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS