text:
function [i,j]=getindxn(R,Idx);
%GETINDXN
%
%[i,j]=getindxn(R,Idx)

% Copyright (C) 1995-2006  Rasmus Bro & Claus Andersson
% Copenhagen University, DK-1958 Frederiksberg, Denmark, rb@life.ku.dk
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
% Street, Fifth Floor, Boston, MA  02110-1301, USA.

l=size(Idx,2);

i=Idx(1);
j=Idx(2);

if l==3,
  j = j + R(2)*(Idx(3)-1);
 else
  for q = 3:l,
    j = j + prod(R(2:(q-1)))*(Idx(q)-1);
  end;
end;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
Wxkpixbg=reshape(Wxk*Hkt,peval.nx,peval.ny,peval.nt)+peval.bg;
resid=(Wxkpixbg-dpixc);
resid_norm=resid./Wxkpixbg;
rnm=mean(resid_norm,3);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% book : Signals and Systems Laboratory with MATLAB
% authors : Alex Palamides & Anastasia Veloni
%
%
%
% problem 8- periodic signal
% in one period x[n]=[0.5,1,1]

 x=[0.5 1 1];

 xp=repmat(x,1,40);

 stem(10:length(xp)+9,xp);

author: Sable

repo: mcbench-benchmarks


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

%%
clear; close all; clc;

%% Syntax
% |[betaAngle] = patchFaceAngles(F,V);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


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

%%
clear; close all; clc;

%% Syntax
% |M_rice=imrician(M,s);|

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


text:
function pass = test_get( pref )
% Test GET.

if ( nargin < 1 )
    pref = chebfunpref;
end
tol = 1e2 * pref.cheb2Prefs.chebfun2eps;
j = 1;

f = spherefun(@(x,y,z) 1 + sin(pi*x.*y) + sin(pi*x.*z));
[C, D, R] = cdr( f );
pass(j) = norm( [-pi pi 0 pi] - f.domain ) < tol; j = j + 1;
pass(j) = norm( C - f.cols ) < tol; j = j + 1;
pass(j) = norm( R - f.rows ) < tol; j = j + 1;
pass(j) = norm( 1./diag(D) - f.pivotValues ) < tol; j = j + 1;
pass(j) = all( size(f.pivotLocations) == [length(f), 2] ); j = j + 1;
pass(j) = ( f.nonZeroPoles ) & ( abs(f(0,0,1)) > tol );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


%==============================================================================
% This code is part of the Matlab-based toolbox
% LagLDDDM - A Lagrangian Gauss--Newton--Krylov Solver for Mass- and
%                        Intensity-Preserving Diffeomorphic Image Registration
%
% For details and license info see
% - https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDMM
%
% function [Sc,dS,d2S] = diffusionST(vc,omega,m,varargin)
%
% Matrix-free spatio-temporal diffusion regularization energy for vc, where
% vc is cell-centered
%
% S(v) = 0.5 * \int_{\omega} alpha(1)*v(x)'*A*v(x)+ alpha(2)*v(x)'*B*v(x) dx,
%
% where A is the spatial gradient operator and B the time derivative operator.
%
% Input:
%
%   vc          instationary velocity field (cell-centered)
%   omega       spatial domain
%   m           number of discretization points
%   varargin    optional parameters (see below)
%
% Optional Input:
%
%   tspan       time span (default: [0 1])
%   nt          number of time points for velocity (default: computed from input)
%
%
% Output:
%
%   Sc          current value   (0.5 * hd * vc'* A *vc + 0.5*dt*vc'*B*vc)
%   dS          derivative      (hd * vc'*A )
%   d2S         Hessian, struct A
% ==================================================================================

function [Sc,dS,d2S] = mfDiffusionST(vc,omega,m,varargin)
if nargin == 0
    help(mfilename);
    return;
end

if strcmp(vc,'para')
    Sc = 'cell-centered';       % grid
    dS = 1;                     % matrixFree
    d2S = @spectralPrecondPCG;  % solver
    return;
end

alpha       = [1 1e-3];
tspan       = [0 1];
nt          = [];
for k=1:2:length(varargin) % overwrites default parameter
    eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

dim = numel(omega)/2;

if isempty(nt) % roughly estimate nt
    nt = round(numel(vc)/(prod(m)*dim))-1;
end

d2S.regularizer = regularizer;
d2S.alpha  = alpha;
d2S.B      = @(omega,m)     getDiffusionMatrixST(omega,tspan,m,nt,alpha);
d2S.d2S    = @(u,omega,m)   diffusionOperatorST(u,omega,tspan,m,nt,alpha);
d2S.diag   = @(omega,m)     diag(omega,tspan,m,nt,alpha);
d2S.solver = @spectralPrecondPCG;
d2S.res    = vc;
dS         = d2S.d2S(vc,omega,m)';
Sc         = .5*dS*vc;

function B = getDiffusionMatrixST(omega,tspan,m,nt,alpha)

h   = (omega(2:2:end)-omega(1:2:end))./m;
hd  = prod(h);
% compute time-stepsize
dt = abs(tspan(2)-tspan(1))/nt;

% build gradient matrix for one transformation
Bx =  getSpaceGradientMatrix(omega,m);
% a = sqrt(alpha(1).*hd.*dt)*ones(nt+1,1);
a = sqrt(alpha(1).*hd.*dt*[1/2;ones(nt-1,1);1/2]);
% apply spatial regularization to all transformations and sum up in
% time
Bx =  kron(sdiag(a(:)),Bx);

% get time regularization matrix
b = sqrt(alpha(2)*hd.*dt);
Bt = b * getTimeDiffusionMatrix(nt,dt,prod(m),length(omega)/2);

% build B
B = [Bx;Bt];

function D = sdiag(v)
	D = diag(sparse(v(:)));

% get diagonal of d2S (interesting in matrix free mode)
function D = diag(omega,tspan,m,nt,alpha)
dim   = numel(omega)/2;
one   = @(i,j) One(omega,m,i,j);
hd    = prod((omega(2:2:end)-omega(1:2:end))./m);
% compute time stepsize
dt = abs(tspan(1)-tspan(2))/nt;

if dim == 2
    Dx = [ one(1,1) + one(1,2);
           one(2,1) + one(2,2) ];
    Dx = kron(hd*dt*ones(nt+1,1),Dx);

    Dt = [1/2;ones(nt-1,1);1/2]/dt^2;

    Dt = kron(hd*dt*Dt , ones(prod(m)*dim,1));

    D  = alpha(1)*Dx + alpha(2)*Dt;

else
    Dx   = [ ...
        one(1,1)+one(1,2)+one(1,3);
        one(2,2)+one(2,1)+one(2,3);
        one(3,3)+one(3,1)+one(3,2)];
    Dx   = kron(hd*dts(:),Dx);

    Dt          = zeros(mTime+1,1);
    Dt(1:end-1) = (1./dt(:));
    Dt(2:end  ) = Dt(2:end) + (1./dt(:));

    Dt = kron(hd*Dt, ones(prod(mSpace)*dim,1));

    D  = alpha(1)*Dx + alpha(2)*Dt;

end;

% helper for computation of diag(d2S)
function o = One(omega,m,i,j)
h = (omega(2:2:end)-omega(1:2:end))./m;
o = ones(m)/h(j)^2;
switch j
    case 1, o(2:end-1,:,:) = 2*o(2:end-1,:,:);
    case 2, o(:,2:end-1,:) = 2*o(:,2:end-1,:);
    case 3, o(:,:,2:end-1) = 2*o(:,:,2:end-1);
end;
o = o(:);

% matrix free implementation of spatio-temporal diffusion operator
function Ay = diffusionOperatorST(vc,omega,tspan,m,nt,alpha)
dim = numel(omega)/2;
h   = (omega(2:2:end)-omega(1:2:end))./m;
hd  = prod(h(1:dim));
dt  = abs(tspan(1)-tspan(2))/nt;
w   = dt*[1/2;ones(nt-1,1);1/2];

switch dim
    case 2
        d1 = @(Y) (Y(2:end,:)-Y(1:end-1,:))/h(1);
        d2 = @(Y) (Y(:,2:end)-Y(:,1:end-1))/h(2);

        d1T = @(Y) reshape([-Y(1,:);Y(1:end-1,:)-Y(2:end,:);Y(end,:)],[],1)/h(1);
        d2T = @(Y) reshape([-Y(:,1),Y(:,1:end-1)-Y(:,2:end),Y(:,end)],[],1)/h(2);

        vc   = reshape(vc,[],nt+1);
        Ay = zeros(prod(m)*dim,nt+1);
        % spatial diffusion
        for k=1:nt+1
            vct = reshape(vc(:,k),[m dim]);
            Ay(:,k) = w(k) * hd*  alpha(1) * ....
                [ (d1T(d1(vct(:,:,1))) + d2T(d2(vct(:,:,1)))); ...
                (d1T(d1(vct(:,:,2))) + d2T(d2(vct(:,:,2))))];
        end
        Ay = Ay(:);

    case 3
        d1 = @(Y) (Y(2:end,:,:)-Y(1:end-1,:,:))/h(1);
        d2 = @(Y) (Y(:,2:end,:)-Y(:,1:end-1,:))/h(2);
        d3 = @(Y) (Y(:,:,2:end)-Y(:,:,1:end-1))/h(3);

        d1T = @(Y) reshape(d1t(Y),[],1)/h(1);
        d2T = @(Y) reshape(d2t(Y),[],1)/h(2);
        d3T = @(Y) reshape(d3t(Y),[],1)/h(3);

        vc   = reshape(vc,[],nt+1);
        Ay = zeros(prod(m)*dim,nt+1);
        % spatial diffusion
        for k=1:nt+1
            vct = reshape(vc(:,k),[m dim]);
            Ay(:,k) = w(k) * hd*  alpha(1) * ....
                [ ...
                (d1T(d1(vct(:,:,:,1))) + d2T(d2(vct(:,:,:,1))) + d3T(d3(vct(:,:,:,1)))); ...
                (d1T(d1(vct(:,:,:,2))) + d2T(d2(vct(:,:,:,2))) + d3T(d3(vct(:,:,:,2)))); ...
                (d1T(d1(vct(:,:,:,3))) + d2T(d2(vct(:,:,:,3))) + d3T(d3(vct(:,:,:,3)))); ...
                ];
        end
        Ay = Ay(:);
    otherwise
        error('%s - dimension %d not supported.',mfilename,dim);
end
% time diffusion
if alpha(2)>0
    Dt  = @(Y) (Y(:,2:end)-Y(:,1:end-1))/dt;
    DtT = @(Y) ([-Y(:,1),Y(:,1:end-1)-Y(:,2:end),Y(:,end)])/dt;
    At  = DtT(Dt(vc));
    Ay  = Ay + (alpha(2)* hd*dt)*At(:);
end

% partial derivative operator for x1 (mf)
function y = d1t(Y)
m = size(Y);
y = zeros(m+[1,0,0]);
y(1,:,:) = -Y(1,:,:);
y(2:end-1,:,:) = Y(1:end-1,:,:)-Y(2:end,:,:);
y(end,:,:) = Y(end,:,:);
% partial derivative operator for x2 (mf)
function y = d2t(Y)
m = size(Y);
y = zeros(m+[0,1,0]);
y(:,1,:) = -Y(:,1,:);
y(:,2:end-1,:) = Y(:,1:end-1,:)-Y(:,2:end,:);
y(:,end,:) = Y(:,end,:);
% partial derivative operator for x3 (mf)
function y = d3t(Y)
m = size(Y); if length(m) == 2, m = [m,1]; end;
y = zeros(m+[0,0,1]);
y(:,:,1) = -Y(:,:,1);
y(:,:,2:end-1) = Y(:,:,1:end-1)-Y(:,:,2:end);
y(:,:,end) = Y(:,:,end);

function A = getSpaceGradientMatrix(omega,m)
dim = length(omega)/2;
h   = (omega(2:2:end)- omega(1:2:end)) ./m ;

I = @(i)  speye(m(i));

% setup regularizer for y
%
% S(y) =  INT |Dy(.,t)|^2 + INT |d/dt y(x,.)|^2 dx
switch dim
    case 2
        % build discrete derivative operators
        D1=spdiags(ones(m(1),1)*[-1 1],0:1,m(1)-1,m(1)); D1=D1/(h(1));
        D2=spdiags(ones(m(2),1)*[-1 1],0:1,m(2)-1,m(2)); D2=D2/(h(2));

        % spatial regularization
        A = [kron(I(2),D1); kron(D2,I(1))];
    case 3
        % build discrete derivative operators
        D1=spdiags(ones(m(1),1)*[-1 1],0:1,m(1)-1,m(1)); D1=D1/(h(1));
        D2=spdiags(ones(m(2),1)*[-1 1],0:1,m(2)-1,m(2)); D2=D2/(h(2));
        D3=spdiags(ones(m(3),1)*[-1 1],0:1,m(3)-1,m(3)); D3=D3/(h(3));
        % build gradient (scalar)
        A = [kron(I(3),kron(I(2),D1)); ...
             kron(I(3),kron(D2,I(1))); ...
             kron(D3,kron(I(2),I(1)))];

end
A = kron(speye(dim), A);                  % nD gradient

function A = getTimeDiffusionMatrix(nt,dt,n,dim)
Dt=spdiags(ones(nt+1,1)*[-1/dt 1/dt],0:1,nt,nt+1);
A = kron(Dt, speye(n*dim)); % 2nd deriveative over time for two components


text:
function GUI = sessionGUI_selectROIType(type);
% Select the type of ROI to browse using the mrVista Session GUI.
%
% GUI = sessionGUI_selectROIType(<type name, or uicontrol handle>, <GUI>);
%
% Can provide the name or flag of an ROI type (described below), or else a
% handle to a uicontrol used to  select the ROI type.
%
% ROI types include:
%	'inplane', 1:	load an ROI from the local Inplane/ROIs directory.
%
%	'volume', 'volume shared', 'gray shared', 2: load an ROI from the
%	shared volume/gray directory. Usually this is 3DAnatomy/ROIs.
%
%	'gray local', 3: load an ROI from the local Gray/ROIs directory.
%
%	'volume local', 4: load an ROI from the local Volume/ROIs directory.
%
%	'browse', 'browse...', 5: browse for an ROI file. The ROI will be
%	assumed to be in the volume (I|P|R) coordinate system, and will be
%	aligned to the current data based on the saved transform into this
%	space.
%
%  All specifications are case-insensitive.
%
% If the handle to a uicontrol (popup) is provided, will select the
% session based on the currently-selected value of the uicontrol.
%
% If no argument is provided, will get from the GUI.controls.roiType
% listbox.
%
% The GUI input/output args do not generally need to be called, as this
% function will update the global variable GUI, as well as the mrGlobals
% variables. However, it will formally allow keeping this information in
% other, independent, structures.
%
% ras, 07/03/06.
mrGlobals2;

if notDefined('type'), type = GUI.controls.roiType; end

types = {'inplane' 'volume' 'gray local' 'volume local' 'browse...'};

% parse how the argument was passed
if ischar(type)
    type = cellfind(types, lower(type));
    if isempty(type)
        error('Invalid string value for ''type'' argument.')
    end

elseif isnumeric(type) & mod(type, 1)==0
    % ok as long as it's in [1-5]
    if (type < 1) | (type > 5)
        error('Invalid numeric value for type argument.')
    end

elseif ishandle(type)
    type = get(type, 'Value');

end

% update the ROI listbox to have a list of relevant ROIs for this type;
% select the first one by default:
switch type
	case {'inplane' 1}
		w = what( roiDir(INPLANE{1}) );

	case {'volume' 'volume shared' 'gray shared' 2}
		w = what( roiDir(VOLUME{1}, 0) );

	case {'gray local' 3}
		w = what( roiDir(VOLUME{1}, 1) );

	case {'volume local' 4}
		w = what( fullfile(HOMEDIR, 'Volume', 'ROIs') );

	case {'browse' 'browse...' 5}
		% for browsing, we ask the user to select the directory to browse
		% (but only if one hasn't already been selected)
		browseDir = get(GUI.controls.roiType, 'String');
		browseDir = browseDir{end};
		if strncmp( lower(browseDir), 'browse', 5 )
			txt = 'Select an ROI directory to browse';
			browseDir = uigetdir( fileparts(HOMEDIR), txt );
		end
		w = what(browseDir);

	otherwise
		error('Invalid ROI type.')
end

% if the user selected the 'browse' option, set the popup to reflect the
% browsing directory. Otherwise, reset the text to 'Browse...'.
str = get(GUI.controls.roiType, 'String');
if strncmp(type, 'browse', 5) | type==5
	str{end} = browseDir;
else
	str{end} = 'Browse...';
end
set(GUI.controls.roiType, 'String', str);

rois = {};
for i = 1:length(w.mat), [p rois{i} ext] = fileparts(w.mat{i}); end
rois = sort(rois);

set(GUI.controls.roi, 'String', rois, 'Value', 1);

% keep values in settings
GUI.settings.roiType = type;

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


function isBelow = isVersionBelow(versionA, versionB)
% Checks if versionA is smaller than versionB
    vA         = versionArray(versionA);
    vB         = versionArray(versionB);
    n          = min(length(vA), length(vB));
    deltaAB    = vA(1:n) - vB(1:n);
    difference = find(deltaAB, 1, 'first');
    if isempty(difference)
        isBelow = false; % equal versions
    else
        isBelow = (deltaAB(difference) < 0);
    end
end


text:
function magic3 ( n )

%*****************************************************************************80
%
%% MAGIC3 is a "main program" for a MATLAB compiler example.
%
%  Discussion:
%
%    This example is similar to magicsquare_print, but now it
%    includes calls to the timestamp() function, which is a separate
%    .M file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer/string N, the order of the magic square.
%
  if ( ischar ( n ) )
    n = str2num ( n );
  end

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MAGIC3\n' );
  fprintf ( 1, '  MATLAB version.\n' );
%
%  Compute the magic square.
%
  m = magic ( n );
%
%  Print the magic square.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Magic square of order %d:\n', n );
  fprintf ( 1, '\n' );
  for i = 1 : n
    for j = 1 : n
      fprintf ( 1, '  %2d', m(i,j) );
    end
    fprintf ( 1, '\n' );
  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MAGIC3\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%SingleCellVariation
%
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 3/23/2011
classdef SingleCellVariation
    methods (Static)
        function run(sim, fileName)
            import edu.stanford.covert.cell.sim.analysis.SingleCellVariation;
            import edu.stanford.covert.cell.sim.util.PlotUtil;
            import edu.stanford.covert.cell.sim.util.PrintUtil;

            %sample growth rate distribution
            initialGrowthFilterWidth = sim.state('MetabolicReaction').initialGrowthFilterWidth;
            sim.state('MetabolicReaction').initialGrowthFilterWidth = Inf;
            growthRates = SingleCellVariation.sampleSingleCellVariationDistribution(sim, 100);
            sim.state('MetabolicReaction').initialGrowthFilterWidth = initialGrowthFilterWidth;
            assertElementsAlmostEqual(sim.state('MetabolicReaction').meanInitialGrowthRate, mean(growthRates), 'relative', 0.50, 0);

            %% excel file
            [content, colLabels, indentation] = SingleCellVariation.printSingleCellVariationDistribution(sim, growthRates);
            if nargin == 1
                PrintUtil.printToStdIO(content, colLabels, struct('indentation', indentation));
            else
                PrintUtil.printToFile(content, colLabels, [fileName '.xls'], 'GrowthRate', struct('indentation', indentation));
            end

            %% plots
            if nargin == 1
                SingleCellVariation.plotSingleCellVariationGrowthRateDistribution(sim, growthRates, PlotUtil.newAxesHandle());
                SingleCellVariation.plotSingleCellVariationDoublingTimeDistribution(sim, growthRates, PlotUtil.newAxesHandle());
                SingleCellVariation.plotExpected_Vs_Simulated_RNA_Weight_Fractions(sim, PlotUtil.newAxesHandle());
                SingleCellVariation.plotExpected_Vs_Simulated_RNA_Expression(sim, PlotUtil.newAxesHandle());
                SingleCellVariation.plotExpected_Vs_Simulated_rRNAExpression(sim, PlotUtil.newAxesHandle());
                SingleCellVariation.plotExpected_Vs_Simulated_sRNA_Expression(sim, PlotUtil.newAxesHandle());
                SingleCellVariation.plotExpected_Vs_Simulated_tRNA_Expression(sim, PlotUtil.newAxesHandle());
                SingleCellVariation.plotExpected_Vs_Simulated_Monomer_Expression(sim, PlotUtil.newAxesHandle());
                SingleCellVariation.plotExpected_Vs_Simulated_NTP_Incorporation(sim, PlotUtil.newAxesHandle());
                SingleCellVariation.plotExpected_Vs_Simulated_AA_Incorporation(sim, PlotUtil.newAxesHandle());
            else
                [axesHandle, figHandle] = PlotUtil.newAxesHandle();

                cla(axesHandle);
                SingleCellVariation.plotSingleCellVariationGrowthRateDistribution(sim, growthRates, axesHandle);
                saveas(figHandle, [fileName '-GrowthRate.pdf']);

                cla(axesHandle);
                SingleCellVariation.plotSingleCellVariationDoublingTimeDistribution(sim, growthRates, axesHandle);
                saveas(figHandle, [fileName '-DoublingTime.pdf']);

                cla(axesHandle);
                SingleCellVariation.plotExpected_Vs_Simulated_RNA_Weight_Fractions(sim, axesHandle);
                saveas(figHandle, [fileName '-Expected_Vs_Simulated_RNA_Weight_Fractions.pdf']);

                cla(axesHandle);
                SingleCellVariation.plotExpected_Vs_Simulated_RNA_Expression(sim, axesHandle);
                saveas(figHandle, [fileName '-Expected_Vs_Simulated_RNA_Expression.pdf']);

                cla(axesHandle);
                SingleCellVariation.plotExpected_Vs_Simulated_rRNAExpression(sim, axesHandle);
                saveas(figHandle, [fileName '-Expected_Vs_Simulated_rRNA_Expression.pdf']);

                cla(axesHandle);
                SingleCellVariation.plotExpected_Vs_Simulated_sRNA_Expression(sim, axesHandle);
                saveas(figHandle, [fileName '-Expected_Vs_Simulated_sRNA_Expression.pdf']);

                cla(axesHandle);
                SingleCellVariation.plotExpected_Vs_Simulated_tRNA_Expression(sim, axesHandle);
                saveas(figHandle, [fileName '-Expected_Vs_Simulated_tRNA_Expression.pdf']);

                cla(axesHandle);
                SingleCellVariation.plotExpected_Vs_Simulated_Monomer_Expression(sim, axesHandle);
                saveas(figHandle, [fileName '-Expected_Vs_Simulated_Monomer_Expression.pdf']);

                cla(axesHandle);
                SingleCellVariation.plotExpected_Vs_Simulated_NTP_Incorporation(sim, axesHandle);
                saveas(figHandle, [fileName '-Expected_Vs_Simulated_NTP_Incorporation.pdf']);

                cla(axesHandle);
                SingleCellVariation.plotExpected_Vs_Simulated_AA_Incorporation(sim, axesHandle);
                saveas(figHandle, [fileName '-Expected_Vs_Simulated_AA_Incorporation.pdf']);

                close(figHandle);
            end
        end
    end

    methods (Static = true)
        function growthRates = sampleSingleCellVariationDistribution(sim, nTrials)
            warnStatus = warning('query', 'WholeCell:warning');
            warning('off', 'WholeCell:warning');

            growthRates = zeros(nTrials, 1);
            r = sim.state('MetabolicReaction');
            for i = 1:nTrials
                %seed rand stream
                sim.applyOptions(struct('seed', i));
                for j = 1:numel(sim.states)
                    o = sim.states{j};
                    o.seed = i;
                end
                for j = 1:numel(sim.processes)
                    o = sim.processes{j};
                    o.seed = i;
                end

                %calculate growth rate
                sim.initializeState();
                growthRates(i) = r.growth;
            end

            warning(warnStatus.state, 'WholeCell:warning');
        end
    end

    %printing
    methods (Static = true)
        function [content, colLabels, indentation] = printSingleCellVariationDistribution(sim, growthRates)
            %import classes
            import edu.stanford.covert.util.ConstantUtil;
            import edu.stanford.covert.cell.sim.Simulation;

            content = cell(0, 4);
            colLabels = {'Trial', 'Growth Rate (cell/s)', 'Doubling Time (hr)'};

            %time
            content = [content;
                num2cell(zeros(size(growthRates)))  cellfun(@(x) num2str(x), num2cell((1:numel(growthRates))'), 'UniformOutput', false)  num2cell(growthRates)  num2cell(1./(growthRates*3600)/sim.state('Mass').timeAveragedCellWeight)
                ];

            content = [content;{
                0 'Mean'          mean(growthRates)  1/(mean(growthRates)*3600)/sim.state('Mass').timeAveragedCellWeight
                0 'Min'           min(growthRates)   1/(max(growthRates)*3600)/sim.state('Mass').timeAveragedCellWeight
                0 'Max'           max(growthRates)   1/(min(growthRates)*3600)/sim.state('Mass').timeAveragedCellWeight
                0 'Experimental'  1/(sim.state('Time').cellCycleLength)/sim.state('Mass').timeAveragedCellWeight sim.state('Time').cellCycleLength/3600
                }];

            %format output
            indentation = cell2mat(content(:, 1));
            content = content(:, 2:end);
        end
    end

    %plotting
    methods (Static = true)
        function plotSingleCellVariationGrowthRateDistribution(sim, growthRates, axesHandle)
            hist(axesHandle, growthRates, 10);
            xlims = xlim(axesHandle);
            ylims = ylim(axesHandle);
            xlabel('Growth (cell/s)', 'FontSize', 12);
            ylabel('Frequency', 'FontSize', 12);

            h2 = line(mean(growthRates) * [1 1], ylims);
            set(h2, 'Color', 'g');

            h3 = line(median(growthRates) * [1 1], ylims);
            set(h3, 'Color', 'c');

            mu = mean(growthRates);
            sigma = std(growthRates);
            h4 = line([mu-sigma; mu-sigma], ylims);
            h5 = line([mu+sigma; mu+sigma], ylims);
            set(h4, 'Color', 'r');
            set(h5, 'Color', 'r');

            mr = sim.state('MetabolicReaction');
            mu = mr.meanInitialGrowthRate;
            wd = mr.initialGrowthFilterWidth;
            h6 = line(mu*(1-wd)*[1 1], ylims);
            h7 = line(mu*(1+wd)*[1 1], ylims);
            set(h6, 'Color', 'y');
            set(h7, 'Color', 'y');

            h = legend([h2 h3 h4 h6], 'Mean', 'Median', 'Mean \pm 1 Std', 'Filter');
            set(h, 'Location', 'NorthWest');

            xlim(axesHandle, xlims);
            ylim(axesHandle, ylims);
        end

        function plotSingleCellVariationDoublingTimeDistribution(sim, growthRates, axesHandle)
            thresh = sim.state('Time').cellCycleLength/3600 * 2;

            doublingTimes = 1./(growthRates * 3600)/sim.state('Mass').timeAveragedCellWeight;
            hist(axesHandle, doublingTimes(doublingTimes <= thresh), 10);
            h1 = findobj(gca, 'Type', 'patch');
            xlabel('Doubling Time (h)', 'FontSize', 12);
            ylabel('Frequency', 'FontSize', 12);
            xlim([floor(min(doublingTimes(doublingTimes <= thresh))) ceil(max(doublingTimes(doublingTimes <= thresh)))]);

            h2 = line(1/(mean(growthRates) * 3600) / sim.state('Mass').timeAveragedCellWeight * [1 1], [0 max(ylim)]);
            set(h2, 'Color', 'g');

            h3 = line(1/(median(growthRates) * 3600) / sim.state('Mass').timeAveragedCellWeight * [1 1], [0 max(ylim)]);
            set(h3, 'Color', 'c');

            h4 = line(sim.state('Time').cellCycleLength * [1 1] / 3600, [0 max(ylim)]);
            set(h4, 'Color', 'r');

            h = legend([h1(1) h2 h3 h4], sprintf('Simulation, growth (%.1f%%)', 100 * sum(doublingTimes <= thresh) / numel(doublingTimes)), 'Sim-Mean', 'Sim-Median', 'Exp-Mean');
            set(h, 'Location', 'NorthEast');
        end

        function plotExpected_Vs_Simulated_RNA_Weight_Fractions(sim, axesHandle)
            import edu.stanford.covert.cell.sim.analysis.Constants;
            import edu.stanford.covert.util.ConstantUtil;

            r = sim.state('Rna');

            weightFractions = [...
                r.molecularWeights(r.matureIndexs(r.matureMRNAIndexs))' * ...
                sum(sum(r.counts(r.matureIndexs(r.matureMRNAIndexs), :, :), 3), 2);...
                r.molecularWeights(r.matureIndexs(r.matureRibosomalRRNAIndexs)) .* ...
                sum(sum(r.counts(r.matureIndexs(r.matureRibosomalRRNAIndexs), :, :), 3), 2);...
                r.molecularWeights(r.matureIndexs(r.matureSRNAIndexs))' * ...
                sum(sum(r.counts(r.matureIndexs(r.matureSRNAIndexs), :, :), 3), 2)
                r.molecularWeights(r.matureIndexs(r.matureTRNAIndexs))' * ...
                sum(sum(r.counts(r.matureIndexs(r.matureTRNAIndexs), :, :), 3), 2)] / ...
                ConstantUtil.nAvogadro;

            plot(axesHandle, r.expectedWeightFractions, weightFractions, '.', 'MarkerSize', 20);
            line([0 1], [0 sum(weightFractions)], 'Parent', axesHandle, 'Color', 'r');
            xlim([min(r.expectedGeneDecayRates) max(r.weightFractions)]);
            ylim([min(weightFractions) max(weightFractions)]);
            xlabel(axesHandle, 'Expected Weight Fraction', 'fontSize', 16);
            ylabel(axesHandle, 'Simulated Weight (g)', 'fontSize', 16);
        end

        function plotExpected_Vs_Simulated_RNA_Expression(sim, axesHandle)
            import edu.stanford.covert.cell.sim.analysis.Constants;

            r = sim.state('Rna');

            expression = r.expression(r.matureIndexs);

            RNAs = sum(sum(...
                r.counts(r.processedIndexs,     :, :) + ...
                r.counts(r.matureIndexs,        :, :) + ...
                r.counts(r.boundIndexs,         :, :) + ...
                r.counts(r.misfoldedIndexs,     :, :) + ...
                r.counts(r.damagedIndexs,       :, :) + ...
                r.counts(r.aminoacylatedIndexs, :, :), ...
                3), 2);

            plot(axesHandle, expression, RNAs, '.', 'MarkerSize', 20);
            line([0 1], [0 sum(RNAs)], 'Parent', axesHandle, 'Color', 'r');
            xlim([min(expression) max(expression)]);
            ylim([min(RNAs) max(RNAs)]);
            xlabel(axesHandle, 'Expected Expression', 'fontSize', 16);
            ylabel(axesHandle, 'Simulated Counts', 'fontSize', 16);
        end

        function plotExpected_Vs_Simulated_rRNAExpression(sim, axesHandle)
            import edu.stanford.covert.cell.sim.analysis.Constants;

            r = sim.state('Rna');
            c = sim.state('ProteinComplex');

            expression = r.geneExpression(sim.gene.rRNAIndexs);
            expression = expression / sum(expression);

            RNAs = sum(sum(...
                r.counts(r.processedIndexs(    r.matureRRNAIndexs), :, :) + ...
                r.counts(r.matureIndexs(       r.matureRRNAIndexs), :, :) + ...
                r.counts(r.boundIndexs(        r.matureRRNAIndexs), :, :) + ...
                r.counts(r.misfoldedIndexs(    r.matureRRNAIndexs), :, :) + ...
                r.counts(r.damagedIndexs(      r.matureRRNAIndexs), :, :) + ...
                r.counts(r.aminoacylatedIndexs(r.matureRRNAIndexs), :, :), ...
                3), 2) + ...
                sum(c.proteinComplexComposition(sim.gene.rRNAIndexs, :, :), 3) * ...
                sum(sum(...
                c.counts(c.matureIndexs,      :, :) + ...
                c.counts(c.inactivatedIndexs, :, :) + ...
                c.counts(c.boundIndexs,       :, :) + ...
                c.counts(c.misfoldedIndexs,   :, :) + ...
                c.counts(c.damagedIndexs,     :, :), 3), 2);

            plot(axesHandle, expression, RNAs, '.', 'MarkerSize', 20);
            line([0 1], [0 sum(RNAs)], 'Parent', axesHandle, 'Color', 'r');
            xlim([min(0.3, min(expression)) max(0.4, max(expression))]);
            ylim([min(RNAs)-1 max(RNAs)+1]);
            xlabel(axesHandle, 'Expected Expression', 'fontSize', 16);
            ylabel(axesHandle, 'Simulated Counts', 'fontSize', 16);
        end

        function plotExpected_Vs_Simulated_sRNA_Expression(sim, axesHandle)
            import edu.stanford.covert.cell.sim.analysis.Constants;

            r = sim.state('Rna');

            expression = r.geneExpression(sim.gene.sRNAIndexs);
            expression = expression / sum(expression);

            RNAs = sum(sum(...
                r.counts(r.processedIndexs(    r.matureSRNAIndexs), :, :) + ...
                r.counts(r.matureIndexs(       r.matureSRNAIndexs), :, :) + ...
                r.counts(r.boundIndexs(        r.matureSRNAIndexs), :, :) + ...
                r.counts(r.misfoldedIndexs(    r.matureSRNAIndexs), :, :) + ...
                r.counts(r.damagedIndexs(      r.matureSRNAIndexs), :, :) + ...
                r.counts(r.aminoacylatedIndexs(r.matureSRNAIndexs), :, :), ...
                3), 2);

            plot(axesHandle, expression, RNAs, '.', 'MarkerSize', 20);
            line([0 1], [0 sum(RNAs)], 'Parent', axesHandle, 'Color', 'r');
            xlim([min(expression) max(expression)]);
            ylim([min(RNAs) max(RNAs)]);
            xlabel(axesHandle, 'Expected Expression', 'fontSize', 16);
            ylabel(axesHandle, 'Simulated Counts', 'fontSize', 16);
        end

        function plotExpected_Vs_Simulated_tRNA_Expression(sim, axesHandle)
            import edu.stanford.covert.cell.sim.analysis.Constants;

            r = sim.state('Rna');

            expression = r.geneExpression(sim.gene.tRNAIndexs);
            expression = expression / sum(expression);

            RNAs = sum(sum(...
                r.counts(r.processedIndexs(    r.matureTRNAIndexs), :, :) + ...
                r.counts(r.matureIndexs(       r.matureTRNAIndexs), :, :) + ...
                r.counts(r.boundIndexs(        r.matureTRNAIndexs), :, :) + ...
                r.counts(r.misfoldedIndexs(    r.matureTRNAIndexs), :, :) + ...
                r.counts(r.damagedIndexs(      r.matureTRNAIndexs), :, :) + ...
                r.counts(r.aminoacylatedIndexs(r.matureTRNAIndexs), :, :), ...
                3), 2);

            plot(axesHandle, expression, RNAs, '.', 'MarkerSize', 20);
            line([0 1], [0 sum(RNAs)], 'Parent', axesHandle, 'Color', 'r');
            xlim([min(expression) max(expression)]);
            ylim([min(RNAs) max(RNAs)]);
            xlabel(axesHandle, 'Expected Expression', 'fontSize', 16);
            ylabel(axesHandle, 'Copy Number', 'fontSize', 16);
        end

        function plotExpected_Vs_Simulated_Monomer_Expression(sim, axesHandle)
            import edu.stanford.covert.cell.sim.analysis.Constants;

            r = sim.state('Rna');
            m = sim.state('ProteinMonomer');
            c = sim.state('ProteinComplex');

            expression = r.geneExpression(sim.gene.mRNAIndexs) ./ m.halfLives(m.matureIndexs);
            expression = expression / sum(expression);

            monomers = ...
                sum(sum(...
                m.counts(m.matureIndexs,      :, :) + ...
                m.counts(m.inactivatedIndexs, :, :) + ...
                m.counts(m.boundIndexs,       :, :) + ...
                m.counts(m.misfoldedIndexs,   :, :) + ...
                m.counts(m.damagedIndexs,     :, :), 3), 2) + ...
                sum(c.proteinComplexComposition(sim.gene.mRNAIndexs, :, :), 3) * ...
                sum(sum(...
                c.counts(c.matureIndexs,      :, :) + ...
                c.counts(c.inactivatedIndexs, :, :) + ...
                c.counts(c.boundIndexs,       :, :) + ...
                c.counts(c.misfoldedIndexs,   :, :) + ...
                c.counts(c.damagedIndexs,     :, :), 3), 2);

            plot(axesHandle, expression, monomers, '.', 'MarkerSize', 10);
            line([0 1], [0 sum(monomers)], 'Parent', axesHandle, 'Color', 'r');
            xlim([min(expression) max(expression)]);
            ylim([min(monomers) max(monomers)]);
            xlabel(axesHandle, 'Expected Expression', 'fontSize', 16);
            ylabel(axesHandle, 'Simulated Counts', 'fontSize', 16);
        end

        function plotExpected_Vs_Simulated_NTP_Incorporation(sim, axesHandle)
            import edu.stanford.covert.cell.sim.analysis.Constants;

            r = sim.state('Rna');
            m = sim.state('Metabolite');

            nmpComposition = sum(m.nmpComposition, 2);

            ntps = sum(sum(multiprod(r.baseCounts(:, sim.state('Metabolite').nmpIndexs)', r.counts, [1 2], [1 2]), 3), 2);

            plot(axesHandle, nmpComposition, ntps, '.', 'MarkerSize', 20);
            line([0 1], [0 sum(ntps)], 'Parent', axesHandle, 'Color', 'r');
            xlim([min(nmpComposition) max(nmpComposition)]);
            ylim([min(ntps) max(ntps)]);
            xlabel(axesHandle, 'Expected Incorporation', 'fontSize', 16);
            ylabel(axesHandle, 'Simulated Counts', 'fontSize', 16);
        end

        function plotExpected_Vs_Simulated_AA_Incorporation(sim, axesHandle)
            import edu.stanford.covert.cell.sim.analysis.Constants;

            met = sim.state('Metabolite');
            m = sim.state('ProteinMonomer');
            c = sim.state('ProteinComplex');

            aaComposition = sum(met.aaComposition, 2);

            aas = multiprod(m.baseCounts(:, sim.state('Metabolite').aminoAcidIndexs)', m.counts, [1 2], [1 2]) + ...
                multiprod(c.baseCounts(:, sim.state('Metabolite').aminoAcidIndexs)', c.counts, [1 2], [1 2]);
            aas = sum(sum(aas, 3), 2);

            plot(axesHandle, aaComposition, aas, '.', 'MarkerSize', 20);
            line([0 1], [0 sum(aas)],'Parent', axesHandle, 'Color', 'r');
            xlim([min(aaComposition) max(aaComposition)]);
            ylim([min(aas) max(aas)]);
            xlabel(axesHandle, 'Expected Incorporation', 'fontSize', 16);
            ylabel(axesHandle, 'Simulated Counts', 'fontSize', 16);
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function h = heaviside(f)
%HEAVISIDE   Heaviside function of a CHEBFUN.
%   HEAVISIDE(F) returns a CHEBFUN which is 0 when F < 0, +1 when F > 0, and
%   0.5 when F == 0.
%
% See also DIRAC, SIGN.

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

% HEAVISIDE() is basically a wrapper for SIGN().
h = .5*(sign(f) + 1);

end


text:
function [U,Q] = lscm(V,F,b,bc,Aeq,Beq,varargin)
  % LSCM Compute Least Squares Conformal Mapping for mesh
  %
  % U = lscm(V,F,b,bc)
  % U = lscm(V,F,b,bc,Aeq,Beq)
  %
  % Inputs:
  %   V  #V by dim list of rest domain positions
  %   F  #F by 3 list of triangle indices into V
  %   b  #b list of indices of constraint (boundary) vertices
  %   bc  #b by 2 list of constraint positions for b
  %   Aeq   #Aeq by 2*#V matrix of linear equality constraints {[]}
  %   Beq   #Aeq vector of linear equality constraint right-hand sides {[]}
  %   Optional:
  %     'Method' followed by one of the following:
  %        'desbrun'  "Intrinsic Parameterizations of Surface Meshes" [Desbrun
  %          et al. 2002]
  %        'levy'  "Least Squares Conformal Maps for Automatic Texture Atlas
  %          Generation" [Lévy et al. 2002]
  %        {'mullen'}  "Spectral Conformal Parameterization" [Mullen et al. 2008]
  % Outputs:
  %   U  #V by 2 list of new positions
  %   Q  #V*2 by #*2  quadratic coefficients matrix
  %
  % Note: This is the same system as takeo_asap up to a factor of 2.5
  %
  % See also: arap, takeo_arap, takeo_asap
  %

  method = 'mullen';
  % default values
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Method'}, ...
    {'method'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  if nargin<=4
    Aeq = [];
    Beq = [];
  end

  % number of vertices
  n = size(V,1);
  % number of triangles
  nt = size(F,1);
  % number of original dimensions
  dim = size(V,2);

  switch method
  case 'levy'
    % first need to convert each triangle to its orthonormal basis, if coming
    % from 3D
    assert(dim == 2);

    %% Indices of each triangle vertex, I, and its corresponding two neighbors, J
    %% and K
    %I = [F(:,1)];
    %J = [F(:,2)];
    %K = [F(:,3)];

    %X = [V(I,1) V(J,1) V(K,1)];
    %Y = [V(I,2) V(J,2) V(K,2)];

    %WRe = [X(:,3)-X(:,2) X(:,1)-X(:,3) X(:,2)-X(:,1)];
    %WIm = [Y(:,3)-Y(:,2) Y(:,1)-Y(:,3) Y(:,2)-Y(:,1)];

    %% sqrt root of twice the area of each triangle
    %dT = sqrt(doublearea(V,F));

    %% build M matrix, real and imaginary parts
    %II = [1:nt 1:nt 1:nt];
    %JJ = [I;J;K]';
    %VVRe = [WRe(:,1)./dT WRe(:,2)./dT WRe(:,3)./dT];
    %VVIm = [WIm(:,1)./dT WIm(:,2)./dT WIm(:,3)./dT];

    %WWRe = sparse(II,JJ,WRe,nt,n);
    %WWIm = sparse(II,JJ,WIm,nt,n);
    %% These look like blocks in the gradient matrix
    %MRe = sparse(II,JJ,VVRe,nt,n);
    %MIm = sparse(II,JJ,VVIm,nt,n);

    %% build A matrix
    %A = [MRe -MIm; MIm MRe];

    %% quadratic system matrix
    %Q = A'*A;

    % Or equivalently

    % compute gradient matrix
    G = grad(V,F);

    % Extract each coordinate's block
    Gx = G(1:nt,:);
    Gy = G(nt+(1:nt),:);

    % Triangle areas
    TA = repdiag(diag(sparse(doublearea(V,F))/2),2);

    % Build quadratic coefficients matrix
    Q = [Gx -Gy;Gy Gx]'*TA*[Gx -Gy;Gy Gx];

    % solve
    U = min_quad_with_fixed(Q,zeros(2*n,1),[b(:);b(:)+n],bc(:));
    % reshape into columns
    U = reshape(U,n,2);
  case 'mullen'
    A = vector_area_matrix(F);
    L = repdiag(cotmatrix(V,F),2);
    Q = -L - 2*A;
    if ~exist('b','var') || isempty(b)
      M = repdiag(massmatrix(V,F),2);
      [EV,ED] = eigs(Q,M,2,'sm');
      U = EV(:,2);
    else
      U = min_quad_with_fixed(Q,zeros(2*n,1),[b(:);b(:)+n],bc(:),Aeq,Beq);
    end
    % reshape into columns
    U = reshape(U,n,2);
  case 'desbrun'
    error('not implemented.');
    % This implements the Dirichlet + Chi energies but usually when people
    % refer to the [Desbrun et al. 2002] paper they mean the [Mullen et
    % al. 2008] interpretation: zero chi energy + ""natural"" boundary
    % conditions.
    lambda = 1;
    mu = 1;
    L = cotmatrix(V,F);
    % chi energy
    C = cotangent(V,F);
    l = edge_lengths(V,F);
    X = sparse( ...
      F(:,[1 2 3 1 2 3]), ...
      F(:,[2 3 1 3 1 2]), ...
      C(:,[2 3 1 3 1 2])./ ...
      l(:,[3 1 2 2 3 1]), ...
      size(V,1),size(V,1));
    Q = repdiag(lambda*L+mu*X,2);
  otherwise
    error('unsupported method: %s',method);
  end
end

author: alecjacobson

repo: gptoolbox


text:
classdef MimTemplateCallback < CoreBaseClass
    % MimTemplateCallback. Used by plugins to fetch empty template images
    %
    %     Template images are empty images with metadata matching the current
    %     image. They are used to create images that will be compatible with
    %     other images from the same dataset. The metadata image size and origin
    %     depends on the chosen context.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Access = private)
        LinkedDatasetChooser  % Sends the API calls to the correct dataset
        DatasetStack       % Handle to the current call stack for the primary dataset
        Reporting
    end

    methods
        function obj = MimTemplateCallback(linked_dataset_chooser, dataset_call_stack, reporting)
            obj.DatasetStack = dataset_call_stack;
            obj.LinkedDatasetChooser = linked_dataset_chooser;
            obj.Reporting = reporting;
        end

        % Returns an image template for the requested context
        function template = GetTemplateImage(obj, context, dataset_name)
            if nargin < 2
                context = [];
            end
            if nargin < 3
                dataset_name = [];
            end
            template = obj.LinkedDatasetChooser.GetDataset(obj.Reporting, dataset_name).GetTemplateImage(context, obj.DatasetStack, obj.Reporting);
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function d = gcdfactor(A)
d = A(1);
for n = 2:numel(A)
    if d == 1
        return
    else
        d = gcd(d,A(n));
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%% DEMO 17: Detector Rotation
%
%
%
%  Some systems have a slight rotation of the detector due to mechanicall
%  inacuracies.
%
%  According to the article "A geometric calibration method for cone beam
%  CT systems" (DOI: 10.1118/1.2198187), only Roll needs to be corrected
%  for in the algorithmic part, as the other 2 can be easily ignored if
%  "sufficiently small". In TIGRE we decided to leave that to the users
%  discretion and implemented the 3 possible detector rotation, per angle
%  if needed.
%
%  This demo shows how to use it.
%
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri
%--------------------------------------------------------------------------
%% Initialize

clear;
close all;
%% Define Geometry
%
% VARIABLE                                   DESCRIPTION                    UNITS
%-------------------------------------------------------------------------------------
geo.DSD = 1536;                             % Distance Source Detector      (mm)
geo.DSO = 1000;                             % Distance Source Origin        (mm)
% Detector parameters
geo.nDetector=[512; 512];					% number of pixels              (px)
geo.dDetector=[0.8; 0.8]; 					% size of each pixel            (mm)
geo.sDetector=geo.nDetector.*geo.dDetector; % total size of the detector    (mm)
% Image parameters
geo.nVoxel=[128;128;128];                   % number of voxels              (vx)

% a bit smaller than usual because the demo includes a very big detector
% angle for showcase
geo.sVoxel=[256;256;256]/2;                 % total size of the image       (mm)

geo.dVoxel=geo.sVoxel./geo.nVoxel;          % size of each voxel            (mm)
% Offsets
geo.offOrigin =[0;0;0];                     % Offset of image from origin   (mm)
geo.offDetector=[0; 0];                     % Offset of Detector            (mm)

% Auxiliary
geo.accuracy=0.5;                           % Accuracy of FWD proj          (vx/sample)

%% Geometry has an additional field for the rotation.
%
%
angles=linspace(0,2*pi,100);

% lets define 3 angles, with big variability in each direction

roll=angles;
pitch=0.7*linspace(0,1,size(angles,2));
yaw=0.7*linspace(0,1,size(angles,2));

% Fill rotDetector. It can also be a 3x1, if its not angle dependent. e.g.
% geo.rotDetector=[pi/5; 0; pi/10];

geo.rotDetector=[roll;pitch;yaw];

%% Load data and generate projections
% define angles
% Load thorax phatom data
head=headPhantom(geo.nVoxel);
% generate projections
projections=Ax(head,geo,angles,'interpolated');

%% Lets plot the projections with a rotated detector
plotProj(projections,angles);

%% lets reconstruct with and without detector rotation
%
imgRotDet=OS_SART(projections,geo,angles,50);

% No rotation
geo.rotDetector=[0;0;0];
projections2=Ax(head,geo,angles,'interpolated');

imgnoRot=OS_SART(projections2,geo,angles,50);

% %% Plot to show that indeed the reconstruction is right
%
%
%
plotImg([imgRotDet imgnoRot] ,'dim',3)

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
%  demo_RoPS_FeatureMatching_Mesh.m
%  Author: Yulan Guo {yulan.guo@nudt.edu.cn}
%  NUDT, China & CSSE, UWA, Australia
% This function performs feature matching on two input meshes to obtain feature
% correspondences

close all;
clc;
clear all;

keypntNum = 1000;

load ../../StanfordData/Bunny
SelIdx = [ 1 3];
pcData0 = bunny{SelIdx(1)}';
tmp = pcdownsample( pointCloud(pcData0'), 'gridAverage', 0.0015 );
pcData0 = tmp.Location';
pcData1 = bunny{SelIdx(2)}';
tmp = pcdownsample( pointCloud(pcData1'), 'gridAverage', 0.0015 );
pcData1 = tmp.Location';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load data\pointcloud_view1;
pointcloud = pcData0';
% %============================transform a pointcloud into a triangular mesh============================%
mesh = pointCloud2mesh(pointcloud,[0 0 1],0.4);                                         %other methods can also be used to perform triangulation
%============================preprocessing============================%
out = preprocessingFunc(mesh);
mesh.faceCenter = out.centroid;
mesh.faceArea = out.area;
mesh.res = out.res ;
%============================detect keypoints============================%
%keypoints are randomly seleted in this demo, any other 3D keypoint detection methods can be used
temp = randperm(length(mesh.vertices));
mesh.keypntIdx = temp(1:keypntNum);
%============================extract RoPS features at the keypoints on a mesh============================%
para.RoPS_nbSize = 15*mesh.res;
para.RoPS_binSize = 5;
para.RoPS_rotaSize = 3;
mesh.LRF =  LRFforMeshFunc(mesh, mesh.keypntIdx, para.RoPS_nbSize);
disp('LRFs calculated');
RoPS = RoPSFunc(mesh, para.RoPS_nbSize, para.RoPS_binSize, para.RoPS_rotaSize,mesh.LRF);
mesh.RoPS = RoPS;
disp(['RoPS features generated']);
mesh1 = mesh;
mesh1Features = [];
for keypntIdx = 1:keypntNum
    temp = trans2Dto1DFunc(mesh.RoPS{keypntIdx});
    mesh1Features = [mesh1Features; temp];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load data\pointcloud_view2;
pointcloud = pcData1';
% %============================transform a pointcloud into a triangular mesh============================%
mesh = pointCloud2mesh(pointcloud,[0 0 1],0.4);                                         %other methods can also be used to perform triangulation
%============================preprocessing============================%
out = preprocessingFunc(mesh);
mesh.faceCenter = out.centroid;
mesh.faceArea = out.area;
mesh.res = out.res ;
%============================detect keypoints============================%
%keypoints are randomly seleted in this demo, any other 3D keypoint detection methods can be used
temp = randperm(length(mesh.vertices));
mesh.keypntIdx = temp(1:keypntNum);
%============================extract RoPS features at the keypoints on a mesh============================%
para.RoPS_nbSize = 15*mesh.res;
para.RoPS_binSize = 5;
para.RoPS_rotaSize = 3;
mesh.LRF =  LRFforMeshFunc(mesh, mesh.keypntIdx, para.RoPS_nbSize);
disp('LRFs calculated');
RoPS = RoPSFunc(mesh, para.RoPS_nbSize, para.RoPS_binSize, para.RoPS_rotaSize,mesh.LRF);
mesh.RoPS = RoPS;
disp(['RoPS features generated']);
mesh2 = mesh;
mesh2Features = [];
for keypntIdx = 1:keypntNum
    temp = trans2Dto1DFunc(mesh.RoPS{keypntIdx});
    mesh2Features = [mesh2Features; temp];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %============================feature matching============================%
NNDRthreshold = 0.9;
corNum = 0;
kdtreeMesh1Features = KDTreeSearcher(mesh1Features,'Distance','euclidean');
for keypntIdx1 = 1:size(mesh2Features,1)
    [idxSort,distSort] = knnsearch(kdtreeMesh1Features, mesh2Features(keypntIdx1,:),'k',2,'Distance','euclidean');
    IDX = idxSort(1);
    if distSort(1)/distSort(2)<=NNDRthreshold
        corNum = corNum+1;
        corPntIdx(corNum,:) = [IDX, keypntIdx1];
        featureDis(corNum) = distSort(1);
    end
end
showCorresFunc(mesh1, mesh2, mesh1.keypntIdx(corPntIdx(:,1)), mesh2.keypntIdx(corPntIdx(:,2)), [0,200,0]);

%============================links============================%
%we may find more test datasets via the following links
% url = 'https://sites.google.com/site/yulanguo66/research-resources/3d-object-recognition-datasets';
% web(url,'-browser')

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR


text:
function size = pivparpool(cmd,core)
% parpool size: pivparpool('size')
% parpool close: pivparpool('close')
% parpool open: pivparpool('open', 2)

if nargin < 1
	error('Requires at least one input argument')
end
if nargin == 1
    if verLessThan('matlab', '8.2')
        if strcmp(cmd,'size')
            size = matlabpool('size');
        elseif strcmp(cmd,'close')
            matlabpool close;
        else
            error('Wrong input arguments.')
        end
    else
        poolobj = gcp('nocreate');
        if strcmp(cmd,'size')
            if isempty(poolobj)
                size = 0;
            else
                size = poolobj.NumWorkers;
            end
        elseif strcmp(cmd,'close')
            delete(poolobj);
        else
            error('Wrong input arguments.')
        end
    end
elseif nargin == 2
    if verLessThan('matlab', '8.2')
        if strcmp(cmd,'open')
            matlabpool('open','local',core);
        else
            error('Wrong input arguments.')
        end
    else
        if strcmp(cmd,'open')
            parpool('local',core);
        else
            error('Wrong input arguments.')
        end
    end
elseif nargin > 2
	error('Too many input arguments.')
end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function b = GetB_Values(protocol)
% Computes the b value for each measurement in the protocol
%
% b = GetB_Values(protocol)
% returns an array of b values, one for each measurement defined
% in the list of measurements in the protocol.
%
% author: Daniel C Alexander (d.alexander@ucl.ac.uk)
%

GAMMA = 2.675987E8;
if(strcmp(protocol.pulseseq, 'PGSE') || strcmp(protocol.pulseseq, 'STEAM'))
    modQ = GAMMA*protocol.smalldel.*protocol.G;
    diffTime = protocol.delta - protocol.smalldel/3;
    b = diffTime.*modQ.^2;

elseif(strcmp(protocol.pulseseq, 'DSE'))
    b = getB_ValuesDSE(protocol.G, protocol.delta1, protocol.delta2, protocol.delta3, protocol.t1, protocol.t2, protocol.t3);

elseif(strcmp(protocol.pulseseq, 'OGSE'))
    b = GetB_ValuesOGSE(protocol.G, protocol.delta, protocol.smalldel, protocol.omega);
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function cordic_test001 ( )

%*****************************************************************************80
%
%% TEST001 demonstrates the use of COSSIN_CORDIC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 June 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST001:\n' );
  fprintf ( 1, '  COSSIN_CORDIC computes the cosine and sine of an angle\n' );
  fprintf ( 1, '  using the CORDIC algorithm.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      A    N           Cos(A)     Cos(A)      Difference\n' );
  fprintf ( 1, '                     Tabulated   CORDIC\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, a, c1 ] = cos_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '\n' );

    for n = 0 : 5 : 25

      v = cossin_cordic ( a, n );
      c2 = v(1);
      d = c1 - c2;

      fprintf ( 1, '  %12f  %4d  %16.8f  %16.8f  %9e\n', a, n, c1, c2, d );

    end

  end

  return
end


text:
function [J, S] = GALA_invert(BF,tS)

YTY   = sparse(0);

Nl = tS.Nl;

for p=1:Nl

    BFp = load(BF{p});

    % load data
    D{p} = spm_eeg_load(BFp.data.D);

    % load gain matrix
    modality='MEG'; % temporary - should be replaced by tS.modality
    L1 = BFp.sources.L.(modality);
    Nd = length(L1);
    Lp{p} = [];
    for i = 1:Nd
        lf = L1{i};
        Lp{p} = cat(2, Lp{p}, lf);
    end
    % eliminate low SNR spatial modes
    U = spm_svd((Lp{p}*Lp{p}'),exp(-16));
    ULp{p} = U'*Lp{p};
    Nm(p) = size(ULp{p},1);
    Up{p} = U;

    trial = D{p}.condlist;      % get conditions list
    Nt = 2;                     % temporary invert only first one
    Nb = size(D{p},2);          % Number of time bins

    % get (spatially aligned) data
    YY = 0; N=0;
    for j = 1:Nt                                % loop over conditions
        Y{p,j}  = Up{p}'*D{p}(:,:,j);
        YY    = YY + Y{p,j}'*Y{p,j};
        N     = N + Nb;
    end

    YTY        = YTY + YY;

end

% Initialize matrices comtrolling covariance reparametrization
Is = 1:Nd;
for p=1:Nl
    Isp{p} = 1:Nd;
end
lIs = 1:Nd*Nl;
Vs = speye(Nd*Nl);

mNm = mean(Nm);
scale = 1/sqrt(trace(YTY)/(N*mNm));
YTY = YTY*(scale^2);

% temporal projector
[V E]  = spm_svd(YTY,exp(-8));              % get temporal modes
E      = diag(E)/trace(YTY);                % normalise variance
Nr  = length(E);                            % number of temporal modes
S   = V(:,1:Nr);                            % temporal modes
VE  = sum(E(1:Nr));                         % variance explained

fprintf('Using %i temporal modes, ',Nr)
fprintf('accounting for %0.2f percent average variance\n',full(100*VE))

% scale and align data (spatial and temporary modes)
for p=1:Nl
    MY=[];
    for j = 1:Nt                % loop over Nt conditions
        MY{j}   = Y{p,j}*S*scale;
    end
    UYp{p} = spm_cat(MY);
end

% prepare long data and leadfields and channel noise components
UL=[]; Qe=[];
for p=1:Nl
    UL = blkdiag(UL,ULp{p});
    Qep{p} = Up{p}'*Up{p};
    Qe = full(blkdiag(Qe,Qep{p}));
end
UY = spm_cat(UYp');

Scale = sqrt(trace(UL*UL')/mNm);
UL = UL/Scale;

% prepare smoothing kernel - it could be done as in spm but this scheme
% applicable for FreeSurfer meshes
% -----------------------------------------------------------------------
radius = 4;

vert  = D{1}.inv{D{1}.val}.mesh.tess_mni.vert;
face  = D{1}.inv{D{1}.val}.mesh.tess_mni.face;

AA=[];
AAA     = spm_mesh_distmtx(struct('vertices',vert,'faces',face),0);
% find subsets of certain distance
AA{1}=speye(length(vert));
AA{2}=AAA;
B = AA{1}+AA{2};
A = B;
for i=2:radius-1
    in = find(AAA^i);
    Bn=sparse(length(vert),length(vert));
    Bn(in)=1;
    AA{i+1}=full(Bn-B);
    B=Bn;
end

y=[1 0.8 0.37 0.2];

QG = sparse(length(vert),length(vert));
for i=1:radius
    QG=QG + AA{i}*y(i);
end
QG =sparse(QG(:,:));
K    = sparse(QG.*(QG > exp(-8)));
clear AA AAA B Bn QG
% -----------------------------------------------------------------------

% covariance of the data
YY = UY*UY';

% noise covariance component
Qn{1} = Qe; % it's the first one for channel noise

% template for source-noise covariance
sQt = kron(speye(Nl),speye(Nd));

thresh = 0.5;       % threshold

for it=1:tS.iter

    % starting from 2nd iteration add source-noise covariance matrices
    if it>1
        sQn = Vn*sQt*Vn; %Vn - noise-vertices matrix, defined below
        Qn{it} = UL*sQn*UL';
    end

    % first ROI covariance matrix - strong correlation between subjects
    psQ1 = kron(ones(Nl),K);
%     psQ1 = kron(ones(Nl),speye(Nd));
    sQ1 = Vs*psQ1*Vs;
    Q1 = full(UL*sQ1*UL');

    % second ROI covariance matrix - subjects specific activity
    psQ2 = kron(speye(Nl),K);
%     psQ2 = kron(speye(Nl),speye(Nd));
    sQ2 = Vs*psQ2*Vs;
    Q2 = full(UL*sQ2*UL');

    Qs = {Q1 Q2};
%     Qs = {Q2};

    % ML estimation of hyperparameters
    [Cy,h,Ph,F] = spm_reml_sc(YY,[],[Qn Qs],1,-4,16);
    % take ROI components and discard noise components
    DD = h(end-1)*sQ1 + h(end)*sQ2;
%     DD = h(end)*sQ2;
    Cy = full(Cy);

    M = DD*UL'/Cy;

    % long (all subjects) source activity
    J = M*UY;

    % goodness of fit
    SSR  = sum(var((UY - UL*J),0,2));
    SST  = sum(var( UY,0,2));
    R2   = 100*(SST - SSR)/SST;
    fprintf('\nPercent variance explained GALA %.2f\n',full(R2));
    fprintf('F = %.10f\n',full(F));
%     fprintf('L = %.10f\n',full(-spm_logdet(Cy)- trace(UY'*inv(Cy)*UY)));
%     fprintf('ML fit = %.10f\n',-trace(UY'*inv(Cy)*UY));
%     fprintf('ML det = %.10f\n\n',full(-spm_logdet(Cy)));

    % idea - get vertices with max covariance both within and between
    % subjects. As the result they should not be the same (indices) for
    % different subjects, so correlated between subjects part of covarince
    % prior is overlapping but not identical patches

    % pay attention - it's to exclude disconected vertices from Jcov
    % it seems ssQ1 = spones(sQ1) works better than simple sQ1
    % may be because in Jcov.*sQ1 there is double attenuation of tails

    if it<tS.iter

        ssQ1 = spones(sQ1);

        % real calculation for big matrices
        J = full(J);

        lJcov = zeros(1,Nl*Nd);
        for i=1:Nd*Nl
            if any(J(i,:))
                Jcovi = J*J(i,:)';
                sJcovi = Jcovi.*ssQ1(:,i);
                lJcov(i) = squeeze(sum(sJcovi));
            end
        end

        % reparametrization of covariance matrices
        lIso = lIs;
        lIs = find(lJcov>quantile(lJcov,thresh));   % ROI vertices update
        lnIs = setdiff(lIso,lIs);                   % noise vertices update

        Vs = zeros(Nd*Nl,1);
        Vn = zeros(Nd*Nl,1);
        Vs(lIs) = 1;
        Vn(lnIs) = 1;
        Vs = spdiags(Vs,0,Nd*Nl,Nd*Nl);     % ROI vertices matrix update
        Vn = spdiags(Vn,0,Nd*Nl,Nd*Nl);     % noise vertices matrix update

        thresh = thresh+(1-thresh)/2;               % threshold update

    end

    % display

%     MM.vertices = vert;
%     MM.faces = face;
%
%     ahp=[];
%     figure;
%     for p=1:Nl
%         lJcovp{p} = lJcov(1+(p-1)*Nd:Nd+(p-1)*Nd);
%
%         ahp(p)=subplot(2,2,p);
%         srcs_disp1 = lJcovp{p}/max(abs(lJcovp{p}));
%
%         cla; axis off
%
%         fig1 = patch('vertices',MM.vertices,'faces',MM.faces,'FaceVertexCData',srcs_disp1');
%
%         set(fig1,'FaceColor',[.5 .5 .5],'EdgeColor','none');
%         shading interp
%         lighting gouraud
%         %camlight
%         zoom off
%         lightangle(0,270);lightangle(270,0),lightangle(90,0),lightangle(0,45),lightangle(0,135);
%         material([.1 .1 .4 .5 .4]);
%         %view(140,15);
%         caxis([0 1])
%         colormap(jet);
%
%     end
%
%     hlink = linkprop(ahp, {'CameraPosition','CameraUpVector'});
%     key = 'graphics_linkprop';
%     setappdata(ahp(1),key,hlink);

end

end

author: spm


text:
function [x,ebuf,Zi] = celpsyn(cb,kappa,k,theta0,P,b,ebuf,Zi)
%  celpsyn --> CELP synthesizer (decoder).
%
%    [x,ebuf,Zi] = celpsyn(cb,kappa,k,theta0,P,b,ebuf,Zi)
%
%    The CELP synthesizer consists of the cascade of the pitch filter
%    and the vocal-tract filter with excitation signal taken from the
%    codebook.
%
%      ------------  Gain, theta0 ----------------      ---------
%      | Gaussian |          |    |      1       |      |   1   |
%      | codebook |--------->X--->| ------------ |----->| ----- |--->
%      |   cb     | rho_k(n)      | 1 - b*z^(-P) | e(n) |  A(z) | x(n)
%      ------------               ----------------      ---------
%
%    First, the Gaussian codebook given by the L-by-K matrix cb, and
%    the excitation parameters k, theta0, P, and b, are used to generate
%    the excitation sequence, e(n). This is done in blocks of length L,
%    so if the length of the parameter vector k is J, then L*J samples
%    of e(n) are generated. The vector ebuf contains previous excitation
%    samples due to the memory hangover in the pitch filter, and the
%    length of this vector must at least be the max possible pitch period.
%
%    Then, the excitation sequence, e(n), is filtered by the vocal-tract
%    filter, 1/A(z), where the coefficients are obtained from the
%    reflection coefficients, kappa. Zi is the memory hangover in this
%    filter.
%

[L,K] = size(cb);                       % Block length and codebook size.
F = length(ebuf);                       % No. of previous excitation samples.
J = length(k);                          % No. blocks per frame.
N = L*J;                                % Frame length.

e = zeros(N,1);

for (j=1:J)
  n = (j-1)*L+1:j*L;

 % Find the signal e(n) based on the parameters b, P, theta0, and k.
 if (P(j) < L)
    Zp   = b(j)*ebuf(F-P(j)+1:F);
    e(n) = filter(1,[1 zeros(1,P(j)-1) -b(j)],theta0(j)*cb(:,k(j)),Zp);
  else
    e(n) = theta0(j)*cb(:,k(j)) + b(j)*ebuf(F-P(j)+1:F-P(j)+L);
  end
  ebuf = [ebuf(L+1:F); e(n)];           % Update e(n) buffer.
end

a = rf2lpc(kappa);                      % Convert kappa to a parameters.
[x,Zi] = filter(1,[1; -a],e,Zi);        % Vocal-tract filter.


text:
function [ x, w ] = laguerre_ek_compute ( n )

%*****************************************************************************80
%
%% LAGUERRE_EK_COMPUTE: Laguerre quadrature rule by the Elhay-Kautsky method.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 April 2011
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Parameters:
%
%    Input, integer N, the order.
%
%    Output, real X(N), the abscissas.
%
%    Output, real W(N), the weights.
%

%
%  Define the zero-th moment.
%
  zemu = 1.0;
%
%  Define the Jacobi matrix.
%
  bj = zeros ( n, 1 );

  for i = 1 : n
    bj(i) = i;
  end

  x = zeros ( n, 1 );
  for i = 1 : n
    x(i) = 2 * i - 1;
  end

  w = zeros ( n, 1 );
  w(1) = sqrt ( zemu );
%
%  Diagonalize the Jacobi matrix.
%
  [ x, w ] = imtqlx ( n, x, bj, w );

  w(1:n) = w(1:n).^2;

  return
end


text:
function butterfly
% BUTTERFLY Visualize payoff schedule of option strategy
% PURPOSE : The function aims to assist students and instructors of
%           finance in visualizing payoffs of simple option strate-
%           gies. Portfolios of n < 9 securities including a stock,
%           a zero-coupon bond,  a forward contract, and a European
%           call or  put option, are allowed. Terminal  payoff  and
%           terminal payoff adjusted for accrued time-0 expense are
%           plotted. Discounting and valuation follow Black-Scholes.
% EXAMPLE : buttefly (Oh, the FEX code metrics..)
% AUTHOR  : Dimitri Shvorob, dimitri.shvorob@vanderbilt.edu, 11/5/07

fw = 560;     fh = 350;
bx = .62*fw;  by = .07*fh;
ew = .08*fw;  eh = .5*ew;
tx = .1*ew;   ty = .25*eh;
bw = 2*ew+tx; bh = .25*bw;
wy = 1.5*bh;
x2 = bx+bw+tx;
x3 = x2+ew+tx;

n = 8;  % max number of securities
st = {' Stock', ...
      ' Bond', ...
      ' Forward', ...
      ' Call option', ...
      ' Put option'};
pT = {'S', ...
      'K', ...
      '(S-K)', ...
      'max(S-K,0)', ...
      'max(K-S,0)'};
p0 = {'(-S0)', ...
      '(-exp(-R*T)*K)', ...
      '(exp(-R*T)*K-S0)', ...
      '(-blscall(S0,K,R,T,sigma))', ...
      '(-blsput(S0,K,R,T,sigma))'};

ax = setupAxisAndButtons;
hp = setupParameterSelection;
[hs,hk,hn] = setupSecuritySelection;
clearSecuritySelection

function[ax] = setupAxisAndButtons
ax = .05; ay = (by + bh + wy)/fh;
ah = .7;  aw = (bx - ew)/fw;
figure('Name','Butterfly: visualize payoffs of an option strategy', ...
       'NumberTitle','off', 'Position',[232 246 fw fh], ...
       'Color','white','Resize','off','MenuBar','none');
ax = axes('Position' ,[ax ay aw ah],'FontSize',8);
guiinput('pushbutton',[bx by bw bh],'Plot' ,@evalSecuritySelection);
guiinput('pushbutton',[x2 by bw bh],'Reset',@clearSecuritySelection);
end

function[h] = setupParameterSelection
tw = .75*bw;
iw = .75*ew;
tt = .10*ew;
v = {'90' ,'100','0.05';
     '110','0.1','1'};
l = {'xmin','Stock price','Interest rate';
     'xmax','Volatility','Time to expiry'};
h = nan(2,3);
for i = 1:2
   y = .4*by + (i-1)*(eh + ty);
   for j = 1:3
       x = (j-1)*(tw + tx + iw + tt);
       guilabel([x y-4 tw eh],l{i,j});
       h(i,j) = guiinput('edit',[x+tw+tx y iw eh],v{i,j},@checkIfPositive);
   end
end
end

function[hs,hk,hn] = setupSecuritySelection
bx = .62*fw;
by = .07*fh;
[hs,hk,hn] = deal(nan(n,1));
for i = 1:n
    y = by + bh + wy + (i-1)*(eh + ty);
    hs(i) = guiinput('popupmenu',[bx y-3 bw eh+2],st,'');
    hk(i) = guiinput('edit',[x2 y ew eh],'',@checkIfPositive);
    hn(i) = guiinput('edit',[x3 y ew eh],'',@checkIfNumber);
end
y = y + ty + 1.2*eh;
guilabel([bx y bw eh],'Security type')
guilabel([x2 y-10 ew 2*eh],{'Strike /','Face value'})
guilabel([x3 y ew eh],'Units')
end

function[h] = guiinput(style,position,string,callback)
h = uicontrol('Style',style,'Position',position,'BackgroundColor','white',...
              'FontSize',8,'String',string,'Callback',callback);
end

function guilabel(position,string)
guiinput('text',position,string,'');
end

function checkIfPositive(hObject,eventdata)      %#ok
i = str2double(get(hObject,'String'));
if isnan(i) || i < 0, beep, end
end

function checkIfNumber(hObject,eventdata)        %#ok
i = str2double(get(hObject,'String'));
if isnan(i), beep, end
end

function evalSecuritySelection(hObject,eventdata)%#ok
cla(ax)
[p,e] = evalSchedule;
plotSchedule(p,e)
end

    function[p,e] = evalSchedule
    p = ''; e = '';
    for i = 1:n
        m = get(hn(i),'String');
        if ~isempty(m)
           if strcmp(m(1),'-')
              m = ['(' m ')'];
           end
           s = get(hs(i),'Value');
           k = get(hk(i),'String');
           p = [p '+(' m '*'           strrep(pT{s},'K',k)  ')'];
           e = [e '+(exp(R*T)*(' m '*' strrep(p0{s},'K',k) '))'];
        end
    end
    v = {'S0','R'; 'sigma','T'};
    for i = 1:2
        for j = 1:2
            m = get(hp(i,j+1),'String');
            p = strrep(p,v{i,j},m);
            e = strrep(e,v{i,j},m);
        end
    end
    end

    function plotSchedule(p,e)
    try
       a = str2double(get(hp(1,1),'String'));
       b = str2double(get(hp(2,1),'String'));
       x = linspace(a,b,100);
       p = feval(eval(['@(S) ' p]),x);
       e = feval(eval(['@(S) ' e]),x);
       plot(ax,x,p  ,'Color',[.1 .1 .1],'LineWidth',1.5); hold on
       plot(ax,x,p+e,'Color',[.8 .8 .8],'LineWidth',1.5); hold on
       legend(ax,'\pi_{\itT}','\pi_{\itT} + \ite^{rT}\rm\pi_{0}','Location','Best')
       legend('boxoff')
       plot(ax,[a b],[0 0],'k:')
    catch
       beep
    end

end

function clearSecuritySelection(hObject,eventdata) %#ok
% Clear/initialize security selections
m = fliplr([1:5 4 4 5]);
for i = 1:n
    set(hn(i),'String','')
    set(hk(i),'String','')
    set(hs(i),'Value',m(i))
end
cla(ax)
end

function[c] = blscall(S,K,r,T,sigma)
% Black-Scholes call price
d = (log(S/K) + (r + [1 -1]*sigma^2/2)*T)/(sigma*sqrt(T));
c = S*normcdf(d(1)) - K*exp(-r*T)*normcdf(d(2));
end

function[p] = blsput(S,K,r,T,sigma)
% Black-Scholes put price
d = -(log(S/K) + (r + [1 -1]*sigma^2/2)*T)/(sigma*sqrt(T));
p = K*exp(-r*T)*normcdf(d(2)) - S*normcdf(d(1));
end

end


text:
classdef prtClassAdaBoostFastAuc < prtClassAdaBoost
    %prtClassAdaBoostFastAuc AdaBoost classifier (fast training)
    %
    %   prtClassAdaBoostFastAuc is a version of prtClassAdaBoost that can
    %   be trained significantly more quickly than prtClassAdaBoost.
    %   prtClassAdaBoostFastAuc acheives this by assuming a linear
    %   classifier when picking the feature to be used in the weak learner.
    %   Unlike regular adaBoost, where the weak learner is trained and
    %   evaluated on each feature, in prtClassAdaBoostFastAuc, the feature
    %   is selected using simple ROC metrics, and this feature is used to
    %   train the weak learner for the current iteration.  This can be
    %   significantly faster than prtClassAdaBoost when the base learner is
    %   slow, or there are a very large number of features.
    %
    %   a = prtClassAdaBoostFastAuc;
    %   a = a.train(prtDataGenBimodal);
    %   plot(a)
    %
    %    See also: prtClassAdaBoost

    methods
        % Constructor
        function self = prtClassAdaBoostFastAuc(varargin)
            self = prtUtilAssignStringValuePairs(self,varargin{:});
        end
    end

    methods (Access = protected, Hidden = true)

        function self = trainAction(self,dataSet)

            d = ones(dataSet.nObservations,1)./dataSet.nObservations;

            classifier = self.baseClassifier + prtDecisionBinaryMinPe;
            classifier.verboseStorage = false;
            %classifier.verboseFeatureNames = false;

            y = double(dataSet.getTargetsAsBinaryMatrix);
            y = y(:,2);
            y(y == 0) = -1;

            globalX = dataSet.getX;
            globalY = dataSet.getY;
            rocFlip = ones(1,dataSet.nFeatures);

            for t = 1:self.maxIters
                if t == 1
                    dataSetBootstrap = dataSet;
                    if self.downSampleBootstrap
                        dataSetBootstrap = dataSetBootstrap.bootstrap(self.downSampleBootstrap);
                    end
                else
                    dataSetBootstrap = dataSet.bootstrap(dataSet.nObservations,d);
                    if self.downSampleBootstrap
                        dataSetBootstrap = dataSetBootstrap.bootstrapByClass(self.downSampleBootstrap);
                    end
                end

                localX = dataSetBootstrap.getX;
                localY = dataSetBootstrap.getY;
                if length(unique(localY)) < 2
                    %d told us not to re-sample
                    return;
                end

                pe = nan(1,size(localX,2));
                 if self.verbose
                    fprintf('\nIter: %d ',t);
                end
                for feature = 1:dataSet.nFeatures
                    if ~mod(feature,1000) & self.verbose
                        %disp(feature./dataSet.nFeatures);
                        if feature == 1000
                            fprintf('\n');
                        end
                        fprintf('.');
                    end

                    %Evaluate the ROC curve for a feature:
                    [pf,pd,tau,auc] = prtClassAdaBoost.fastRoc(localX(:,feature)*rocFlip(feature),localY);
                    if auc < .5
                        rocFlip(feature) = -1*rocFlip(feature);
                        %try again, flip the ROC
                        [pf,pd,tau] = prtClassAdaBoost.fastRoc(localX(:,feature)*rocFlip(feature),localY);
                    end

                    %Figure out the optimal expected threshold:
                    tempPe = prtUtilPfPd2Pe(pf,pd);
                    [~,ind] = min(tempPe);
                    threshold = tau(ind);

                    %Use the threshold to estimate pe(feature) from all the
                    %data
                    decisions = globalX(:,feature)*rocFlip(feature) >= threshold;
                    correctLogical = decisions == (globalY == 1);
                    pe(feature) = sum(double(~correctLogical).*d);
                end
                [minDeltaPe,minInd] = max(abs(pe-.5));

                if minDeltaPe < self.deltaPeThreshold
                    return;
                else

                    baseClass = self.baseClassifier;
                    baseClass.verboseStorage = false;
                    algo = prtFeatSelStatic('selectedFeatures',minInd) + baseClass + prtDecisionBinaryMinPe;
                    algo.verboseStorage = false;
                    %algo.verboseFeatureNames = false;

                    if t == 1
                        self.classifierArray = train(algo,dataSetBootstrap);
                    else
                        self.classifierArray(t) = train(algo,dataSetBootstrap);
                    end

                    yOut = run(self.classifierArray(t),dataSet);
                    wrong = yOut.getX ~= yOut.getY;
                    eps_t = sum(d.*wrong);

                    if eps_t == 0 %nothing wrong!
                        self.alpha(t) = 10;
                    else
                        self.alpha(t) = 1/2*log((1-eps_t)/eps_t);
                    end

                    h = double(yOut.getObservations);
                    h(h == 0) = -1;

                    d = d.*exp(-self.alpha(t).*y.*h);

                    if sum(d) == 0
                        return;
                    end

                    d = d./sum(d);

                end

            end
        end

    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function xval = r8vec_even_select ( xlo, xhi, n, ival )

%*****************************************************************************80
%
%% R8VEC_EVEN_SELECT returns the I-th of N evenly spaced values in [ XLO, XHI ].
%
%  Formula:
%
%    XVAL = ( (N-IVAL) * XLO + (IVAL-1) * XHI ) / dble ( N - 1 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 November 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real XLO, XHI, the low and high values.
%
%    Input, integer N, the number of values.
%
%    Input, integer IVAL, the index of the desired point.
%    IVAL is normally between 1 and N, but may be any
%    integer value.
%
%    Output, real XVAL, the IVAL-th of N evenly spaced values
%    between XLO and XHI.
%
%    Unless N = 1, X(1) = XLO and X(N) = XHI.
%
%    If N = 1, then X(1) = 0.5*(XLO+XHI).
%
  if ( n == 1 )

    xval = 0.5 * ( xlo + xhi );

  else

    xval = ( ( n - ival     ) * xlo   ...
           + (     ival - 1 ) * xhi ) ...
           / ( n        - 1 );

  end

  return
end


text:
function mpc = apply_changes(label, mpc, chgtab)
%APPLY_CHANGES  Applies a set of changes to a MATPOWER case
%   mpc_modified = apply_changes(label, mpc_original, chgtab)
%
%   Applies the set of changes identified by LABEL to the case in MPC, where
%   the change sets are specified in CHGTAB.
%
%   LABEL is an integer which identifies the set of changes of interest
%
%   MPC is a MATPOWER case struct with at least fields bus, gen and branch
%
%   CHGTAB is the table of changes that lists the individual changes required
%   by each change set, one "change" per row (multiple rows or changes
%   allowed for each change set). Type "help idx_ct" for more complete
%   information about the format.
%   1st column: change set label, integer > 0
%   2nd column: change set probability
%   3rd column: table to be modified (1:bus, 2:gen, 3:branch) or
%               4: bus area changes, apply to all gens/buses/branches in
%               a given area; 5: gen table area changes apply to all
%               generators in a given area; or 6: branch area changes apply
%               to all branches connected to buses in a given area.
%   4th column: row of table to be modified (if 3rd col is 1-3),
%               or area number for overall changes (if 3rd column is 4-6).
%   5th column: column of table to be modified. It is best to use the
%               named column index defined in the corresponding idx_bus,
%               idx_gen, idx_branch or idx_cost M-files in MATPOWER.
%   6th column: type of change: 1: absolute (replace)
%                               2: relative (multiply by factor)
%                               3: additive (add to value)
%   7th column: new value or multiplicative or additive factor
%
%   Examples:
%
%   chgtab = [ ...
%       1   0.1   CT_TGEN       2  GEN_STATUS     CT_REP  0;
%       2   0.05  CT_TGEN       3  PMAX           CT_REP  100;
%       3   0.2   CT_TBRCH      2  BR_STATUS      CT_REP  0;
%       4   0.1   CT_TAREALOAD  2  CT_LOAD_ALL_P  CT_REL  1.1;
%    ];
%
%   Description of each change set:
%   1. Turn off generator 2, 10% probability.
%   2. Set generator 3's max output to 100 MW, 5% probability.
%   3. Take branch 2 out of service, 20% probability.
%   4. Scale all loads in area 2 (real & reactive, fixed and dispatchable)
%      by a factor of 1.1, 10% probability.
%
%   See IDX_CT.

%   To do:
%       - check for valid row number

%   MATPOWER
%   Copyright (c) 2000-2016, Power Systems Engineering Research Center (PSERC)
%   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
%   and Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%% define named indices into data matrices
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;
[CT_LABEL, CT_PROB, CT_TABLE, CT_TBUS, CT_TGEN, CT_TBRCH, CT_TAREABUS, ...
    CT_TAREAGEN, CT_TAREABRCH, CT_ROW, CT_COL, CT_CHGTYPE, CT_REP, ...
    CT_REL, CT_ADD, CT_NEWVAL, CT_TLOAD, CT_TAREALOAD, CT_LOAD_ALL_PQ, ...
    CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ, CT_LOAD_ALL_P, CT_LOAD_FIX_P, ...
    CT_LOAD_DIS_P, CT_TGENCOST, CT_TAREAGENCOST, CT_MODCOST_F, ...
    CT_MODCOST_X] = idx_ct;

%% find the change set to be applied
if numel(label) > 1
    error('apply_changes: LABEL must be a scalar');
end
kk = find(label == chgtab(:, CT_LABEL));
if isempty(kk)
    error('apply_changes: LABEL %d not found in CHGTAB', label);
end

%% create map of external bus numbers to bus indices
i2e = mpc.bus(:, BUS_I);
e2i = sparse(max(i2e), 1);
e2i(i2e) = (1:size(mpc.bus, 1))';

for c = 1:length(kk)
    %% extract data for individual change to be applied
    %% (c-th change in specified change set)
    change = num2cell(chgtab(kk(c), :));
    [tbl, row, col, typ, val] = ...
        deal(change{[CT_TABLE, CT_ROW, CT_COL, CT_CHGTYPE, CT_NEWVAL]});

    %% apply the change
    if tbl == CT_TBUS                               %% modify bus table
      if ~any(col == [PD QD GS BS VMAX VMIN])
        error('apply_changes: modification to column %d of bus table not supported', col);
      end
      if row == 0                                       %% modify all rows
        if typ == CT_REP                                    %% replace
          mpc.bus(:, col) = val * ones(size(mpc.bus, 1), 1);
        elseif typ == CT_REL                                %% scale
          mpc.bus(:, col) = val * mpc.bus(:, col);
        elseif typ == CT_ADD                                %% shift
          mpc.bus(:, col) = val + mpc.bus(:, col);
        else
          error('apply_changes: unsupported modification type %d for bus table', typ);
        end
      else                                              %% modify single row
        if typ == CT_REP                                    %% replace
          mpc.bus(row, col) = val;
        elseif typ == CT_REL                                %% scale
          mpc.bus(row, col) = val * mpc.bus(row, col);
        elseif typ == CT_ADD                                %% shift
          mpc.bus(row, col) = val + mpc.bus(row, col);
        else
          error('apply_changes: unsupported modification type %d for bus table', typ);
        end
      end
    elseif tbl == CT_TBRCH                          %% modify branch table
      if ~any(col == ...
            [BR_R BR_X BR_B RATE_A RATE_B RATE_C TAP SHIFT BR_STATUS ANGMIN ANGMAX])
        error('apply_changes: modification to column %d of branch table not supported', col);
      end
      if row == 0                                       %% modify all rows
        if typ == CT_REP                                    %% replace
          mpc.branch(:, col) = val * ones(size(mpc.branch, 1), 1);
        elseif typ == CT_REL                                %% scale
          mpc.branch(:, col) = val * mpc.branch(:, col);
        elseif typ == CT_ADD                                %% shift
          mpc.branch(:, col) = val + mpc.branch(:, col);
        else
          error('apply_changes: unsupported modification type %d for branch table', typ);
        end
      else                                              %% modify single row
        if typ == CT_REP                                    %% replace
          mpc.branch(row, col) = val;
        elseif typ == CT_REL                                %% scale
          mpc.branch(row, col) = val * mpc.branch(row, col);
        elseif typ == CT_ADD                                %% shift
          mpc.branch(row, col) = val + mpc.branch(row, col);
        else
          error('apply_changes: unsupported modification type %d for branch table', typ);
        end
      end
    elseif tbl == CT_TGEN                           %% modify gen table
      if ~any(col == ...
            [QMAX QMIN GEN_STATUS PMAX PMIN PC1 PC2 QC1MIN QC1MAX QC2MIN  QC2MAX ...
            RAMP_AGC RAMP_10 RAMP_30 RAMP_Q APF])
        error('apply_changes: modification to column %d of gen table not supported', col);
      end
      if row == 0                                       %% modify all rows
        if typ == CT_REP                                    %% replace
          mpc.gen(:, col) = val * ones(size(mpc.gen, 1), 1);
        elseif typ == CT_REL                                %% scale
          mpc.gen(:, col) = val * mpc.gen(:, col);
        elseif typ == CT_ADD                                %% shift
          mpc.gen(:, col) = val + mpc.gen(:, col);
        else
          error('apply_changes: unsupported modification type %d for gen table', typ);
        end
      else                                              %% modify single row
        if typ == CT_REP                                    %% replace
          mpc.gen(row, col) = val;
        elseif typ == CT_REL                                %% scale
          mpc.gen(row, col) = val * mpc.gen(row, col);
        elseif typ == CT_ADD                                %% shift
          mpc.gen(row, col) = val + mpc.gen(row, col);
        else
          error('apply_changes: unsupported modification type %d for gen table', typ);
        end
      end
    elseif tbl == CT_TGENCOST                       %% modify gencost table
      if col == CT_MODCOST_F || col == CT_MODCOST_X     %% use modcost to scale/shift cost
        if typ == CT_REL
          if col == CT_MODCOST_F
            modcost_type = 'SCALE_F';
          else  %% col == CT_MODCOST_X
            modcost_type = 'SCALE_X';
          end
        elseif typ == CT_ADD
          if col == CT_MODCOST_F
            modcost_type = 'SHIFT_F';
          else  %% col == CT_MODCOST_X
            modcost_type = 'SHIFT_X';
          end
        else
            error('apply_changes: unsupported modification type %d for gencost table CT_MODCOST_F/X modification', typ);
        end
        if row == 0
          mpc.gencost = modcost(mpc.gencost, val, modcost_type);
        else
          mpc.gencost(row, :) = modcost(mpc.gencost(row, :), val, modcost_type);
        end
      else                                              %% normal individual column mod
        if col < 1 || fix(col) ~= col   %% needs to be positive integer
          error('apply_changes: modification to column %d of gencost table not supported', col);
        end
        if row == 0                                       %% modify all rows
          if typ == CT_REP                                    %% replace
            mpc.gencost(:, col) = val * ones(size(mpc.gencost, 1), 1);
          elseif typ == CT_REL                                %% scale
            mpc.gencost(:, col) = val * mpc.gencost(:, col);
          elseif typ == CT_ADD                                %% shift
            mpc.gencost(:, col) = val + mpc.gencost(:, col);
          else
            error('apply_changes: unsupported modification type %d for gencost table', typ);
          end
        else                                              %% modify single row
          if typ == CT_REP                                    %% replace
            mpc.gencost(row, col) = val;
          elseif typ == CT_REL                                %% scale
            mpc.gencost(row, col) = val * mpc.gencost(row, col);
          elseif typ == CT_ADD                                %% shift
            mpc.gencost(row, col) = val + mpc.gencost(row, col);
          else
            error('apply_changes: unsupported modification type %d for gencost table', typ);
          end
        end
      end
    elseif tbl == CT_TAREABUS                   %% area-wide mod to bus table
      if ~any(col == [PD QD GS BS VMAX VMIN])
        error('apply_changes: area-wide modification to column %d of bus table not supported', col);
      end
      jj = find(mpc.bus(:, BUS_AREA) == row);
      if typ == CT_REP                                      %% replace
        mpc.bus(jj, col) = val * ones(size(jj));
      elseif typ == CT_REL                                  %% scale
        mpc.bus(jj, col) = val * mpc.bus(jj, col);
      elseif typ == CT_ADD                                  %% shift
        mpc.bus(jj, col) = val + mpc.bus(jj, col);
      else
        error('apply_changes: unsupported area-wide modification type %d for bus table', typ);
      end
    elseif tbl == CT_TAREABRCH                  %% area-wide mod to branch table
      if ~any(col == ...
            [BR_R BR_X BR_B RATE_A RATE_B RATE_C TAP SHIFT BR_STATUS ANGMIN ANGMAX])
        error('apply_changes: area-wide modification to column %d of branch table not supported', col);
      end
      jj = find((row == mpc.bus(e2i(mpc.branch(:, F_BUS)), BUS_AREA)) | ...
                (row == mpc.bus(e2i(mpc.branch(:, T_BUS)), BUS_AREA)) );
      if typ == CT_REP                                      %% replace
        mpc.branch(jj, col) = val * ones(size(jj));
      elseif typ == CT_REL                                  %% scale
        mpc.branch(jj, col) = val * mpc.branch(jj, col);
      elseif typ == CT_ADD                                  %% shift
        mpc.branch(jj, col) = val + mpc.branch(jj, col);
      else
        error('apply_changes: unsupported area-wide modification type %d for branch table', typ);
      end
    elseif tbl == CT_TAREAGEN                   %% area-wide mod to gen table
      if ~any(col == ...
            [QMAX QMIN GEN_STATUS PMAX PMIN PC1 PC2 QC1MIN QC1MAX QC2MIN  QC2MAX ...
            RAMP_AGC RAMP_10 RAMP_30 RAMP_Q APF])
        error('apply_changes: area-wide modification to column %d of gen table not supported', col);
      end
      jj = find(row == mpc.bus(e2i(mpc.gen(:, GEN_BUS)), BUS_AREA));
      if typ == CT_REP                                      %% replace
        mpc.gen(jj, col) = val * ones(size(jj));
      elseif typ == CT_REL                                  %% scale
        mpc.gen(jj, col) = val * mpc.gen(jj, col);
      elseif typ == CT_ADD                                  %% shift
        mpc.gen(jj, col) = val + mpc.gen(jj, col);
      else
        error('apply_changes: unsupported area-wide modification type %d for gen table', typ);
      end
    elseif tbl == CT_TAREAGENCOST               %% area-wide mod to gencost table
      jj = find(row == mpc.bus(e2i(mpc.gen(:, GEN_BUS)), BUS_AREA));
      if col == CT_MODCOST_F || col == CT_MODCOST_X
        if typ == CT_REL
          if col == CT_MODCOST_F
            modcost_type = 'SCALE_F';
          else  %% col == CT_MODCOST_X
            modcost_type = 'SCALE_X';
          end
        elseif typ == CT_ADD
          if col == CT_MODCOST_F
            modcost_type = 'SHIFT_F';
          else  %% col == CT_MODCOST_X
            modcost_type = 'SHIFT_X';
          end
        else
            error('apply_changes: unsupported area-wide modification type %d for gencost table CT_MODCOST_F/X modification', typ);
        end
        mpc.gencost(jj, :) = modcost(mpc.gencost(jj, :), val, modcost_type);
      else
        if col < 1 || fix(col) ~= col   %% needs to be positive integer
          error('apply_changes: area-wide modification to column %d of gencost table not supported', col);
        end
        if typ == CT_REP                                      %% replace
          mpc.gencost(jj, col) = val * ones(size(jj));
        elseif typ == CT_REL                                  %% scale
          mpc.gencost(jj, col) = val * mpc.gencost(jj, col);
        elseif typ == CT_ADD                                  %% shift
          mpc.gencost(jj, col) = val + mpc.gencost(jj, col);
        else
          error('apply_changes: unsupported area-wide modification type %d for gencost table', typ);
        end
      end
    elseif tbl == CT_TLOAD || tbl == CT_TAREALOAD           %% modify loads
      if ~any(abs(col) == (CT_LOAD_ALL_PQ:CT_LOAD_DIS_P))
        error('apply_changes: column=%d for load modifications is not supported', col);
      end
      switch abs(col)
        case {CT_LOAD_ALL_PQ, CT_LOAD_FIX_PQ, CT_LOAD_DIS_PQ}
            opt.pq = 'PQ';
        case {CT_LOAD_ALL_P, CT_LOAD_FIX_P, CT_LOAD_DIS_P}
            opt.pq = 'P';
      end
      switch abs(col)
        case {CT_LOAD_ALL_PQ, CT_LOAD_ALL_P}
            opt.which = 'BOTH';
        case {CT_LOAD_FIX_PQ, CT_LOAD_FIX_P}
            opt.which = 'FIXED';
        case {CT_LOAD_DIS_PQ, CT_LOAD_DIS_P}
            opt.which = 'DISPATCHABLE';
      end

      %% define load_zone
      if tbl == CT_TLOAD
        nb = size(mpc.bus, 1);
        if row == 0                                     %% modify load at all buses
          load_zone = ones(nb, 1);
        else                                            %% modify load at single bus
          load_zone = zeros(nb, 1);
          load_zone(row) = 1;
        end
      elseif tbl == CT_TAREALOAD
        load_zone = double(mpc.bus(:, BUS_AREA) == row);  %% modify load in single area
      end

      switch typ
        case CT_REP                                     %% replace
            opt.scale = 'QUANTITY';
            dmd = val;
        case CT_REL                                     %% scale
            opt.scale = 'FACTOR';
            dmd = val;
        case CT_ADD                                     %% shift
            opt.scale = 'QUANTITY';
            old_val = total_load(mpc, load_zone);
            dmd = old_val + val;
        otherwise
          error('apply_changes: unsupported modification type %d for loads', typ);
      end
      %% scale the loads ...
      if col < 0    %% ... including dispatchable load costs
          [mpc.bus, mpc.gen, mpc.gencost] = scale_load(dmd, mpc.bus, mpc.gen, load_zone, opt, mpc.gencost);
      else          %% ... not including dispatchable load costs
          [mpc.bus, mpc.gen] = scale_load(dmd, mpc.bus, mpc.gen, load_zone, opt);
      end
    else
      error('apply_changes: CHGTAB attempts to modify unsupported table type')
    end
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower