function Spath = Simulate_SLV_func( N_sim, M, T, S_0, r, q, model, params)
% Simulates Paths of Stochastic Local Volatility Models (basic Euler Scheme for most cases)
% N_sim = # paths
% M = #time steps on [0,T], ie dt =T/M
% Note: returns paths of dimension (N_sim,M+1), since they include S_0
%
% SVLModel:                         (parameters)
%        1 = Heston:                alpha, v0, rho, theta, eta
%        2 = SABR:                  alpha, v0, rho, beta
%        3 = Shifted SABR:          alpha, v0, rho, beta, shift
%        4 = Quadratic SLV:         alpha, v0, rho, theta, eta, a, b, c
%        5 = TanHyp-Heston:         alpha, v0, rho, theta, eta, beta
%        6 = Heston-SABR:           alpha, v0, rho, theta, eta, beta
%
%==============================
% Initialize Common Params/Vectors
%==============================
alpha  = params.alpha;   % (vol of vol)
v0     = params.v0;  % initial vol/var
rho    = params.rho; % covar bewteen asset and vol innovations

Sigmav = alpha; %NOTE: alpha is same as Sigmav (vol of vol)

if model == 1 || model == 4 || model == 5 || model == 6
    theta = params.theta;
    eta   = params.eta;
end

Spath = zeros(N_sim,M+1);
Spath(:,1) = S_0;

dt = T/M;
sqdt = sqrt(dt);
sqdtrho1 = sqdt*rho;
sqdtrho2 = sqdt*sqrt(1-rho^2);

Zeta = r - q;
%==============================
% Simulate Based on Specific StochVol Model
%==============================

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if model == 1 %HESTON MODEL

    expEta = exp(-eta*dt);
    driftv = theta*(1 - expEta);  %analytical drift for variance process
    vOld   = v0*ones(N_sim,1); %used to store variance process
    sqvOld = sqrt(v0)*ones(N_sim,1);

    for m = 1:M
        W1 = randn(N_sim,1); W2 = randn(N_sim,1);  %Generate two Brownian motions
        Spath(:,m+1) = Spath(:,m).*exp((Zeta - .5*vOld)*dt + sqdt*sqvOld.*W1);  %log scheme
        %Spath(:,m+1) = Spath(:,m)*( (1+r)*dt + sqrt(vOld)*(sqrho1*W1 + sqrho2*W2));  %level scheme
        vNew = driftv + vOld*expEta + Sigmav*sqvOld.*(sqdtrho1*W1 + sqdtrho2*W2);
        %vNew = (sqvOld +.5*Sigmav*(sqdtrho1*W1 + sqdtrho2*W2)).^2 + eta*(theta - vOld)*dt -.25*Sigmav*dt; %Milstein
        vOld = abs(vNew); %Always stays positive, the "reflection" scheme
        %vOld = max(0,vNew); %least biased scheme
        sqvOld = sqrt(vOld);
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif model == 2  % SABR
    beta  = params.beta;
    cons1 = -.5*(alpha)^2*dt;
    vOld  = v0*ones(N_sim,1); %used to store variance process
    for m = 1:M
        W1 = randn(N_sim,1); W2 = randn(N_sim,1);  %Generate two Brownian motions
        Spath(:,m+1) = max(0, Spath(:,m) + Spath(:,m).^beta.* vOld .*sqdt.*W1);  %level scheme
        vOld = vOld.*exp(cons1 + alpha*(sqdtrho1*W1 + sqdtrho2*W2));
        %vOld = vOld + alpha*vOld.*(sqdtrho1*W1 + sqdtrho2*W2);
    end

%     %%% LOG SCHEME
%     cons2 = 2*(beta - 1);
%     for m = 1:M
%         W1 = randn(N_sim,1); W2 = randn(N_sim,1);  %Generate two Brownian motions
%         Spath(:,m+1) = Spath(:,m).*exp(-.5*vOld.^2.*Spath(:,m).^cons2*dt + Spath(:,m).^(beta-1).* vOld .*sqdt.*W1);  %level scheme
%         vOld = vOld.*exp(cons1 + alpha*(sqdtrho1*W1 + sqdtrho2*W2));
%     end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif model == 3  % SHIFTED SABR
    beta  = params.beta;
    shift = params.shift;
    cons1 = -.5*(alpha)^2*dt;
    vOld  = v0*ones(N_sim,1); %used to store variance process
    for m = 1:M
        W1 = randn(N_sim,1); W2 = randn(N_sim,1);  %Generate two Brownian motions
        %%% Note: because of shift, -shift <= S_t
        Spath(:,m+1) = max(-shift, Spath(:,m) + (Spath(:,m) + shift).^beta.* vOld .*sqdt.*W1);  %level scheme
        vOld = vOld.*exp(cons1 + alpha*(sqdtrho1*W1 + sqdtrho2*W2));
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif model == 4  % QUADRATIC SLV
    a = params.a; b = params.b; c = params.c;

    expEta = exp(-eta*dt);
    driftv = theta*(1 - expEta);  %analytical drift for variance process
    vOld   = v0*ones(N_sim,1); %used to store variance process
    sqvOld = sqrt(v0)*ones(N_sim,1);

    drift = 1+ Zeta*dt;
    for m = 1:M
        W1 = randn(N_sim,1); W2 = randn(N_sim,1);  %Generate two Brownian motions

        Spath(:,m+1) = max(0,Spath(:,m)*drift + sqdt*sqvOld.*W1.*(Spath(:,m).*(a*Spath(:,m) +b) + c));

        %Spath(:,m+1) = Spath(:,m)*( (1+r)*dt + sqrt(vOld)*(sqrho1*W1 + sqrho2*W2));  %level scheme
        vNew = driftv + vOld*expEta + Sigmav*sqvOld.*(sqdtrho1*W1 + sqdtrho2*W2);
        %vNew = (sqvOld +.5*Sigmav*(sqdtrho1*W1 + sqdtrho2*W2)).^2 + eta*(theta - vOld)*dt -.25*Sigmav*dt; %Milstein
        %vOld = abs(vNew); %Always stays positive, the "reflection" scheme
        vOld = max(0,vNew); %least biased scheme
        sqvOld = sqrt(vOld);
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif model == 5  % TAN-HYP-HESTON
    beta = params.beta;

    expEta = exp(-eta*dt);
    driftv = theta*(1 - expEta);  %analytical drift for variance process
    vOld   = v0*ones(N_sim,1); %used to store variance process
    sqvOld = sqrt(v0)*ones(N_sim,1);

    drift = 1+ Zeta*dt;
    for m = 1:M
        W1 = randn(N_sim,1); W2 = randn(N_sim,1);  %Generate two Brownian motions

        Spath(:,m+1) = max(0,Spath(:,m)*drift + sqdt*sqvOld.*W1.*tanh(beta*Spath(:,m)));
        vNew = driftv + vOld*expEta + Sigmav*sqvOld.*(sqdtrho1*W1 + sqdtrho2*W2);
        %vNew = (sqvOld +.5*Sigmav*(sqdtrho1*W1 + sqdtrho2*W2)).^2 + eta*(theta - vOld)*dt -.25*Sigmav*dt; %Milstein
        %vOld = abs(vNew); %Always stays positive, the "reflection" scheme
        vOld = max(0,vNew); %least biased scheme
        sqvOld = sqrt(vOld);
    end

elseif model == 6   % HESTON-SABR
    drift = 1+ Zeta*dt;
    expEta = exp(-eta*dt);
    driftv = theta*(1 - expEta);  %analytical drift for variance process
    vOld   = v0*ones(N_sim,1); %used to store variance process
    sqvOld = sqrt(v0)*ones(N_sim,1);

    beta  = params.beta;

    for m = 1:M
        W1 = randn(N_sim,1); W2 = randn(N_sim,1);  %Generate two Brownian motions

        %%%% LEVEL scheme
        Spath(:,m+1) = max(0, Spath(:,m)*drift + Spath(:,m).^beta.* sqvOld .*sqdt.*W1);

%         %%%% LOG scheme
%         Spath(:,m+1) = Spath(:,m).*exp(-.5*vOld.^2.*Spath(:,m).^cons2*dt + Spath(:,m).^(beta-1).* sqvOld .*sqdt.*W1);

        vNew = driftv + vOld*expEta + Sigmav*sqvOld.*(sqdtrho1*W1 + sqdtrho2*W2);
        vOld = abs(vNew); %Always stays positive, the "reflection" scheme
        %vOld = max(0,vNew); %least biased scheme
        sqvOld = sqrt(vOld);
    end

%     beta  = params.beta;
%     cons1 = -.5*(alpha)^2*dt;
%     vOld  = v0*ones(N_sim,1); %used to store variance process
%     for m = 1:M
%         W1 = randn(N_sim,1); W2 = randn(N_sim,1);  %Generate two Brownian motions
%         Spath(:,m+1) = max(0, Spath(:,m) + Spath(:,m).^beta.* vOld .*sqdt.*W1);  %level scheme
%         vOld = vOld.*exp(cons1 + alpha*(sqdtrho1*W1 + sqdtrho2*W2));
%         %vOld = vOld + alpha*vOld.*(sqdtrho1*W1 + sqdtrho2*W2);
%     end
%
%     expEta = exp(-eta*dt);
%     driftv = theta*(1 - expEta);  %analytical drift for variance process
%     vOld   = v0*ones(N_sim,1); %used to store variance process
%     sqvOld = sqrt(v0)*ones(N_sim,1);
%
%     for m = 1:M
%         W1 = randn(N_sim,1); W2 = randn(N_sim,1);  %Generate two Brownian motions
%         Spath(:,m+1) = Spath(:,m).*exp((Zeta - .5*vOld)*dt + sqdt*sqvOld.*W1);  %log scheme
%         %Spath(:,m+1) = Spath(:,m)*( (1+r)*dt + sqrt(vOld)*(sqrho1*W1 + sqrho2*W2));  %level scheme
%         vNew = driftv + vOld*expEta + Sigmav*sqvOld.*(sqdtrho1*W1 + sqdtrho2*W2);
%         %vNew = (sqvOld +.5*Sigmav*(sqdtrho1*W1 + sqdtrho2*W2)).^2 + eta*(theta - vOld)*dt -.25*Sigmav*dt; %Milstein
%         vOld = abs(vNew); %Always stays positive, the "reflection" scheme
%         %vOld = max(0,vNew); %least biased scheme
%         sqvOld = sqrt(vOld);
%     end

end

% %%%% REFORMULATED VERSION
% etahat = eta*theta; thetahat = (eta + Sigmav^2)/etahat; Sigmavhat = -Sigmav; v0hat = 1/v0;
% vOld = v0hat*ones(N_sim,1); %used to store variance process
% sqvOld = sqrt(v0hat)*ones(N_sim,1);
%
% expkt = exp(-etahat*dt);
% driftv = thetahat*(1 - expkt);  %analytical drift for variance process
% sqtSv = sqdt*Sigmavhat;  %analytical drift for variance process
%
% for m = 1:M
%     W1 = randn(N_sim,1); W2 = randn(N_sim,1);  %Generate two Brownian motions
%
%     Spath(:,m+1) = Spath(:,m).*exp((r - .5*vOld)*dt + (sqrho1*W1 + sqrho2*W2)./sqvOld);  %log scheme
%     %Spath(:,m+1) = Spath(:,m)*( (1+r)*dt + sqrt(vOld)*(sqrho1*W1 + sqrho2*W2));  %level scheme
%
%     vNew = driftv + vOld*expkt + sqtSv*sqvOld.*W1;
%
%     vOld = abs(vNew); %Always stays positive, the "reflection" scheme
%     sqvOld = sqrt(vOld);
% end


text:
%% This function is used for moving the endeffector on a circle, for the KUKA iiwa 7 R 800.
function [ state ] = movePTPCirc1OrintationInterCheck( t , f1,f2, relVel)
% This function is used to move the end-effector on a circle,

%% Arreguments
% t: is the TCP/IP connection
% f1: intermediate frame, to specify a point from the circle, it is 1x6 cell array.
% f2: final frame, to specify the end point of the circle, it is 1x6 cell array.
% the first three elements of cell array represent the X,Y and Z position of
% the frame
% the second three elements of cell array represent the alpha,beta
% and gamma rotaion angles (rads) that represent the frame orientation
% relVel: is a double, the over-ride relative velocity.

%% Return value:
% state: returns (true) if the circular motion is valid, returns false
% otherwise

% Copyright, Mohammad SAFEEA, 9th of May 2017

    theCommand=['jRelVel_',num2str(relVel),'_']; % set over ride.
    fprintf(t, theCommand);
    message=fgets(t);
    % The new position of end-effector, described in robot
    newPos={0,0,0,0,0,0};

    newPos{1}=f1{1};
    newPos{2}=f1{2};
    newPos{3}=f1{3};
    newPos{4}=f1{4};
    newPos{5}=f1{5};
    newPos{6}=f1{6};

    sendCirc1FramePos( t ,newPos); % send first frame of circle to server on controller.

    newPos{1}=f2{1};
    newPos{2}=f2{2};
    newPos{3}=f2{3};
    newPos{4}=f2{4};
    newPos{5}=f2{5};
    newPos{6}=f2{6};

    sendCirc2FramePos( t ,newPos); % send second frame of circle to server on controller.

    theCommand='doPTPinCSCircle1_';
    fprintf(t, theCommand); % start the point to point motion.
    message=fgets(t);

    readingFlag=false;

    message='';

    while readingFlag==false
        message=fgets(t);

        if checkAcknowledgment(message)
            state=true;
            break;
        end

        if checkErrorMessage(message)
            state=false;
            break;
         end
        pause(0.1);  %% This is to enable breaking the program from the outside, using Ctrl+C for example
    end

end

function [ output_args ] = sendCirc1FramePos( t ,jPos)
%% sendCircFramePos
% This function is used to send the first frame of the circle
% to the robot

% Pos: is 6 cells array of doubles
% t: is the TCP/IP connection object
% Copy right, Mohammad SAFEEA, 3rd of May 2017

theCommand='cArtixanPositionCirc1_';

for i=1:6
    x=sprintf('%0.2f',jPos{i});
    theCommand=[theCommand,x,'_'];
end

fprintf(t, theCommand);
message=fgets(t);
end

function [ output_args ] = sendCirc2FramePos( t ,jPos)
%% sendCircFramePos
% This function is used to send the first frame of the circle
% to the robot

% Pos: is 6 cells array of doubles
% t: is the TCP/IP connection object
% Copy right, Mohammad SAFEEA, 3rd of May 2017

theCommand='cArtixanPositionCirc2_';

for i=1:6
    x=sprintf('%0.2f',jPos{i});
    theCommand=[theCommand,x,'_'];
end

fprintf(t, theCommand);
message=fgets(t);
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function display(f,varargin)
% standard output

refSystems = [char(f.CS,'compact') ' ' getMTEXpref('arrowChar') ' ' char(f.SS,'compact')];

displayClass(f,inputname(1),varargin{:},'moreInfo',refSystems);

if length(f)~=1, disp([' size: ' size2str(f)]); end

% display symmetry
%dispLine(f.CS);
%dispLine(f.SS);

if f.antipodal, disp(' antipodal: true'); end
disp(' ');

if length(f)~=1, return; end

disp(['  h || r: ' char(round(f.h)) ' || (' char(round(f.r)) ')']);

% display starting and end orientation
if angle(f.o2,f.o1,'noSymmetry')>0
  disp([' o1 ' getMTEXpref('arrowChar') ' o2: ' char(f.o1) ' ' getMTEXpref('arrowChar') ' ' char(f.o2)]);
end

end

author: mtex-toolbox

repo: mtex

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


text:
function cheby_u_poly_values_test ( )

%*****************************************************************************80
%
%% CHEBY_U_POLY_VALUES_TEST demonstrates the use of CHEBY_U_POLY_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHEBY_U_POLY_VALUES_TEST:\n' );
  fprintf ( 1, '  CHEBY_U_POLY_VALUES returns values of\n' );
  fprintf ( 1, '  the Chebyshev U polynomials.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     N       X      U(N)(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, n, x, fx ] = cheby_u_poly_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %4d  %8f  %24.16f\n', n, x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function c = oned_binary_imdilate(a,b)

sza = size(a);
szb = size(b);
if nnz(sza>1) > 1 || nnz(szb>1) > 1,
  error('oned_binary_imdilate is only meant for one-dimensional arrays');
end
nb = length(b);
if mod(nb,2) == 0,
  b = [b,false];
end

c = myconv(a,b,'same',false) > 0;

author: kristinbranson

repo: JAABA


text:
function outputStr = formatRW(inputStr)
%adds whitespace after \reservedwordplot(1,1);

inputStrCell = cell(1,length(inputStr));

for c = 1:length(inputStr)
    inputStrCell{c} = inputStr(c);
end

rW = {'\\alpha','\\upsilon','\\sim','\\angle','\\phi','\\leq',...
    '\\ast','\\chi','\\infty','\\beta','\\psi','\\clubsuit',...
    '\\gamma','\\omega','\\diamondsuit','\\delta',...
    '\\Gamma','\\heartsuit','\\epsilon','\\Delta',...
    '\\spadesuit','\\zeta','\\Theta','\\leftrightarrow',...
    '\\eta','\\Lambda','\\leftarrow','\\theta','\\Xi',...
    '\\Leftarrow','\\vartheta','\\Pi','\\uparrow',...
    '\\iota','\\Sigma','\\rightarrow','\\kappa',...
    '\\Upsilon','\\Rightarrow','\\lambda','\\Phi',...
    '\\downarrow','\\mu','\\Psi','\\circ','\\nu',...
    '\\Omega','\\pm','\\xi','\\forall','\\geq','\\pi',...
    '\\exists','\\propto','\\rho','\\ni','\\partial',...
    '\\sigma','\\cong','\\bullet','\\varsigma',...
    '\\approx','\\div','\\tau','\\Re','\\neq','\\equiv',...
    '\\oplus','\\aleph','\\Im','\\cup','\\wp','\\otimes',...
    '\\subseteq','\\oslash','\\cap','\\in','\\supseteq',...
    '\\supset','\\lceil','\\subset','\\int','\\cdot','\\o',...
    '\\rfloor','\\neg','\\nabla','\\lfloor','\\times','\\lots',...
    '\\perp','\\surd','\\prime','\\wedge','\\varpi','\\0',...
    '\\rceil','\\rangle','\\mid','\\vee','\\copyright','\\langle'};

for w = 1:length(rW)
    [startInd endInd] = regexp(inputStr,rW{w});
    for ind = 1:length(endInd)
        %add space at end of reserved words
        inputStrCell{endInd(ind)}  = [inputStrCell{endInd(ind)} ' '];
    end
end

outputStr = [inputStrCell{:}];

end

author: plotly


text:
%
% FILE NAME:    applyJSON2RxnAssoc.m
%
% PURPOSE: Sort out and clean up the previous .mat files generated during
%          reaction association/curation of human-GEM:
%          1. identify redundant rxn-association .mat files and remove them;
%          2. convert the remained .mat files to JSON format as discussed
%          in #75, while retaining the content unchanged.
%

%% sort out and clean the reaction association files to HMR2

% There are two Matlab model structures (ihumanRxns2BiGG.mat and
% ihumanRxns2MNX.mat) that were generated from the HMR2 model during
% reaction-association curation. Since all the fields in ihumanRxns2BiGG
% are shared in ihumanRxns2MNX and with the same content, and the latter
% includes addtional information. The ihumanRxns2BiGG.mat file is moved to
% `deprecated` subfolder for reducing repetition.

% load the reaction association file and assign its info to a temp variable
load('ihumanRxns2MNX.mat');

r.rxns                = ihuman.rxns;                    % HMR id
r.rxnKEGGID           = ihuman.rxnKEGGID;               % KEGG
r.rxnEHMNID           = ihuman.rxnEHMNID;               % EHMN
r.rxnBiGGID           = ihuman.rxnBiGGID;               % BiGG
r.rxnHepatoNET1ID     = ihuman.rxnHepatoNET1ID;         % HepatoNET1
r.rxnREACTOMEID       = ihuman.rxnREACTOMEID;           % REACTOME
r.BiGG2BiGG           = ihuman.BiGG2BiGG;               % map to BiGG db id
r.HepatoNet12BiGG     = ihuman.HepatoNet12BiGG;         % HepatoNet1 map to BiGG id
r.EHMN2BiGG           = ihuman.EHMN2BiGG;               % EHMN map to BiGG id
r.HMR2BiGG            = ihuman.HMR2BiGG;                % combined BiGG ids from above
r.rxnBiGGDB2MNX       = ihuman.rxnBiGGDB2MNX;           % BiGG id to MNX id
r.rxnKEGG2MNX         = ihuman.rxnKEGG2MNX;             % KEGG id to MNX id
r.rxnREACTOMEStableID = ihuman.rxnREACTOMEStableID;     % REACTOME Stable ID
r.rxnReactome2MNX     = ihuman.rxnReactome2MNX;         % Reactome id 2 MNX id
r.rxnMNXID            = reformatElements(ihuman.rxnMNXID,'cell2str');  % MetaNetX, non-unique association
r.rxnCompIdx          = ihuman.rxnCompIdx;              % index of involved compartments

% encode to JSON format and pretty the layout before saving as plaintext file
jsonStr = jsonencode(r);
fid = fopen('ihumanRxns2MNX.JSON', 'w');
fwrite(fid, prettyJson(jsonStr));
fclose(fid);

% check if the content of JSON file is identical to that of .mat file
check = jsondecode(fileread('ihumanRxns2MNX.JSON'));
if isequal(r, check)
    fprintf('\nThe file ihumanRxns2MNX.JSON is confirmed with the same content to ihumanRxns2MNX.mat, which therefore can be removed!\n\n');
end

%% sort and clean up reaction association files to Recon3D

% There are two Matlab model structures (Recon3Rxns2MNX.mat and Recon3Rxns2HMR.mat)
% that were generated for Recon3D reaction-association curation. Since
% all the fields of Recon3Rxns2MNX.mat have already been included in
% Recon3Rxns2HMR.mat, which also includes some addtional fields. Hence the
% Recon3Rxns2MNX.mat is thus moved to `deprecated` subfolder for avoidance
% of repetition.

% load the reaction association file and assign its info to a temp variable
load('Recon3Rxns2HMR.mat');
clear r;   % clean and reuse the temp variable

r.rxns         = Recon3D.rxns;          % Recon3D id
r.rxnKEGGID    = Recon3D.rxnKEGGID;     % KEGG
r.rxnBiGGID    = Recon3D.rxnBiGGID;     % BiGG
r.rxnMNXID     = Recon3D.rxnMNXID;      % MetaNetX
r.rxnHMRID     = Recon3D.rxnHMRID;      % HMR id
r.BiGG2HMR     = Recon3D.BiGG2HMR;      % BiGG map to HMR id
r.rxnCompIdx   = Recon3D.rxnCompIdx;    % index of involved compartments
%r.comps       = Recon3D.comps;         % identical field is in Recon3DRAVEN
%r.compNames   = Recon3D.compNames;     % identical field is in Recon3DRAVEN
r.withMNXnoHMR = Recon3D.withMNXnoHMR;  % can be associated to MNX but not HMR

% Note that BiGG2HMR and rxnHMRID are different, the latter is used for the
% generation of rxnAssoc.mat

% encode to JSON format and pretty the layout before saving as plaintext file
jsonStr = jsonencode(r);
fid = fopen('Recon3Rxns2HMR.JSON', 'w');
fwrite(fid, prettyJson(jsonStr));
fclose(fid);

% check the content of JSON file, make sure they are identical to that of .mat file
check = jsondecode(fileread('Recon3Rxns2HMR.JSON'));
if isequal(r, check)
    fprintf('\nThe file Recon3Rxns2HMR.JSON is confirmed with the same content to Recon3Rxns2HMR.mat, which therefore can be removed!\n\n');
end

author: SysBioChalmers

repo: Human-GEM


text:
function transition_net(varargin)

opts.expDir = fullfile('../data','DCFNet-net-7') ;
opts = vl_argparse(opts, varargin);

modelPath = dir(fullfile(opts.expDir, 'net-epoch-*0.mat'));
modelPath = sort({modelPath.name});

[~,out_name,~] = fileparts(opts.expDir);

for i = numel(modelPath)
    load(fullfile(opts.expDir, modelPath{i}));
    net = dagnn.DagNN.loadobj(net) ;
    net = deployDCFNet(net);
    if ~exist('../model', 'dir')
        mkdir('../model');
    end
    save(fullfile('../model', [out_name '.mat']), 'net') ;
end
end

author: foolwood

repo: DCFNet

save_path: github-repos/MATLAB/foolwood-DCFNet


text:
function e1_values_test ( )

%*****************************************************************************80
%
%% E1_VALUES_TEST demonstrates the use of E1_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'E1_VALUES_TEST:\n' );
  fprintf ( 1, '  E1_VALUES stores values of\n' );
  fprintf ( 1, '  the exponential integral function E1(X).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X          E1(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = e1_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end


text:
%% Double Pendulum
% Clear work space and initaliza figure
clear; clc; close all; %figure;

%% Initial and Constant Values

% Initial values for tuned state (1):
% gamma0    = pi/12;
% dtgamma0  = pi;
% alpha0    = 0;
% dtalpha0  = 0;

% Initial values for tuned state (2):
% gamma0    = 0;
% dtgamma0  = 0;
% alpha0    = pi/6;
% dtalpha0  = pi/36;

% Initial values for tuned state (3):
% gamma0    = 0;
% dtgamma0  = 0;
% alpha0    = pi/24;
% dtalpha0  = pi/2;

% Initial values for chaotic solution:
% gamma0    = 0;
% dtgamma0  = 0;
% alpha0    = pi/2;
% dtalpha0  = 5;

% Initial values for chaotic solution2:
gamma0    = 0;
dtgamma0  = 0;
alpha0    = pi/2;
dtalpha0  = 5.025;

% Constant Values
beta0     = pi/2;
lambda    = 2.74; %rad/s
omega     = 5.48; %rad/s
psi       = 0.96;
eta       = 0.24;

fps       = 10;
%movie     = true;

%% Using solver ode45

duration  = 100;
ivp       = [gamma0; dtgamma0; alpha0; dtalpha0; beta0; lambda; omega; psi; eta];
[t,x]     = ode45(@dpendulum,[0 duration], ivp);

%options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
%[t,x]=ode45(@dpendulum,[0 duration], ivp, options);

%% Plot

nframes=duration*fps;
%t = linspace(0,duration,nframes);
%x = deval(sol,t);

gamma    = x(:,1);
dtgamma  = x(:,2);
alpha    = x(:,3);
dtalpha  = x(:,4);
% we assume:
l1=1; % b distance
l2=2; % c distance

%% Plot alpha vs Gamma
figure;
plot (t,x(:,3),'-.',t,x(:,1),'-')

%% Plot d_alpha vs d_Gamma
figure;
plot (t,x(:,4),'-.',t,x(:,2),'-')

repo: Aero-matlab


text:
function writeMesh(fileName, vertices, faces, varargin)
%WRITEMESH Write 3D mesh data by inferring format from file name.
%
%   writeMesh(FNAME, V, F)
%
%   writeMesh(FNAME, MESH)
%
%   Example
%   writeMesh
%
%   See also
%     meshes3d, readMesh, writeMesh_off, writeMesh_ply, writeMesh_stl
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2020-11-20, using Matlab 9.8.0.1323502 (R2020a)
% Copyright 2020-2022 INRAE - BIA Research Unit - BIBS Platform (Nantes)

% check inputs
if ~ischar(fileName)
    error('First argument must contain the name of the file');
end

% optionnaly parses data
if isstruct(vertices)
    if nargin > 2
        varargin = [{faces} varargin{:}];
    end
    faces = vertices.faces;
    vertices = vertices.vertices;
end

[~, ~, ext] = fileparts(fileName);
switch lower(ext)
    case '.off'
        writeMesh_off(fileName, vertices, faces);
    case '.ply'
        writeMesh_ply(fileName, vertices, faces, varargin{:});
    case '.stl'
        writeMesh_stl(fileName, vertices, faces, varargin{:});
    otherwise
        error('Unrecognized file format for rezading mesh: %s', ext);
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
% problem 8 - Computation of circular convolution

% c)
x1=[1,2,3,4];
x2=[5,4,3,2,1];
y1=circonv2(x1,x2)
y2=circonv3(x1,x2)

