text:
function out = correl_compare_dep(y1,y2,varargin)
% Compare dependent correlations between pairs of vectors in y1 and y2.
%
% :Usage:
% ::
%
%     out = correl_compare_dep(y1,y2,['alpha',myalpha],['rank'],['table'])
%
% Each of y1 and y2 would contain at least 2 columns, which would be
% correlated and saved in r1 and r2 matrices in output.
% Then, the r1 and r2 matrices are subtracted, and P-values are
% returned for the differences.
%   - In the simplest case, y1 would contain vectors [a b] and y2 would
%     contain vectors [a c].  tests are provided on the a-b vs. a-c
%     difference in correlations.
%
% Repeats dep. correl. analysis for each pair of columns
%
% Returns results in correlation matrix form, where number of rows and
% cols. are the number of pairs [y1(:,i) y2(:,i)]
%
% myalpha is 2-tailed alpha value; p-values are 2-tailed
% FDR correction is at .01, 2-tailed
%
% Based on Steiger, 1980, tests for comparing dependent correlations.
%
% :Examples:
% ::
%
%    for i = 1:length(cl), y1(:,i) = cl.CONTRAST.data(:,2); y2(:,i) = cl.CONTRAST.data(:,1); end
%    for i = 1:length(cl), y1(:,i) = cl(i).CONTRAST.data(:,2); y2(:,i) = cl(i).CONTRAST.data(:,1); end
%
%    % y1 is matrix of obs x data vectors for condition 1
%    % y2 is matrix of obs x data vectors for condition 2
%    out = correl_compare_dep(y1,y2)
%
%    figure('Color','w');nmdsfig(c.GroupSpace,c.ClusterSolution.classes, ...
%    c.names,out.sig,1,{'Pos' 'Neg'});
%    nmdsfig_legend(c.ClusterSolution.X,c.r)
%
% :Examples:
% ::
%
%    c_compare = correl_compare_dep(y1avg,y2avg,'alpha',.05,'rank','table','names',c.APPLY_CLUSTER.names);
%
%    out = correl_compare_dep([ypred pain],[ypred temp], 'alpha', .06, 'table', 'names', {'biomarker resp' 'pain or temp'});

    myalpha = .05;
    dorankdata = 0;
    dotable = 0;
    names = [];

    for i = 1:length(varargin)
        if isstr(varargin{i})
            switch varargin{i}
                % functional commands
                case 'alpha', myalpha = varargin{i+1};
                case 'rank', dorankdata = 1;
                case 'table', dotable = 1;
                case 'names', names = varargin{i+1};
                otherwise, warning(['Unknown input string option:' varargin{i}]);
            end
        end
    end

    [N,npairs] = size(y1);

    [rows,cols,ncorr] = corrcoef_indices(npairs);

    if nargin < 4, dorankdata = 0; end
    if dorankdata
        str = sprintf('Ranking data: Nonparametric correlations'); fprintf(1,str);

        for i = 1:npairs
            y1(:,i) = rankdata(y1(:,i));
            y2(:,i) = rankdata(y2(:,i));
        end
    else
        str = sprintf('Assuming continuous data (no ranks).'); fprintf(1,str);
    end

    erase_string(str);

    i = 1; j = 2; k = 3; h = 4;

    str = sprintf('Computing differences among correlations %04d',0); fprintf(1,str);

    diffr = zeros(ncorr,1);
    Zstar2 = zeros(ncorr,1);
    rr = zeros(ncorr,2);

    for cc = 1:ncorr

        fprintf(1,'\b\b\b\b%04d',cc);

        % get full correlation matrix for the 4 variables involved
        dat = [y1(:,[rows(cc) cols(cc)]) y2(:,[rows(cc) cols(cc)])];

        % get differences between correlation z-values (estimate)
        [diffr(cc),r,rr(cc,:),diffrz,z] = correlation_diffs(dat,i,j,k,h);

        s = covcorr(r,i,j,k,h);         % covariance of corr coeffs

        Zstar2(cc) = diffrz * sqrt( (N-3) ./ (2-(2*s)) );

        % bootstrap
        %vals = bootstrp(5000,@correlation_diffs,dat,i,j,k,h);
        %Zboot(cc) = mean(vals) ./ std(vals);

        %if Zstar2(cc) > 2, keyboard, end
    end

    erase_string(str);

    pvec = 2 * (1 - normcdf(abs(Zstar2)));

    r1 = reconstruct(rr(:,1),npairs,ncorr,rows,cols);
    r2 = reconstruct(rr(:,2),npairs,ncorr,rows,cols);

    Z = reconstruct(Zstar2,npairs,ncorr,rows,cols);
    p = reconstruct(pvec,npairs,ncorr,rows,cols);

    dat = [rr diffr Zstar2 pvec];
    dat = [rows cols dat];
    dat = dat(pvec <= myalpha,:);

    diffr = reconstruct(diffr,npairs,ncorr,rows,cols);

    sig = (p <= myalpha - eye(size(p))) .* sign(diffr);

    % FDR corrected
    pthr = FDR(pvec,.05);
    if isempty(pthr), pthr = 0; end

    sigfdr = (p <= pthr) .* sign(diffr);

    out = struct('alpha',myalpha,'r1',r1,'r2',r2,'diffr',diffr, ...
        'Z',Z,'p',p,'sig',sig,'pthr',pthr,'sigfdr',sigfdr,'sigstats',dat);

    % output table...
    if dotable
        if isempty(names)
            disp(['No names entered; try ''names'' keyword to add them.']);
            for i = 1:npairs, names{i} = ['R' num2str(i)]; end
        end

        out.names = names;

        disp(['Uncorrected, p < ' num2str(myalpha)])
        maketable(out,'sig',names);
        disp('')

        disp(['FDR corrected, p < ' num2str(myalpha)])
        maketable(out,'sigfdr',names);
        disp('')
    end

    return

function [diffr,r,rr,diffrz,z] = correlation_diffs(dat,i,j,k,h)
    % get differences between correlation z-values
    r = corrcoef(dat);

    rr = [r(i,j) r(k,h)];                % correls to be compared
    diffr = -diff(rr);              % negative sign means we get (1) - (2)

    if nargout > 3
        z = .5 .* log( (1+rr) ./ (1-rr) );
        diffrz = -diff(z);               % diff btwn correl z values
    end

    return

function [rows,cols,ncorr] = corrcoef_indices(npairs)
    % upper triangle only
    tmp = triu(ones(npairs));
    tmp = tmp - eye(npairs);
    [rows,cols] = find(tmp);
    ncorr = length(rows);
    return

function s = covcorr(R,i,j,k,h)

    % pool correl. coeff for more stable var est.
    % as they are equal under Ho.  Steiger, 1980, eq. 14 for Z*
    pooledr = mean([R(i,j) R(k,h)]);
    s = pearsonf(R,i,j,k,h) ./  ( corvar(pooledr) );

    return

function cv = corvar(r)
    % variance of a correlation coefficient
    % simplified form (special case) of pearsonf
    cv = (1 - r^2) ^ 2;
    return

    % % function cv = co(R,x,y)
    % % % x and y are 2-vector indices of matrix R to compute covariance for
    % % % cv is covariance
    % % cv = pearsonf(R,x(1),x(2),y(1),y(2));
    % % return

function pf = pearsonf(R,i,j,k,h)

    % Pearson-Filon: covariance (or var) of element i,j with element k,h
    % depending on correlation values

    % for variance of 1 correl, works as well, and reduces to:
    % (1 - rr(1)^2)^2

    pf = (1/2) .* R(i,j) .* R(k,h) .* ...
        ( R(i,k).^2 + R(i,h).^2 + R(j,k).^2 + R(j,h).^2 ) + ...
        R(i,k) .* R(j,h) + R(i,h) .* R(j,k) - ...
        R(i,j) .* ( R(j,k) .* R(j,h) + R(i,k) .* R(i,h) ) - ...
        R(k,h) .* ( R(j,k) .* R(i,k) + R(j,h) .* R(i,h) );

    return

    % function pf = pearsonf2(r,j,k,h,m)
    %
    % pf = .5 * ( (r(j,h) - r(j,k)*r(k,h)) * (r(k,m) - r(k,h)*r(h,m))  + ...
    %     (r(j,m) - r(j,h)*r(h,m)) * (r(k,h) - r(k,j)*r(j,h)) + ...
    %     (r(j,h) - r(j,m)*r(m,h)) * (r(k,m) - r(k,j)*r(j,m)) + ...
    %     (r(j,m) - r(j,k)*r(k,m)) * (r(k,h) - r(k,m)*r(m,h))    );
    %
    % return

function valmat = reconstruct(vals,npairs,ncorr,rows,cols)

    valmat = zeros(npairs);
    for i = 1:ncorr
        valmat(rows(i),cols(i)) = vals(i);
    end
    valmat = valmat + valmat';

    return

function erase_string(str1)
    fprintf(1,repmat('\b',1,length(str1))); % erase string
    return

function maketable(c_compare,whfield,names)

    str = {'-' '+'};
    [rows,cols] = find(triu(c_compare.(whfield)));
    if isempty(rows)
        disp('No significant results.')
    else
        fprintf(1,'Name1\tName2\trow\tcol.\t+ or -\tr1\tr2\tZ\tp\n');
        for i = 1:length(rows)
            fprintf(1,'%s\t%s\t%3.0f\t%3.0f\t%s\t%3.3f\t%3.3f\t%3.2f\t%3.4f\n', ...
                names{rows(i)}, names{cols(i)},rows(i),cols(i), ...
                str{1.5 + (.5.*c_compare.(whfield)(rows(i),cols(i)))}, ...          % + or - sign
                c_compare.r1(rows(i),cols(i)),c_compare.r2(rows(i),cols(i)), ...   % correlations
                c_compare.Z(rows(i),cols(i)),c_compare.p(rows(i),cols(i)));        % Z and p
        end
    end
    fprintf(1,'\n');
    return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [X,m]=getscreen(varargin)
%GETSCREEN replacement for GETFRAME for systems with multiple monitors. Use it just like you use GETFRAME.
%Can also be used to take screenshots of non-MATLAB windows. Requires JAVA.
%
%Calling options:
% Output X is m*n*3 cdata array. Output m is corresponding colormap (if applicable.)
% [X,m]=getscreen         Gets screenshot of current figure
% [X,m]=getscreen(h)      Gets screenshot of axes or figure with handle h
% [X,m]=getscreen(h,rect) Gets screenshot of at position rect relative to axes or figure with handle h.
% [X,m]=getscreen(rect) Gets screenshot of whatever is on your screen at position vector rect. (Relative to bottom left corner.)
%   rect is a 4 element position vector in pixel units, referenced from bottom left as usual for MATLAB.
%
% By Matt Daskilewicz <mattdaskil@gatech.edu> 11/2008
% www.asdl.gatech.edu
%

%process inputs:
if nargin==0
    h=gcf;
    rect=get(h,'position');
    rect(1:2)=[0 0];

elseif nargin==1
    if numel(varargin{1})==1 %it should be a handle
        h=varargin{1};
        if ~strcmp(get(h,'type'),'axes') && ~strcmp(get(h,'type'),'figure')
            error('First input must be a handle to an axes or figure, or a position vector');
        end

        rect = hgconvertunits(h, get(h,'position'), get(h,'units'), 'pixels', get(h,'parent')); %x/y coords of fighandle in pixels
        rect(1:2)=[1 1];

    elseif numel(varargin{1})==4 %it is a position vector
        h=[];
        rect=varargin{1};
    else % isnt a handle or a position vector
        error('First input must be a handle to an axes or figure, or a position vector');
    end

elseif nargin==2 %should be a handle and a position vector
    h=varargin{1};
    if ~ismember(get(h,'type'),{'figure';'axes'})
        error('First input must be a handle to an axes or figure, or a position vector');
    end

    rect=varargin{2};
    if length(rect)~=4
        error('Second input must be a 4 element position vector.')
    end

else
    error('Too many input arguments')
end

%get monitor resolutions... because java references coords from top of screen, need to know how big
%screens is.
monitorz=get(0,'monitorpositions');
maxheight=max(monitorz(:,4)); %max vertical resolution of all monitors

if ~isempty(h)

    %h is the handle to figure or axes whose origin (first 2 elemetns of position vector) we're taking
    %rect with respect to.
    origin=getpixelposition(h,true);

    %if axes are in a uipanel, nudge rect slightly to account for panel border.
    if strcmp(get(get(h,'parent'),'type'),'uipanel')
        origin=origin+[1 1 0 0];
    end

    %also need position of the figure containing h, since origin is w.r.t this figure:
    fighandle = ancestor(h,'figure'); %handle of parent figure to h

    figorigin = hgconvertunits(fighandle, get(fighandle,'position'), get(fighandle,'units'), 'pixels', 0); %x/y coords of fighandle in pixels

    if strcmp(get(h,'type'),'figure')
        origin=[1 1 0 0]; %avoid double counting figorigin and origin if fighandle==h
    end

    %% calculate dimensions of rectangle to take screenshot of: (in java coordinates)
    %java coordinates start at top left instead of bottom left, and start at [0,0]. Also "top" is defined as the top of your
    %largest monitor, even if that's not the one you're taking a screenshot from.

    x=figorigin(1)+origin(1)+rect(1)-3;
    y=maxheight-figorigin(2)-origin(2)-rect(2)-rect(4)+3;
    w=floor(rect(3));
    h=floor(rect(4));

    figure(fighandle); %make sure figure is visible on screen
    drawnow;

else %there is no figure... just take a screenshot of whatever's on the screen at rect

    rect=floor(rect);

    x=rect(1);
    y=maxheight-rect(2)-rect(4); %java coordinates start at top of screen.
    w=rect(3);
    h=rect(4);
end

%% do java:
robo = java.awt.Robot;
target=java.awt.Rectangle(x,y,w,h);
image = robo.createScreenCapture(target); %take screenshot at rect
rasta=image.getRGB(0,0,w,h,[],0,w); %get RGB data from bufferedimage

%convert java color integers to matlab RGB format:
rasta=256^3+rasta;
B=uint8(mod(rasta,256));
G=uint8(mod((rasta-int32(B))./256,256));
R=uint8(mod((rasta-256*int32(G))./65536,256));

X.cdata=uint8(zeros(h,w,3));
X.cdata(:,:,1)=reshape(R,[w h])';
X.cdata(:,:,2)=reshape(G,[w h])';
X.cdata(:,:,3)=reshape(B,[w h])';
X.colormap=[];

if (nargout == 2)
    m=X.colormap;
    X=X.cdata;
end

author: brainstorm-tools


text:
function [g1, g2, g3] = sdlfmaXsdlfmaKernGradientBlockIGJ(lfmKern1, lfmKern2, ...
    t1, t2,  i, j, generalConst, generalConstGrad, covGrad, g1Mean, ...
    g2Mean, gsp1Mean, gsp2Mean, typeMeanParam, typeMeanSwitching, ...
    typeKernParam1, typeKernParam2, typeKernSwitching1, typeKernSwitching2)

% SDLFMAXSDLFMAKERNGRADIENTBLOCKIGJ
% FORMAT
% DESC computes the gradients of the parameters for system 1 and system 2
% when i is greater than j. It assummes both systems represent
% accelerations.
% ARG lfmKern1 : structure containing parameters for the system 1
% ARG lfmKern2 : structure containing parameters for the system 2
% ARG t1 : times at which the system 1 is evaluated
% ARG t2 : times at which the system 2 is evaluated
% ARG covGrad : partial derivatives of the objective function wrt portion
% of the corresponding kernel matrix
% ARG g1Mean : gradients of the parameter of the system 1 obtained from the
% part of the function that uses the funcitons accommpanying the initial
% conditions.
% ARG g2Mean : gradients of the parameter of the system 2 obtained from the
% part of the function that uses the funcitons accommpanying the initial
% conditions.
% ARG gsp1Mean : gradient of the switching point of the system 1 obtained
% from the part of the function that uses the funcitons accommpanying the
% initial conditions.
% ARG gsp2Mean : gradient of the switching point of the system 2 obtained
% from the part of the function that uses the funcitons accommpanying the
% initial conditions.
% ARG typeMeanParam : specify the mean functions used to compute this part
% of the kernel
% ARG typeMeanSwitching : specify the functions used to compute the
% gradients of the swicthing points in part thst uses mean functions
% ARG typeKernParam1 : specify the first kernel function used to compute
% this part of the kernel
% ARG typeKernParam2 : specify the second kernel function used to compute
% this part of the kernel
% ARG typeKernSwitching1 : specify the functions used to compute the
% gradients of the swicthing points in both sides of the kernel function 1
% ARG typeKernSwitching2 : specify the functions used to compute the
% gradients of the swicthing points in both sides of the kernel function 2
% RETURN g1 : gradients of parameters for the system 1
% RETURN g2 : gradients of parameters for the system 2
% RETURN g3 : gradients of switching points
%
% COPYRIGHT : Mauricio A. Alvarez, 2010.

% KERN

if nargin < 14
    typeMeanParam = 'sdlfm';
    typeMeanSwitching = 'sdlfmv';
    typeKernParam1 = 'lfmXlfm';
    typeKernParam2 = 'lfmvXlfm';
    typeKernSwitching1= {'lfmvXlfm', 'lfmvXlfm'};
    typeKernSwitching2 ={'lfmaXlfm', 'lfmvXlfmv'};
end

g3 = [];

fhandleMeanParam = str2func([typeMeanParam 'MeanCompute']);
fhandleMeanGradParam = str2func([typeMeanParam 'MeanGradient']);
fhandleMeanGradSwitching = str2func([typeMeanSwitching 'MeanCompute']);
fhandleKernPosPos = str2func([typeKernParam1 'KernCompute']);
fhandleKernGradPosPos = str2func([typeKernParam1 'KernGradient']);
fhandleKernGradSwitchingPosPos1 = str2func([typeKernSwitching1{1} 'KernCompute']);
fhandleKernGradSwitchingPosPos2 = str2func([typeKernSwitching1{2} 'KernCompute']);
fhandleKernVelPos = str2func([typeKernParam2 'KernCompute']);
fhandleKernGradVelPos = str2func([typeKernParam2 'KernGradient']);
fhandleKernGradSwitchingVelPos1 = str2func([typeKernSwitching2{1} 'KernCompute']);
fhandleKernGradSwitchingVelPos2 = str2func([typeKernSwitching2{2} 'KernCompute']);

c1 = fhandleMeanParam(lfmKern1(1), t1, 'Pos');
e1 = fhandleMeanParam(lfmKern1(1), t1, 'Vel');

if isempty(generalConst{i,j})
    coeffPosPos = c1;
    coeffVelPos = e1;
else
    coeffPosPos = generalConst{i,j}(1,1)*c1 + generalConst{i,j}(2,1)*e1;
    coeffVelPos = generalConst{i,j}(1,2)*c1 + generalConst{i,j}(2,2)*e1;
end
g1Kern = zeros(length(lfmKern1), 5);
g2Kern = zeros(length(lfmKern1), 5);
PosPos = zeros(1,length(t2));
VelPos = zeros(1,length(t2));

for k=1:length(lfmKern1)
    PosPos = PosPos + fhandleKernPosPos(lfmKern2(k), lfmKern1(k), t2, ...
        lfmKern2(k).limit).';
    [g2KLocal, g1KLocal] = fhandleKernGradPosPos(lfmKern2(k), lfmKern1(k), ...
        t2, lfmKern2(k).limit, covGrad.', coeffPosPos.');
    g1Kern(k,:) = g1Kern(k, :) + g1KLocal;
    g2Kern(k,:) = g2Kern(k, :) + g2KLocal;
    VelPos = VelPos + fhandleKernVelPos(lfmKern2(k), lfmKern1(k), ...
        t2, lfmKern2(k).limit).';
    [g2KLocal, g1KLocal] = fhandleKernGradVelPos(lfmKern2(k), lfmKern1(k), ...
        t2, lfmKern2(k).limit, covGrad.', coeffVelPos.');
    g1Kern(k,:) = g1Kern(k, :) + g1KLocal;
    g2Kern(k,:) = g2Kern(k, :) + g2KLocal;
end
[gcAlpha, gcOmega] = fhandleMeanGradParam(lfmKern1(1), t1, 'Pos');
[geAlpha, geOmega] = fhandleMeanGradParam(lfmKern1(1), t1, 'Vel');
[gradAlpha, gradOmega] = getLocalGradAlphaOmega(lfmKern1);
g1Pos = fhandleMeanGradSwitching(lfmKern1(1), t1, 'Pos');
h1Vel = fhandleMeanGradSwitching(lfmKern1(1), t1, 'Vel');
gsp1Kern = zeros(1, length(lfmKern1));
gsp2Kern = zeros(1, length(lfmKern1));
% This means that are only two switching points involved, t1
% and t0 appear in k, like k_{ff}(t1-t0,t'-t0) and
% k_{mf}(t1-t0,t'-t0). The other points appear in c1(t - t1)
% and e1(t - t1)
for k=1:length(lfmKern1)
    temp = fhandleKernGradSwitchingPosPos1(lfmKern2(k), lfmKern1(k), t2, lfmKern2(k).limit).';
    gsp2Kern(k) = gsp2Kern(k) - sum(sum((coeffPosPos*temp).*covGrad));
    gsp1Kern(k) = gsp1Kern(k) + sum(sum((coeffPosPos*temp).*covGrad));
    temp = fhandleKernGradSwitchingPosPos2(lfmKern2(k), lfmKern1(k), t2, lfmKern2(k).limit).';
    gsp2Kern(k) = gsp2Kern(k) - sum(sum((coeffPosPos*temp).*covGrad));
    temp = fhandleKernGradSwitchingVelPos1(lfmKern1(k), lfmKern2(k), lfmKern2(k).limit, t2);
    gsp2Kern(k) = gsp2Kern(k) - sum(sum((coeffVelPos*temp).*covGrad));
    gsp1Kern(k) = gsp1Kern(k) + sum(sum((coeffVelPos*temp).*covGrad));
    temp = fhandleKernGradSwitchingVelPos2(lfmKern2(k), lfmKern1(k), t2, lfmKern2(k).limit).';
    gsp2Kern(k) = gsp2Kern(k) - sum(sum((coeffVelPos*temp).*covGrad));
end
if isempty(generalConst{i,j})
    matGradAlpha = gcAlpha*PosPos + geAlpha*VelPos;
    matGradOmega = gcOmega*PosPos + geOmega*VelPos;
    gCoeff = gradAlpha*sum(sum(matGradAlpha.*covGrad)) + ...
        gradOmega*(sum(sum(matGradOmega.*covGrad)));
    matGradSp1 = g1Pos*PosPos + h1Vel*VelPos;
    gsp1 = - sum(sum(matGradSp1.*covGrad));
    g3(i) = gsp1Mean + sum(gsp1Kern) + gsp1; % switching point 1
    g3(j) = gsp2Mean + sum(gsp2Kern);        % switching point 2
else
    constGradAlpha = generalConstGrad{1}{i,j};
    constGradOmega = generalConstGrad{2}{i,j};
    constVal = generalConst{i,j};
    % Derivative wrt parameters
    matGradAlpha = (constVal(1,1)*gcAlpha + constGradAlpha(1,1)*c1 ...
        + constVal(2,1)*geAlpha + constGradAlpha(2,1)*e1)*PosPos ...
        + (constVal(1,2)*gcAlpha + constGradAlpha(1,2)*c1 ...
        + constVal(2,2)*geAlpha + constGradAlpha(2,2)*e1)*VelPos;
    matGradOmega = (constVal(1,1)*gcOmega + constGradOmega(1,1)*c1 ...
        + constVal(2,1)*geOmega + constGradOmega(2,1)*e1)*PosPos ...
        + (constVal(1,2)*gcOmega + constGradOmega(1,2)*c1 ...
        + constVal(2,2)*geOmega + constGradOmega(2,2)*e1)*VelPos;
    gCoeff = gradAlpha*sum(sum(matGradAlpha.*covGrad)) + ...
        gradOmega*(sum(sum(matGradOmega.*covGrad)));
    % Derivative wrt switching points
    % Firts, notice that gsp1Kern doesn't correspond to the switching point
    % of the interval i (or say 1). It's just the derivative of the inner
    % switching point that lead to the actual switching point for the
    % interval 1. So we rename it first.
    gspInt = sum(gsp1Kern);
    % Compute the derivative of the swicthing point i, not appearing in the
    % constant
    matGradSp1 = (constVal(1,1)*g1Pos + constVal(2,1)*h1Vel)*PosPos + ...
        (constVal(1,2)*g1Pos + constVal(2,2)*h1Vel)*VelPos;
    gsp1 = - sum(sum(matGradSp1.*covGrad));
    % Compute the derivatives for all the switching points apearing in the
    % constant
    constGradSPoint = generalConstGrad{3}{i,j};
    numberSP = size(constGradSPoint,2);
    gspInBetween = zeros(1, numberSP);
    for k=1:numberSP
        temp = (constGradSPoint(1,k)*c1 + constGradSPoint(3,k)*e1)*PosPos + ...
            (constGradSPoint(2,k)*c1 + constGradSPoint(4,k)*e1)*VelPos;
        gspInBetween(k) = sum(sum(temp.*covGrad));
    end
    % Assign derivatives wrt all other switching points, with a correction
    % for the derivative of the innermost switching point in the constant
    gspInBetween(end) =  gspInBetween(end) + gspInt;
    g3(i) = gsp1Mean + gsp1 + gspInBetween(1);
    g3(j) = gsp2Mean + sum(gsp2Kern);
    g3(j+1:i-1) = fliplr(gspInBetween(2:end));
end
% Assign derivatives wrt first system
g1{1} = g1Mean + sum(g1Kern(:,1:3), 1) + gCoeff; % mass 1, spring 1, damper 1
g1{2} = g1Kern(:,4).';                            % inverse widths
g1{3} = g1Kern(:,5).';                            % seinsitivities 1
% Assign derivatives wrt second system
g2{1} = g2Mean + sum(g2Kern(:,1:3), 1); % mass 2, spring 2, damper 2
g2{2} = g2Kern(:,4).';                   % inverse widths
g2{3} = g2Kern(:,5).';

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%function S = sampleFromMatrixNormal(M,V,K,nSamples=1)
function S = sampleFromMatrixNormal(M,sqrtV,sqrtinvK,nSamples)

if ~exist('nSamples','var'), nSamples = 1; end

[mu,sqrtsigma] = matrixNormalToNormal(M,sqrtV,sqrtinvK);

S = mu + sqrtsigma'*randn(length(mu),1);
S = reshape(S,size(M));

author: michaelchughes

repo: NPBayesHMM


text:
function Y=sparse(varargin)
%SPARSE (overloaded)

if nargin < 3
    error('At-least 3 arguments needed');
end

data = varargin{3};
ns = varargin{1}(:);
ms = varargin{2}(:);

if length(ns)~=length(ms) | length(ms)~=length(data)
    error('Length of first 3 arguments must be equal');
end

if min(size(data))>1
    error('Third argument should be a vector');
end

if nargin < 4
    n = max(ns);
else
    n = varargin{4};
end
if nargin < 5
    m = max(ms);
else
    m = varargin{5};
end

if any(ms>m)
    error('Dimension mismatch')
end

if any(ns>n)
    error('Dimension mismatch')
end

Y = data;
Y.dim(1) = n;
Y.dim(2) = m;
[i1,j1,s1] = find(data.basis);
ind = ns+(ms-1)*n;
Y.basis = sparse(ind(i1),j1,s1,n*m,1+length(Y.lmi_variables));
% Reset info about conic terms
Y.conicinfo = [0 0];

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
classdef TestTempfile
    %TestTempfile

    methods (Static)
        function test_1
            fname = cv.tempfile();
            validateattributes(fname, {'char'}, {'vector'});
        end

        function test_2
            fname = cv.tempfile('Suffix','.txt');
            validateattributes(fname, {'char'}, {'vector'});
            assert(strncmp(fname(end-3:end), '.txt', 4));
        end

        function test_error_unrecognized_option
            try
                cv.tempfile('foo','bar');
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function setM(OCP,M)
    setM@DynamicSystem(OCP,M);
end

author: deng-haoyang

repo: ParNMPC


text:
function x = subsasgn( x, S, y )
error( 'Subscripted assignment not allowed for cvx problem objects.' );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function X = MvnRndMatchCrossCov(S, varargin)
% X = MvnRndMatchCrossCov(S, J) simulates a panel of J scenarios of a
% multivariate normal random vector with zero sample mean and sample
% covariance exactly equal to S. Each row of X is a scenario and each
% column is a variable.
%
% X = MvnRndMatchCrossCov(S, P, F, CovF) generates a multivariate normal
% panel X with zero sample mean, sample covariance S and sample cross
% covariance P with another zero-mean panel F. CovF, the sample covariance
% of F, is computed if not provided as input.

% IMPORTANT: The order of computations in the following expressions is
% important for efficiency. Please use the same order, or have a good reason
% to change it!

% Code by S. Gollamudi. This version December 2009.

% check input
N = size(S,1);
if nargin==2,
    J = varargin{1};
    K = 0;
else
    P = varargin{1};
    F = varargin{2};
    [J,K] = size(F);
    assert(all(size(P)==[K,N]), 'Incorrect dimensions of cross-covariance P.');
    if length(varargin)<=2,
        CovF = cov(F,1);
    else
        CovF = varargin{3};
    end
    sqrtJ = sqrt(J);
end

% tolerance level
eps = 1e-9;

% Algorithm: X = F*B + V*R
% where F'*V=0, V'V=J*I and R is upper triangular

if K>0,

    % Compute B
    B = CovF \ P;

    % Compute R
    Sigma_orth = S - P'*B;
    var_orth = diag(Sigma_orth);
    var_x = diag(S);
    if any(var_orth < -eps*var_x),
        error('MvnRndMatchCrossCov: S and P are incompatible.')
    end
    i_positive = (var_orth > eps*var_x);
    num_positive = sum(i_positive);
    R = zeros(num_positive,N);
    R(:,i_positive) = chol(S(i_positive,i_positive) - P(:,i_positive)'*B(:,i_positive));

    % Generate V
    Z = randn(J/2,N);
    Z=[Z
        -Z];

    % using QR decomposition
    %[V,dummy] = qr([F Z],0);
    %V = V(:,K+1:K+num_positive)*sqrtJ;

    % using modified Gram-Schmidt
    Z = Z - (F/CovF)*((F'*Z)/J);
    V = zeros(J,0);
    for n = 1:num_positive,
        u = Z(:,n) - V*((V'*Z(:,n))/J);
        V = [V u*(sqrtJ/norm(u))];
    end

    % generate panel X
    X = F*B + V*R;

else

    % generate X to have zero sample mean and sample covariance S
    Z = randn(J/2,N);
    Z = [Z
        -Z];
    X = Z * (chol(cov(Z,1)) \ chol(S));

end

author: Sable

repo: mcbench-benchmarks


text:
function size = varsize(X, units)
%VARSIZE Returns the size of a variable in bytes.
% Usage:
%   size = varsize(X)
%   size = varsize(X, units) % 'bytes' (default), 'KB', 'MB', 'GB'
%
% See also: whos
narginchk(1, 2)
if nargin < 2
    units = 'bytes';
end

% Get variable size in bytes
S = whos('X');
size = S.bytes;

if nargout < 1
    % Just output to console if not storing variable
    disp(bytes2str(size, 4))
    clear size
    return
end

% Convert to units
switch lower(units)
    case {'kb', 'k', 'kilo', 'kilobyte', 'kib'}
        size = size / 1024;
    case {'mb', 'm', 'mega', 'megabyte', 'mib'}
        size = size / 1024 / 1024;
    case {'gb', 'g', 'giga', 'gigabyte', 'gib'}
        size = size / 1024 / 1024 / 1024;
    case {'auto', 'str'}
        size = bytes2str(size, 4);
    otherwise
        return
end

end

author: talmo

repo: leap

save_path: github-repos/MATLAB/talmo-leap


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

%%
clear; close all; clc;

%% Syntax
% |Cmapped=cmaperise(C,cmap,clim);|

%% Description
% This function creates RGB colors for the data C using the colormap cmap
% and the limits clim.

%% Examples
%

[X,Y,Z]=peaks(50);
[F,V,C]=surf2patch(X,Y,Z,Z);

clim=[min(C(:)) max(C(:))];

%RGB color set with map 1
cmap1=viridis(25);
C_rgb_1=cmaperise(C,cmap1,clim);

%RGB color set with map 2
cmap2=gray(25);
C_rgb_2=cmaperise(C,cmap2,clim);

%RGB color set with map 3
cmap3=gjet(25);
C_rgb_3=cmaperise(C,cmap3,clim);

%% Visualize

% Offset coordinates so plots are side by side
V2=V;
V2(:,1)=V2(:,1)-min(V2(:,1))+max(V(:,1));
V3=V2;
V3(:,1)=V3(:,1)-min(V3(:,1))+max(V2(:,1));

cFigure;

subplot(1,2,1);
title('Colormapped')
hp=gpatch(F,V,C,'none');
legend(hp,'Colormapped data');
axisGeom;
camlight headlight; colormap(gca,cmap1); colorbar;

subplot(1,2,2);
title('RGB (red-green-blue) painted')
hp1=gpatch(F,V,C_rgb_1,'none');
hp2=gpatch(F,V2,C_rgb_2,'none');
hp3=gpatch(F,V3,C_rgb_3,'none');
legend([hp1 hp2 hp3],{'RGB colored with map 1','RGB colored with map 2','RGB colored with map 3'});
axisGeom;
camlight headlight;

drawnow;

%%
%
% <<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 [ulf, vlf, lf, hf, lfhf, ttlpwr] = CalcLfHfParams(PSD, F, limits,plot_on)
% [ulf, vlf, lf, hf, lfhf, ttlpwr] = CalcLfHfParams(PSD, F, limits,plot_on)
%
%   OVERVIEW: Compute the frequency domain features for a given PSD and
%             frequency bans limits
%
%   INPUT:
%        PSD     - power spectral density
%        F       - frequency vector
%        limits  - frequency domain analysis limits
%        plot_on -
%
%   OUTPUT:
%	- ulf     : (ms^2) Power in the ultra low frequency range (default < 0.003 Hz)
%	- vlf     : (ms^2) Power in very low frequency range (default 0.003 <= vlf < 0.04 Hz)
%	- lf      : (ms^2) Power in low frequency range (default 0.04Hz  <= lf < 0.15 Hz)
%	- hf      : (ms^2) Power in high frequency range (default 0.15 <= hf < 0.4 Hz)
%	- lfhf    : Ratio LF [ms^2]/HF [ms^2]
%	- ttlpwr  : (ms^2) Total spectral power (approximately <0.4 Hz)
%
%	REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Main script written by Adriana N. Vest
%       Dependent scripts written by various authors
%       (see functions for details)
%	COPYRIGHT (C) 2016
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information
%%
if nargin <3
    ULF = [0 .003];
    VLF = [0.003 .04];
    LF = [.04 .15];
    HF = [0.15 0.4];
    limits = [ULF; VLF; LF; HF];
end
if nargin < 4
    plot_on =1;
end

Indx_ULF = find( (limits(1,1) <= F) & (F <= limits(1,2)) );
Indx_VLF = find( (limits(2,1) <= F) & (F <= limits(2,2)) );
Indx_LF = find( (limits(3,1) <= F) & (F <= limits(3,2)) );
Indx_HF = find( (limits(4,1) <= F) & (F <= limits(4,2)) );
space = F(2)-F(1);

ulf = sum(PSD(Indx_ULF)*space) * 1e6; % convert to ms^2
vlf = sum(PSD(Indx_VLF)*space) * 1e6; % convert to ms^2
lf = sum(PSD(Indx_LF)*space) * 1e6;   % convert to ms^2
hf = sum(PSD(Indx_HF)*space) * 1e6;   % convert to ms^2

ttlpwr = sum([ulf vlf lf hf]);

lf_n = lf/ttlpwr; % normalized
hf_n = hf/ttlpwr;
lfhf = round(lf_n/hf_n*100)/100; % lf/hf ratio

if plot_on
    figure
    % plot PSD
    plot(F,10*log10(PSD),'b','linewidth',2)
    hold on
    % plot limits on graph for lf and hf
    plot([F(Indx_LF(1)) F(Indx_LF(1))],[-80 40],'k:')
    hold on
    plot([F(Indx_LF(end)) F(Indx_LF(end))],[-80 40],'k:')
    hold on
    plot([F(Indx_HF(end)) F(Indx_HF(end))],[-80 40],'k:')

    % labelsc
    text(0.07,30,'LF','Fontname','Times New Roman','Fontsize',10)
    text(0.25,30,'HF','Fontname','Times New Roman','Fontsize',10)
    %text(0.15, 35, 'Power Spectral Density','Fontname','Times New Roman','Fontsize',10)
    text(0.3, -60, strcat('LF/HF=',num2str(lfhf)),'Fontname','Times New Roman','Fontsize',10)
    ylabel('Normalized PSD (db/Hz)','Fontname','Times New Roman','fontsize',10)
    xlabel('Frequency (Hz)','Fontname','Times New Roman','fontsize',10)
    axis([0 .45 -80 40]);
    box off
end % end plot

end % end function

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function [xsecx,xsecy] = lc_xsect(eqlat,eqlon,depth,width,length,...
        lat1,lon1,lat2,lon2)

    %LC_XSECTION
    %
    %	[xsecx, xsecy] = LC_xsection(eqlat,eqlon,depth,width,length,...
    %                                        lat1,lon1,lat2,lon2)        (1)
    %
    %	[xsecx, xsecy] = LC_xsection(eqlat,eqlon,depth,width,length,...
    %                                        lat0,lon0,azimuth)          (2)
    %
    %	[xsecx, xsecy] = LC_xsection(eqlat,eqlon,depth,width)         (3)
    %
    %	Function to make a cross section of data points on a map
    %	created by LC_MAP (Lambert Conformal).
    %	The WIDTH of the zone from which the data points is given
    %	in "km" and represent the total width (1/2 on one side, 1/2 on
    %	the other side).
    %	The LENGTH of the xsection in "km" is only used in method (2),
    %	but this argument is still neccessary in argument list of
    %	method (1) in order to keep the argument list in the right order.
    %	The data points location are given by EQLAT, EQLON and DEPTH.
    %	The cross section location can be given in any one of three ways:
    %
    %	  (1) given latitudes and longitudes of two points on the map,
    %             using the arguments as described above.
    %
    %	  (2) given latitude and longitude of a center point and an azimuth.
    %
    %	  (3) using the cursor to select two points on the map by clicking
    %	      a mouse button above the desired points.
    %
    %	If the output argument is used, the distance-depth data is kept
    %	in variables for other use; otherwise the xsection will be plotted
    %	on a new figure window.
    %
    %	It is possible to set the symbol type, size and line width by
    %	setting the following global variables:
    %	"symb_type", "symb_size" and "symb_width" respectively. Otherwise
    %	it will use the defaults.
    %
    %	It is also possible to set the minimum and maximum depth of the
    %	cross-section by setting the following global variables:
    %	"mindepth" and "maxdepth".  If either or both are not set, it
    %	will use 0 km as the minimum depth and/or the depth of the deepest
    %	data point as the maximum depth.
    %
    %	NOTE:
    %	It is assumed that LC_MAP was used before using this function!
    %	This is neccessary to set global variables used by this function.

    report_this_filefun(mfilename('fullpath'));

    global torad Re scale
    global sine_phi0 phi0 lambda0 phi1 phi2
    global maxlatg minlatg maxlong minlong
    global symb_type symb_size symb_width
    global label1 label2
    global mindepth maxdepth

    todeg = 180 / pi;

    if nargin < 9

        if nargin == 8	% method 2: given lat & lon of center point and angle

            lat0 = lat1;
            lon0 = lon1;
            [x0, y0] = lc_tocart(lat0,lon0);
            azimuth = lat2;

            if azimuth >= 180, azimuth = azimuth - 180; end
            theta0 = ((lon0*torad - lambda0) * sine_phi0) * todeg;
            alpha = azimuth - theta0;
            beta = (90 - azimuth) + theta0;

            x2 = ((length / 2) * cos(beta*torad));
            y2 = ((length / 2) * sin(beta*torad));
            x1 = x0 - x2;
            y1 = y0 - y2;
            x2 = x0 + x2;
            y2 = y0 + y2;

            [lat1, lon1] = lc_froca(x1,y1);
            [lat2, lon2] = lc_froca(x2,y2);

        elseif nargin == 4	% method 3: selection of the end points by mouse

            limits = ginput(2);
            x1 = limits(1,1);
            y1 = limits(1,2);
            x2 = limits(2,1);
            y2 = limits(2,2);

            if x1 > x2
                xtemp = x1; ytemp = y1;
                x1 = x2; y1 = y2;
                x2 = xtemp; y2 = ytemp;
            end

            [lat1, lon1] = lc_froca(x1,y1);
            [lat2, lon2] = lc_froca(x2,y2);

            x0 = (x1 + x2) / 2;
            y0 = (y1 + y2) / 2;
            [lat0, lon0] = lc_froca(x0,y0);
            dx = x2 - x1;
            dy = y2 - y1;

            alpha = 90 - (atan(dy/dx)*todeg);
            length = sqrt(dx^2 + dy^2);

        else
            disp('ERROR: incompatible number of arguments')
            help lc_xsection
            return
        end

    elseif nargin == 9	% method 1: given lat & lon of the two end points

        [x1, y1] = lc_tocart(lat1,lon1);
        [x2, y2] = lc_tocart(lat2,lon2);

        if x1 > x2
            xtemp = x1; ytemp = y1;
            x1 = x2; y1 = y2;
            x2 = xtemp; y2 = ytemp;
        end

        x0 = (x1 + x2) / 2;
        y0 = (y1 + y2) / 2;
        [lat0, lon0] = lc_froca(x0,y0);
        dx = x2 - x1;
        dy = y2 - y1;

        alpha = 90 - (atan(dy/dx)*todeg);
        length = sqrt(dx^2 + dy^2);

    else

        disp('ERROR: incompatible number of arguments')
        help lc_xsection
        return
    end

    % correction factor to correct for longitude away from the center meridian
    theta0 = ((lon0*torad - lambda0) * sine_phi0) * todeg;

    % correct the XY azimuth of the Xsection line with the above factor to obtain
    % the true azimuth
    azimuth = alpha + theta0;
    if azimuth < 0, azimuth = azimuth + 180; end

    % convert XY coordinate azimuth to a normal angle like we used to deal with
    sigma = 90 - alpha;

    % transformation matrix to rotate the data coordinate w.r.t the Xsection line
    transf = [cos(sigma*torad) sin(sigma*torad)
        -sin(sigma*torad) cos(sigma*torad)];

    % inverse transformation matrix to rotate the data coordinate back
    invtransf = [cos(-sigma*torad) sin(-sigma*torad)
        -sin(-sigma*torad) cos(-sigma*torad)];

    % convert the map coordinate of the events to cartesian coordinates
    idx_map = find(minlatg < eqlat & eqlat < maxlatg & ...
        minlong < eqlon & eqlon < maxlong);
    [eq(1,:) eq(2,:)] = lc_tocart(eqlat,eqlon);

    % create new coordinate system at center of Xsection line
    eq0(1,:) = eq(1,:) - x0;
    eq0(2,:) = eq(2,:) - y0;

    % rotate this last coordinate system so that X-axis correspond to Xsection line
    eq0p = transf * eq0;

    % project the event data to the Xsection line
    eq1(1,:) = eq0p(1,:);
    eq1(2,:) = eq0p(2,:) * 0;

    % convert back to the original coordinate system
    eq1p = invtransf * eq1;
    eq2(1,:) = eq1p(1,:) + x0;
    eq2(2,:) = eq1p(2,:) + y0;

    % plot the Xsection line on the map
    plot([x1 x2],[y1 y2],'--','LineWidth',1.5)

    % label the Xsection end points
    xlim = get(gca,'XLim');
    ylim = get(gca,'YLim');
    label_dist = (2 / 100) * sqrt((2*xlim(2))^2 + (2*ylim(2))^2);
    label_pt(1,1) = -(length/2 + label_dist);
    label_pt(2,1) = 0;
    label_pt(1,2) = length/2 + label_dist;
    label_pt(2,2) = 0;
    rlabel_pt = invtransf * label_pt;
    label_pt(1,:) = rlabel_pt(1,:) + x0;
    label_pt(2,:) = rlabel_pt(2,:) + y0;
    if isempty(label1), label1 = 'A'; end
    if isempty(label2), label2 = 'A'''; end
    lbl1_h = text(label_pt(1,1),label_pt(2,1),label1,'FontSize',14,...
        'Vertical','middle','Horizontal','center','FontWeight','bold');
    lbl2_h = text(label_pt(1,2),label_pt(2,2),label2,'FontSize',14,...
        'Vertical','middle','Horizontal','center','FontWeight','bold');

    % create a box of width "width" around the Xsection line and plot it
    box(1,1) = -length/2; box(2,1) = width/2;
    box(1,2) = length/2; box(2,2) = width/2;
    box(1,3) = length/2; box(2,3) = -width/2;
    box(1,4) = -length/2; box(2,4) = -width/2;
    xbox = [box(1,:) box(1,1)];
    ybox = [box(2,:) box(2,1)];
    rbox = invtransf * [xbox ; ybox];
    rbox(1,:) = rbox(1,:) + x0;
    rbox(2,:) = rbox(2,:) + y0;
    plot(rbox(1,:),rbox(2,:),'-.k','LineWidth',0.8)

    % check if symbol parameters global variables are set, if not --> defaults
    if isempty(symb_type), symb_type = '+'; end
    if isempty(symb_size), symb_size = 6; end
    if isempty(symb_width), symb_width = [0.5]; end

    % plot the events on the map
    plot(eq(1,idx_map),eq(2,idx_map),symb_type,'MarkerSize',symb_size,...
        'LineWidth',symb_width)

    % find index of all events which are within the given box width
    idx_box = find(abs(eq0p(2,:)) <= width/2 & abs(eq0p(1,:)) <= length/2);

    % Open another graphic window for the cross section
    map_fig = gcf;
    xsec_fig = map_fig + 1;
    figure_w_normalized_uicontrolunits(xsec_fig)
    set(xsec_fig,'PaperPosition',[1 .5 9 6.9545])

    % Plot events on cross section figure
    xdist = eq1(1,idx_box) + (length / 2);

    global Xwbz Ywbz
    Xwbz = xdist;
    Ywbz = depth(idx_box);

    xsecx = xdist;
    xsecy = depth(idx_box);

    plot(xdist,depth(idx_box),symb_type,'MarkerSize',symb_size,...
        'LineWidth',symb_width,'YDir','reverse')

    if isempty(maxdepth)
        maxZ = max(depth(idx_box));
    else
        maxZ = maxdepth;
    end

    if isempty(mindepth)
        minZ = 0;
    else
        minZ = mindepth;
    end

    if length > (maxZ - minZ)*11/8.5
        position = [.1 .1 .7 ((maxZ-minZ)/length)*0.7*11/8.5];
    else
        position = [.1 .1 (length/(maxZ-minZ))*0.7*8.5/11 .7];
    end
    set(gca,'Position',position,'XLim',[0 length],'Ylim',[-maxZ -minZ],...
        'LineWidth',2)

    % Plot labels
    Xstring = ['Distance from ' label1 ' (km)'];
    set(gca,'XLabel',text(0,0,Xstring),'YLabel',text(0,0,'Depth (km)'))
    label_base1 = 1 + .04;
    label_base2 = 1 + .06;
    label_base3 = 1 + .08;
    lbl3_h = text(0,label_base2,label1,'FontSize',14,'Horizontal','center',...
        'FontWeight','bold','Vertical','middle','Units','norm');
    lat1_dm = sprintf('    %2.2i N %4.2f''',fix(lat1),(frac(lat1)*60));
    lon1_dm = sprintf('   %3.3i W %4.2f''',abs(fix(lon1)),abs(frac(lon1)*60));
    lbl5_h = text(0,label_base1,lat1_dm,'FontSize',12,'Horizontal','left',...
        'Vertical','bottom','Units','norm');
    lbl6_h = text(0,label_base3,lon1_dm,'FontSize',12,'Horizontal','left',...
        'Vertical','bottom','Units','norm');

    lbl4_h = text(1,label_base2,label2,'FontSize',14,'Horizontal','center',...
        'FontWeight','bold','Vertical','middle','Units','norm');
    lat2_dm = sprintf('%2.2i N %4.2f''    ',fix(lat2),(frac(lat2)*60));
    lon2_dm = sprintf('%3.3i W %4.2f''    ',abs(fix(lon2)),abs(frac(lon2)*60));
    lbl7_h = text(1,label_base1,lat2_dm,'FontSize',12,...
        'Horizontal','right','Vertical','bottom','Units','norm');
    lbl8_h = text(1,label_base3,lon2_dm,'FontSize',12,...
        'Horizontal','right','Vertical','bottom','Units','norm');

    % Go back to map figure
    %figure_w_normalized_uicontrolunits(map_fig)

author: CelsoReyes


text:
function [a,b,c] = FCE(L,t)
% finite strain ellipsoid
%
% Syntax
%   [a,b,c] = FCE(L,t)
%   r = FCE(l,t)
%
% Input
%  L -
%  t - time
%
% Output
%

Gamma = L.vorticity;

if nargout == 1

  % r = log(a/c)
  a = 2 * asinh(sinh(sqrt(1-Gamma.^2)) .* L.strainRate .* t ./ sqrt(1-Gamma.^2));

else

end

author: mtex-toolbox

repo: mtex

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


text:
%%*************************************************************************
%% symqmr: symmetric QMR with left (symmetric) preconditioner.
%%         The preconditioner used is based on the analytical
%%         expression of inv(A).
%%
%% [x,resnrm,solve_ok] = symqmr(A,b,L,tol,maxit)
%%
%% child function: linsysolvefun.m
%%
%% A = [mat11 mat12; mat12' mat22].
%% b = rhs vector.
%% if matfct_options = 'chol' or 'spchol'
%%    L = Cholesky factorization of (1,1) block.
%%    M = Cholesky factorization of
%%        Schur complement of A ( = mat12'*inv(mat11)*mat12-mat22).
%% else
%%    L = triangular factors of A.
%%    M = not relevant.
%% end
%% resnrm = norm of qmr-generated residual vector b-Ax.
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*************************************************************************

   function  [xx,resnrm,solve_ok] = symqmr(A,b,L,tol,maxit,printlevel)

   N = length(b);
   if (nargin < 6); printlevel = 1; end
   if (nargin < 5) | isempty(maxit); maxit = max(30,length(A.mat22)); end;
   if (nargin < 4) | isempty(tol); tol = 1e-10; end;
   tolb = min(1e-4,tol*norm(b));

   solve_ok = 1;
   x = zeros(N,1);
   if (norm(x))
      if isstruct(A); Aq = matvec(A,x); else; Aq=A*x; end;
      r = b-Aq;
   else
      r = b;
   end
   err = norm(r); resnrm(1) = err; minres = err; xx = x;
   if (err < 1e-3*tolb); return; end

   q = precond(A,L,r);
   tau_old   = norm(q);
   rho_old   = r'*q;
   theta_old = 0;
   d = zeros(N,1);
   res = r; Ad = zeros(N,1);
%%
%% main loop
%%
   tiny = 1e-30;
   for iter = 1:maxit

       if isstruct(A); Aq = matvec(A,q); else; Aq=A*q; end;
       sigma = q'*Aq;
       if (abs(sigma) < tiny)
          solve_ok = 2;
          if (printlevel); fprintf('*'); end;
          break;
       else
          alpha = rho_old/sigma;
          r = r - alpha*Aq;
       end
       u = precond(A,L,r);

       theta = norm(u)/tau_old; c = 1/sqrt(1+theta^2);
       tau = tau_old*theta*c;
       gam = (c^2*theta_old^2); eta = (c^2*alpha);
       d = gam*d + eta*q;
       x = x + d;
%%
       Ad = gam*Ad + eta*Aq;
       res = res - Ad;
       err = norm(res); resnrm(iter+1) = err;
       if (err < minres); xx = x; minres = err; end
       if (err < tolb); break; end
       if (iter > 10)
          if (err > 0.98*mean(resnrm(iter-10:iter)))
             solve_ok = 0.5; break;
          end
       end
%%
       if (abs(rho_old) < tiny)
          solve_ok = 2;
          if (printlevel); fprintf('*'); end;
          break;
       else
          rho  = r'*u;
          beta = rho/rho_old;
          q = u + beta*q;
       end
       rho_old = rho;
       tau_old = tau;
       theta_old = theta;
   end
   if (iter == maxit); solve_ok = 0.3; end;
%%
%%*************************************************************************
%% precond:
%%*************************************************************************

   function Mx = precond(A,L,x)

   m = length(L.perm); m2 = length(x)-m;
   Mx = zeros(length(x),1);

   for iter = 1:1
      if norm(Mx); r = full(x - matvec(A,Mx)); else; r = full(x); end
      r1 = r(1:m);
      if (m2 > 0)
         r2 = r(m+[1:m2]);
         w = linsysolvefun(L,r1);
         z = mexMatvec(A.mat12,w,1) - r2;
         z = L.Mu \ (L.Ml \ (L.Mp*z));
         r1 = r1 - mexMatvec(A.mat12,z);
      end
      d = linsysolvefun(L,r1);
      if (m2 > 0)
         d = [d; z];
      end
      Mx = Mx + d;
   end
%%*************************************************************************
%% matvec: matrix-vector multiply.
%% matrix = [A.mat11 A.mat12; A.mat12' A.mat22]
%%*************************************************************************

   function Ax = matvec(A,x);

   m = length(A.mat11); m2 = length(x)-m;

   if issparse(x); x = full(x); end
   if (m2 > 0)
      x1 = x(1:m);
   else
      x1 = x;
   end
   Ax = mexMatvec(A.mat11,x1);
   if (m2 > 0)
      x2 = x(m+[1:m2]);
      Ax = Ax + mexMatvec(A.mat12,x2);
      Ax2 = mexMatvec(A.mat12,x1,1) + mexMatvec(A.mat22,x2);
      Ax = [full(Ax); full(Ax2)];
   end
%%*************************************************************************


text:
function view = plotResidualError(view,scan,baseFrame)
%
% view = plotResidualError(view,scan,[baseFrame])
%
% Plots the rmse between each frame and a baseFrame
%
% If you change this function make parallel changes in:
%   plotMaxTSErr
%
% djh, 2/2001. rmse instead of var, conver to 3.0

if ~exist('scan','var')
   scan = getCurScan(view);
end
if ~exist('baseFrame','var')
   baseFrame=1;
   %baseFrame=nFrames;
end

slices = sliceList(view,scan);
nSlices = length(slices);
nFrames = numFrames(view,scan);

% Load tSerises and compute RMSE slice by slice
waitHandle = mrvWaitbar(0,'Loading tSeries and computing RMSE.  Please wait...');
vres = zeros(nSlices,nFrames);
for slice=slices
   mrvWaitbar(slice/nSlices);
   tSeries = loadtSeries(view,scan,slice);
   for frame=1:nFrames
      vres(slice,frame) = sqrt(mse(tSeries(frame,:),tSeries(baseFrame,:)));
   end
end
close(waitHandle)

% mean across slices
vres = mean(vres);

% dont plot anything for the base frame (otherwise it would be 0)
vres(baseFrame)=NaN;

% plot it
selectGraphWin;
fontSize = 14;
set(gcf,'Name',['RMSE, scan: ',num2str(scan),', baseframe:',int2str(baseFrame)]);
x = 1:nFrames;
plot(x,vres,'-b','LineWidth',2)
set(gca,'FontSize',fontSize)
set(gca,'XLim',[0,nFrames]);
xlabel('Frame number','FontSize',fontSize)
ylabel('RMSE (raw intensity units)','FontSize',fontSize)

% Save the data in gca('UserData')
data.frameNumbers = x;
data.tSeries  =  vres;
set(gca,'UserData',data);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [X,names]=elipcyl
% [X,names]=elipcyl defines elliptic
% cylinder coordinates
syms et ps z real; names=[et ps z];
X=[cosh(et)*cos(ps);
   sinh(et)*sin(ps); z];


text:
%% file example_rMTFL.m
% this file shows the usage of Least_rMTFL.m function
% and study how to detect outlier tasks.
%
%% OBJECTIVE
%  argmin_W ||X(P+Q) - Y||_F^2 + lambda1*||P||_{1,2} + lambda2*||Q^T||_{1,2}
%   s.t. W = P + Q
%
%% Copyright (C) 2012 Jiayu Zhou, and Jieping Ye
%
% You are suggested to first read the Manual.
% For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
% Last modified on April 16, 2012.
%
%% Related papers
%
% [1] Gong, P. and Ye, J. and Zhang, C. Robust Multi-Task Feature Learning,
% Submitted, 2012
%

clear;
clc;
close all;

addpath('../MALSAR/functions/rMTFL/'); % load function
addpath('../MALSAR/utils/'); % load utilities

%rng('default');     % reset random generator. Available from Matlab 2011.

%generate synthetic data.
dimension = 500;
sample_size = 50;
task = 50;
X = cell(task ,1);
Y = cell(task ,1);
for i = 1: task
    X{i} = rand(sample_size, dimension);
    Y{i} = rand(sample_size, 1);
end

opts.init = 0;      % guess start point from data.
opts.tFlag = 1;     % terminate after relative objective value does not changes much.
opts.tol = 10^-6;   % tolerance.
opts.maxIter = 500; % maximum iteration number of optimization.

rho_1 = 90;%   rho1: P
rho_2 = 280; %   rho2: Q

[W funcVal P Q] = Least_rMTFL(X, Y, rho_1, rho_2, opts);

% draw figure
close;
figure();
subplot(3,1,1);
%imshow(1- (abs(S')~=0), 'InitialMagnification', 'fit');
imshow(P'==0, 'InitialMagnification', 'fit')
ylabel('P^T (feature)');
title('Visualization of Robust Multi-Task Feature Learning Model');
subplot(3,1,2);
%imshow(1- (zscore(L')), 'InitialMagnification', 'fit')
imshow(Q'==0, 'InitialMagnification', 'fit')
ylabel('Q^T (outliers)');
subplot(3,1,3);
%imshow(1- (zscore(W')), 'InitialMagnification', 'fit')
imshow(W'==0, 'InitialMagnification', 'fit')
ylabel('W^T');
xlabel('Dimension')
print('-dpdf', '-r600', 'LeastrMTFLExp');

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function test_suite = test_minimumCaliperDiameter
%TESTMINIMUMCALIPERDIAMETER  One-line description here, please.
%
%   output = testMinimumCaliperDiameter(input)
%
%   Example
%   testMinimumCaliperDiameter
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-04-14,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testSquare(testCase) %#ok<*DEFNU>

p1 = [10 10];
p2 = [20 10];
p3 = [20 20];
p4 = [10 20];
square = [p1;p2;p3;p4];

width = minimumCaliperDiameter(square);
testCase.assertEqual(10, width);

function testRectangle(testCase)

p1 = [10 10];
p2 = [20 10];
p3 = [20 50];
p4 = [10 50];
square = [p1;p2;p3;p4];

width = minimumCaliperDiameter(square);
testCase.assertEqual(10, width);

function testCross(testCase)

pts = [...
    10 40; ...
    50 40; ...
    50 20; ...
    70 20; ...
    70 40; ...
    110 40; ...
    110 60; ...
    70 60; ...
    70 80; ...
    50 80; ...
    50 60; ...
    10 60];
width = minimumCaliperDiameter(pts);
testCase.assertEqual(60, width);

% try again by shuffling vertices
pts = pts([ 4 10 6 12 2 11 8 1 5 3 9 7], :);
width = minimumCaliperDiameter(pts);
testCase.assertEqual(60, width);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% Exact solution of Riemann problem

% Copyright 2001 P. Wesseling
% This program and its subprograms may be freely used, modified and distributed
% under the GNU General Public License: http://www.gnu.org/copyleft/gpl.html

% Theory is given in Section 10.2 of
% 	Principles of Computational Fluid Dynamics, by P. Wesseling
% 	Springer-Verlag, Berlin etc., 2001. ISBN 3-540-67853-0
% See http://dutita0.twi.tudelft.nl/nw/users/wesseling/

% This program generates Figs. 10.2 -- 10.4 in the book

% Functions called: f, problem_specification

global  PRL  CRL MACHLEFT  gamma  pleft  pright  rholeft  rhoright  uleft...
	uright  tend  lambda		% lambda = dt/dx

problem_specification

gamma = 1.4; gammab = 1/(gamma - 1); gam1 = gamma-1;

% Assumed structure of exact solution
%
%    \         /      |con |       |s|
%     \   f   /       |tact|       |h|
% left \  a  /  state |disc| state |o| right
% state \ n /    2    |cont|   3   |c| state
%   1    \ /          |tinu|       |k|   4
%         |           |ity |       | |

PRL = pright/pleft;
cright = sqrt(gamma*pright/rhoright); cleft = sqrt(gamma*pleft/rholeft);
CRL = cright/cleft;
MACHLEFT = (uleft - uright)/cleft;

p34 = fzero('f',3);		% p34 = p3/p4
p3 = p34*pright; 	alpha = (gamma+1)/(gamma-1);
rho3 = rhoright*(1+alpha*p34)/(alpha+p34);
rho2 = rholeft*(p34*pright/pleft)^(1/gamma);
u2 = uleft-uright+(2/(gamma-1))*cleft*...
	(1-(p34*pright/pleft)^((gamma-1)/(2*gamma)));
c2 = sqrt(gamma*p3/rho2);
spos = 0.5 + ...		% Shock position
	tend*cright*sqrt((gamma-1)/(2*gamma) + (gamma+1)/(2*gamma)*p34)+...
	tend*uright;

conpos = 0.5 + u2*tend + tend*uright;	% Position of contact discontinuity
pos1 = 0.5 + (uleft - cleft)*tend;	% Start of expansion fan
pos2 = 0.5 + (u2+uright-c2)*tend;	% End of expansion fan
xx = 0:0.002:1;
pexact = zeros(size(xx)); uexact= zeros(size(xx)); rhoexact = zeros(size(xx));
machexact = zeros(size(xx));  cexact = zeros(size(xx));
for i = 1:length(xx)
  if xx(i) <= pos1
    pexact(i) = pleft;    rhoexact(i) = rholeft;
    uexact(i) = uleft;    cexact(i)   = sqrt(gamma*pexact(i)/rhoexact(i));
    machexact(i) = uexact(i)/cexact(i);
  elseif xx(i) <= pos2
    pexact(i) = pleft*(1+(pos1-xx(i))/(cleft*alpha*tend))^(2*gamma/(gamma-1));
    rhoexact(i) = rholeft*(1+(pos1-xx(i))/(cleft*alpha*tend))^(2/(gamma-1));
    uexact(i) = uleft + (2/(gamma+1))*(xx(i)-pos1)/tend;
    cexact(i) = sqrt(gamma*pexact(i)/rhoexact(i));
    machexact(i) = uexact(i)/cexact(i);
  elseif xx(i) <= conpos
    pexact(i) = p3;    	      rhoexact(i) = rho2;
    uexact(i) = u2+uright;    cexact(i)   = sqrt(gamma*pexact(i)/rhoexact(i));
    machexact(i) = uexact(i)/cexact(i);
  elseif xx(i) <= spos
    pexact(i) = p3;    rhoexact(i) = rho3;    uexact(i) = u2+uright;
    cexact(i) = sqrt(gamma*pexact(i)/rhoexact(i));
    machexact(i) = uexact(i)/cexact(i);
  else
    pexact(i) = pright;    rhoexact(i) = rhoright;
    uexact(i) = uright;    cexact(i)   = sqrt(gamma*pexact(i)/rhoexact(i));
    machexact(i) = uexact(i)/cexact(i);
  end
end
entroexact = log(pexact./rhoexact.^gamma);

figure(1), clf
subplot(2,3,1), hold on, title('DENSITY','fontsize',12),    plot(xx,rhoexact)
subplot(2,3,2), hold on, title('VELOCITY','fontsize',12),   plot(xx,uexact)
subplot(2,3,3), hold on, title('PRESSURE','fontsize',12),   plot(xx,pexact)
subplot(2,3,4), hold on, title('MACHNUMBER','fontsize',12), plot(xx,machexact)
subplot(2,3,5), hold on, title('ENTROPY','fontsize',12),    plot(xx,entroexact)

repo: Aero-matlab


text:
function [mm] = um2mm(um)
% Convert length from micrometers (or microns) to millimeters.
% Chad A. Greene 2012
mm = um*0.001;


text:
function b = gray2bi( g )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                %
%%   GRAY2BI converts Gray encoded sequence into the binary       %
%%   sequence. It is assumed that the most significant bit is     %
%%   the left hand side bit.                                      %
%%                                                                %
%%   Comments and suggestions to: batlles@gmail.com               %
%%                                                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % copy the msb:
    b(:,1) = g(:,1);

    for i = 2:size(g,2),
        b(:,i) = xor( b(:,i-1), g(:,i) );
    end

return;


text:
function [x1,err] = perform_debiasing(A,x,y, options)

% perform_debiasing - remove bias by orthogonal projection
%
%   x1 = perform_debiasing(A,x,y, options);
%
%   Compute x1 with same support I=find(abs(x)>Thresh) as x that minimize
%       min | y - A(:,I)*x(I) |
%   Thresh is set in options.Thresh
%
%   Usefull to remove some bias after L1 minimization.
%
%   If A is an implicit callback function or if options.use_pinv=0, then it
%   used an itertive gradient descent.
%
%   Copyright (c) 2008 Gabriel Peyre

options.null = 0;
use_pinv = getoptions(options, 'use_pinv', 1);
Thresh = getoptions(options, 'Thresh', 1e-5);
verb = getoptions(options, 'verb', 1);

if isnumeric(A) && use_pinv
    if size(x,2)>1
        % multiple signals
        x1 = x*0;
        for i=1:size(x,2)
            [x1(:,i),err] = perform_debiasing(A,x(:,i),y(:,i), options);
        end
        return;
    end
    % use explicit pseudo inverse
    I = find(abs(x)>Thresh);
    x1 = x*0;
    x1(I) = A(:,I)\y;
    err = [];
    return;
end

niter = getoptions(options, 'niter_debiasing', 100);
tau = getoptions(options, 'tau', 1);
x1 = getoptions(options, 'xguess', x);

err = [];
x1 = perform_support_projection(x1,x,Thresh);
for i=1:niter
    if verb
        progressbar(i,niter);
    end
    % gradient descent
    % x <- x + 1/tau * A'*( y-A*x )
    r   = y - applyop_single( A, x1,+1,options  );
    err(i) = norm(r, 'fro');
    if not(iscell(x))
        x1  = x1 + 1/tau * applyop_single( A, r,-1,options  );
    else
        xx = applyop_single( A, r,-1,options  );
        for k=1:length(x)
            x1{k}  = x1{k} + 1/tau * xx{k};
        end
    end
    % project on support
    x1 = perform_support_projection(x1,x,Thresh);
end

%%
function x = perform_support_projection(x,x0,Thresh)

if iscell(x)
    for i=1:length(x)
        x{i} = perform_support_projection(x{i},x0{i},Thresh);
    end
    return;
end
% impose same support
x(abs(x0)<Thresh) = 0;

%%
function y = applyop_single(A,x,dir,options, pA)

if isempty(A)
    y = x;
    return;
end
if isnumeric(A)
    if dir==1
        y = A*x;
    elseif dir==-1
        y = A'*x;
    else
        if exist('pA') && not(isempty(pA))
            y = pA*x;
        else
            y = pinv(A)*x;
        end
    end
else
    y = feval( A,  x, dir, options );
end

author: gpeyre

repo: matlab-toolboxes


text:
%% t_gifti_1
%

% Download from RDT
%{
 fullFolderName = fullfile(vistaRootPath,'local');
 rdt = RdtClient('vistasoft');

 rdt.crp('/vistadata/anatomy/anatomyNIFTI');
 leftMeshFile = rdt.readArtifact('leftMesh',...
    'type','mat',...
    'destinationFolder',fullFolderName);

 rdt.crp('/vistadata/gifti/BV_GIFTI/Base64');
 gFile = rdt.readArtifact('sujet01_Lwhite.surf',...
    'type','gii',...
    'destinationFolder',fullFolderName);

%}

%%
chdir(fullfile(vistaRootPath,'local'));

%% Converting  gifti to mrMesh (see above for the other direction)

mshFile = fullfile(vistaRootPath,'local','leftMesh.mat');
load(mshFile)
meshVisualize(msh);

% Load a GIFTI file
gFile = fullfile(vistaRootPath,'local','sujet01_Lwhite.surf.gii');
g = gifti(gFile);
% Blue shaded
% mrvNewGraphWin; plot(g);

% Seen with mrMesh - shading not yet right, but it comes up
msh = meshCreate;
msh.triangles = double(g.faces' - 1);
msh.vertices  = double(g.vertices');
% msh = meshColor(msh);

c = ones(3,size(msh.vertices,2))*120;
c(4,size(msh.vertices,2)) = 255;
msh.colors = c;
meshVisualize(msh);

% Can't smooth without the normals.  Perhaps we can create these using
% isonormals().  Also, we might be able to effectively smooth using the
% Matlab function reducepatch.  I don't know how to meshColor.
%
% msh2 = meshSmooth(msh);
% meshVisualize(msh);
% Also doesn't work just yet.
% foo = mrMeshCurvature(msh);

%% End


text:
function gray = ip2volTSeries(inplane,gray,selectedScans,method)
%
% function gray = ip2volTSeries(inplane,gray,[selectedScans],[method])
%
% Uses point sampling and nearest neighbor interpolation to map
% tSeries from inplane view to gray view. Inplane and
% gray views must already be open. Loads the inplane tSeries as
% it goes.
%
% Output tSeries matrices are (as usual) nFrames x nVoxels in size
% where nVoxels is the number of gray voxels that correspond to the
% inplanes, i.e., size(gray.grayCoords,2).
%
% selectedScans:
%   0 - do all scans
%   number or list of numbers - do only those scans
%   default - prompt user via selectScans dialog
% %
% method : 'nearest' [default], 'linear' interpolation
%
% If you change this function make parallel changes in:
%    ip2volCorAnal, ip2volParMap, ip2volSpatialGradient,
%    vol2flatCorAnal, vol2flatParMap, vol2flatTSeries
%
% djh, 2/2001
% ras, 10/2005, fixed to agree with a concomittant change in
% upSampleFactor.
% sod, 11/2005 added linear interpolation

% Don't do this unless inplane is really an inplane and gray is really a gray
if ~strcmp(inplane.viewType,'Inplane') || ~(strcmp(gray.viewType,'Gray') || strcmp(gray.viewType,'Volume'))
    myErrorDlg('ip2grayTSeries can only be used to transform from inplane to gray.');
end

% Check that both gray & flat are properly initialized
if isempty(inplane)
    myErrorDlg('Inplane view must be open.  Use "Open Inplane Window" from the Window menu.');
end
if isempty(gray)
    myErrorDlg('Gray view must be open.  Use "Open Gray Window" from the Window menu.');
end

% check for compatible data types
checkTypes(inplane,gray);

nScans = viewGet(gray, 'numScans');

% (Re-)set scanList
if ~exist('selectedScans','var') || isempty(selectedScans),
    selectedScans = er_selectScans(inplane);
elseif selectedScans == 0
    selectedScans = 1:nScans;
end
if isempty(selectedScans)
    disp('Analysis aborted')
    return
end

if nargin < 4,
    method = 'nearest';
end;
fprintf('[%s]: using %s interpolation.\n',mfilename,method);

% Size of the output tSeries matrices is: nFrames x nVoxels
% Also need to know number of inplane slices
nVoxels = size(gray.coords,2);

% open mrvWaitbar
verbose = prefsVerboseCheck;
if verbose,
    waitHandle = mrvWaitbar(0, 'Interpolating tSeries.  Please wait...');
end

% Tranform gray coords to inplane functional coords. Previously, the code
% to do this xform was duplicated in many functions, including this one.
% It is now a separate routine. The third argument when set to true returns
% the precise (non-integer) functional coords, which are interpolated
% below.
ipCoords = ip2volXformCoords(gray, inplane, true);

% Loop through the scans
for scan = selectedScans

    % Scale and round the grayCoords
    fprintf('Xforming scan %i ...\n',scan);

    % only round for nearest neighbor interpolation
    switch method,
        case 'nearest', ipCoords=round(ipCoords);
        case 'linear'
        otherwise,
            fprintf('Unknown interpolation method: %s\n',method);
            return
    end

    nFrames = viewGet(gray,'numFrames',scan);

    % Reset to NaNs
    tSeries = repmat(single(NaN), [nFrames nVoxels]);

    [~,nii] = loadtSeries(inplane,scan);
    funcData = niftiGet(nii, 'data');

    for frame = 1:nFrames
        subData = double(funcData(:,:,:,frame));
        tSeries(frame,:) = interp3(subData, ...
            ipCoords(2,:), ...
            ipCoords(1,:), ...
            ipCoords(3,:), ...
            method);
    end

    % Save tSeries
    % This should not need to be changed to new version because it is only
    % for the gray view
    savetSeries(tSeries, gray, scan, 1);

    % update the mrvWaitbar
    if verbose,
        mrvWaitbar(find(selectedScans==scan)/nScans, waitHandle);
    end
end %for

% close mrvWaitbar
if verbose, close(waitHandle); end

fprintf('Done xforming tSeries.\n');

return


text:
function [ a, ipvt, rcond ] = r8mat_geco ( a, n )

%*****************************************************************************80
%
%% R8MAT_GECO factors a real matrix and estimates its condition number.
%
%  Discussion:
%
%    For the system A * X = B, relative perturbations in A and B
%    of size EPSILON may cause relative perturbations in X of size
%    EPSILON/RCOND.
%
%    If RCOND is so small that the logical expression
%      1.0 + RCOND == 1.0
%    is true, then A may be singular to working precision.  In particular,
%    RCOND is zero if exact singularity is detected or the estimate
%    underflows.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2005
%
%  Author:
%
%    Original FORTRAN77 version by Cleve Moler.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, real A(N,N), a matrix to be factored.
%
%    Input, integer N, the order of the matrix A.
%
%    Output, real A(LDA,N), the LU factorization of the matrix.
%
%    Output, integer IPVT(N), the pivot indices.
%
%    Output, real RCOND, an estimate of the reciprocal
%    condition number of A.
%

%
%  Compute the L1 norm of A.
%
  anorm = 0.0;
  for j = 1 : n
    anorm = max ( anorm, sum ( abs ( a(1:n,j) ) ) );
  end
%
%  Compute the LU factorization.
%
  [ a, ipvt, info ] = r8mat_gefa ( a, n );
%
%  RCOND = 1 / ( norm(A) * (estimate of norm(inverse(A))) )
%
%  estimate of norm(inverse(A)) = norm(Z) / norm(Y)
%
%  where
%    A * Z = Y
%  and
%    A' * Y = E
%
%  The components of E are chosen to cause maximum local growth in the
%  elements of W, where U'*W = E.  The vectors are frequently rescaled
%  to avoid overflow.
%
%  Solve U' * W = E.
%
  ek = 1.0;
  z(1:n) = 0.0;

  for k = 1 : n

    if ( z(k) ~= 0.0 )
      ek = - abs ( ek ) * r8_sign ( z(k) );
    end

    if ( abs ( a(k,k) ) < abs ( ek - z(k) ) )
      s = abs ( a(k,k) ) / abs ( ek - z(k) );
      z(1:n) = s * z(1:n);
      ek = s * ek;
    end

    wk = ek - z(k);
    wkm = -ek - z(k);
    s = abs ( wk );
    sm = abs ( wkm );

    if ( a(k,k) ~= 0.0 )
      wk = wk / a(k,k);
      wkm = wkm / a(k,k);
    else
      wk = 1.0;
      wkm = 1.0;
    end

    if ( k+1 <= n )

      for j = k+1 : n
        sm = sm + abs ( z(j) + wkm * a(k,j) );
        z(j) = z(j) + wk * a(k,j);
        s = s + abs ( z(j) );
      end

      if ( s < sm )
        t = wkm - wk;
        wk = wkm;
        z(k+1:n) = z(k+1:n) + t * a(k,k+1:n);
      end

    end

    z(k) = wk;

  end

  z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) );
%
%  Solve L' * Y = W
%
  for k = n : -1 : 1

    z(k) = z(k) + z(k+1:n) * a(k+1:n,k);

    if ( 1.0 < abs ( z(k) ) )
      z(1:n) = z(1:n) / abs ( z(k) );
    end

    l = ipvt(k);

    t = z(l);
    z(l) = z(k);
    z(k) = t;

  end

  z(1:n) = z(1:n) / sum ( abs ( z(1:n) ) );

  ynorm = 1.0;
%
%  Solve L * V = Y.
%
  for k = 1 : n

    l = ipvt(k);

    t = z(l);
    z(l) = z(k);
    z(k) = t;

    z(k+1:n) = z(k+1:n) + t * a(k+1:n,k)';

    if ( 1.0 < abs ( z(k) ) )
      ynorm = ynorm / abs ( z(k) );
      z(1:n) = z(1:n) / abs ( z(k) );
    end

  end

  s = sum ( abs ( z(1:n) ) );
  z(1:n) = z(1:n) / s;
  ynorm = ynorm / s;
%
%  Solve U * Z = V.
%
  for k = n : -1 : 1

    if ( abs ( a(k,k) ) < abs ( z(k) ) )
      s = abs ( a(k,k) ) / abs ( z(k) );
      z(1:n) = s * z(1:n);
      ynorm = s * ynorm;
    end

    if ( a(k,k) ~= 0.0 )
      z(k) = z(k) / a(k,k);
    else
      z(k) = 1.0;
    end

    z(1:k-1) = z(1:k-1) - z(k) * a(1:k-1,k)';

  end
%
%  Normalize Z in the L1 norm.
%
  s = 1.0 / sum ( abs ( z(1:n) ) );
  z(1:n) = s * z(1:n);
  ynorm = s * ynorm;

  if ( anorm ~= 0.0 )
    rcond = ynorm / anorm;
  else
    rcond = 0.0;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [pt degen] = line_intersect(p1, p2, p3, p4)
% function pt = line_intersect(p1, p2, p3, p4)
% intersection point of line defined by p1 & p2 and line defined by p3 & p4
% http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/
x1 = p1(1); y1 = p1(2);
x2 = p2(1); y2 = p2(2);
x3 = p3(1); y3 = p3(2);
x4 = p4(1); y4 = p4(2);

if (y4-y3)*(x2-x1)-(x4-x3)*(y2-y1) == 0
% 	warning('line_intersect.m degenerate --dclee');
% 	pt = (p1 + p2 + p3 + p4)/4;
    pt = [];
    degen = 1;
	return;
end

pt = [ ...
	x1 + (x2-x1) * ((x4-x3)*(y1-y3)-(y4-y3)*(x1-x3))/((y4-y3)*(x2-x1)-(x4-x3)*(y2-y1)) ...
	y1 + (y2-y1) * ((x4-x3)*(y1-y3)-(y4-y3)*(x1-x3))/((y4-y3)*(x2-x1)-(x4-x3)*(y2-y1)) ];
degen = 0;

author: zouchuhang

repo: LayoutNet


text:
function CPDpot = convert_dbn_CPDs_to_pots(bnet, evidence, pot_type, softCPDpot)
% CONVERT_DBN_CPDS_TO_POTS Convert CPDs of (possibly instantiated) DBN nodes to potentials
% CPDpot = convert_dbn_CPDs_to_pots(bnet, evidence, pot_type, softCPDpot)
%
% CPDpot{n,t} is a potential containing P(n,t|pa(n,t), ev)
% softCPDpot{n,t} is a potential containing P(n,t|pa(n,t), ev) insted of using n's CPD

[ss T] = size(evidence);

if nargin < 4, softCPDpot = cell(ss,T); end
CPDpot = softCPDpot;

% Convert CPDs of instantiated nodes to potential form
t = 1;
for n=1:ss
  fam = family(bnet.dag, n);
  e = bnet.equiv_class(n, 1);
  if isempty(softCPDpot{n,t})
    CPDpot{n,t} = convert_to_pot(bnet.CPD{e}, pot_type, fam(:), evidence(:,1));
  end
end
for n=1:ss
  fam = family(bnet.dag, n, 2);
  e = bnet.equiv_class(n, 2);
  for t=2:T
    if isempty(softCPDpot{n,t})
      CPDpot{n,t} = convert_to_pot(bnet.CPD{e}, pot_type, fam(:), evidence(:,t-1:t));
    end
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute joint speed according to speeddata variable
% vmax means maximum speed for each joint
% the other cases are scaled according to a maximum value of v6000
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function vel=obtain_joint_speed(robot, speeddata)

if strncmp(speeddata, 'vmax',4)
    vel=robot.velmax;
else
    [tag,remain] = strtok(speeddata, 'v');
    vel = robot.velmax.*str2num(tag)/6000;
end

repo: ARTE


text:
function [w, dwdx, dwdxx] = Weight(type, para, di, dmI)
% EVALUATE WEIGHT FUNCTION
%
% SYNTAX: [w, dwdr, dwdrr] = GaussWeight(type, para, di, dmi)
%
% INPUT PARAMETERS
%    type - Type of weight function
%    para - Weight function parameter
%    di   - Distance
%    dmi  - Support size
% OUTPUT PARAMETERS
%    w    - Value of weight function at r
%    dwdx - Value of first order derivative of weight function with respect to x at r
%    dwdxx- Value of Second order derivative of weight function with respect to x at r
%
r  = abs(di) / dmI;

if (di >= 0.0)
   drdx = 1.0/dmI;
else
   drdx = -1.0/dmI;
end

% EVALUATE WEIGHT FUNCTION AND ITS FIRST AND SECOND ORDER OF DERIVATIVES WITH RESPECT r AT r
if (type == 'GAUSS')
   [w,dwdr,dwdrr] = Gauss(para,r);
elseif (type == 'CUBIC')
   [w,dwdr,dwdrr] = Cubic(r);
elseif (type == 'SPLIN')
   [w,dwdr,dwdrr] = Spline(r);
elseif (type == 'power')
   [w,dwdr,dwdrr] = power_function(para,r);
elseif (type == 'CSRBF')
   [w,dwdr,dwdrr] = CSRBF2(r);
else
   error('Invalid type of weight function.');
end

dwdx  = dwdr * drdx;
dwdxx = dwdrr * drdx * drdx;

function [w,dwdr,dwdrr] = Gauss(beta,r)
if (r>1.0)
   w     = 0.0;
   dwdr  = 0.0;
   dwdrr = 0.0;
else
   b2 = beta*beta;
   r2 = r*r;
   eb2 = exp(-b2);

   w     = (exp(-b2*r2) - eb2) / (1.0 - eb2);
   dwdr  = -2*b2*r*exp(-b2*r2) / (1.0 - eb2);
   dwdrr = -2*b2*exp(-b2*r2)*(1-2*b2*r2) / (1.0 - eb2);
end

function [w,dwdr,dwdrr] = Cubic(r)
if (r>1.0)
   w     = 0.0;
   dwdr  = 0.0;
   dwdrr = 0.0;
else
   w     = 1-6*r^2+8*r^3-3*r^4;
   dwdr  = -12*r+24*r^2-12*r^3;
   dwdrr = -12+48*r-36*r^2;
end

function [w,dwdr,dwdrr] = power_function(arfa,r)
if (r>1.0)
   w     = 0.0;
   dwdr  = 0.0;
   dwdrr = 0.0;
else
    a2 = arfa*arfa;
   r2 = r*r;
    w     = exp(-r2/a2);
   dwdr  = (-2*r/a2)*exp(-r2/a2);
   dwdrr  = (-2/a2+(-2*r/a2).^2)*exp(-r2/a2);
end

function [w,dwdr,dwdrr] = Spline(r)
if (r>1.0)
   w     = 0.0;
   dwdr  = 0.0;
   dwdrr = 0.0;
elseif (r<=0.5)
   w     = 2/3 - 4*r^2 + 4*r^3;
   dwdr  = -8*r + 12*r^2;
   dwdrr = -8 + 24*r;
else
   w     = 4/3 - 4*r + 4*r^2 - 4*r^3/3;
   dwdr  = -4 + 8*r -4*r^2;
   dwdrr = 8 - 8*r;
end

function [w,dwdr,dwdrr] = CSRBF2(r)
if (r>1.0)
   w     = 0.0;
   dwdr  = 0.0;
   dwdrr = 0.0;
else
	w     = (1-r)^6*(6+36*r+82*r^2+72*r^3+30*r^4+5*r^5);
	dwdr  = 11*r*(r+2)*(5*r^3+15*r^2+18*r+4)*(r-1)^5;
   dwdrr = 22*(25*r^5+100*r^4+142*r^3+68*r^2-16*r-4)*(r-1)^4;
end


text:
 function [xs, info] = pwls_sps_os(x, yi, wi, Ab, R, niter, pixmax, ...
		denom, aai, relax0, chat)
%function [xs, info] = pwls_sps_os(x, yi, wi, Ab, R, niter, pixmax, ...
%|		denom, aai, relax0, chat)
%|
%| penalized weighted least squares estimation/reconstruction
%| using separable paraboloidal surrogates algorithm with
%| (optionally relaxed) ordered subsets.  (relaxation ensures convergence.)
%|
%| cost(x) = (y-Gx)' W (y-Gx) / 2 + R(x)
%|
%| in
%|	x	[np 1]		initial estimate
%|	yi	[nb na]		measurements (noisy sinogram)
%|	wi	[nb na]		weighting sinogram (or [] for uniform)
%|	Ab	[nd np]		Gblock object, aij >= 0 required!
%|					or sparse matrix (implies nsubset=1)
%|	R			penalty object (see Reg1.m)
%|	niter			# of iterations (including 0)
%|
%| optional
%|	pixmax	[1] or [2]	max pixel value, or [min max] (default [0 inf])
%|	denom	[np 1]		precomputed denominator
%|	aai	[nb na]		precomputed row sums of |Ab|
%|	relax0	[1] or [2]	relax0 or (relax0, relax_rate)
%|	chat
%|
%| out
%|	xs	[np niter]	iterates
%|	info	[niter 1]	time
%|
%| Copyright 2002-2-12, Jeff Fessler, University of Michigan

if nargin < 4, ir_usage, end

cpu etic
info = zeros(niter,1);

Ab = block_op(Ab, 'ensure'); % make it a block object (if not already)
nblock = block_op(Ab, 'n');
starts = subset_start(nblock);

if ~isvar('niter')	|| isempty(niter),	niter = 1;	end
if ~isvar('pixmax')	|| isempty(pixmax),	pixmax = inf;	end
if ~isvar('chat')	|| isempty(chat),	chat = false;	end
if isempty(wi)
	wi = ones(size(yi));
end
if ~isvar('aai') || isempty(aai)
	aai = reshape(sum(Ab'), size(yi)); % a_i = sum_j |a_ij|
					% requires real a_ij and a_ij >= 0
end

if ~isvar('relax0') || isempty(relax0)
	relax0 = 1;
end
if length(relax0) == 1
	relax_rate = 0;
elseif length(relax0) == 2
	relax_rate = relax0(2);
	relax0 = relax0(1);
else
	error relax
end

if length(pixmax) == 2
	pixmin = pixmax(1);
	pixmax = pixmax(2);
elseif length(pixmax) == 1
	pixmin = 0;
else
	error pixmax
end

%
% likelihood denom, if not provided
%
if ~isvar('denom') || isempty(denom)
	denom = Ab' * col(aai .* wi);	% requires real a_ij and a_ij >= 0
end, clear aai

if ~isvar('R') || isempty(R)
	pgrad = 0;		% unregularized default
	Rdenom = 0;
end

[nb na] = size(yi);

%
% loop over iterations
%

xs = zeros(numel(x), niter, class(x));
x = max(x,pixmin);
x = min(x,pixmax);
xs(:,1) = x;

for iter = 2:niter
	ticker(mfilename, iter, niter)

	relax = relax0 / (1 + relax_rate * (iter-2));

	%
	% loop over subsets
	%
	for iset=1:nblock
		iblock = starts(iset);
		ia = iblock:nblock:na;

		li = Ab{iblock} * x;
		li = reshape(li, nb, length(ia));
		resid = wi(:,ia) .* (yi(:,ia) - li);
		grad = Ab{iblock}' * resid(:); % G' * W * (y - G*x)

		if ~isempty(R)
			pgrad = R.cgrad(R, x);
			Rdenom = R.denom(R, x);
		end

		num = nblock * grad - pgrad;
		den = denom + Rdenom;

		x = x + relax * num ./ den;	% relaxed update
		x = max(x,pixmin);		% lower bound
		x = min(x,pixmax);		% upper bound
	end

	if chat, printm('range(x) = %g %g', min(x), max(x)), end
	xs(:,iter) = x;
	info(iter,1) = cpu('etoc');
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [J, grad] = costFunction(theta, X, y)
%COSTFUNCTION Compute cost and gradient for logistic regression
%   J = COSTFUNCTION(theta, X, y) computes the cost of using theta as the
%   parameter for logistic regression and the gradient of the cost
%   w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Note: grad should have the same dimensions as theta
%

H_theta = sigmoid(X * theta);

J = -1/m * sum(y .* log(H_theta) + (1 - y) .* log(1 - H_theta));

grad = 1/m .* X' * (H_theta - y);

% =============================================================

end

author: rmarquis

repo: coursera-machinelearning


text:
% example: edge-preserving smoothing
% figure 1 in our paper

close all;

I = double(imread('.\img_smoothing\cat.bmp')) / 255;
p = I;
r = 4; % try r=2, 4, or 8
eps = 0.2^2; % try eps=0.1^2, 0.2^2, 0.4^2

q = guidedfilter(I, p, r, eps);

figure();
imshow([I, q], [0, 1]);

author: drakeguan


text:
function [s] = ADEM_sample_image(V,o,R)
% samples a (memory mapped) image at displacement o
% FORMAT [s] = ADEM_sample_image(V,o,R)
% FORMAT [s] = ADEM_sample_image(o,h)
%
% V - a structure array containing image volume information
% o - coordinates of foveal sampling:
%   o(1) - oculomotor angle
%   o(2) - oculomotor angle
% R - retinal modulation (n x n)
%
% or
%
% o - coordinates of foveal sampling
% h - vector of coefficients weighting images in STIM.H{:}
%
% s - sensory sample (n x n)
%
% requires a global variable with the following fields:
% STIM.R = contrast modulation matrix that defines resolution
% STIM.W = width of foveal sampling of an image   (default: 1/6)
% STIM.P = image position in retinal  coordinates (default: [0;0])
% STIM.B = basis functions or receptive fields    (default: 1)
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: ADEM_sample_image.m 6932 2016-11-16 12:11:01Z karl $

% retinotopic predictions
%--------------------------------------------------------------------------
global STIM
if nargin < 3
    s     = 0;
    for i = 1:numel(o)
        s = s + o(i)*ADEM_sample_image(STIM.H{i},V,STIM.R);
    end
    return
end

% preliminaries
%--------------------------------------------------------------------------
if ~isfield(STIM,'R'), STIM.R = ones(64,64); end
if ~isfield(STIM,'W'), STIM.W = 1/6;         end
if ~isfield(STIM,'P'), STIM.P = [0;0];       end
if ~isfield(STIM,'B'), STIM.B = 1;           end
if ~isfield(STIM,'A'), STIM.A = 512;         end

% retinotopic sampling
%--------------------------------------------------------------------------
dim = size(R);
dx  = V.dim(1)/dim(1)*STIM.W;

i   = dx*((1:dim(1)) - dim(1)/2) + V.dim(1)/2  + (o(1) + STIM.P(1))*16;
j   = dx*((1:dim(2)) - dim(2)/2) + V.dim(2)/2  + (o(2) + STIM.P(2))*16;
x   = kron(ones(1,dim(2)),i);
y   = kron(j,ones(1,dim(1)));
z   = ones(1,dim(1)*dim(2));

x   = min(max(x,1),V.dim(1));
y   = min(max(y,1),V.dim(2));

s   = spm_sample_vol(V,x,y,z,-2); s(~s) = 1;
s   = reshape(s,dim(1),dim(2)).*R;
s   = STIM.B'*s*STIM.B;

% eccentricity attenuation
%--------------------------------------------------------------------------
s   = s*exp(-o'*o/STIM.A);

author: spm


text:
function double_pendulum(ivp, duration, fps, movie)
% DOUBLE_PENDULUM Animates the double pendulum's (mostly) chaotic behavior.
%
%   author:  Alexander Erlich (alexander.erlich@gmail.com)
%
%   parameters:
%
%   ivp=[phi1; dtphi1; phi2; dtphi2; g; m1; m2; l1; l2]
%
%                               Initial value problem. phi1 and dtphi1 are
%                               the initial angle and anglular velocity. g
%                               is gravity, m1 and l1 mass and rod length.
%                               For an explaining picture, see
%                               documentation file in same folder.
%
%   duration                    The time interval on which the ode is
%                               solved spans from 0 to duration (in sec).
%
%   fps                         Frames Per Second. The framerate is
%                               relevant both for normal (realtime)
%                               animation and movie recording.
%
%   movie                       If false, a normal realtime animation of
%                               the motion of the double pendulum (the
%                               framerate being fps) is shown.
%                               If true, a movie (.avi) is recorded. The
%                               filename is 'doublePendulumAnimation.avi'
%                               and the folder into which it is saved is
%                               the current working directory.
%
%   This function calls double_pendulum_ODE and is, in turn, called by
%   double_pendulum_init.
%
%   Example call:    >> double_pendulum([pi;0;pi;5;9.81;1;1;2;1],100,10,false)
%   Or, simply call  >> double_pendulum_init
%
%   ---------------------------------------------------------------------

clear All; clf;

nframes=duration*fps;
sol=ode45(@double_pendulum_ODE,[0 duration], ivp);
t = linspace(0,duration,nframes);
y=deval(sol,t);

phi1=y(1,:)'; dtphi1=y(2,:)';
phi2=y(3,:)'; dtphi2=y(4,:)';
l1=ivp(8); l2=ivp(9);
% phi1=x(:,1); dtphi1=x(:,2);
% phi2=x(:,3); dtphi2=x(:,4);
% l1=ivp(8); l2=ivp(9);

h=plot(0,0,'MarkerSize',30,'Marker','.','LineWidth',2);
range=1.1*(l1+l2); axis([-range range -range range]); axis square;
set(gca,'nextplot','replacechildren');

    for i=1:length(phi1)-1
        if (ishandle(h)==1)
            Xcoord=[0,l1*sin(phi1(i)),l1*sin(phi1(i))+l2*sin(phi2(i))];
            Ycoord=[0,-l1*cos(phi1(i)),-l1*cos(phi1(i))-l2*cos(phi2(i))];
            set(h,'XData',Xcoord,'YData',Ycoord);
            drawnow;
            F(i) = getframe;
            if movie==false
                pause(t(i+1)-t(i));
            end
        end
    end
    if movie==true
        movie2avi(F,'doublePendulumAnimation.avi','compression','Cinepak','fps',fps)
    end


text:
function [dist, t] = distancePointEdge3d(point, edge)
%DISTANCEPOINTEDGE3D Minimum distance between a 3D point and a 3D edge
%
%   DIST = distancePointEdge3d(POINT, EDGE);
%   Return the euclidean distance between edge EDGE and point POINT.
%   EDGE has the form: [x1 y1 z1 x2 y2 z2], and POINT is [x y z].
%
%   If EDGE is N-by-6 array, result is N-by-1 array computed for each edge.
%   If POINT is a N-by-3 array, the result is computed for each point.
%   If both POINT and EDGE are array, they must have the same number of
%   rows, and the result is computed for each couple point(i,:);edge(i,:).
%
%   [DIST POS] = distancePointEdge3d(POINT, EDGE);
%   Also returns the position of closest point on the edge. POS is
%   comprised between 0 (first point) and 1 (last point).
%
%   See also:
%   edges3d, points3d, distancePoints3d, distancePointLine3d
%

%   ---------
%   author : David Legland
%   INRA - CEPIA URPOI - MIA MathCell
%   created the 07/04/2004.
%

%   HISTORY
%   2005-06-24 rename, and change arguments sequence
%   2009-04-30 add possibility to return position of closest point
%   2011-04-14 add checkup for degenerate edges, improve speed, update doc

% direction vector of each edge
vl = edge(:, 4:6) - edge(:, 1:3);

% compute position of points projected on the supporting line
% (Size of t is the max number of edges or points)
t = linePosition3d(point, [edge(:,1:3) vl]);

% change position to ensure projected point is located on the edge
t(t < 0) = 0;
t(t > 1) = 1;

% difference of coordinates between projected point and base point
p0 = bsxfun(@plus, edge(:,1:3), [t .* vl(:,1) t .* vl(:,2) t .* vl(:,3)]);
p0 = bsxfun(@minus, point, p0);

% compute distance between point and its projection on the edge
dist = sqrt(sum(p0 .* p0, 2));

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
%% mri_sense_demo1.m
%|
%| Example illustrating regularized iterative reconstruction for parallel MRI
%| (sensitivity encoding imaging or SENSE), from nonuniform k-space samples.
%| (This example does not include field inhomogeneity or relaxation.)
%|
%| Copyright 2006-4-18, Jeff Fessler, University of Michigan

%% make coil sensitivity maps
if ~isvar('smap'), printm 'sense maps'
	ig = image_geom('nx', 64, 'fov', 250); % 250 mm FOV
	ig.mask = ig.circ(ig.dx * (ig.nx/2-2), ig.dy * (ig.ny/2-1)) > 0;
	f.ncoil = 4;
	smap = mri_sensemap_sim('nx', ig.nx, 'ny', ig.ny, 'dx', ig.dx, ...
		'rcoil', 120, 'ncoil', f.ncoil, 'chat', 1);
prompt
end

%% true object (discrete because of discrete sense map)
if ~isvar('xtrue'), printm 'true object'
	xtrue = ellipse_im(ig, 'shepplogan-emis', 'oversample', 2);
	clim = [0 8];
	im plc 2 2
	im(1, ig.x, ig.y, xtrue, '\x true', clim), cbar
	im(4, ig.mask, 'mask')
prompt
end

%% trajectory and system matrix
if ~isvar('Am'), printm 'system matrix A'
	f.traj = 'spiral1';
	f.dens = 'voronoi';

	N = [ig.nx ig.ny];
	[kspace, omega, wi_traj] = mri_trajectory(f.traj, {}, ...
		N, ig.fov, {f.dens});

	% create Gnufft class object
	J = [6 6];
	nufft_args = {N, J, 2*N, N/2, 'table', 2^10, 'minmax:kb'};
	Am = Gmri(kspace, ig.mask, ...
		'fov', ig.fov, 'basis', {'rect'}, 'nufft', nufft_args);

	if im
		im subplot 2
		plot(omega(1:5:end,1), omega(1:5:end,2), '.')
		axis_pipi, axis square
		titlef('%s: %d', f.traj, size(omega,1))
	end
prompt
end

wi_basis = wi_traj ./ Am.arg.basis.transform;

%% SENSE system object
if ~isvar('Ab'), printm 'Ab object with sense maps within'
	for ic=1:f.ncoil
		tmp = smap(:,:,ic);
		tmp = Gdiag(tmp(ig.mask), 'mask', ig.mask);
		Ac{ic} = Am * tmp; % cascade
	end
	Ab = block_fatrix(Ac, 'type', 'col'); % [A1; A2; ... ]
end

%% noisy data
if ~isvar('yi'), printm 'data yi'
	ytrue = Ab * xtrue(ig.mask);

	% add noise
	rng(0)
	yi = ytrue + 0 * randn(size(ytrue));
	% todo: visualize data...
end

%% CP recon
if ~isvar('xcp'), printm 'conj. phase reconstruction'
	xcp = zeros(ig.nx, ig.ny, f.ncoil, 'single');
	y4 = reshape(yi, [], f.ncoil);
	for ic=1:f.ncoil
		tmp = Am' * (wi_basis .* y4(:,ic));
		xcp(:,:,ic) = ig.embed(tmp);
	end

	im(2, abs(xcp), 'Conj. Phase Recon for each coil'), cbar
prompt
end

%% SSoS
if ~isvar('xssos'), printm 'sqrt sum-of-squares reconstruction'
	xssos = sqrt(sum(abs(xcp).^2, 3));
	xssos = ir_wls_init_scale(1, xtrue, xssos); % cheat trick for scaling
	im(3, abs(xssos), 'sum-of-squares recon'), cbar
	xlabelf('NRMSE %.1f%%', 100*nrms(xssos(:), xtrue(:)))
prompt
end

%% regularizer
if ~isvar('R'), printm 'regularizer'
	f.beta = 2^11;
	R = Reg1(ig.mask, 'beta', f.beta, 'type_penal', 'mat'); % complex
	if 1
		psf = qpwls_psf(Ab, R, 1, ig.mask, 1, 'offset', [0 0]);
		im(4, psf)
	end
prompt
end

%% PCG
if ~isvar('xpcg'), printm 'PCG with quadratic penalty'
	f.niter = 10;
	xpcg = qpwls_pcg(xssos(ig.mask), Ab, 1, yi(:), 0, R.C, 1, f.niter);
	xpcg = ig.embed(xpcg(:,end)); % convert last vector to image for display

	im(4, abs(xpcg), '|\x| PCG quad', clim), cbar
	xlabelf('NRMSE %.1f%%', 100*nrms(xpcg(:), xtrue(:)))
end