author: Sable

repo: mcbench-benchmarks


text:
function box = rectToBox(rect)
%RECTTOBOX Convert rectangle data to box data.
%
%   RECT = rectToBox(BOX)
%   Converts from rectangle representation to box representation.
%   BOX is given by [XMIN XMAX YMIN YMAX].
%   RECT is given by [X0 Y0 WIDTH HEIGHT], with WIDTH and HEIGHT > 0
%
%   Example
%   rectToBox
%
%   See also
%   boxToRect, drawBox, drawRect

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2012-08-23, using Matlab 7.9.0.529 (R2009b)
% Copyright 2012-2022 INRA - Cepia Software Platform

box = [rect(:,1) rect(:,1)+rect(:,3) rect(:,2) rect(:,2)+rect(:,4)];


text:
% PSWARM  Solve a Global NLP using PSWARM
%
% THIS IS A LOW LEVEL FUNCTION - USE opti_pswarm() INSTEAD!
%
% pswarm uses the Pattern and Particle Swarm Optimization library.
%
%   [x,fval,exitflag,iter,feval] = pswarm(fun,x0,lb,ub,A,b,opts)
%
%   Input arguments:
%       fun - nonlinear function handle
%       x0 - initial solution guess
%       lb - decision variable lower bounds
%       ub - decision variable upper bounds
%       A - linear inequality matrix (dense)
%       b - linear inequality rhs
%       opts - solver options (see below)
%
%   Return arguments:
%       x - solution vector
%       fval - objective value at the solution
%       exitflag - exit status (see below)
%       iter - number of iterations taken by the solver
%       feval - number of function evaluations taken by the solver
%
%   Option Fields (all optional - see pswarmset):
%       display - solver display level [0,1,2]
%       tolfun - function tolerance
%       maxiter - maximum solver iterations
%       maxfeval - Maximum number of function evaluations
%       maxtime - Maximum solver execution time
%       swarm_size - Swarm Size
%       vectorized - Objective function is vectorized
%       mu - Cognitial Parameter
%       nu - Social Parameter
%       iweight - Initial Weight
%       fweight - Final Weight
%       delta - Initial Delta
%       idelta - Increase Delta
%       ddelta - Decrease Delta
%
%   Return Status:
%       1 - converged
%       0 - maximum iterations / function evaluations exceeded
%      -1 - abnormal exit
%      -2 - memory error
%      -3 - population error
%
%
%   Copyright (C) 2012 Jonathan Currie (IPL)

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function[varargout]=arrayify(varargin)
%ARRAYIFY  Converts a set of scalars or arrays into column arrays.
%
%   [O1,O2,O3,...,ON]=ARRAYIFY(I1,I2,I3,...,IN) where each of the IN is
%   either a scalar or an array with M elements, converts each input
%   variable into a length M column vector.
%
%   The input variables that are M-element arrays are resized to column
%   vectors, while scalars are replicated to have M identical elements.
%
%   ARRAYIFY will return an error if any of the input variables that are
%   not scalars have different numbers of elements from each other.
%
%   ARRAYIFY(I1,I2,...IN); with no output arguments overwrites the original
%   input variables.
%
%   'arrayify --t' runs a test.
%
%   Usage: [o1,o2,o3]=arrayify(i1,i2,i3);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2013--2016 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--t')
    arrayify_test,return
end

M=zeros(length(varargin),1);
for i=1:length(varargin)
    varargin{i}=varargin{i}(:);
    M(i)=length(varargin{i});
end

if any(M>1)
    if any(M(M~=1)~=max(M))
        error('Input parameters must either be scalars or arrays having the same number of elements.')
    end
    M=max(M);
    for i=1:length(varargin)
        varargin{i}=varargin{i}+0*ones(M,1);
    end
    varargout=varargin;
    eval(to_overwrite(nargin));
else
    %Do nothing.  No need to overwrite in this case.
    varargout=varargin;
end

function[]=arrayify_test

[o1,o2,o3]=arrayify(4,1:10,[1:10]');

reporttest('ARRAYIFY',aresame(o1,4+0*[1:10]')&aresame(o2,[1:10]')&aresame(o3,[1:10]'))

i1=4;
i2=1:10;
i3=[1:10]';

arrayify(i1,i2,i3);

reporttest('ARRAYIFY with overwriting',aresame(i1,4+0*[1:10]')&aresame(i2,[1:10]')&aresame(i3,[1:10]'))

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
% Fig. 6.59a   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%
% response comparison of continuous and digital control step.

clear all;
%close all;
clf

sim('fig6_58step')
r=[1 1];   %reference input
t=[0 5];
plot(t,r,'r')
hold on
plot(ycd(:,1),ycd(:,2))
plot(ycd(:,1),ycd(:,3),'m:')
title('Figure 6.59(a) Step Responses of Digital and Continuous Controllers')
ylabel('y')
xlabel('Time (sec)')
text(1.1,1.06, '\leftarrow continuous  controller')
text(.7,1.23, '\leftarrow digital  controller')
nicegrid
hold off


text:
function [nFiles, aSize, anat] = CheckAnatomy(dirName)

% [nFiles, size] = CheckAnatomy(dirName);
%
% Looks at the files in the specified directory and finds all files
% with names of form I.*. Counts the number of files that continuously
% run in the sequence I.001, I.002, etc., and have the same non-zero
% size. Outputs the number of files found and their size.
%
% DBR  4/99

nFiles = 0;
aSize = [0 0];
anat = [];

if ~exist(dirName, 'dir'); return; end

dS = dir(fullfile(dirName, 'I.*'));
nList = length(dS);
if nList == 0; return; end

% Create a subset of files that match our filename criterion: I.nnn
nFList = 0;
for iList=1:nList
  fName = dS(iList).name;
  if length(fName) == 5 & strcmp('I.', upper(fName(1:2)))
    seqNo = str2num(fName(3:5));
    if length(seqNo)
      nFList = nFList + 1;
      fileList{nFList} = fName;
      seqNos(nFList) = seqNo;
    end
  end
end

% Sort the matching files in ascending numerical order.
[seqNos, sortInds] = sort(seqNos);
fileList = fileList(sortInds);

% Check the files for unbroken sequence and matching image size.
%
% Put up a status bar if there is more than a single file:
if nFList > 1
  hBar = mrvWaitbar(0, 'Scanning MRI files');
end
% Begin the second pass. Break if sequence or size isn't right.
for iList=1:nFList
  if seqNos(iList) ~= iList; break; end
  img = ReadMRImage(fullfile(dirName, fileList{iList}));
  imSize = size(img);
  if exist('firstSize', 'var')
    if ~all(firstSize == imSize); break; end
    nFiles = nFiles + 1;
    anat = cat(3, anat, img);
  else
    firstSize = imSize;
    nFiles = 1;
    anat = img;
  end
  if nFList > 1
    mrvWaitbar(iList/nFList)
  end
end
% Get rid of any status bar
if nFList > 1
  close(hBar)
end

aSize = firstSize;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ts,sr]=readnist(filename,sr);

% READNIST
%
%  [ts,sr]= READNIST (filename,[sr])
%
%  Reads NIST WAV/SPHERE files, such as those in N/C/TIMIT or
%  Switchboard or Mandarin-VOA TDT files.
%
%  filename can be provided with or without the .wav/.sph/.WAV/.SPH extension.
%
%  Example usage: to read '/blah/blah1/timit/dr1/fcjf0/sa1.wav' just say
%    [ts,sr] = readnist ('/blah/blah1/timit/dr1/fcjf0/sa1');
%  or
%    [ts,sr] = readnist ('/blah/blah1/timit/dr1/fcjf0/sa1.wav');
%
%  To read the ntimit file '/blah/blah1/ntimit/dr1/fcjf0/SA1.SPH'
%  just say
%    [ts,sr] = readnist ('/blah/blah1/ntimit/dr1/fcjf0/SA1');
%  or
%    [ts,sr] = readnist ('/blah/blah1/ntimit/dr1/fcjf0/SA1.SPH');
%
%  ts is a Nx1 vector of entries representing the sound sample in filename
%  sr represents the sample rate, and you DON'T have to supply it - this
%     function works it out from the file headers automatically.
%     If you do supply it, any sample rate info contained in headers is ignored.
%  The sr that is returned is
%    : the one input, or (if none is input)
%    : the sample rate in the file, or
%    : 0 (if no sample rate is supplied by the user or the file headers)
%
%  The first line of filename should have NIST_1A in its header
%    otherwise ts = [] and sr = 0
%
% Assumptions:
%  - first line of file has NIST_1A
%  - second line has an integer, usually 512 or 1024, that has the number of
%     bytes in the header
%  - the header has a line of the form "sample_rate -i D", where D is an
%     integer like 16000 that is the sampling rate in Hertz.
%  - the header has a line with "sample_byte_format" (if not, little endian
%     is assumed) and if so, that this line is either
%     "sample_byte_format -s2 01" (little endian) or
%     "sample_byte_format -s2 10" (big endian)
%
% You should be able to play the sound with the Matlab command
%   soundsc (ts,sr)
% afterwards, assuming your sound card exists and works. (If it doesnt,
% pray hard, and try again.)
%
% This file has only been tested on TIMIT, NTIMIT, CTIMIT and Switchboard (ICSI-WS97) files.
%
%	Author : D.Surendran 04/24/04 (dinojatcsdotuchicagodotedu)
%
% NOTE ON SHORTENED FILES
% Update: if you are on a unix-based system, you can also read
%   'shorten'ed files. This requires that you change the current
%   Matlab directory to one where you have write permission and
%   enough space to store an uncompressed version of the sound
%   file. The uncompressed file will be deleted afterwards, dont
%   worry!
%     The current directory should also have a copy of shorten
%   (If not, modify the command in this file that reads
%
%     cmd = ['shorten -x -d' num2str(headerlen) ' ' filename ' ' tmpfilename ';'];
%

% Acknowledgements to documentation http://ftp.cwi.nl/audio/NIST-SPHERE
% and the readtimit.m file by C.E.Ho 8/20/97 (which this file is meant to replace)

ts = [];

if nargin < 2
  sr = 0;
end
if nargin < 1
  error ('You need to have at least one input!');
end

exts = {'wav','WAV','sph','SPH'};
extsor = '';
if ~exist (filename)
  for i = 1 : length (exts)
    extsor = [extsor '.' exts{i}];
    if (i ~= length(exts))
      extsor = [extsor ' or '];
    end
    fn = [filename '.' exts{i}];
    if exist (fn)
      filename = fn;
      break;
    end
  end
end
if ~exist (filename)
  error (sprintf ('%s not found (even with %s added)',filename, extsor));
end

s=textread(filename,'%s',2);

if length(s) >= 1
  a= findstr(s{1},'NIST_1A');
  if ~(length(a)) | ~a
    error (sprintf('Not a NIST Sphere file : first line of %s doesnt contain NIST_1A',filename));
  end
else
  error (sprintf('Not a NIST Sphere file : first line of %s isnt a header',filename));
end

if length(s) >= 2
  headerlen = floor (str2num (s{2}));
  if 0 == headerlen
    error (sprintf('Second line of %s should have the header length in bytes',filename));
  end
end

fid = fopen (filename,'rt','ieee-le');             % doesnt matter which endian it is when reading ascii headers
[a,count] = fread (fid,headerlen,'uchar');         % uchar is default anyway
header = char(a');                                 % a is a vector of ascii values,
                                                   % header is an ascii string

headersplit = strread (header, '%s', length(header),'delimiter','\n');

% headersplit{i} has the i-th line of the file header

endian = 0;
useshorten = 0;
for i = 1 : length(headersplit)
  if findstr (headersplit{i},'sample_rate') & (0 == sr)
    sr = strread (headersplit{i}, 'sample_rate -i %d');
  end
  if findstr (headersplit{i},'sample_byte_format')
    endian = strread (headersplit{i}, 'sample_byte_format -s2 %d');
  end
  if findstr (headersplit{i},'sample_coding') & findstr (headersplit{i},'shorten')
    useshorten = 1;
  end
end

if (0 == sr)
  error ('What is the sample rate for this file? Cant find the sample_rate line');
end
if (0 == endian)
  warning (sprintf ('%s does not have byte format coded, assuming little endian'));
  endian = 'le';
elseif (1 == endian)
    endian = 'le';
elseif (10 == endian)
  endian = 'be';
else
  error ('error parsing line with sample_byte_format - cant work out whether this file is little or big endian');
end
fclose (fid);

if 1 == useshorten
  if ~isunix
     error ([filename 'needs to be uncompressed with Tony Robinsons' ...
	     'program - this can be done, but only when' ...
	     'this Matlab code is run on a unix-based system.' ...
	     'I cant get the ! command to behave like the' ...
	     'unix command in Matlab. Sorry! Note that' ...
	     'readnist works ok on windows for unshortened' ...
	     'NIST_1A sound files.']);
  end
  if 0 == exist ('shorten')
     error ([filename 'needs to be uncompressed with Tony Robinsons' ...
	     ' shorten program. Please installit from' ...
	     ' http://www.hornig.net/shorten/ and add its'...
	     ' directory to the Matlab path.']);
  end

  tmpfilename = ['uncompressed_' filename];

  cmd = ['./shorten -x -d' num2str(headerlen) ' ' filename ' ' tmpfilename ';'];
  unix (cmd);

  fid=fopen(tmpfilename,'r',['ieee-' endian] );
else
  fid=fopen(filename,'r',['ieee-' endian] );
end;

junk=fread(fid,headerlen,'short');
ts=fread(fid,Inf, 'short');
fclose(fid);

if (1 == useshorten) & isunix
  unix (['rm ' tmpfilename]);
else
  warning (['temporary file ' tmpfilename ' still exists!']);
end

author: posenhuang

repo: deeplearningsourceseparation


text:
% load data
% load('../results/LightenedCNN_A_lfw.mat');      % model A
% load('../results/LightenedCNN_B_lfw.mat');      % model B
load('../results/LightenedCNN_C_lfw.mat');      % model C
load('lfw_pairs.mat');

% pos
for i = 1: length(pos_pair)
    feat1 = features(pos_pair(1, i), :)';
    feat2 = features(pos_pair(2, i), :)';
    pos_scores(i) = distance.compute_cosine_score(feat1, feat2);
%     pos_scores(i) = -distance.compute_L2_score(feat1, feat2);
end
pos_label = ones(1, length(pos_pair));

%neg
for i = 1: length(neg_pair)
    feat1 = features(neg_pair(1, i), :)';
    feat2 = features(neg_pair(2, i), :)';
    neg_scores(i) = distance.compute_cosine_score(feat1, feat2);
%     neg_scores(i) = -distance.compute_L2_score(feat1, feat2);
end
neg_label = -ones(1, length(neg_pair));

scores = [pos_scores, neg_scores];
label = [pos_label neg_label];

% ap
ap = evaluation.evaluate('ap', scores, label);

% roc
roc = evaluation.evaluate('roc', scores, label);

%% output
fprintf('ap:           %f\n', ap.measure);
fprintf('eer:          %f\n', roc.measure);
fprintf('tpr001:       %f\n', roc.extra.tpr001*100);
fprintf('tpr0001:      %f\n', roc.extra.tpr0001*100);
fprintf('tpr00001:     %f\n', roc.extra.tpr00001*100);
fprintf('tpr000001:    %f\n', roc.extra.tpr000001*100);
fprintf('tpr0:         %f\n', roc.extra.tpr0*100);
result = [ap.measure/100 roc.measure/100  roc.extra.tpr001 roc.extra.tpr0001 roc.extra.tpr00001 roc.extra.tpr000001 roc.extra.tpr0];

author: AlfredXiangWu


text:
function timemetric
tic; [x,n,bco,bcn,gco,gcn,cs1,cs2]=metr(@elipsod,1); toc

function [x,n,bco,bcn,gco,gcn,cs1,cs2]=metr(func,method)
% [x,n,bco,bcn,gco,gcn,cs1,cs2]=metr(func,method)

% This function computes base vectors, metric tensor
% components, and Christoffel symbols for a general
% curvilinear coordinate system
% func     - handle of a function which returns the
%            cartesian radius vector x as a function
%            of the curvilinear coordinate variables,
%            and a vector containing the names of the
%            coordinate variables used to define x. In
%            spherical coordinates, for example, names
%            might be 'r, theta, phi'.
% method   - Use 1 to compute the Christoffel symbols by
%            differentiating the metric tensor components.
%            Use 2 to compute the Christoffel symbols by
%            by differentiating vector x.
% x        - cartesian components of the radius vector
%            expressed in terms of the curvilinear
%            coordinate variables
% n          a vector containing the names of the
%            curvilinear coordinate variables
% bco, bcn - cartesian components of the covariant and
%            contravariant base vectors. The i'th column
%            of each array contains components of the
%            i'th base vector.
% gco, gcn - covariant and contravariant components of
%            the metric tensor.
% cs1, cs2 - Christoffel symbols or the first and second
%            kinds. The symbol for index i,j,k is in
%            row i, column j, and layer k of each array.

if nargin<2, method=1; end
if nargin==0, func=@sphr; end
[x,n]=feval(func); x=simple(x(:));

% Differentiate the radius vector to get the
% contravariant base vectors.
% bco=[diff(x,n(1)),diff(x,n(2)),diff(x,n(3))]
bco=jacobian(x,n(:).');

% Use orthogonality to compute the contravariant
% base vectors.
bcn=simple(inv(bco.'));

% Obtain the metric tensor components as dot products of
% the base vectors
gco=simple(bco.'*bco); gcn=simple(bcn.'*bcn);

% If Christoffel symbols are not required, then exit
if nargout<6, return, end

% Compute the Christoffel symbols.
cs1=sym(zeros(3,3,3)); cs2=cs1;
if method==1
  % Obtain symbols of the first kind by differentiating
  % the covariant metric tensor components.
  for k=1:3
    for i=1:3, for j=1:i
      cs1(i,j,k)=1/2*(diff(gco(j,k),n(i))+...
             diff(gco(i,k),n(j))-diff(gco(i,j),n(k)));
      if j~=i, cs1(j,i,k)=cs1(i,j,k); end
    end, end
  end
else % method==2
  % Obtain symbols of the first kind using derivatives
  % of vector x.
  h=sym(zeros(3,3,3));
  for k=1:3, h(:,:,k)=simple(diff(bco(:,:),n(k))); end

  for k=1:3
    cs1(:,:,k)=squeeze(h(1,:,:)*bco(1,k)+h(2,:,:)*bco(2,k)+...
               h(3,:,:)*bco(3,k));
  end
%   for k=1:3, for i=1:3, for j=1:3
%     cs1(i,j,k)=h(1,i,j)*bco(1,k)+h(2,i,j)*bco(2,k)+...
%                h(3,i,j)*bco(3,k);
%   end,end,end
  cs1=simple(cs1);
end

% Obtain Christoffel symbols of the second kind by
% using the contravariant metric tensor components to
% raise the third index.
for k=1:3
  cs2(:,:,k)=cs1(:,:,1)*gcn(1,k)+cs1(:,:,2)*gcn(2,k)+...
             cs1(:,:,3)*gcn(3,k);
end
x=simple(x); bco=simple(bco); bcn=simple(bcn);
gco=simple(gco); gcn=simple(gcn);
cs1=simple(cs1); cs2=simple(cs2);

function [X,names]=elipsod
% [X,names]=elipsod defines ellipsoidal
% coordinates with 0 < lm < b,
% b < t < c, and c < e <infinity
syms e t lm real
names=[lm,t,e];
b=1; c=2*b; b2=b^2; c2=c^2; cb=c2-b2;
x=e*t*lm/(b*c);
y=sqrt((e^2-b2)*(t^2-b2)*...
    (b2-lm^2)/(b2*cb));
z=sqrt((e^2-c2)*(c2-t^2)*...
     (c2-lm^2)/(c2*cb));
X=[x;y;z];


text:
function [lb, ub] = bounds(var_struct, bounds)

lb_var = ocl.Variable.create(var_struct, -inf);
ub_var = ocl.Variable.create(var_struct, inf);

for k=1:length(bounds.data)

  d = bounds.data{k};
  id = d.id;
  lower = d.lower;
  upper = d.upper;

  lb_var.get(id).set(lower);
  ub_var.get(id).set(upper);

end

lb = lb_var.value;
ub = ub_var.value;

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
% Stop playback
PsychPortAudio('Stop', audios.pahandle);
% Close the audio device
PsychPortAudio('Close', audios.pahandle);
% Clear the screen
sca;

author: translationalneuromodeling

repo: tapas


text:
function UNew = diffusionNeumann2D(varargin);
% diffusionNeumann2D: solve diffusion registraion in 2D with Neumann
%        boundary conditions
%
%
% author: Nathan D. Cahill
% email: nathan.cahill@rit.edu
% affiliation: Rochester Institute of Technology
% date: January 2014
% licence: GNU GPL v3
%
% Copyright Nathan D. Cahill
% Code available from https://github.com/tomdoel/npReg
%
%

% parse input arguments
[U,F,PixSize,NumPix,RegularizerFactor] = parse_inputs(varargin{:});

% divide by regularizer factor
Fnew = F/RegularizerFactor;

% compute dct of new force field
FnewF1 = real(fft(real(fft(Fnew(:,:,1),2*NumPix(1)-2,1)),2*NumPix(2)-2,2));
FnewF2 = real(fft(real(fft(Fnew(:,:,2),2*NumPix(1)-2,1)),2*NumPix(2)-2,2));
FnewF1 = FnewF1(1:NumPix(1),1:NumPix(2));
FnewF2 = FnewF2(1:NumPix(1),1:NumPix(2));

% construct images of coordinates scaled by pi/(N or M)
[alpha,beta] = ndgrid(pi*(0:(NumPix(1)-1))/(NumPix(1)-1),pi*(0:(NumPix(2)-1))/(NumPix(2)-1));

% construct LHS factor
LHSfactor = 2*cos(alpha) + 2*cos(beta) - 4;

% set origin term to 1, as DC term does not matter
LHSfactor(1,1) = 1;

% solve for FFT of U
UF1 = FnewF1./LHSfactor;
UF2 = FnewF2./LHSfactor;

% if gamma is zero, set DC term to 0
UF1(1,1) = 0;
UF2(1,1) = 0;

% perform inverse dct
U1 = real(ifft(real(ifft(UF1,2*NumPix(1)-2,1)),2*NumPix(2)-2,2));
U2 = real(ifft(real(ifft(UF2,2*NumPix(1)-2,1)),2*NumPix(2)-2,2));

% crop and concatenate
UNew = cat(3,U1(1:NumPix(1),1:NumPix(2)),U2(1:NumPix(1),1:NumPix(2)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [U,F,PixSize,NumPix,RegularizerFactor] = parse_inputs(varargin);

% get displacement field and check size
U = varargin{1};
F = varargin{2};
PixSize = varargin{4}(1:2);
NumPix = [varargin{5} varargin{6}];
RegularizerFactor = varargin{10};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

author: tomdoel

repo: pulmonarytoolkit


text:
%% Gaussian Process Regression example from the gpml toolbox
clear all;
%% Generate some data

nbSamples   = 100;
epsilon     = 0.2;
X           = linspace(0,50,nbSamples);
y           = sin(X*0.2) + normrnd(0,0.2,1,nbSamples);

X           = X(:);
y           = y(:);

% Make hole in Data

id    = (X > 20 & X < 30);
X(id) = [];
y(id) = [];

% Plot data
options             = [];
options.points_size = 10;
options.title       = 'noisy sinusoidal data';

if exist('h1','var') && isvalid(h1), delete(h1);end
h1      = ml_plot_data([X(:),y(:)],options);

%% Train GP [gpml toolbox]

meanfunc = {@meanZero};
covfunc  = {@covSEiso};

ell      = 10;  % kernel width of RBF covariance function.
sf       = 1;   % signal variance (not measurement noise)
sn       = 0.0001; % measurement noise

hyp      = [];
hyp.cov  = log([ell; sf]);
hyp.lik  = log(sn);

% Test GP [gpml toolbox]

X_test = linspace(0,50,200)';

[m,s2] = gp(hyp, @infExact, meanfunc, covfunc, @likGauss, X, y, X_test);

if exist('h2','var') && isvalid(h2), delete(h2);end
h2 = figure;

f = [m+2*sqrt(s2); flipdim(m-2*sqrt(s2),1)];
fill([X_test; flipdim(X_test,1)], f, [7 7 7]/8);
hold on;

options             = [];
options.no_figure   = true;
options.points_size = 10;
ml_plot_data([X(:),y(:)],options);

author: epfl-lasa


text:
function undockfig(fig)
% Copyright 2008 The MathWorks, Inc.
% Programmatically undock one or all open figures.
%
% SYNTAX:
% undockfig
%    Undocks the current figure
%
% undockfig(fig)
%    Undocks figure with handle FIG, and brings it to the front.
%
% undockfig('all')
%    Undocks all open figures
%
% Written by Brett Shoelson, PhD
% brett.shoelson@mathworks.com
% 01/01/08

% Revision 1
% 01/03/08: Providing default behavior when no input argument is provided
% (operates on current figure).

if nargin == 0
    fig = gcf;
end
if ischar(fig) && strcmpi(fig,'all')
    fig = findall(0,'type','figure');
end
for ii = 1:numel(fig)
    set(fig(ii),'windowstyle','modal');
    set(fig(ii),'windowstyle','normal');
end
if numel(fig) == 1
    figure(fig);
end


text:
function [ lines ] = lineFromTwoPoint( pt1, pt2 )
%LINEFROMTWOPOINT Generate line segment based on two points on panorama
%   pt1, pt2: two points on panorama
%   lines:
%       1~3-th dim: normal of the line
%       4-th dim: the projection dimension ID
%       5~6-th dim: the u of line segment endpoints in projection plane
%   use paintParameterLine to visualize

numLine = size( pt1, 1);
lines = zeros(numLine, 6);
n = cross( pt1, pt2, 2);
n = n./repmat( sqrt(sum(n.^2,2)), 1, 3);
lines(:,1:3) = n;

areaXY = abs(sum(n.*repmat([0 0 1], [numLine 1]),2));
areaYZ = abs(sum(n.*repmat([1 0 0], [numLine 1]),2));
areaZX = abs(sum(n.*repmat([0 1 0], [numLine 1]),2));
[~, planeIDs] = max([areaXY areaYZ areaZX], [], 2); % 1:XY 2:YZ 3:ZX
lines(:,4) = planeIDs;

for i = 1:numLine
    uv = xyz2uvN([pt1(i,:); pt2(i,:)], lines(i,4));
    umax = max(uv(:,1))+pi;
    umin = min(uv(:,1))+pi;
    if umax-umin>pi
        lines(i,5:6) = [umax umin]/2/pi;
    else
        lines(i,5:6) = [umin umax]/2/pi;
    end
end

end

author: zouchuhang

repo: LayoutNet


text:
function varargout = mldivide(N, rhs, varargin)
%\    MLDIVIDE   Solve CHEBOP BVP or IVP system.
%
% MLDIVIDE is a convenient wrapper for CHEBOP/SOLVEBVP and CHEBOP/SOLVEIVP. It
% calls the appropriate method, depending on whether the problem being solved is
% a boundary-value problem, or an initial/final-value problem. Problems,
% specified by a CHEBOP N are determined to be a boundary-value or
% initial/final-value problems as follows:
%   * If N.LBC is non-empty, but both N.RBC and N.BC are empty, it's considered
%     to be an initial-value problem. In this case, CHEBOP/SOLVEIVP is called.
%   * If N.RBC is non-empty, but both N.LBC and N.BC are empty, it's considered
%     to be a final-value problem. In this case, CHEBOP/SOLVEIVP is called.
%   * Otherwise, the problem is considered to be a boundary-value problem. In
%     this case, CHEBOP/SOLVEBVP is called.
%
% See the documentation of the corresponding methods for further details.
%
% See also CHEBOP/SOLVEBVP, CHEBOP/SOLVEIVP.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Are we dealing with an initial or a final value problem. In that case, either
% we have
%   - N.LBC is nonempty, but both N.RBC and N.BC are empty. Here, we are dealing
%     with an initial value problem, where all conditions are imposed via N.LBC.
%   - N.RBC is nonempty, but both N.LBC and N.BC are empty. Here, we are dealing
%     with a final value problem, where all conditions are imposed via N.RBC.
isIVP =  ( ~isempty(N.lbc) && isempty(N.rbc) && isempty(N.bc) ) || ...
     ( isempty(N.lbc) && ~isempty(N.rbc) && isempty(N.bc) );

% If we did not get preferences passed in, we need to create a cheboppref so
% that we know wheter we want to solve the IVP via time-stepping or globally:
if ( nargin < 3 )
    pref = cheboppref();
else
    pref = varargin{1};
end

% Look at what solver we're using for IVPs:
ivpSolver = func2str(pref.ivpSolver);

% If we are solving an IVP, and we've specified one of the MATLAB solvers
% (ode113, ode15s, ode45) as our choice, call the solveivp method. Otherwise,
% solve problems globally:
if ( isIVP && ~isempty(strfind(ivpSolver, 'chebfun.ode')) )
    % Pass PREF here, as we'd have to start by constructing a CHEBOPPREF anyway
    % that start of solveivp otherwise.
    [varargout{1:nargout}] = solveivp(N, rhs, pref, varargin{:});
else
    % We have conditions in other fields, or we want to solve IVPs globablly,
    % call CHEBOP/SOLVEBVP. However, here, we don't want to pass PREF if it
    % wasn't included in VARARGIN, as that'd mean we couldn't to an automatic
    % switch to Fourier methods in the periodic case (i.e. change the
    % discretization to TRIGTECH).
    [varargout{1:nargout}] = solvebvp(N, rhs, varargin{:});
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function ds_plumb=cosmo_fmri_deoblique(ds)
% de-oblique a dataset
%
% Input:
%     ds                fmri dataset struct
%
% Output:
%
% Example:
%     % start with a simple dataset
%     x=cosmo_synthetic_dataset('size','huge','ntargets',1,'nchunks',1);
%     % make dataset oblique (manually)
%     x.a.vol.mat(1,1)=.8;
%     x.a.vol.mat(2,1)=.6;
%     y=cosmo_fmri_deoblique(x);
%     cosmo_disp(x.a.vol)
%     %|| .mat
%     %||   [ 0.8         0         0        -3
%     %||     0.6         2         0        -3
%     %||       0         0         2        -3
%     %||       0         0         0         1 ]
%     %|| .dim
%     %||   [ 20        17        19 ]
%     %|| .xform
%     %||   'scanner_anat'
%     cosmo_disp(y.a.vol)
%     %|| .mat
%     %||   [ 1         0         0      -3.2
%     %||     0         2         0      -2.4
%     %||     0         0         2        -3
%     %||     0         0         0         1 ]
%     %|| .dim
%     %||   [ 20        17        19 ]
%     %|| .xform
%     %||   'scanner_anat'
%     %
%     % other attributes are unchanged:
%     assert(isequal({x.samples x.fa x.sa x.a.fdim},...
%                    {y.samples y.fa y.sa y.a.fdim}));
%     %
%     % a plump dataset does not change after de-obliqueing
%     z=cosmo_fmri_deoblique(y);
%     isequal(y,z)
%     %|| true
%
% Notes:
%   - Using this function changes the location of the voxels in
%     world-space, that is world coordinates (x, y, z).
%   - This function is intended for AFNI and BrainVoyager, as these
%     programs prefer 'plump' (non-oblique) volumes.
%   - When using this function, it is recommended to inspect the result
%     visually.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

% get canonical orthogonal matrix
[unused,rot_ortho]=cosmo_fmri_orientation(ds);
mat=ds.a.vol.mat;

% convert base1 -> base0
mat(1:3,4)=mat(1:3,4)+mat(1:3,1:3)*[1 1 1]';

% use pixel dimension to set non-zero elements
mat_ortho=mat;
pixdim=sqrt(sum(mat(1:3,1:3).^2,1));
mat_ortho(1:3,1:3)=bsxfun(@times,rot_ortho(1:3,1:3),pixdim);

% convert base0 -> base1
mat_ortho(1:3,4)=mat_ortho(1:3,1:3)*-[1 1 1]'+mat_ortho(1:3,4);

% ensure single element in rotation part
assert(all(sum(mat_ortho(:,1:3)~=0,1)==1));
assert(all(sum(mat_ortho(1:3,1:3)~=0,2)==1));

ds_plumb=ds;
ds_plumb.a.vol.mat=mat_ortho;

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [p] = mesh_fit_elec(p)

% mesh_fit_elec - find mesh vertices nearest to electrodes
%
% [p] = mesh_fit_elec(p)
%
% This function is in development.  A transformation matrix,
% with rotations and translations, is in development to
% facilitate coregistration of elec vertices to the scalp
% vertices.  This needs to be refined by a minimisation
% function for the difference between the electrode
% locations and the nearest scalp vertices. The scalp
% mesh should be sufficiently refined to provide unique
% vertices for each electrode vertex.
%
% At present, it works with the BrainStorm meshes and the
% elec_124_cart.txt file in the eeg_example_data folder.
% The BrainStorm meshes contain a 'scalp' mesh with vertices
% in meters and the electrode file contains position
% vectors in centimeters.  The electrodes are converted
% to meters by ELEC_LOAD.  In the example data, they were
% previously coregistered and fitted to the scalp mesh.
%
% The return vertices of the scalp that are nearest to
% the electrodes are returned in the following:
%
% p.mesh.data.meshtype{p.mesh.current}
% p.mesh.data.vertices{p.mesh.current}
% p.mesh.data.faces{p.mesh.current}
%
% where p.mesh.current is also returned.
%
% See also: fiducial_coregister
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:57 $

% Licence:  GNU GPL, no express or implied warranties
% History:  04/2002, Darren.Weber_at_radiology.ucsf.edu
%           09/2002, Darren.Weber_at_radiology.ucsf.edu
%                    - bug fixing location of scalp mesh vertices
%                    that correspond to electrode vertices.
%                    - Still working with example dataset, not
%                    a genuine coregistration algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% debugging on/off
global DB;
DB = 0;
test = 0;  % coregistration testing

fprintf('MESH_FIT_ELEC: In development, not yet a true coregistration\n');

p.mesh.current = mesh_check(p,'scalp');

if isempty(p.mesh.current),
    error('MESH_FIT_ELEC: No ''scalp'' mesh in p struct.\n');
end

% Extract scalp triangulation
FV.vertices = p.mesh.data.vertices{p.mesh.current};
FV.faces    = p.mesh.data.faces{p.mesh.current};

if isempty(p.elec.data),
    error('MESH_FIT_ELEC: ''p.elec.data'' is empty.\n');
else
    x = p.elec.data.x;
    y = p.elec.data.y;
    z = p.elec.data.z;
end

fprintf('MESH_FIT_ELEC: ''Coregistering'' Electrodes...');
tic

if test,
	% Electrode fiducial points
	En = p.elec.data.nasion;
	Er = p.elec.data.rpa;
	El = p.elec.data.lpa;
	reg.fiducials.head = [ En; Er; El ];

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% NEED AN INTERACTIVE PROCESS TO SELECT THESE
	% Mesh fiducial points
	if isempty(p.mri.fiducials),
        if ~isempty(p.mri.data),
            fprintf('\n\nMESH_FIT_ELEC: Select Fiducial Points in MRI...\n\n');
            avw_view(p.mri.data);
            waitfor(gcf);
            p.mri.fiducials = evalin('base',p.mri.fiducials);
		else
            fprintf('\n\nMESH_FIT_ELEC: Select Fiducial Points in MRI...\n\n');
           [p] = mri_open(p);
            avw_view(p.mri.data);
            waitfor(gcf);
            p.mri.fiducials = evalin('base',p.mri.fiducials);
		end
	end
	reg.fiducials.mri = p.mri.fiducials;

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	T = fiducial_coregister(reg.fiducials.head,reg.fiducials.mri);

    E = [p.elec.data.x p.elec.data.y p.elec.data.z];
    E2MRI = [E,ones(size(E,1),1)] * T;
end

% perform least squares fit of elec to scalp???
%options = optimset('fminsearch');
%[sumd,fval,exitflag,output] = fminsearch('mesh_fit_elec_optim',...
%                              0, options, vertices, [x,y,z]);
%fprintf('\n%s%f\n', 'Iterations = ', output.iterations);
%fprintf('\n%s%d\n', 'Exit = ', exitflag);

% Try to adjust the height of the electrodes to
% those of the scalp mesh - really need a full
% translation/rotation/scaling least squares solution
elecmax = max(z);
vertmax = max(FV.vertices(:,3));
z = z + (vertmax - elecmax);

if DB,
    figure; plot3(x,y,z,'ro'); hold on;
    patch('vertices',FV.vertices,'faces',FV.faces,...
          'FaceColor',[0.0 0.0 0.6],'Edgecolor',[.6 .6 .6]);
    set(gca,'Projection','perspective');
    set(gca,'DataAspectRatio',[1 1 1]);
    axis off tight vis3d
end

% Find nearest scalp vertices to electrode positions.
% In the process, reorder the scalp mesh vertices
% and correct the scalp mesh faces so that all scalp
% vertices that lie near the electrodes
% are numbered 1-N, in the same order as the
% electrodes.  This is done here so that the
% mesh_laplacian_interp function will work
% correctly.

% This mess might be better handled by adding
% electrode vertices into the scalp mesh, rather
% than finding the scalp mesh vertices nearest
% to those of the electrodes.  However, this would
% require retesselation of the scalp mesh, which
% is a tricky process

[FV,k] = mesh_nearest_reorder(FV,[x,y,z]);

% Assign the adjusted triangulation back
% into the original structure
p.mesh.data.vertices{p.mesh.current} = FV.vertices;
p.mesh.data.faces{p.mesh.current}    = FV.faces;

p.mesh.data.elecindex{p.mesh.current} = k;

x = FV.vertices(k,1);
y = FV.vertices(k,2);
z = FV.vertices(k,3);

% A quick triangulation of the scalp vertices
% that lie nearest to the electrode vertices
Efaces = convhulln(FV.vertices(k,:));

% Now store the vertices
[p.mesh.current,meshExists] = mesh_check(p,'elec');

p.mesh.data.meshtype{p.mesh.current} = 'elec';
p.mesh.data.vertices{p.mesh.current} = [x,y,z];
p.mesh.data.faces{p.mesh.current} = Efaces;

t = toc; fprintf('done (%6.2f sec).\n',t);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [FV,k] = mesh_nearest_reorder(FV,elecverts)

    % Find the indices of the nearest vertices in
    % FV to those in vertices
    k = dsearchn(FV.vertices,elecverts);

    if size(k,1) ~= size(elecverts,1),
        msg = sprintf('Sorry, duplicate vertices found in electrode/mesh coregistration');
        error(msg);
    end

    % Debug visualization
    global DB;
    if DB,
        figure; hold on;
        patch('vertices',FV.vertices,'faces',FV.faces,...
              'FaceColor',[0.8 0.8 0.8],'Edgecolor',[.6 .6 .6]);
        set(gca,'Projection','perspective');
        set(gca,'DataAspectRatio',[1 1 1]);
        axis off tight vis3d
    end

    % OK, now reorder the triangulation in FV
    % so that all the mesh vertices are ordered
    % 1:N

    Nelec = size(elecverts,1);

    for Eindex = 1:Nelec,

        % It's important in this process that the
        % mesh indices are found over again, in case
        % they are modified for electrode Eindex-X
        elecpos = elecverts(Eindex,:);
        Kindex = dsearchn(FV.vertices,elecpos);

        if Eindex == Kindex,
            continue;  % great, do nothing
        end

        % Get vertex coordinates at Eindex
        Evertex = FV.vertices(Eindex,:);
        % Get vertex coordinates at Kindex
        Kvertex = FV.vertices(Kindex,:);

        if DB,
            fprintf('Fitting Electrode: %3d\n',Eindex);
            H1 = plot3(Evertex(1),Evertex(2),Evertex(3),'ro');
            H2 = plot3(elecpos(1),elecpos(2),elecpos(3),'rd');
            H3 = plot3(Kvertex(1),Kvertex(2),Kvertex(3),'bo');
            pause
            delete(H1)
            delete(H2)
            delete(H3)
        end

        % Now swap them
        FV.vertices(Kindex,:) = Evertex;
        FV.vertices(Eindex,:) = Kvertex;

        % Swapping the vertices requires
        % rearranging the faces

        % find all faces that contain Vindex
        EfaceIndices = find(FV.faces == Eindex);
        % find all faces that contain Kindex
        KfaceIndices = find(FV.faces == Kindex);

        % Where faces refer to Kvertex, make them
        % refer to Vvertex and vice versa
        FV.faces(KfaceIndices) = Eindex;
        FV.faces(EfaceIndices) = Kindex;

        % Debug testing....
        if DB,
            % Get vertex coordinates at Vindex
            Evertex = FV.vertices(Eindex,:);
            % Get vertex coordinates at Kindex
            Kvertex = FV.vertices(Kindex,:);
            H1 = plot3(Evertex(1),Evertex(2),Evertex(3),'ro');
            H2 = plot3(elecpos(1),elecpos(2),elecpos(3),'rd');
            H3 = plot3(Kvertex(1),Kvertex(2),Kvertex(3),'bo');
            pause

            delete(H1)
            delete(H2)
            delete(H3)
        end

    end

    % Now the nearest vertices in FV to vertices
    % should be numbered 1:Nelec
    k = dsearchn(FV.vertices,elecverts);

    if ~isequal(k,[1:Nelec]'),
        msg = sprintf('MESH_FIT_ELEC: Mesh vertices not in order of electrodes!');
        warning(msg);
    end

return

author: PatternRecognition

repo: OpenBMI


text:
function hh = errorbar2(x, y, l,u,symbol)
%ERRORBAR Error bar plot.
%   ERRORBAR(X,Y,L,U) plots the graph of vector X vs. vector Y with
%   error bars specified by the vectors L and U.  L and U contain the
%   lower and upper error ranges for each point in Y.  Each error bar
%   is L(i) + U(i) long and is drawn a distance of U(i) above and L(i)
%   below the points in (X,Y).  The vectors X,Y,L and U must all be
%   the same length.  If X,Y,L and U are matrices then each column
%   produces a separate line.
%
%   ERRORBAR(X,Y,E) or ERRORBAR(Y,E) plots Y with error bars [Y-E Y+E].
%   ERRORBAR(...,'LineSpec') uses the color and linestyle specified by
%   the string 'LineSpec'.  See PLOT for possibilities.
%
%   H = ERRORBAR(...) returns a vector of line handles.
%
%   For example,
%      x = 1:10;
%      y = sin(x);
%      e = std(y)*ones(size(x));
%      errorbar(x,y,e)
%   draws symmetric error bars of unit standard deviation.

%   L. Shure 5-17-88, 10-1-91 B.A. Jones 4-5-93
%   Copyright 1984-2000 The MathWorks, Inc.
%   $Revision: 1.1 $  $Date: 2003/05/12 22:36:15 $

if min(size(x))==1,
  npt = length(x);
  x = x(:);
  y = y(:);
    if nargin > 2,
        if ~isstr(l),
            l = l(:);
        end
        if nargin > 3
            if ~isstr(u)
                u = u(:);
            end
        end
    end
else
  [npt,n] = size(x);
end

if nargin == 3
    if ~isstr(l)
        u = l;
        symbol = '-';
    else
        symbol = l;
        l = y;
        u = y;
        y = x;
        [m,n] = size(y);
        x(:) = (1:npt)'*ones(1,n);;
    end
end

if nargin == 4
    if isstr(u),
        symbol = u;
        u = l;
    else
        symbol = '-';
    end
end

if nargin == 2
    l = y;
    u = y;
    y = x;
    [m,n] = size(y);
    x(:) = (1:npt)'*ones(1,n);;
    symbol = '-';
end

u = abs(u);
l = abs(l);

if isstr(x) | isstr(y) | isstr(u) | isstr(l)
    error('Arguments must be numeric.')
end

if ~isequal(size(x),size(y)) | ~isequal(size(x),size(l)) | ~isequal(size(x),size(u)),
  error('The sizes of X, Y, L and U must be the same.');
end

tee = (max(x(:))-min(x(:)))/100;  % make tee .02 x-distance for error bars
xl = x - tee;
xr = x + tee;
ytop = y + u;
ybot = y - l;
n = size(y,2);

% Plot graph and bars
hold_state = ishold;
cax = newplot;
next = lower(get(cax,'NextPlot'));

% build up nan-separated vector for bars
xb = zeros(npt*9,n);
xb(1:9:end,:) = x;
xb(2:9:end,:) = x;
xb(3:9:end,:) = NaN;
xb(4:9:end,:) = xl;
xb(5:9:end,:) = xr;
xb(6:9:end,:) = NaN;
xb(7:9:end,:) = xl;
xb(8:9:end,:) = xr;
xb(9:9:end,:) = NaN;

yb = zeros(npt*9,n);
yb(1:9:end,:) = ytop;
yb(2:9:end,:) = ybot;
yb(3:9:end,:) = NaN;
yb(4:9:end,:) = ytop;
yb(5:9:end,:) = ytop;
yb(6:9:end,:) = NaN;
yb(7:9:end,:) = ybot;
yb(8:9:end,:) = ybot;
yb(9:9:end,:) = NaN;

[ls,col,mark,msg] = colstyle(symbol); if ~isempty(msg), error(msg); end
symbol = [ls mark col]; % Use marker only on data part
esymbol = [ '-' col]; % Make sure bars are solid

h = plot(xb,yb,esymbol); hold on
h = [h;plot(x,y,symbol)];

if ~hold_state, hold off; end

if nargout>0, hh = h; end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%% Misc Geometry Functions
%
% This section of geometry tools especially contains methods to convert
% directions and rotations from one parametrization into another.
% Additionally, some basic geometrical objects are predefined.
%%
%

author: mtex-toolbox

repo: mtex


text:
function [fJointProbability, fSigFirst, fSigSecond] = st_probability(fBValueFirst, vBValuesFirst, fBValueSecond, vBValuesSecond)

    vBValuesFirst = vBValuesFirst(~isnan(vBValuesFirst),:);
    vBValuesSecond = vBValuesSecond(~isnan(vBValuesSecond),:);

    vBValuesFirstHi_ = vBValuesFirst(vBValuesFirst >= fBValueFirst);
    vBValuesFirstLo_ = vBValuesFirst(vBValuesFirst <= fBValueFirst);
    if fBValueSecond < fBValueFirst
        fSigFirst = calc_GetPercentile(fBValueSecond, vBValuesFirstLo_, 1)/100;
    elseif fBValueSecond > fBValueFirst
        fSigFirst = calc_GetPercentile(fBValueSecond, vBValuesFirstHi_, 0)/100;
    else
        fSigFirst = 1;
    end

    vBValuesSecondHi_ = vBValuesSecond(vBValuesSecond >= fBValueSecond);
    vBValuesSecondLo_ = vBValuesSecond(vBValuesSecond <= fBValueSecond);
    if fBValueFirst < fBValueSecond
        fSigSecond = calc_GetPercentile(fBValueFirst, vBValuesSecondLo_, 1)/100;
    elseif fBValueFirst > fBValueSecond
        fSigSecond = calc_GetPercentile(fBValueFirst, vBValuesSecondHi_, 0)/100;
    else
        fSigSecond = 1;
    end

    fJointProbability = fSigFirst * fSigSecond;
end

author: CelsoReyes


text:
function [p1,p2,dp1,dp2] = acrobotKinematics(z,p)
% [p1,p2,dp1,dp2] = acrobotKinematics(z,p)
%
% This function computes the mechanical energy of the acrobot.
%
% INPUTS:
%   z = [4,n] = state vector
%   p = parameter struct:
%       .m1 = elbow mass
%       .m2 = wrist mass
%       .g = gravitational acceleration
%       .l1 = length shoulder to elbow
%       .l2 = length elbow to wrist
%
% OUTPUTS:ls
%   p1 = [2,n] = position of the elbow joint
%   p2 = [2,n] = position of the wrist
%   dp1 = [2,n] = velocity of the elbow joint
%   dp2 = [2,n] = velocity of the wrist
%
% NOTES:
%
%   states:
%       1 = q1 = first link angle
%       2 = q2 = second link angle
%       3 = dq1 = first link angular rate
%       4 = dq2 = second link angular rate
%
%   angles: measured from negative j axis with positive convention
%

q1 = z(1,:);
q2 = z(2,:);
dq1 = z(3,:);
dq2 = z(4,:);

[p1,p2,dp1,dp2] = autoGen_acrobotKinematics(q1,q2,dq1,dq2,p.l1,p.l2);

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function [new_center,cl]=update_centers(SD,II,E,K,num,center)
center;
for iter1=2:num,
    com=SD(:,iter1);
    i1=II(:,iter1);
    comp=[com,II(:,1),i1];
    if size(find(center),1)==num,
       aa=sortrows(comp);
       bb=flipud(aa);
    else
        for i=1:size(find(center),1)
            a=center(find(center));
            b(i,:)=comp(a(i),:);
        end
       aa=sortrows(b);
       bb=flipud(aa);
    end
    for iter2=1:size(bb,1),
        row=bb(iter2,2);
        column=bb(iter2,3);
        sum1=E(row,1);
        sum2=E(column,1);
        if sum1>sum2,
           center(column,:)=[0];
           comp(column,3)=[0];
        else
           center(row,:)=[0];
           comp(row,2)=[0];
        end
        new_center=find(center);
        kk=size(new_center,1);
        if kk==K,
            break;
        end
    end
    if kk==K
       break;
    end
end
for i=1:num,
    for j=1:K,
        ind2=find(II(i,:)==new_center(j));
        temp(j,:)=SD(i,ind2);
    end
    [a,b]=max(temp);
    cl(i,:)=new_center(b);
end

author: Sable


text:
function [valid, saveflg]= isUniformized(scanNum, planC)
%"isUniformized"
%   Checks to see if the uniform data for the specified scanNum exists and
%   that all fields appear valid.  Returns true if they are, false if they
%   arent (therefore suggesting the data should be regenerated using
%   setUniformizedData).
%
%JRA 05/24/05
%DK  04/20/06 check uniformization for old plans
%Usage:
%   function valid = isUniformized(scanNum, planC);
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};
pathStr = getCERRPath;
optName = [pathStr 'CERROptions.json'];
optS = opts4Exe(optName);
saveflg = 0;
%Assume valid to start.
valid = 1;

%Check for required planC fields.
if ~isfield(indexS, 'scan') || ~isfield(indexS, 'structureArray')
    valid = 0;
    return;
end

%Check the two major planC components.
scanStruct = planC{indexS.scan}(scanNum);
if length(planC{indexS.structureArray})>=scanNum
    uniStruct  = planC{indexS.structureArray}(scanNum);
else
    uniStruct = [];
end

if isempty(scanStruct.uniformScanInfo)
    valid = 0;
    return;
end

%if isfield(uniStruct,'bitsArray') & isempty(uniStruct.bitsArray)& length(planC{indexS.structures})~= 0
if isfield(uniStruct,'bitsArray') && isempty(uniStruct.bitsArray) && ismember(scanNum,getAssociatedScan({planC{indexS.structures}.assocScanUID},planC))
    valid = 0;
    return;
elseif isempty(planC{indexS.structures})
    return;
end

if ~isfield(scanStruct, 'uniformScanInfo') ||...
        ~isfield(scanStruct, 'scanArrayInferior') || ~isfield(scanStruct, 'scanArraySuperior')
    valid = 0;
    return;
end

if isempty(uniStruct) || isempty(uniStruct.indicesArray) || isempty(uniStruct.bitsArray)
    %The uniformized structs are empty.  This should only happen if no
    %structures associated with this scan have no rasterSegments or if the
    %structures are excluded from uniformization in the options file.

    nStructs = length(planC{indexS.structures});
    assocScansV = getStructureAssociatedScan(1:nStructs, planC);

    %Check only structs that are associated with this scan.
    structsToCheck = find(assocScansV == scanNum);

    for i=1:length(structsToCheck)
        structNum = structsToCheck(i);
        structName = planC{indexS.structures}(structNum).structureName;

        %Ignore structures that are on the exclude list.
        if ismember(lower(structName), lower(optS.uniformizeExcludeStructs))
            continue;
        end

        [rS, planC] = getRasterSegments(structNum, planC);
        if ~isempty(rS) && ~(isfield(planC{indexS.structures},'rasterized') && planC{indexS.structures}(1).rasterized == 1)
            valid = 0;
            return;
        end
    end
else
    [~, ~, zV] = getUniformScanXYZVals(planC{indexS.scan}(scanNum));
    indicesM = planC{indexS.structureArray}(scanNum).indicesArray;
    askQuestion = 1;
%     for cellNum = 1:length(indicesC)
%         indicesM = indicesC{cellNum};
        if ~isempty(indicesM) && max(indicesM(:,3))>length(zV) && askQuestion
            valid = 0;
            querry = {'This is an old archived plan which was uniformized incorrectly'...
                'CERR will uniformized the plan again'...
                ''...
                'Would you like to save changes to original plan?'};
            ButtonName = questdlg(querry,'Uniformization Check','Yes','No','Yes');
            if strcmpi(ButtonName,'yes')
                saveflg = 1;
            end
        end
%     end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
filename='Gripping_tetrahedra_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'nonadjoint_compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
constraint_case = 'INEQUALITY';
optimizer = 'PROJECTED GRADIENT'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.2;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
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:
function [X12,x1,x2] = DigitalOscillator(f1,f2,Fs,t1)

% Function to simulate Digital Sinusodial Oscillator
% Input Variables
%    f1 = first input analog frequency    Suggested Range(697-1477)
%    f2 = second input analog frequency   Suggested Range(697-1477
%    Fs = Sampling Frequency              Default = 8000 cycles/sec
%    t1 = Tone Length                     Default = .1 sec
%
% Output Variables
%    Y  = Sine Signal for two given frequencies
%    x1 = Sine Signal for first frequency
%    x2 = Sine Signal for second frequency
%
%    Rajiv Singla        DSP Final Project            Fall 2005
%=====================================================================

% Checking for minimum number of arguments
if nargin < 2
    error('Not enough input arguments');
end

% Setting Default values
if nargin == 2
    Fs = 8000;
    t1 = 0.1;
end

Ts=1/Fs;   %Sampling Time
samples=t1/Ts;  %Number of samples in the tone

%% Generating first Signal
w0=2*3.1416*f1/Fs;
A=1;
a1=-2*cos(w0);
Y1=0;
Y2=-A*sin(w0);

for n=1:samples
    x1(n)=-a1*Y1-Y2;
    Y2=Y1;
    Y1=x1(n);
end
x1=[0 x1]; %appending zero in the beginning

%% Generating second Signal
w1=2*3.1416*f2/Fs;
A=1;
a2=-2*cos(w1);
P1=0;
P2=-A*sin(w1);

for n=1:samples
    x2(n)=-a2*P1-P2;
    P2=P1;
    P1=x2(n);
end
x2=[0 x2]; %appending zero in the beginning

%% Combining both signal
X12 = .5*(x1 + x2);

repo: mcbench-benchmarks


text:
function [m,Lambda,Cm] = spm_nwrnd (M,N)
% Generate N samples from Normal-Wishart density
% FORMAT [m,Lambda,Cm] = spm_nwrnd (M,N)
%
% Parameters M
%           .a,.B,.beta,.m
% N         number of samples
%
% m         Means
% Lambda    precisions
% Cm        covariances
%
% See J. Bernardo and A. Smith (2000)
% Bayesian Theory, Wiley (page 435)
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_nwrnd.m 6548 2015-09-11 12:39:47Z will $

%iB=inv(M.B);
for n=1:N,
    % Sample precisions
    %L=wishrnd(iB,M.a);
    L=spm_wishrnd(M.B,M.a);
    C=pinv(L);

    if nargout>1
        Lambda(:,:,n)=L;
        Cm(:,:,n)=C;
    end

    % Sample means
    %C=pinv(M.beta*L);
    %m(:,n)=spm_normrnd(M.m,C,1);

    m(:,n)=spm_normrnd(M.m,C/M.beta,1);
end

if N==1
    Lambda=squeeze(Lambda);
    Cm=squeeze(Cm);
end

author: spm


text:
function cdf = empirical_discrete_cdf ( x, a, b, c )

%*****************************************************************************80
%
%% EMPIRICAL_DISCRETE_CDF evaluates the Empirical Discrete CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the CDF.
%
%    Input, integer A, the number of values.
%    0 < A.
%
%    Input, real B(A), the weights of each value.
%    0 <= B(1:A) and at least one value is nonzero.
%
%    Input, real C(A), the values.
%    The values must be distinct and in ascending order.
%
%    Output, real CDF, the value of the CDF.
%
  cdf = 0.0;

  bsum = sum ( b(1:a) );

  for i = 1 : a

    if ( x < c(i) )
      return
    end

    cdf = cdf + b(i) / bsum;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [pcaicaAnalysisOutput] = runPcaIca(inputMovie,nPCs,nICs,varargin)
	% Wrapper to run PCA-ICA (Mukamel, 2009) cell extraction using two existing versions on CIAPKG.
	% Biafra Ahanonu
	% started: 2019.10.31 [10:01:18]
	% inputs
		% inputMovie - Matrix ([x y frames]) or char string pointing to movie path.
		% nPCs - Int: number of principal components to request from PCA.
		% nICs - Int: number of independent components to request from ICA. Should be less than nPCs.
	% outputs
		% pcaicaAnalysisOutput - structure containing output filters and traces for the requested number of nICs.
		% pcaicaAnalysisOutput.IcaFilters - Matrix with dimensions [x y nICs].
		% pcaicaAnalysisOutput.IcaTraces - Matrix with dimensions [nICs frames].

	% changelog
		% 2020.10.14 [13:52:47] - Updated to complete conversion from class to independent function.
		% 2020.10.17 [19:05:12] - Additional updates to deal with matrix inputs and give additional information for each PCA-ICA version output.
		% 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
	% TODO
		%

	import ciapkg.api.* % Import CIAtah package API

	%========================
	% Int:
		% 2 = Hakan/Tony version. (Preferred)
		% 1 = SpikeE version (Maggie Carr, Eran Mukamel, Jerome Lecoq, and Lacey Kitch, Biafra Ahanonu)
	options.version = 2;
	% Termination tolerance, e.g. 1e-5.
	options.TermTolICs = 10^(-5);
	% Str: hierarchy name in hdf5 where movie data is located
	options.inputDatasetName = '/1';
	% Float: parameter (between 0 and 1) specifying weight of temporal information in spatio-temporal ICA
	options.mu = 0.1;
	% Float: termination tolerance, e.g. 1e-5.
	options.term_tol = 5e-6;
	% Int: max iterations of FastICA, e.g. 750.
	options.max_iter = 1e3;
	% Str: output units, options: string of fluorescence ('fl'), standard deviation ('std'), '2norm', or variance ('var').
	options.outputUnits = 'fl';
	% Str: character string with regular expression if inputMovie is a folder path. In general ignore and give direct path to movie.
	options.fileFilterRegexp = 'concat';
	% Int vector: list of specific frames to load.
	options.frameList = [];
	% DEPRECIATED
		options.selectPCs = 0;
		options.selectICs = 0;
		options.displayIcImgs = 0;
	% get options
	options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	try
		if options.version==1
			disp('running PCA-ICA, old version...')
			startTime = tic;
			[PcaFilters, PcaTraces] = ciapkg.signal_extraction.pca_ica.runPCA(inputMovie, '', nPCs, options.fileFilterRegexp,'inputDatasetName',options.inputDatasetName,'frameList',options.frameList);
			if isempty(PcaFilters)
				disp('PCs are empty, skipping...')
				return;
			end
			[IcaFilters, IcaTraces, IcaInfo] = ciapkg.signal_extraction.pca_ica.runICA(PcaFilters, PcaTraces, '', nICs, '');
			pcaicaAnalysisOutput.IcaInfo = IcaInfo;
			traceSaveDimOrder = '[nComponents frames]';
			% reorder if needed
			options.IcaSaveDimOrder = 'xyz';
			if strcmp(options.IcaSaveDimOrder,'xyz')
				IcaFilters = permute(IcaFilters,[2 3 1]);
				imageSaveDimOrder = 'xyz';
			else
				imageSaveDimOrder = 'zxy';
			end
		elseif options.version==2
			disp('running PCA-ICA, new version...')
			startTime = tic;
			[PcaOutputSpatial, PcaOutputTemporal, PcaOutputSingularValues, PcaInfo] = ciapkg.signal_extraction.pca_ica_2.run_pca(inputMovie, nPCs, 'movie_dataset_name',options.inputDatasetName,'frameList',options.frameList);

			if isempty(PcaOutputTemporal)
				disp('PCs are empty, skipping...')
				return;
			end

			disp('+++')
			if ischar(inputMovie)==1
				movieDims = loadMovieList(inputMovie,'convertToDouble',0,'frameList',[],'inputDatasetName',options.inputDatasetName,'treatMoviesAsContinuous',1,'getMovieDims',1,'frameList',options.frameList);
			else
				movieDimsTmp = size(inputMovie);
				movieDims.x = movieDimsTmp(1);
				movieDims.y = movieDimsTmp(2);
				movieDims.z = movieDimsTmp(3);
			end

			% output_units = 'fl';
			% output_units = 'std';
			% options.PCAICA.term_tol = 5e-6;
			% options.PCAICA.max_iter = 1e3;
			[IcaFilters, IcaTraces, IcaInfo] = ciapkg.signal_extraction.pca_ica_2.run_ica(...
				PcaOutputSpatial,...
				PcaOutputTemporal,...
				PcaOutputSingularValues,...
				movieDims.x, movieDims.y,...
				nICs,...
				'output_units',options.outputUnits,...
				'mu',options.mu,...
				'term_tol',options.term_tol,...
				'max_iter',options.max_iter);
			IcaTraces = permute(IcaTraces,[2 1]);
			traceSaveDimOrder = '[nComponents frames]';
			% reorder if needed
			options.IcaSaveDimOrder = 'xyz';
			if strcmp(options.IcaSaveDimOrder,'xyz')
				imageSaveDimOrder = 'xyz';
			else
				IcaFilters = permute(IcaFilters,[3 1 2]);
				imageSaveDimOrder = 'zxy';
			end
			pcaicaAnalysisOutput.IcaInfo = IcaInfo;
			pcaicaAnalysisOutput.PcaInfo = PcaInfo;
		else
			disp('Incorrect version requested.')
			pcaicaAnalysisOutput.status = 0;
			return;
		end

		pcaicaAnalysisOutput.PcaIcaVersion = options.version;
		pcaicaAnalysisOutput.IcaFilters = IcaFilters;
		pcaicaAnalysisOutput.IcaTraces = IcaTraces;
		pcaicaAnalysisOutput.imageSaveDimOrder = imageSaveDimOrder;
		pcaicaAnalysisOutput.traceSaveDimOrder = traceSaveDimOrder;
		pcaicaAnalysisOutput.nPCs = nPCs;
		pcaicaAnalysisOutput.nICs = nICs;
		pcaicaAnalysisOutput.time.startTime = startTime;
		pcaicaAnalysisOutput.time.endTime = toc(startTime);
		pcaicaAnalysisOutput.time.dateTime = datestr(now,'yyyymmdd_HHMM','local');
		if ischar(inputMovie)
			pcaicaAnalysisOutput.movieList = inputMovie;
		else
			pcaicaAnalysisOutput.movieList = 'matrixInput';
		end
		pcaicaAnalysisOutput.status = 1;

	catch err
		disp(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		disp(repmat('@',1,7))
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [ i_lo, i_hi ] = r8vec_sorted_range ( n, r, r_lo, r_hi )

%*****************************************************************************80
%
%% R8VEC_SORTED_RANGE searches a sorted vector for elements in a range.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of items in the vector.
%
%    Input, real R(N), the sorted vector.
%
%    Input, real R_LO, R_HI, the limits of the range.
%
%    Output, integer I_LO, I_HI, the range of indices
%    so that I_LO <= I <= I_HI => R_LO <= R(I) <= R_HI.  If no
%    values in R lie in the range, then I_HI < I_LO will be returned.
%

%
%  Cases we can handle immediately.
%
  if ( r(n) < r_lo )
    i_lo = 0;
    i_hi = - 1;
    return
  end

  if ( r_hi < r(1) )
    i_lo = 0;
    i_hi = - 1;
    return
  end
%
%  Are there are least two intervals?
%
  if ( n == 1 )
    if ( r_lo <= r(1) && r(1) <= r_hi )
      i_lo = 1;
      i_hi = 1;
    else
      i_lo = 0;
      i_hi = -1;
    end
    return
  end
%
%  Bracket R_LO.
%
  if ( r_lo <= r(1) )

    i_lo = 1;

  else
%
%  R_LO is in one of the intervals spanned by R(J1) to R(J2).
%  Examine the intermediate interval [R(I1), R(I1+1)].
%  Does R_LO lie here, or below or above?
%
    j1 = 1;
    j2 = n;
    i1 = floor ( ( j1 + j2 - 1 ) / 2 );
    i2 = i1 + 1;

    while ( 1 )

      if ( r_lo < r(i1) )
        j2 = i1;
        i1 = floor ( ( j1 + j2 - 1 ) / 2 );
        i2 = i1 + 1;
      elseif ( r(i2) < r_lo )
        j1 = i2;
        i1 = floor ( ( j1 + j2 - 1 ) / 2 );
        i2 = i1 + 1;
      else
        i_lo = i1;
        break;
      end

    end

  end
%
%  Bracket R_HI
%
  if ( r(n) <= r_hi )

    i_hi = n;

  else

    j1 = i_lo;
    j2 = n;
    i1 = floor ( ( j1 + j2 - 1 ) / 2 );
    i2 = i1 + 1;

    while ( 1 )

      if ( r_hi < r(i1) )
        j2 = i1;
        i1 = floor ( ( j1 + j2 - 1 ) / 2 );
        i2 = i1 + 1;
      elseif ( r(i2) < r_hi )
        j1 = i2;
        i1 = floor ( ( j1 + j2 - 1 ) / 2 );
        i2 = i1 + 1;
      else
        i_hi = i2;
        break
      end

    end

  end
%
%  We expect to have computed the largest I_LO and smallest I_HI such that
%    R(I_LO) <= R_LO <= R_HI <= R(I_HI)
%  but what we want is actually
%    R_LO <= R(I_LO) <= R(I_HI) <= R_HI
%  which we can usually get simply by incrementing I_LO and decrementing I_HI.
%
  if ( r(i_lo) < r_lo )
    i_lo = i_lo + 1;
    if ( n < i_lo )
      i_hi = i_lo - 1;
    end
  end

  if ( r_hi < r(i_hi) )
    i_hi = i_hi - 1;
    if ( i_hi < 1 )
      i_lo = i_hi + 1;
    end
  end

  return
end


text:
function [node,elem,face]=s2m(v,f,keepratio,maxvol,method,regions,holes)
%
% [node,elem,face]=s2m(v,f,keepratio,maxvol,method)
% [node,elem,face]=s2m(v,f,keepratio,maxvol,'tetgen',regions,holes)
%
% volumetric mesh generation from a closed surface, shortcut for surf2mesh
%
% author: Qianqian Fang (fangq <at> nmr.mgh.harvard.edu)
%
% inputs and outputs are similar to those defined in surf2mesh
%
% if method='cgalpoly', s2m will call cgals2m and keepratio should be a
% structure (as the 'opt' input in cgals2m)
%
% input default values:
%       method: if ignored, iso2mesh uses surf2mesh ('tetgen') to do the
%               tetrahedral mesh generation
%       regions,holes: if ignored, iso2mesh assumes both are empty
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

p0=min(v(:,1:3));
p1=max(v(:,1:3));
if(nargin>=5 && strcmp(method,'cgalpoly'))
    [node,elem,face]=cgals2m(v,f,keepratio,maxvol);
    return;
end
if(nargin<=5)
    regions=[];
end
if(nargin<=6)
    holes=[];
end
[node,elem,face]=surf2mesh(v,f,p0,p1,keepratio,maxvol,regions,holes);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function plot_acorr(x, lag)

error('not ready');
if nargin < 2
  lag = size(x,2);
end

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function h=framemul(f,Fa,Fs,s,varargin)
%FRAMEMUL  Frame multiplier
%   Usage:  h=framemul(f,Fa,Fs,s);
%
%   Input parameters:
%          Fa   : Analysis frame
%          Fs   : Synthesis frame
%          s    : Symbol
%          f    : Input signal
%
%   Output parameters:
%          h    : Output signal
%
%   `framemul(f,Fa,Fs,s)` applies the frame multiplier with symbol *s*
%   to the signal *f*. The frame *Fa* is used for analysis and the frame
%   *Fs* for synthesis.
%
%   Examples:
%   ---------
%
%   In the following example Gabor coefficients obtained through the DGT
%   of pink noise are multiplied by the symbol batmask before resynthesis.
%   The result of this operation is an output signal *h* that is constructed
%   through a Gabor expansion of the modified coefficients.:::
%
%      f = pinknoise(400);
%      a = 10;
%      M = 40;
%      [Fa, Fs] = framepair('dgt', 'gauss', 'dual', a, M);
%      s = framenative2coef(Fa, batmask);
%      fhat = framemul(f, Fa, Fs, s);
%      figure(1);
%      plotframe(Fa,frana(Fa,f),'clim',[-100,-20]);
%      figure(2);
%      plotframe(Fa,s,'lin');
%      figure(3);
%      plotframe(Fa,frana(Fa,fhat),'clim',[-100,-20]);
%
%   See also: iframemul, framemuladj

% Author: Peter L. Søndergaard

if nargin < 4
    error('%s: Too few input parameters.',upper(mfilename));
end;

if size(s,2)>1
    error(['%s: Symbol should be a column vecor i.e. ',...
           'in the common Frames framework coefficient format. ',...
           'See FRAMENATIVE2COEF and FRAMECOEF2NATIVE.' ],upper(mfilename));
end

% Check for compatibility
L1=framelength(Fa,size(f,1));
L2=framelengthcoef(Fs,size(s,1));
if L1~=L2
    error(['%s: The symbol and signal lengths are incompatible.'],upper(mfilename));
end;

% This is not *strictly* necessary, but we cannot check that the symbol
% is complex-valued in just the right way.
if Fa.realinput && ~isreal(s)
    error(['%s: For real-valued-input-only frames, the symbol must also ' ...
           'be real.'],upper(mfilename));
end;

h=frsyn(Fs,bsxfun(@times,frana(Fa,f),s));

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function out = cellflat(celllist)
% CELLFLAT is a helper function to flatten nested cell arrays.
%
% CELLFLAT(celllist) searches every cell element in cellist and put them on
% the top most level. Therefore, CELLFLAT linearizes a cell array tree
% structure.
%
% Example: cellflat({[1 2 3], [4 5 6],{[7 8 9 10],[11 12 13 14 15]},{'abc',{'defg','hijk'},'lmnop'}})
%
% Output:
%Columns 1 through 7
%     [1x3 double]    [1x3 double]    [1x4 double]    [1x5 double]    'abc'    'defg'    'hijk'
%   Column 8
%     'lmnop'
%
% cellflat(({{1 {2 3}} 'z' {'y' 'x' 'w'} {4 @iscell 5} 6}) )
% Output:
% [1]    [2]    [3]    'z'    'y'    'x'    'w'    [4]    @iscell    [5]    [6]
%
% Version: 1.0
% Author: Yung-Yeh Chang, Ph.D. (yungyeh@hotmail.com)
% Date: 12/31/2014
% Copyright 2015, Yung-Yeh Chang, Ph.D.
% See Also: cell
if ~iscell(celllist)
    error('CELLFLAT:ImproperInputAugument','Input argument must be a cell array');
end
out = {};
for idx_c = 1:numel(celllist)
    if iscell(celllist{idx_c})
        out = [out cellflat(celllist{idx_c})];
    else
        out = [out celllist(idx_c)];
    end
end

repo: ifem


text:
function value = factorial_log ( n )

%*****************************************************************************80
%
%% FACTORIAL_LOG returns the logarithm of N!.
%
%  Discussion:
%
%    N! = Product ( 1 <= I <= N ) I
%
%    N! = Gamma(N+1).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the argument of the function.
%    0 <= N.
%
%    Output, real FACTORIAL_LOG, the logarithm of N!.
%
  if ( n < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FACTORIAL_LOG - Fatal error!\n' );
    fprintf ( 1, '  N < 0.\n' );
    error ( 'FACTORIAL_LOG - Fatal error!' );
  end

  value = 0.0;

  for i = 2 : n
    value = value + log ( i );
  end

  return
end


text:
function b = renamefields(a, old, new)

% RENAMEFIELDS renames a selection of the fields in a structure
%
% Use as
%   b = renamefields(a, old, new)
% which renames the fields with the old name to the new name. Fields that
% are specified but not present will be silently ignored.
%
% See also COPYFIELDS, KEEPFIELDS, REMOVEFIELDS

% Copyright (C) 2014, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if isempty(a)
  % this prevents problems if a is an empty double, i.e. []
  return
end

% these should be cell-arrays
if ischar(old)
  old = {old};
end
if ischar(new)
  new = {new};
end

if length(old)~=length(new)
  ft_error('the number of field names does not match between old and new');
end

% keep the fields that were not mentioned
b = keepfields(a, setdiff(fieldnames(a), old));
% copy the fields over with their new name
for i=1:length(old)
  if isfield(a, old{i})
    for j=1:numel(b)
      b(j).(new{i}) = a(j).(old{i});
    end
  end
end


text:
function [sys,x0,str,ts]=Plant2(t,x,u,flag)

switch flag,
    case 0,
        [sys,x0,str,ts]=mdlInitializeSizes;
    case 1,
        sys=mdlDerivatives(t,x,u);
    case 3,
        sys=mdlOutputs(x);
    case {2,4,9},
        sys=[];
    otherwise
        error(['Unhandled flag=',num2str(flag)]);
end
function [sys,x0,str,ts]=mdlInitializeSizes
    sizes=simsizes;
    sizes.NumContStates=2;
    sizes.NumDiscStates=0;
    sizes.NumOutputs=1;
    sizes.NumInputs=1;
    sizes.DirFeedthrough=1;
    sizes.NumSampleTimes=1;
    sys=simsizes(sizes);
    x0=[0;0];
    str=[];
    ts=[0 0];
function sys=mdlDerivatives(t,x,u)
    sys(1)=x(2);
    sys(2)=x(1)+x(2)+sign(sin(t/2))+u;
function sys=mdlOutputs(x)
    sys=x(1);

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
load camera_results;

string_global = 'global n_ima';
for kk = 1:n_ima,
   string_global = [string_global ' x_' num2str(kk) ' X_' num2str(kk) ' xproj_' num2str(kk) ' x_proj_' num2str(kk)];
end;
eval(string_global);

%----------------- global optimization: ---------------------

load projector_data; % load the projector corners (previously saved)
load projector_results;

solution_projector = solution;

load camera_results;

solution_camera = solution;

param_cam = solution_camera([1:10 16:end]);
param_proj = solution_projector([1:10 16:end]);

param = [param_cam;param_proj];

% Restart the minimization from here (if need be):
load camera_results;
load calib_cam_proj_optim2;

options = [1 1e-4 1e-4 1e-6  0 0 0 0 0 0 0 0 0 12000 0 1e-8 0.1 0];

%if 0, % use the full distortion model:

%   fprintf(1,'Take the complete distortion model\n');

   % test the global error function:
%   e_global = error_cam_proj(param);

%   param_init = param;

%   param = leastsq('error_cam_proj',param,options);

%else

   % Use a limitd distortion model (no 6th order)
   fprintf(1,'Take the 6th order distortion coefficient out\n');

   param = param([1:9 11:11+6*n_ima-1  11+6*n_ima:11+6*n_ima+9-1  11+6*n_ima+9+1:end]);

   % test the global error function:
   e_global2 = error_cam_proj2(param);

   param_init = param;

   param = leastsq('error_cam_proj2',param,options);

   param = [param(1:9);0;param(10:10+6*n_ima-1);param(10+6*n_ima:10+6*n_ima+9-1);0;param(10+6*n_ima+9:end)];

%end;

% Extract the parameters:

cam_proj_extract_param;

% Relative prosition of camera wrt world:
omc = omc_1;
Rc = Rc_1;
Tc = Tc_1;

% relative position of projector wrt world:
Rp = R*Rc;
omp = rodrigues(Rp);
Tp = T + R*Tc;

eval(['save calib_cam_proj_optim3  R om T fc fp cc cp alpha_c alpha_p kc kp Rc Rp Tc Tp omc omp param param_init']);

no_image = 0;
% Image size: (may or may not be available)
nx = 640;
ny = 480;

comp_error_calib;

% Save the optimal camera parameters:
saving_calib;
copyfile('Calib_Results.mat','camera_results_optim3.mat');
delete('Calib_Results.mat');

% Save the optimal camera parameters:
fc = fp;
cc = cp;
alpha_c = alpha_p;
kc = kp;

n_ima = 1;
X_1 = X_proj;
x_1 = x_proj;
omc_1 = om;
Tc_1 = T;
Rc_1 = R;

% Image size: (may or may not be available)
nx = 1024;
ny = 768;

% No calibration image is available (only the corner coordinates)
no_image = 1;

comp_error_calib;

saving_calib;
copyfile('Calib_Results.mat','projector_results_optim3.mat');
delete('Calib_Results.mat');

repo: ekfmonoslam


text:
function outpoints = mni2tal(inpoints)
% Converts coordinates from MNI brain to best guess
% for equivalent Talairach coordinates
% FORMAT outpoints = mni2tal(inpoints)
% Where inpoints is N by 3 or 3 by N matrix of coordinates
%  (N being the number of points)
% outpoints is the coordinate matrix with Talairach points
% Matthew Brett 10/8/99

dimdim = find(size(inpoints) == 3);
if isempty(dimdim)
  ft_error('input must be a N by 3 or 3 by N matrix')
end
if dimdim == 2
  inpoints = inpoints';
end

% Transformation matrices, different zooms above/below AC
upT = [
    0.9900    0.0000    0.0000    0.0000
    0.0000    0.9688    0.0460    0.0000
    0.0000   -0.0485    0.9189    0.0000
    0.0000    0.0000    0.0000    1.0000
      ];

downT = [
    0.9900    0.0000    0.0000    0.0000
    0.0000    0.9688    0.0420    0.0000
    0.0000   -0.0485    0.8390    0.0000
    0.0000    0.0000    0.0000    1.0000
    ];

tmp = inpoints(3,:)<0;  % 1 if below AC
inpoints = [inpoints; ones(1, size(inpoints, 2))];
inpoints(:, tmp) = downT * inpoints(:, tmp);
inpoints(:, ~tmp) = upT * inpoints(:, ~tmp);
outpoints = inpoints(1:3, :);
if dimdim == 2
  outpoints = outpoints';
end


text:
% Fig. 6.39   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%

clear all;
close all;

num=[1 1 1];
den=[1 10 24 0 0 0];
rlocus(num,den);
grid;
title('Fig. 6.39 Root Locus of a conditionally stable system');


text:
function A = double(T)
%DOUBLE Convert a sptenmat to a sparse matrix.
%
%   A = double(T) converts T stored as a SPTENMAT to a sparse matrix.
%
%   See also SPTENMAT.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

m = prod(T.tsize(T.rdims));
n = prod(T.tsize(T.cdims));
if isempty(T.subs)
    A = sparse(m,n);
else
    A = sparse(T.subs(:,1), T.subs(:,2), T.vals, m, n);
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt