text:
function [MetsAndRxns, pVal_up, pVal_down, Dir_1, Dir_2] = compareConditions(MetConnectivity1, ConnectedMet1, MetConnectivity2, ConnectedMet2)
% This function compare reaction score of both sampling condition
% and identify metabolites that change significantly
%
% USAGE:
%
%    [MetsAndRxns, pVal_up, pVal_down, Dir_1, Dir_2] = compareConditions(MetConnectivity1, ConnectedMet1, MetConnectivity2, ConnectedMet2)
%
% INPUTS:
%    MetConnectivity1:    A structure that present, for each
%                         metabolite present in the model sampled under first
%                         condition the following sets of fields:
%
%                           * ConnRxns - the reactions that are connected
%                             to the metabolite
%                           * Sij - The stoichiometric coefficient for the
%                             metabolite in each reaction in `ConnRxns`
%                           * RxnScore - Score for each reaction in `ConnRxns`
%                           * Direction - The direction of reaction flux
%                             for each sample point
%                           * MetNotUsed - Whether or not the metabolite is
%                             used in the condition
%    ConnectedMet1:       List of metabolites described in `MetConnectivity1`
%    MetConnectivity2:    Same as `MetConnectivity1` but for `model`
%                         sampled under the second condition
%    ConnectedMet2:       List of metabolites described in `MetConnectivity2`
%
% OUTPUTS:
%    MetsAndRxns:         Cell arrays containing in the first column
%                         the list of metabolites that significantly change
%                         under both sampling conditions and in the second
%                         column the reactions that are connected to
%                         these metabolites
%    pVal_up:             p-value associated to upregulated
%                         `MetsAndRxns`
%    pVal_down:           p-value associated to downregulated
%                         `MetsAndRxns`
%    Dir_1:               reaction directionality for model sampled under first
%                         condition (1 producing metabolite, -1
%                         consuming the metabolite)
%    Dir_2:               Same as `Dir_1` but for model
%                         sampled under the second condition
%
% .. Authors:
%       - Nathan E. Lewis, May 2010-May 2011
%       - Anne Richelle, May 2017

Met1=MetConnectivity1.(ConnectedMet1{1});
Met2=MetConnectivity2.(ConnectedMet2{1});
% Get the minimum number of points used to score reactions in the two
% conditions
MinNumPts = min([length(Met1.RxnScore(1,:)) length(Met2.RxnScore(1,:))]);

% only look at nodes shared between the sampled models
ConnectedMet = intersect(ConnectedMet1,ConnectedMet2);
MetsAndRxns={};
pVal_up=[];
pVal_down=[];
Dir_1 = [];
Dir_2 = [];

for i = 1:length(ConnectedMet)

    Met1=MetConnectivity1.(ConnectedMet{i});
    Met2=MetConnectivity2.(ConnectedMet{i});

    FracNotUsed1 = sum(Met1.MetNotUsed(1:MinNumPts))/(MinNumPts); % fraction of the sample points not using this node in model1
    FracNotUsed2 = sum(Met2.MetNotUsed(1:MinNumPts))/(MinNumPts); % fraction of the sample points not using this node in model2

    % ignore nodes that are used in less than 10% of the sample points in one
    % condition
    if and(FracNotUsed1<.9,FracNotUsed2<.9)

        % subtract scores of 2nd condition from the 1st condition
        tmp = Met1.RxnScore(:,1:MinNumPts) - Met2.RxnScore(:,1:MinNumPts);

        % p of increasing flux is obtained by subtracting 2nd from 1st and
        % finding how many are > or < 0 for the pvalue for up and down,
        % respectively. (many >= 0 means 2nd condition was lower, so high
        % p-value for up and low for down, suggesting that the flux dropped
        % in the 2nd condition
        p_up = sum(tmp>=0,2)/length(tmp(1,:));
        p_down = sum(tmp<=0,2)/length(tmp(1,:));

        % This computes the median direction, just to give an idea which
        % direction the reaction tends to go in the conditions
        tmp_Dir1 = median(Met1.Direction(:,1:MinNumPts),2);
        tmp_Dir2 = median(Met2.Direction(:,1:MinNumPts),2);

        % Go through each reaction and build arrays with the reactions, the
        % branch-point metabolites, and p-values, and reaction directionality
        % (i.e., in or out of the central metabolite)
        for j = 1:length(tmp(:,1))
            if ~all(tmp(j,:)==0) % skip reactions that carry no flux or don't change at all
                MetsAndRxns{end+1,1} =ConnectedMet{i};
                MetsAndRxns{end,2}=Met1.ConnRxns{j};
                pVal_up(end+1,1)=p_up(j);
                pVal_down(end+1,1)=p_down(j);
                Dir_1(end+1,1)= tmp_Dir1(j);
                Dir_2(end+1,1)= tmp_Dir2(j);
            end
        end

    end
end

end

author: opencobra

repo: cobratoolbox


text:
function[dir]=whichdir(name)
%   WHICHDIR  Returns directory name containing file in search path.
%
%   WHICHDIR NAME returns the directory containing a file called NAME,
%   i.e. the full path of the file excluding NAME and the trailing '/'.
%
%   If more than one directory contains a file called NAME, then
%   WHICHDIR returns a cell array of directories names.
%
%   If NAME does not include an extension, then it is interpreted as
%   the name of an m-file, i.e. having extension '.m' .
%   _________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2006--2015 J.M. Lilly --- type 'help jlab_license' for details

if isempty(strfind(name,'.'))
    name=[name '.m'];
end

dir=which(name,'-all');

for i=1:length(dir)
   dir{i}=dir{i}(1:end-length(name)-1);
end

if i==1
    dir=dir{1};
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
% TABLE Table.
%    Tables are used to collect heterogeneous data and metadata into a single
%    container.  Tables are suitable for storing column-oriented or tabular data
%    that are often stored as columns in a text file or in a spreadsheet.  Tables
%    can accommodate variables of different types, sizes, units, etc.  They are
%    often used to store experimental data, with rows representing different
%    observations and columns representing different measured variables.
%
%    Use the TABLE constructor to create a table from variables in the MATLAB
%    workspace.  Use the readtable function to create a table by reading data
%    from a text or spreadsheet file.
%
%    The TABLE constructor can also be used to create tables without
%    providing workspace variables, by providing the size and variable
%    types.
%
%    Tables can be subscripted using parentheses much like ordinary numeric
%    arrays, but in addition to numeric and logical indices, you can use a
%    table's variable and row names as indices.  You can access individual
%    variables in a table much like fields in a structure, using dot
%    subscripting.  You can access the contents of one or more variables using
%    brace subscripting.
%
%    Tables can contain different kinds of variables, including numeric, logical,
%    character, categorical, and cell.  However, a table is a different class
%    than the variables that it contains.  For example, even a table that
%    contains only variables that are double arrays cannot be operated on as if
%    it were itself a double array.  However, using dot subscripting, you can
%    operate on a variable in a table as if it were a workspace variable.  Using
%    brace subscripting, you can operate on one or more variables in a table as
%    if they were in a homogeneous array.
%
%    A table T has properties that store metadata such as its variable and row
%    names.  Access or assign to a property using P = T.Properties.PropName or
%    T.Properties.PropName = P, where PropName is one of the following:
%
%    TABLE metadata properties:
%        Description          - A character vector describing the table
%        DimensionNames       - A two-element cell array of character vectors containing names of
%                               the dimensions of the table
%        VariableNames        - A cell array containing names of the variables in the table
%        VariableDescriptions - A cell array of character vectors containing descriptions of the
%                               variables in the table
%        VariableUnits        - A cell array of character vectors containing units for the variables
%                               in table
%        VariableContinuity   - An array containing a matlab.tabular.Continuity value for each table
%                               variable, specifying whether a variable represents continuous or discrete
%                               data values. You can assign 'unset', 'continuous', 'step', or 'event' to
%                               elements of VariableContinuity.
%        RowNames             - A cell array of nonempty, distinct character vectors containing names
%                               of the rows in the table
%        UserData             - A variable containing any additional information associated
%                               with the table.  You can assign any value to this property.
%
%    TABLE methods and functions:
%      Construction and conversion:
%        table              - Create a table from workspace variables.
%        array2table        - Convert homogeneous array to table.
%        cell2table         - Convert cell array to table.
%        struct2table       - Convert structure array to table.
%        table2array        - Convert table to a homogeneous array.
%        table2cell         - Convert table to cell array.
%        table2struct       - Convert table to structure array.
%      Import and export:
%        readtable          - Create a table by reading from a file.
%        writetable         - Write a table to a file.
%        write              - Write a table to a file.
%      Size and shape:
%        istable            - True for tables.
%        size               - Size of a table.
%        width              - Number of variables in a table.
%        height             - Number of rows in a table.
%        ndims              - Number of dimensions of a table.
%        numel              - Number of elements in a table.
%        horzcat            - Horizontal concatenation for tables.
%        vertcat            - Vertical concatenation for tables.
%      Set membership:
%        intersect          - Find rows common to two tables.
%        ismember           - Find rows in one table that occur in another table.
%        setdiff            - Find rows that occur in one table but not in another.
%        setxor             - Find rows that occur in one or the other of two tables, but not both.
%        unique             - Find unique rows in a table.
%        union              - Find rows that occur in either of two tables.
%      Data manipulation and reorganization:
%        summary            - Print summary of a table.
%        addvars            - Insert new variables at a specified location in a table.
%        movevars           - Move table variables to a specified location.
%        removevars         - Delete the specified table variables.
%        splitvars          - Splits multi-column variables into separate variables.
%        mergevars          - Merges multiple variables into one multi-column variable or a nested table.
%        sortrows           - Sort rows of a table.
%        stack              - Stack up data from multiple variables into a single variable.
%        unstack            - Unstack data from a single variable into multiple variables.
%        join               - Merge two tables by matching up rows using key variables.
%        innerjoin          - Inner join between two tables.
%        outerjoin          - Outer join between two tables.
%        rows2vars          - Reorient rows to be variables of output table.
%        inner2outer        - Invert a nested table-in-table hierarchy.
%        ismissing          - Find elements in a table that contain missing values.
%        standardizeMissing - Insert missing data indicators into a table.
%      Computations on tables:
%        varfun             - Apply a function to variables in a table.
%        rowfun             - Apply a function to rows of a table.
%
%    Examples:
%
%       % Create a table from individual workspace variables.
%       load patients
%       patients = table(LastName,Gender,Age,Height,Weight,Smoker,Systolic,Diastolic)
%
%       % Select the rows for patients who smoke, and a subset of the variables.
%       smokers = patients(patients.Smoker == true, {'LastName' 'Gender' 'Systolic' 'Diastolic'})
%
%       % Convert the two blood pressure variables into a single variable.
%       patients.BloodPressure = [patients.Systolic patients.Diastolic];
%       patients(:,{'Systolic' 'Diastolic'}) = []
%
%       % Pick out two specific patients by the LastName variable.
%       patients(ismember(patients.LastName,{'Smith' 'Jones'}), :)
%
%       % Convert the LastName variable into row names.
%       patients.Properties.RowNames = patients.LastName;
%       patients.LastName = []
%
%       % Use the row names to pick out two specific patients.
%       patients({'Smith' 'Jones'},:)
%
%       % Add metadata to the table.
%       patients.Properties.Description = 'Simulated patient data';
%       patients.Properties.VariableUnits =  {''  'Yrs'  'In'  'Lbs'  ''  'mm Hg'};
%       patients.Properties.VariableDescriptions{6} = 'Systolic/Diastolic';
%       summary(patients)
%
%       % Create a new variable in the table from existing variables.
%       patients.BMI = (patients.Weight * 0.453592) ./ (patients.Height * 0.0254).^2
%       patients.Properties.VariableUnits{'BMI'} =  'kg/m^2';
%       patients.Properties.VariableDescriptions{'BMI'} = 'Body Mass Index';
%
%       % Sort the table based on the new variable.
%       sortrows(patients,'BMI')
%
%       % Make a scatter plot of two of the table's variables.
%       plot(patients.Height,patients.Weight,'o')
%
%       % Create tables from text and spreadsheet files
%       patients2 = readtable('patients.dat','ReadRowNames',true)
%       patients3 = readtable('patients.xls','ReadRowNames',true)
%
%       % Create a table from a numeric matrix
%       load tetmesh.mat
%       t = array2table(X,'VariableNames',{'x' 'y' 'z'});
%       plot3(t.x,t.y,t.z,'.')
%
%    See also TABLE, CATEGORICAL
%
%    Reference page in Doc Center
%       doc table
%
%    Other functions named table
%
%       codistributed/table    tall/table
%

author: jmaih


text:
function [A,B,Ypredict] = opls(X,Y,nhidden,algorithm)
% BLS  Bottleneck least squares aka sparse orthogonalized least squares
%
% X: ninput x nsamples input data matrix
% Y: noutput x nsamples output data matrix
% nhidden: dimension of hidden
%
% A: noutput x nhidden weight matrix
% B: ninput x nhidden weight matrix

if nargin < 3,
  nhidden = min(size(Y,1),2);
end

if nargin < 4,
  algorithm = 'direct';
end

[ninput,nsamples] = size(X);

Sxx = symmetric(X*X')/nsamples;

switch algorithm

  case 'direct'
    Syx = Y*X'/nsamples;
    Sxy = Syx';
    Sxxinv = symmetric(pinv(Sxx));
    S = symmetric(Syx * Sxxinv * Sxy);
    opts.disp = 0;
    opts.issym = 'true';
    [A,D] = eigs(S,nhidden,'LM',opts);   % first nhidden eigenvectors of S
    B = Sxxinv*Sxy*A;

  case 'indirect'          % same result and subspace as 'direct', but then different linear combination

    Syx = Y*X'/nsamples;
    Sxy = Syx';
    Sxxinv = symmetric(pinv(Sxx));

    % initialize randomly

    B = randn(ninput,nhidden);

    iter = 0;
    maxiter = 1000;
    tol = nhidden*ninput*(1e-10);
    Bold = B;
    while iter < maxiter,
      Z = B'*X;
      %		Szz = symmetric(Z*Z')/nsamples;
      Szz = symmetric(B'*Sxx*B);
      %		Syz = symmetric(Y*Z')/nsamples;
      Syz = Syx*B;
      A = Syz*pinv(Szz);

      B = Sxxinv*Sxy*A*pinv(A'*A);
      if sumsqr(B - Bold) < tol,
        fprintf('done after %d iterations!!\n',iter+1);
        iter = maxiter;
      else
        iter = iter + 1;
        Bold = B;
      end
    end

  case 'sequential'      % same result as 'direct', up to minus sign per hidden unit

    R = Y;
    noutput = size(Y,1);
    A = zeros(noutput,nhidden);
    B = zeros(ninput,nhidden);
    for i=1:nhidden,
      [A(:,i),B(:,i)] = bls(X,R,1,'direct');
      if i < nhidden,
        R = R - A(:,i)*B(:,i)'*X;
      end
    end

end

% transform to make latent variable orthonormal

Szz = symmetric(B'*Sxx*B);
Szzinv = symmetric(pinv(Szz));
sqrtSzzinv = chol(Szzinv(end:-1:1,end:-1:1));    % such that the rescaled B(i,:) is a linear combi of the unscaled B(1:i,:); this only makes sense for the direct method
sqrtSzzinv = sqrtSzzinv(end:-1:1,end:-1:1);
B = B*sqrtSzzinv;
A = A*pinv(sqrtSzzinv);

if nargout > 2,
  Ypredict = A*B'*X;
end

function A = symmetric(A)

A = (A+A')/2;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [imr,xa,ya]=rectifyhomog(ims,roc,k0,mode)
%RECTIFYHOMOG Apply rectifying homographies to an image set
%
% Usage:    figure(101);     % Initial figure for rectiied image display
%           rectifyhomog(ims,roc,k0,'ga');   % plot in individual figures
%
% Inputs:
%        ims{nc}       cell array of input images (colour or monochrome)
%        roc(3,3,nc)   rotation matrices from world coordinates to camera coordinates
%        k0            camera matrix or focal length in pixels optionally divided by the image width [0.8]
%        mode          mode string
%                         g  show images on separate figures
%                         G  tile images onto a single figure [default if no output arguments]
%                         k  clip to original image dimensions
%                         l  do not link axes
%                         v  k0 contains diffent values for each camera
%                         a  orient to average camera orientation
%                        [z  make optical axis perpendicular to camera plane]
%                        [x  align x axis with camera displacements]
% Outputs:
%        imr{nc}(my,mx,nc)  output images (uint8)
%        xa{nc}(mx)        x axis for each image
%        ya{nc}(my)        y axis for each image

%      Copyright (C) Mike Brookes 2012
%      Version: $Id: rectifyhomog.m 1642 2012-03-16 16:22:26Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~iscell(ims)
    ims={ims};
end
nc=numel(ims); % number of images
if ~((nc==1 && isequal(size(roc),[3 3])) || isequal(size(roc),[3 3 nc]))
    error('roc should have dimensions [3,3,%d]',nc);
end
if nargin<4 || ~numel(mode)
    mode='';
end
if nargin<3 || ~numel(k0)
    k0=0.8;
end
if nargin<2 || ~numel(roc)
    roc=repmat(eye(3),[1 1 nc]);
end

% sort out options

vv=any(mode=='v');  % different k0 for each camera
nk=numel(k0)/(1+vv*(nc-1));  % number of k parameters per camera
switch nk
    case {1,3}
        k0=reshape(k0,[1 nk 1+vv*(nc-1)]);
    case 9
        k0=reshape(k0,[3 3 1+vv*(nc-1)]);
    otherwise
        error('k0 must have 1, 3 or 9 values per camera');
end
ncr=1+(nargout>0)*(nc-1);   % numer of output images
if any(mode=='g')
    gmode=1;
elseif any(mode=='G') || ~nargout
    gmode=2;
else
    gmode=0;
end
if gmode>0
    fig0=gcf; % initialize the figure
end
if any(mode=='k')
    modeh='kxt';
else
    modeh='t';
end

% determine a global camera rotation

if any(mode=='a')
    qrc=zeros(4,nc);  % calculate the mean camera orientation
    for i=1:nc
        qrc(:,i)=rotro2qr(roc(:,:,i));
    end
    rocmean=rotqr2ro(rotqrmean(qrc));
else
    rocmean=eye(3);
end

% now do image transformations

imr=cell(ncr,1);
xa=imr;
ya=imr;
axh=zeros(nc,1);
splx=ceil(sqrt(nc));
sply=ceil(nc/splx);
for i=1:nc
    k0i=k0(:,:,1+vv*(i-1));  % camera parameters
    if numel(k0i)<9
        imsz=size(ims{i});
        fe=k0(1);
        if fe<0.1*imsz(2) % focal length is a fraction of the width
            fe=k0*imsz(2);
        end
        if numel(k0i)<3
            xy0=(imsz(2:-1:1)+1)/2;
        else
            xy0=k0i(2:3);
        end
        k0=eye(3);
        k0(1:4:5)=fe;
        k0(7:8)=xy0;
    end
    j=min(i,ncr);
    rocall=rocmean*roc(:,:,i)';
    titl=sprintf('%d: pan-tilt-roll = %.1f°, %.1f°, %.1f°',i,-rotro2eu('yxz',rocall)*180/pi);
    [imr{j},xa{j},ya{j}]=imagehomog(uint8(ims{i}),k0*rocall/k0,modeh);  % apply inverse of rotation matrix
    if gmode>0
        if gmode>1
            subplot(sply,splx,i);
        else
            figure(fig0+i-1);
        end
        imagesc(xa{j},ya{j},imr{j});
        axis image
        title(titl);
        axh(i)=gca;
    end
end
if gmode>0 && nc>1
    if ~any(mode=='l')
        linkaxes(axh);
    end
    figure(fig0);
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [f] = spm_fx_fnirs(x,u,P,M)
% State equation for a dynamic model of fNIRS responses
% FORMAT [f] = spm_fx_nirs(x,u,P,M)
%
% x      - state vector
%--------------------------------------------------------------------------
%   x(:,1) - excitatory neuronal activity            ue
%   x(:,2) - vasodilatory signal                          s
%   x(:,3) - rCBF                                  ln(f)
%   x(:,4) - venous volume                         ln(v)
%   x(:,5) - deoxyHb                               ln(q)
%   x(:,6) - totalHb                               ln(p)
%   [x(:,7)] - inhibitory neuronal activity            ui
%--------------------------------------------------------------------------
% u     experimental inputs
% P     prior of latent variables
% M    model structure
%
% f      - dx/dt
%___________________________________________________________________________
%
% References for hemodynamic & neuronal state equations:
% 1. Friston KJ, Mechelli A, Turner R, Price CJ. Nonlinear responses in
%    fMRI: the Balloon model, Volterra kernels, and other hemodynamics.
%    Neuroimage 12:466-477, 2000.
% 2. Stephan KE, Kasper L, Harrison LM, Daunizeau J, den Ouden HE,
%    Breakspear M, Friston KJ. Nonlinear dynamic causal models for fMRI.
%    Neuroimage 42:649-662, 2008.
% 3. Marreiros AC, Kiebel SJ, Friston KJ. Dynamic causal modelling for
%    fMRI: a two-state model.
%    Neuroimage. 39(1):269-78, 2008.
% 4. Buxton RB, Uludag, K, Dubowitz, DJ, Liu, TT. Modeling the hemodynamic
%    response to brain activation. Neuroimage. 2004, 23: 220-233.
% 5. X Cui and S Bray and A Reiss. Functional near infrared spectroscopy (NIRS)
%    signal improvement based on negative correlation between oxygenated and
%    deoxygenated hemoglobin dynamics.
%    Neuroimage 49:3039-3046, 2010.
%
% This script is based on spm_fx_fmri.m written by
% Karl Friston & Klaas Enno Stephan.
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny & Sungho Tak
% $Id: spm_fx_fnirs.m 6942 2016-11-21 13:17:44Z guillaume $

% Neuronal motion
%==========================================================================
P.A   = full(P.A);                       %    linear parameters
P.B   = full(P.B);                       % bi-linear parameters
P.C   = P.C/16;                          % exogenous parameters

% implement differential state equation y = dx/dt (neuronal)
%--------------------------------------------------------------------------
f    = x;

% input dependent modulation
%--------------------------------------------------------------------------
for i = 1:size(P.B,3)
    P.A(:,:,1) = P.A(:,:,1) + u(i)*P.B(:,:,i);
end

if size(x, 2) == 6 % one neuronal state per region
    %======================================================================

    % one neuronal state per region: diag(A) is a log self-inhibition
    %----------------------------------------------------------------------
    SI     = diag(P.A);
    P.A    = P.A - diag(exp(SI)/2 + SI);

    % flow
    %----------------------------------------------------------------------
    f(:,1) = P.A*x(:,1) + P.C*u(:);

elseif size(x,2) == 7 % two neuronal states per region
    %======================================================================

    % extrinsic (two neuronal states): enforce positivity
    %----------------------------------------------------------------------
    n     = size(P.A,1);            % number of regions
    EE    = exp(P.A(:,:,1))/8;
    IE    = diag(diag(EE));         % intrinsic inhibitory to excitatory
    EE    = EE - IE;                % extrinsic excitatory to excitatory
    EI    = eye(n,n);               % intrinsic excitatory to inhibitory
    SE    = eye(n,n)/2;             % intrinsic self-inhibition (excitatory)
    SI    = eye(n,n);               % intrinsic self-inhibition (inhibitory)

    % excitatory proportion
    %----------------------------------------------------------------------
    if size(P.A,3) > 1
        phi = spm_phi(P.A(:,:,2)*2);
        EI  = EI + EE.*(1 - phi);
        EE  = EE.*phi - SE;
    else
        EE  = EE - SE;
    end

    % motion - excitatory and inhibitory: f = dx/dt
    %----------------------------------------------------------------------
    f(:,1) = EE*x(:,1) - IE*x(:, 7) + P.C*u(:);
    f(:,6 + 1) = EI*x(:,1) - SI*x(:, 7);
end

% Hemodynamic motion
%==========================================================================

% hemodynamic parameters
%--------------------------------------------------------------------------
%   H(1) - signal decay                                   d(ds/dt)/ds)
%   H(2) - autoregulation                                 d(ds/dt)/df)
%   H(3) - transit time                                      (t0)
%   H(4) - exponent for Fout(v)                           (alpha)
%   H(5) - resting oxygen extraction                      (E0)
%   H(6) - viscoelastic time constant              (tv)
%--------------------------------------------------------------------------
H = [0.64 0.32 2.00 0.32 0.32 2.00];

% exponentiation of hemodynamic state variables
%--------------------------------------------------------------------------
x(:,3:6) = exp(x(:,3:6));

% signal decay
sd = H(1)*exp(P.decay);

% autoregulatory feedback
af = H(2).*exp(P.afback);

% transit time
tt = H(3).*exp(P.transit);
% viscoelastic time constant
tv = H(6).*exp(P.tv);

% outflow f(v) (steady state)
fv_s = x(:,4).^(1/H(4));

% implement differential state equation f = dx/dt (hemodynamic)
%--------------------------------------------------------------------------
f(:,2)   = x(:,1) - sd.*x(:,2) - af.*(x(:,3) - 1);
f(:,3)   = x(:,2)./x(:,3);
f(:,4)   = (x(:,3) - fv_s)./((tt+tv).*x(:,4));

% outflow (dynamic state)
fv_d = fv_s + tv.*x(:,4).*f(:,4);

% oxygen extraction fraction
ff = (1 - (1 - H(5)).^(1./x(:,3)))/H(5);

f(:,5)   = (x(:,3).*ff - fv_d.*x(:,5)./x(:,4))./(tt.*x(:,5));
f(:,6)   = (x(:,3) - fv_d) ./ (tt.*x(:,4));

f        = f(:);

author: spm


text:
function orientationStruct = readLeksellOrientationFile(filename)
%"readLeksellOrientationFile"
%   Reads a Leksell orientation file using decodeLeksellData and places the
%   fields into a datastructure with meaningful names.  The transformation
%   matrices found in these files have no known use in CERR.
%
%JRA 6/13/05
%
%Usage:
%   orientationStruct = readLeksellOrientationFile(filename)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

fid = fopen(filename, 'r', 'b');

data = decodeLeksellData(fid);

fclose(fid);

%If data does not have at least two elements it is effectively empty.
if length(data) < 2
    orientationStruct.transM = [];
    orientationStruct(1) = [];
    return;
end

%Get what appears to be transformation matrices.
transMats = data{1};

%Iterate over them, reshape to 3x3 and save.
%However, the purposes of these orientation matrices have not been found.
for i=1:length(transMats)
    orientationStruct(i).transM = reshape(transMats{i}, [3 3]);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ spins, detectedPts, locations, pts, density ] = create_mesh_spin_features(vertex, faces, bins,scale_rat, cscale,opnormal)
% detects feature locations, and computes spin images for these points
%Input:
% vertex, faces - input mesh
% bins - spin image resolution (default: 5)
% scale_rat - a width constant for the SI support (default: 3)
% cscale - 1: use constant scale. 0 - scale invariant. (default: 1)
%Output:
% spins - the spin images (a row for each detected point)
% detectedPts - the indices of points detected by the DOG detector for each
%               level
% locations(:,1) - list of the detected points
% locations(:,2) - list of the scale levels
% locations(:,3) - the local density of each detected point
% pts - same as locations(:,1)
% density - density of each vertex of the mesh
if size (vertex,2) ~= 3 || size (faces,2) ~= 3
     error('incompetible vertex or faces array size');
end

spins = [];
if ~exist('opnormal','var')
    opnormal = 0;
end
if ~exist('cscale','var')
    cscale = 0;
end
if ~exist('scale_rat','var')
    scale_rat = 3;
end
if ~exist('bins','var')
    bins = 5;
end
% mesh.vertices = vertex;
% mesh.vertexNormals = compute_normal(vertex,faces)';
normals = compute_normal(vertex,faces)';
if (opnormal == 1)
    normals = -normals;
end
[detectedPts, density] = dog(vertex,faces, 20);
%detectedPts = detectedPts(2:end);
len = length(detectedPts);
pts=[];
ImageSize = bins;
for ll = 1:len
    if (~isempty(detectedPts{ll}))
        pts = [pts detectedPts{ll}];
        if cscale
            spinstmp{ll} = SpinImages(vertex, normals, detectedPts{ll}',cscale* scale_rat * ones(size(density)) / ImageSize, ImageSize);
        else
            spinstmp{ll} = SpinImages(vertex, normals, detectedPts{ll}', scale_rat*density * sqrt(ll) / ImageSize, ImageSize);
        end
        %spinstmp{ll} = SpinImages(vertex, normals, detectedPts{ll}', scale_rat*scale  / ImageSize, ImageSize);
    end
end

sp2 = length(spinstmp);
kind = 1;
ImageSize2 = ImageSize*ImageSize;
for ind2 = 1:sp2
    sp3 = length(spinstmp{ind2});
    for ind3 = 1:sp3
        spins(kind,:) = reshape(spinstmp{ind2}{ind3},1,ImageSize2);
        locations(kind,:) = [detectedPts{ind2}(ind3) ind2 density(detectedPts{ind2}(ind3))];
        kind = kind + 1;
    end
end
end

author: Sable

repo: mcbench-benchmarks


text:
function b = iscellnumeric(C)
% Return 1 if all elements of cell array are numeric
%
% Tim Mullen, 2011, SCCN/INC, UCSD

b = all(cellfun(@(x) isnumeric(x),C));

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function orbEls=state2OrbEls(stateVec,elType,GM,epsVal)
%%STATE2ORBELS Convert a state vector consisting of at least position and
%              velocity into a set of orbital elements. Such elements are
%              for a simple two-body problem where the satellite (object)
%              in motion has negligible mass. The elements can be found
%              with respect to a given epoch time.
%
%INPUTS: stateVec A 6XnumVec matrix of numVec state vectors consisting of
%                 position and velocity in a Cartesian (quasi)-inertial
%                 coordinate system where the gravitating body is at the
%                 origin. Extra state elements are ignored. The units are
%                 assumed to be meters and meters per second.
%          elType A value indicating the type of orbital elements. Possible
%                 values are:
%                 0 (The default if omitted) The elements are Gooding's
%                   universal orbital elements.
%                 2 The elements are direct equinoctial orbital elements.
%                 3 The elements are retrograde equnoctial orbital
%                   elements.
%              GM Optionally, the universal gravitation constant times the
%                 mass of the body about which the object with the given
%                 state vector is orbiting. If this parameter is omitted,
%                 the value in Constants.WGS84GMWithAtmosphere is used. The
%                 units are m^3/sec^2.
%          epsVal If universal orbital elements are chosen, this is a
%                 precision bound used for equality comparisons when
%                 determining whether a trajectory is rectilinear,
%                 parabolic or circular. If omitted, a default value of
%                 eps is used.
%
%OUTPUTS: orbEls An 6XnumVec set of vectors of orbital elements, the format
%                of which depends on the elType parameter. The format of
%                the elements is discussed below.
%
%Gooding's universal orbital elements are consist of (in order)
%alpha=GM/a where a is the semi-major axis in meters
%q=a*(1-e)  where e is the eccentricity (unitless). This is the perifocal
%           distance.
%i          inclination in radians
%Omega      longitude (right ascension) of the ascending node in radians
%omega      argument of periapsis/perigee in radians
%tau        the time at pericenter (seconds)
%
%The equinoctial orbital elements (for both direct and retrograde elements)
%are
%a      semi-major axis in meters
%h      first eccentricity vector component (unitless)
%k      second eccentricity vector component (unitless)
%p      first ascending node vector component (unitless)
%q      second ascending node vector component (unitless)
%lambda mean longitude in radians.
%
%The algorithm using Gooding's universal orbital elements is very robust to
%numerical errors. It is implemented using the algorithm of [1], where the
%elements are also described.
%
%Equinoctial orbital elements are discussed in Section 2 of [2]. and in
%[3]. The diffference between direct and retrograde elements is essentially
%a matter of handedness.
%
%Orbital elements are discussed in general in Chapter 2.2 of [4].
%
%The inverse of this function is orbEls2State, which can also predict the
%state forward in time. Adding a time interval in seconds to the tau
%argument also predicts the trajectory forward in terms of universal
%orbital elements.
%
%REFERENCES:
%[1] R. H. Gooding, "On universal elements, and conversion procedures to
%    and from position and velocity," Celestial mechanics, vol. 44, no. 3,
%    pp. 283-298, 1988.
%[2] D. A. Danielson, C. P. Sagovac, B. Neta, and L. W. Early,
%    "Semianalytic satellite theory," Mathematics Department, Naval
%    Postgraduate School, Monterey, CA, Tech. Rep., 1995. [Online].
%    Available: http://oai.dtic.mil/oai/oai?verb=getRecord&metadataPrefix= html&identifier=ADA531136
%[3] R. A. Broucke and P. J. Cefola, "On the equinoctial orbit elements,"
%    Celestial Mechanics, vol. 5, no. 3, pp. 303-310, 1972.
%[4] O. Montenbruck and E. Gill, Satellite Orbits: Models, Methods
%    Applications. Berlin: Springer, 2000.
%
%January 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4)
    epsVal=eps;
end

if(nargin<3)
   GM=Constants.WGS84GMWithAtmosphere;
end

if(nargin<2)
    elType=0;
end

numVec=size(stateVec,2);
orbEls=zeros(6,numVec);

switch(elType)
    case 0
        for curVec=1:numVec
            orbEls(:,curVec)=state2OrbElsUniv(stateVec(1:6,curVec),GM,epsVal);
        end
    case 1
        for curVec=1:numVec
            orbEls(:,curVec)=state2OrbElsEquinoctial(stateVec(1:6,curVec),1,GM);
        end
    case 2
        for curVec=1:numVec
            orbEls(:,curVec)=state2OrbElsEquinoctial(stateVec(1:6,curVec),-1,GM);
        end
    otherwise
        error('Invalid element type provided.')
end

end

function elEqui=state2OrbElsEquinoctial(state,I,GM)
%%STATE2ORBELSEQUINOCTIAL An algorithm to convert from a target state to a
%                         set of equinoctial orbital elements.
%
%The input I is the retrograde factor ond is +1 for direct equinoctial
%elements and -1 for retrograde equinoctial elements.
%
%The algorithm is taken from Section 2.1.5 of
%D. A. Danielson, C. P. Sagovac, B. Neta, and L. W. Early, "Semianalytic
%satellite theory," Mathematics Department, Naval Postgraduate School,
%Monterey, CA, Tech. Rep., 1995. [Online]. Available:
%http://oai.dtic.mil/oai/oai?verb=getRecord&metadataPrefix= html&identifier=ADA531136
%
%January 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

%position
r=state(1:3);
%velocity
rDot=state(4:6);

%Equation 1
a=1/(2/norm(r)-norm(rDot)^2/GM);%Semi-major axis.
%Equation 2
w=cross(r,rDot)/norm(cross(r,rDot));

%Equation 3
p=w(1)/(1+I*w(3));%First ascending node vector component.
q=-w(2)/(1+I*w(3));%Second ascending node vector component.

%Equation 4
e=-r/norm(r)+cross(rDot,cross(r,rDot))/GM;

%Equation 2.1.4-1
coeff=(1/(1+p^2+q^2));
f=coeff*[1-p^2+q^2;
         2*p*q;
         -2*I*p];
g=coeff*[2*I*p*q;
         (1+p^2-q^2)*I;
         2*q];

%Equation 5
h=dot(e,g);%First eccentricity vector component.
k=dot(e,f);%Second eccentricity vector component.

%Equation 6
X=dot(r,f);
Y=dot(r,g);

%Equation 2.1.4-4
b=1/(1+sqrt(1-h^2-k^2));
%Equation 7
denom=a*sqrt(1-h^2-k^2);
sinF=h+((1-h^2*b)*Y-h*k*b*X)/denom;
cosF=k+((1-k^2*b)*X-h*k*b*Y)/denom;
F=atan2(sinF,cosF);

lambda=F+h*cosF-k*sinF;%Mean longitude

elEqui=[a;h;k;p;q;lambda];

end

function elUniv=state2OrbElsUniv(state,GM,epsVal)
%%STATE2ORBELSUNIV An implementation of Gooding's conversion from a target
%                  state to universal orbital elements. This implements the
%                  PV3ELS function.
%
%The implementation is taken from Appendix D of
%R. H. Gooding, "On universal elements, and conversion procedures to
%and from position and velocity," Celestial mechanics, vol. 44, no. 3,
%pp. 283-298, 1988.
%with minor changes so that equality comparisons are replaced with
%comparisons within a region epsVal.
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.

rVec=state(1:3);
vVec=state(4:6);
x=rVec(1);
y=rVec(2);
z=rVec(3);

xyMag=norm(rVec(1:2));
r=norm(rVec);

%The projection of the velocity onto the radial component.
VR=dot(rVec,vVec)/r;%Radial velocity

%The angular momentum vector.
hVec=cross(rVec,vVec);

if(norm(hVec)^2<epsVal)%Rectilinear orbit
    %This works for an axial orbit as well as for a general rectilinear
    %orbit.
    i=pi/2;%inclination
    %For axial or general rectilinear orbits, atan2 will provide the
    %correct result. When axial and both are zero, it will return 0.
    Omega=atan2(y,x);%Longitude of the ascending node.
    %Angle from assumed reference direction (argument of latitude)
    u=atan2(z,xyMag);
    VT=0;%Transverse velocity
else%Non degenerate orbit
    b=cross(hVec,rVec);

    %r2h is r^2*h as explained before the beginning of Section 5.
    r2h=cross(rVec,b);
    W=dot(r2h(1:2),r2h(1:2));

    i=atan2(sqrt(W),r2h(3));%inclination

    if(W<epsVal)%If the orbit is in the reference plane
        Omega=0;%Longitude of the ascending node.
        u=atan2(y*sign(r2h(3)),x);%Angle from assumed reference direction
    else%General orbit
        Omega=atan2(r2h(1),-r2h(2));%Longitude of the ascending node.
        %Angle from assumed reference direction (argument of latitude)
        u=atan2(norm(r2h)*z,r^2*b(3));
    end
    VT=norm(r2h)/r^3;%Transverse velocity
end

[alpha,q,omega,tau]=PV2ELS(VR,VT,r,u,GM,epsVal);

elUniv=[alpha;
            q;
            i;
        Omega;
        omega;
          tau];
end

function [alpha,q,omega,tau]=PV2ELS(VR,VT,r,u,GM,epsVal)
%%PV2ELS  An implementation of the PV2ELS subroutine of Gooding's algorithm
%        for universal orbital element conversion. This is a
%        two-dimensional conversion suroutine.
%
%INPUTS VR, VT Radial and transerver velocity. Note VT>=0.
%       r      Radial distance.
%       u      Angle from assumed reference direction.
%       GM     Universal gravitational constant times the mass of the
%              gravitating body.
%       epsVal A small number for determining equality within finite
%              precision bounds.
%
%The algorithm is taken from Appendix C of
%R. H. Gooding, "On universal elements, and conversion procedures to
%and from position and velocity," Celestial mechanics, vol. 44, no. 3,
%pp. 283-298, 1988.
%with minor changes so that equality comparisons are not performed.
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.

sw=0.25;

VMag2=VR^2+VT^2;
alpha=2*GM/r-VMag2;

d=r*VR;
h=r*VT;
p=h^2;

esq1=p*alpha;
es=d*sqrt(abs(alpha));
ec=r*VMag2-GM;
if(alpha>0)%One formula superior for the ellipse
    e=sqrt(ec^2+es^2);
else%Different formula superior for the hyperbola
    e=sqrt(GM^2-esq1);
end

q=p/(GM+e);
if(abs(alpha)<=epsVal)%Parabola
    tau=d*(2*q+r)/(3*GM);
    v=2*atan2(VR,VT);%The true anomaly
elseif(abs(e)<=epsVal)%Circle
    tau=0;
    v=0;%The true anomaly
else%Ellipse or hyperbola
    e1=alpha*q;

    if(alpha>0)%Ellipse
        eh=atan2(es,ec);
        if(GM*eh^2/6+e1>GM*sw)%General case
            em=GM*eh-es;
            ecesq=GM*ec-e^2;
        else%For e1 and eh both near zero
            em=GM*GoodingSinDiffFun(e1/GM,eh);
            ecesq=(esq1*ec^2-e^2*es^2)/(e^2+GM*ec);
        end
    else%Hyperbola
        eh=asinh(es/e);

        if(GM*eh^2/6-e1>GM*sw)%General case
            em=es-GM*eh;
            ecesq=e^2-GM*ec;
        else%For e1 and eh both near zero
            em=e*GoodingHyperSinDiffFun(-e1/e,es/e);
            ecesq=-(esq1*ec^2+e^2*es^2)/(e^2+GM*ec);
        end
    end
    %Still ellipse or hyperbola
    en=abs(alpha)^(3/2);
    tau=em/en;
    v=atan2(es*h*sqrt(abs(alpha)),ecesq);%The true anomaly
end
%All orbits
omega=u-v;%The argument of periapsis

%Adjust revolutions, if necessary, so that omega remains in the range
%-pi to pi. The second condition in the if-statement takes care of possible
%parabolic case; this is just for the elliptical case.
if(alpha>0&&abs(alpha)>epsVal)
    adj=2*pi*fix(abs(omega/(2*pi))+(1/2))*sign(omega);
    omega=omega-adj;
    tau=tau+adj/en;
end

end

function x=GoodingSinDiffFun(e,EE)
%%GOODINGSINDIFFFUN Evaluate the function EE-(1-e)*sin(EE) using Gooding's
%                   EMKEP procedure for when e and EE are close
%                   to (1,0). This is supposed to be more accurate than
%                   just directly evaluating the functions, unless EE is
%                   large as it is then supposed to worsen rounding errors.
%
%
%The algorithm is the EMKPL algorithm taken from Appendix C of
%A. W. Odell and R. H. Gooding, "Procedure for solving Kepler's
%equation," Celestial Mechanics, vol. 38, no. 4, pp. 307-334, Apr. 1986.
%modified to solve EE-(1-e)*sin(EE) instead of EE-e*sin(EE).
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.

x=e*sin(EE);
EE2=-EE^2;
term=EE;
d=0;
while(1)
    d=d+2;
    term=term*EE2/(d*(d+1));
    x0=x;
    x=x-term;
    if(x==x0)
        break;
    end
end

end

function x=GoodingHyperSinDiffFun(g1,s)
%%GOODINGHYPERSINDIFFFUN Evaluate the function s-(1-g1)*asinh(s) when
%                        (g1,s) is close to (0,0) using Gooding's method.
%                        This is supposed to have a higher precision than
%                        just explicitly evaluating the function.
%
%The algorithm is the SHMKEP function taken from Appendix B of
%R. H. Gooding and A. W. Odell, "The hyperbolic Kepler's equation,
%and the elliptic equations revisited," Royal Aerospace Executive,
%Procurement Executive, Ministry of Defence, Farnborough, Hants, United
%Kingdom, Tech. Rep. 369, Jul. 1989.
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.

g=1-g1;
t=s/(1+sqrt(1+s^2));
x=s*(g1+g*t^2);
term=2*g*t;
twoI1=1;
%Iterate until convergence or until a maximum number of iterations is
%reached.
maxIter=64;
for curIter=1:maxIter
    twoI1=twoI1+2;
    term=term*t^2;
    x0=x;
    x=x-term/twoI1;
    if(x==x0)
        break;
    end
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function model = dnetCreate(inputDim, outputDim, Y, options)

% DNETCREATE Density network model.
% FORMAT
% DESC creates a structure for a density network.
% ARG inputDimension : dimension of latent data.
% ARG outputDim : dimension of observed data.
% ARG Y : the data to be modelled in design matrix format (as many
% rows as there are data points).
% ARG options : options structure as returned by dnetCreate.
% RETURN model : model structure containing the density network
% specified.
%
% COPYRIGHT : Neil D. Lawrence, 2008
%
% SEEALSO : dnetOptions, mlpCreate, rbfCreate, kbrCreate

% MLTOOLS

model.type = 'dnet';

if size(Y, 2) ~= outputDim
  error(['Input matrix Y does not have dimension ' num2str(d)]);
end

if isstr(options.initX)
  initFunc = str2func([options.initX 'Embed']);
  model.X = initFunc(Y, inputDim);
else
  if size(options.initX, 1) == size(Y, 1) ...
        & size(options.initX, 2) == inputDim
    model.X = options.initX;
  else
    error('options.initX not in recognisable form.');
  end
end

model.q = inputDim;
model.d = outputDim;
model.N = size(Y, 1);

model.mapping = modelCreate(options.mappingType, inputDim, ...
                            outputDim, options.mappingOptions);
model.M = options.M;
model.y = Y;
model.w = repmat(1/model.M, model.N, model.M);
model.betaTransform =  optimiDefaultConstraint('positive');

model.beta = exp(-2);
model.grid = options.grid;

if ~isempty(model.grid)
  minVals = 1.2*min(model.X);
  maxVals = 1.2*max(model.X);
  x = linspace(minVals(1), maxVals(1), model.grid(1));
  y = linspace(minVals(2), maxVals(2), model.grid(2));
  [mx, my] = meshgrid(x, y);
  for i = 1:size(model.grid)
    model.X_u = [mx(:) my(:)];
  end
else
  model.X_u = randn(model.M, model.q);
end

logw = -dist2(model.X, model.X_u);
logw = logw -repmat(max(logw, [], 2), 1, model.M);
model.w = exp(logw)./repmat(sum(exp(logw), 2), 1, model.M);

model.basisStored = options.basisStored;

model.alpha = options.alpha;

params = dnetExtractParam(model);
model = dnetExpandParam(model, params);

if model.basisStored
  model = dnetUpdateOutputWeights(model);
  model = dnetUpdateBeta(model);
end

model.numParams = model.mapping.numParams + 1;
model.outputDim = model.d;
model.inputDim = model.q;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%PCA Principal component analysis (PCA or MCA on overall covariance matrix)
%
%   [W,FRAC] = PCA(A,N)
%   [W,N]    = PCA(A,FRAC)
%
% INPUT
%   A           Dataset
%   N  or FRAC  Number of dimensions (>= 1) or fraction of variance (< 1)
%               to retain; if > 0, perform PCA; otherwise MCA. Default: N = inf.
%
% OUTPUT
%   W           Affine PCA mapping
%   FRAC or N   Fraction of variance or number of dimensions retained.
%
% DESCRIPTION
% This routine performs a principal component analysis (PCA) or minor
% component analysis (MCA) on the overall covariance matrix (weighted
% by the class prior probabilities). It finds a rotation of the dataset A to
% an N-dimensional linear subspace such that at least (for PCA) or at most
% (for MCA) a fraction FRAC of the total variance is preserved.
%
% PCA is applied when N (or FRAC) >= 0; MCA when N (or FRAC) < 0. If N is
% given (abs(N) >= 1), FRAC is optimised. If FRAC is given (abs(FRAC) < 1),
% N is optimised.
%
% Objects in a new dataset B can be mapped by B*W, W*B or by A*PCA([],N)*B.
% Default (N = inf): the features are decorrelated and ordered, but no
% feature reduction is performed.
%
% ALTERNATIVE
%
%   V = PCA(A,0)
%
% Returns the cumulative fraction of the explained variance. V(N) is the
% cumulative fraction of the explained variance by using N eigenvectors.
%
% Use KLM for a principal component analysis on the mean class covariance.
% Use FISHERM for optimizing the linear class separability (LDA).
%
% SEE ALSO
% MAPPINGS, DATASETS, PCLDC, KLLDC, KLM, FISHERM

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: pca.m,v 1.2 2006/03/08 22:06:58 duin Exp $

function [w,truefrac] = pca (varargin)

	prtrace(mfilename);

	[w,truefrac] = pcaklm(mfilename,varargin{:});

	return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [cx, cy, w, h] = getAxisAlignedBB(region)
% GETAXISALIGNEDBB extracts an axis aligned bbox from the ground truth REGION with same area as the rotated one
    cx = mean(region(1:2:end));
    cy = mean(region(2:2:end));
    x1 = min(region(1:2:end));
    x2 = max(region(1:2:end));
    y1 = min(region(2:2:end));
    y2 = max(region(2:2:end));
    A1 = norm(region(1:2) - region(3:4)) * norm(region(3:4) - region(5:6));
    A2 = (x2 - x1) * (y2 - y1);
    s = sqrt(A1/A2);
    w = s * (x2 - x1) + 1;
    h = s * (y2 - y1) + 1;
end

author: bertinetto

repo: staple

save_path: github-repos/MATLAB/bertinetto-staple


text:
classdef LevelSetSmoothRectangleInclusion < SettingsLevelSetCreator

    properties (Access = public)
        widthH
        widthV
        pnorm
    end

    methods (Access = public)

        function obj = LevelSetSmoothRectangleInclusion(varargin)
            obj.loadParams('paramsLevelSetCreator_SmoothRectangle.json')
            if nargin == 1
                obj.loadParams(varargin{1})
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [E,v] = km_kpca_icd(X,m,ktype,kpar,precision)
% KM_KPCA_ICD performs kernel principal component analysis (KPCA) on a data
% set X, using Incomplete Cholesky Decomposition to approximate the large
% kernel matrix.
%
% Input:	- X: data matrix in column format (each data point is a row)
%			- m: the number of principal components to return. If m is
%			smaller than 1, it is interpreted as the fraction of the signal
%			energy that is to be contained within the returned principal
%			components.
%			- ktype: string representing kernel type.
%			- kpar: vector containing the kernel parameters.
% Output:	- E: matrix containing the principal components.
%			- v: array containing the eigenvalues.
% USAGE: [E,v] = km_kpca_icd(X,m,ktype,kpar,precision)
%
% Author: Steven Van Vaerenbergh (steven *at* gtas.dicom.unican.es), 2010.
%
% This file is part of the Kernel Methods Toolbox for MATLAB.
% https://github.com/steven2358/kmbox

if nargin<5
	precision = 10^-6;
end
n = size(X,1);

G = km_kernel_icd(X,ktype,kpar,m,precision);
m1 = min(m,size(G,2));
R = G'*G;
[Er,Vr] = eig(R);
v = diag(Vr);

E = G*Er*diag(1./sqrt(v));
[v,ind] = sort(v,'descend');
E = E(:,ind(1:m1));	% principal components
v = v(1:m1);
for i=1:m1
	E(:,i) = E(:,i)/sqrt(n*v(i));	% normalization
end

author: epfl-lasa


text:
function laguerre_polynomial_test06 ( )

%*****************************************************************************80
%
%% LAGUERRE_POLYNOMIAL_TEST06 tests LM_POLYNOMIAL_COEFFICIENTS.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 March 2012
%
%  Author:
%
%    John Burkardt
%
  n = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAGUERRE_POLYNOMIAL_TEST06\n' );
  fprintf ( 1, '  LM_POLYNOMIAL_COEFFICIENTS determines polynomial coefficients of Lm(n,m,x).\n' );

  for m = 0 : 4

    c = lm_polynomial_coefficients ( n, m );

    for i = 0 : n
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Lm(%d,%d) = \n', i, m );
      fprintf ( 1, '\n' );
      for j = i : -1 : 0
        if ( c(i+1,j+1) == 0.0 )

        elseif ( j == 0 )
          fprintf ( 1, '  %g\n', c(i+1,j+1) );
        elseif ( j == 1 )
          fprintf ( 1, '  %g * x\n', c(i+1,j+1) );
        else
          fprintf ( 1, '  %g * x^%d\n', c(i+1,j+1), j );
        end
      end
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef testPrintingPostProcessComputer < handle

    properties (Access = private)
        iter
        fileName
        postProcessCase
        postProcessCreateData
        postProcessPrintData
    end

    methods (Access = public)

        function obj = testPrintingPostProcessComputer(d)
            obj.init(d);
            obj.loadPostProcessCreateData();
            obj.loadPostProcessPrintData();
            obj.print();
        end

        function print(obj)
            it = obj.iter;
            pC = obj.postProcessCase;
            cD = obj.postProcessCreateData;
            pD = obj.postProcessPrintData;
            postProcess = Postprocess(pC,cD);
            postProcess.print(it,pD);
        end

    end

    methods (Access = private)

        function init(obj,d)
            obj.iter            = d.iter;
            obj.fileName        = d.file;
            obj.postProcessCase = d.postCase;
        end

        function loadPostProcessCreateData(obj)
            loadData = load(obj.fileName);
            cD = loadData.postProcessData.createData;
            obj.postProcessCreateData = cD;
        end

        function loadPostProcessPrintData(obj)
            loadData = load(obj.fileName);
            pD = loadData.postProcessData.printData;
            obj.postProcessPrintData = pD;
        end

    end

end


text:
function f = resolve_function(expr)
% resolve a scoped function, which is normally not possible in MATLAB

% thus, we have to use a crazy trick: arg_report gives us access to nested functions
ofs = strfind(expr,'::@');
f = arg_report('handle',str2func(expr(1:ofs-1)),{expr(ofs+3:end)});
if iscell(f)
    if isempty(f)
        error(['No function corresponding to the expression ' expr ' was found.']);
    else
        if length(f) > 1
            warnonce(['There are multiple functions with the name ' char(f{1}) '; returning the first one.']); end
        f = f{1};
    end
end

repo: REST


text:
function cae = caetrain(cae, x, opts)
    n = cae.inputkernel(1);
    cae.rL = [];
    for m = 1 : opts.rounds
        tic;
        disp([num2str(m) '/' num2str(opts.rounds) ' rounds']);
        i1 = randi(numel(x));
        l  = randi(size(x{i1}{1},1) - opts.batchsize - n + 1);
        x1{1} = double(x{i1}{1}(l : l + opts.batchsize - 1, :, :)) / 255;

        if n == 1   %Auto Encoder
            x2{1} = x1{1};
        else        %Predictive Encoder
            x2{1} = double(x{i1}{1}(l + n : l + n + opts.batchsize - 1, :, :)) / 255;
        end
        %  Add noise to input, for denoising stacked autoenoder
        x1{1} = x1{1} .* (rand(size(x1{1})) > cae.noise);

        cae = caeup(cae, x1);
        cae = caedown(cae);
        cae = caebp(cae, x2);
        cae = caesdlm(cae, opts, m);
%         caenumgradcheck(cae,x1,x2);
        cae = caeapplygrads(cae);

        if m == 1
            cae.rL(1) = cae.L;
        end
%         cae.rL(m + 1) = 0.99 * cae.rL(m) + 0.01 * cae.L;
        cae.rL(m + 1) = cae.L;
%         if cae.sv < 1e-10
%             disp('Converged');
%             break;
%         end
        toc;
    end

end

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function [TrainData, TestData, TrainLabel, TestLabel] = createData_CCR()
    A = load(strcat('data/12_select.data'));
    B = spconvert(A);
    n1 = size(B,2);
    for i = 1:n1
        TrainY1(1,i) = 1;
    end
    A1 = load(strcat('data/17_select.data'));
    B1 = spconvert(A1);
    n2 = size(B1,2);
    for i = 1:n2
        TrainY1(1,i+n1) = -1;
    end
    TrainX1 = [B,B1];
    clear A B A1 B1;

    A = load(strcat('data/13_select.data'));
    B = spconvert(A);
    n1 = size(B,2);
    for i = 1:n1
        TrainY2(1,i) = 1;
    end
    A1 = load(strcat('data/18_select.data'));
    B1 = spconvert(A1);
    n2 = size(B1,2);
    for i = 1:n2
        TrainY2(1,i+n1) = -1;
    end
    TrainX2 = [B,B1];
    clear A B A1 B1;

    A = load(strcat('data/14_select.data'));
    B = spconvert(A);
    n1 = size(B,2);
    for i = 1:n1
        TrainY3(1,i) = 1;
    end
    A1 = load(strcat('data/19_select.data'));
    B1 = spconvert(A1);
    n2 = size(B1,2);
    for i = 1:n2
        TrainY3(1,i+n1) = -1;
    end
    TrainX3 = [B,B1];
    clear A B A1 B1;

    A = load(strcat('data/15_select.data'));
    B = spconvert(A);
    n1 = size(B,2);
    for i = 1:n1
        TestY(1,i) = 1;
    end
    A1 = load(strcat('data/20_select.data'));
    B1 = spconvert(A1);
    n2 = size(B1,2);
    for i = 1:n2
        TestY(1,i+n1) = -1;
    end
    TestX = [B,B1];
    clear A B A1 B1;

    TrainData{1,1} = TrainX1;
    TrainLabel{1,1} = TrainY1;
    TrainData{1,2} = TrainX2;
    TrainLabel{1,2} = TrainY2;
    TrainData{1,3} = TrainX3;
    TrainLabel{1,3} = TrainY3;
    TestData{1,1} = TestX;
    TestLabel{1,1} = TestY;
    save inputData.mat TrainData TrainLabel TestData TestLabel
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function complainif_notenoughargs(fnargin,limit,callfun)

if nargin<3
    callfun = mfilename;
end

if fnargin<limit
   error('%s: Not enough input arguments.',upper(callfun));
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function imout=padimage(im, sizeimout)
% imout=padimage(im, sizeimout)

sizeim = size(im);

sizedifhalf = floor((sizeimout-sizeim)/2);

imout = padarray(im, sizedifhalf);

sizedif = sizeimout-size(imout);

if sizedif(1)>0
    imout = [imout; zeros(1,size(imout,2))];
end
if sizedif(2)>0
    imout = [imout, zeros(size(imout,1),1)];
end

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
 function y = fractional_delay(x, delay)
%function y = fractional_delay(x, delay)
%|
%| given N samples x[n] of a real, periodic, band-limited signal x(t),
%| compute sinc interpolated samples of delayed signal y(t) = x(t - delay)
%| each column of x can be shifted by a different amount if delay is a vector.
%| in
%|	x	[N L]
%|	delay	[L]
%| out
%|	y	[N L]
%|
%| see laakso:96:stu for more ideas.
%|
%| Copyright 2003-11-1, Jeff Fessler, The University of Michigan
%| Extend to allow x to have multiple columns, 2003-11-2, Yingying Zhang.

if nargin < 1, ir_usage, end
if streq(x, 'test'), fractional_delay_test, return, end
if size(x,2) ~= length(delay)
	error 'Need size(x) = [N L]; length(delay) = L'
end

dims = size(x);
N = dims(1);
L = dims(2);
if length(dims) > 2, error 'x must be 1d or 2d', end
X = fft(x); % fft of each column

% it is important to choose the k indices appropriately!
if rem(N,2) % odd
	k = [-(N-1)/2:(N-1)/2]';
else
	k = [-N/2:(N/2-1)]';
end

c = exp(-1i * 2*pi/N * k * delay(:)'); % [N L] outer product

if ~rem(N,2) % even
	mid = 1;
	c(mid,:) = real(c(mid,:)); % this is the other key trick!
end

c = ifftshift(c, 1); % ifftshift differs from fftshift for odd N!

Y = X .* c;
y = ifft(Y);

% fractional_delay_test
% self test
function fractional_delay_test

Nlist = [5 6];
im clf, pl = @(i) subplot(240+i);
for ii=1:2
	N = Nlist(ii);
	n = [0:(N-1)]';
	xt = @(t, N) sinc_periodic(t, N);
	x = xt(n, N);
	xx = [x x];
	delay = [3.7; -2.2];
	y = fractional_delay(xx, delay);

	t = linspace(0,2*N,401);
	yt1 = xt(t-delay(1),N);
	yt2 = xt(t-delay(2),N);

	if im
	pl(0+ii), plot(t, xt(t,N), '-', n, xx(:,1), 'o')
	axis([0 2*N -0.4 1.1]), titlef('N=%d 1st input', N)
	pl(2+ii), plot(t, xt(t,N), '-', n, xx(:,2), 'o')
	axis([0 2*N -0.4 1.1]), titlef('N=%d 2nd input', N)
	pl(4+ii), plot(t, yt1, '-', n, real(y(:,1)), 's', n, imag(y(:,1)), '.')
	axis([0 2*N -0.4 1.1]), titlef('delay=%g 1st input',delay(1))
	pl(6+ii), plot(t, yt2, '-', n, real(y(:,2)), 's', n, imag(y(:,2)), '.')
	axis([0 2*N -0.4 1.1]), titlef('delay=%g 2nd input',delay(2))
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
filename='waferdata.bin';

%% Open
fid=fopen(filename,'rb');
NumRows=1e3;

%% Write Header
a=fread(fid,20*NumRows,'*uint8' );
a=reshape(a,20,NumRows);

%% Close
fclose(fid);
disp('Finished');


text:
function M = unroll_dbn_topology(intra, inter, T, intra1)
% UNROLL_DBN_TOPOLOGY Make the block diagonal adjacency matrix for a DBN consisting of T slices
% M = unroll_dbn_topology(intra, inter, T, intra1)
%
% intra is the connectivity within a slice, inter between two slices.
% M will have intra along the diagonal, and inter one above the diagonal.
% intra1 is an optional argumnet, in case the intra is different for the first slice.

if nargin < 4, intra1 = intra; end

ss = length(intra); % slice size
M = sparse(ss*T, ss*T);

b = 1:ss;
M(b,b) = intra1;
M(b,b+ss) = inter;

for t=2:T-1
  b = (1:ss) + (t-1)*ss;
  M(b,b) = intra;
  M(b,b+ss) = inter;
end

t = T;
b = (1:ss) + (t-1)*ss;
M(b,b) = intra;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [ geom, iner, cpmo ] = polygeom( x, y )
%POLYGEOM Geometry of a planar polygon
%
%   POLYGEOM( X, Y ) returns area, X centroid,
%   Y centroid and perimeter for the planar polygon
%   specified by vertices in vectors X and Y.
%
%   [ GEOM, INER, CPMO ] = POLYGEOM( X, Y ) returns
%   area, centroid, perimeter and area moments of
%   inertia for the polygon.
%   GEOM = [ area   X_cen  Y_cen  perimeter ]
%   INER = [ Ixx    Iyy    Ixy    Iuu    Ivv    Iuv ]
%     u,v are centroidal axes parallel to x,y axes.
%   CPMO = [ I1     ang1   I2     ang2   J ]
%     I1,I2 are centroidal principal moments about axes
%         at angles ang1,ang2.
%     ang1 and ang2 are in radians.
%     J is centroidal polar moment.  J = I1 + I2 = Iuu + Ivv

% H.J. Sommer III - 02.05.14 - tested under MATLAB v5.2
%
% sample data
% x = [ 2.000  0.500  4.830  6.330 ]';
% y = [ 4.000  6.598  9.098  6.500 ]';
% 3x5 test rectangle with long axis at 30 degrees
% area=15, x_cen=3.415, y_cen=6.549, perimeter=16
% Ixx=659.561, Iyy=201.173, Ixy=344.117
% Iuu=16.249, Ivv=26.247, Iuv=8.660
% I1=11.249, ang1=30deg, I2=31.247, ang2=120deg, J=42.496
%
% H.J. Sommer III, Ph.D., Professor of Mechanical Engineering, 337 Leonhard Bldg
% The Pennsylvania State University, University Park, PA  16802
% (814)863-8997  FAX (814)865-9693  hjs1@psu.edu  www.me.psu.edu/sommer/

% begin function POLYGEOM

% check if inputs are same size
if ~isequal( size(x), size(y) ),
  error( 'X and Y must be the same size');
end

% number of vertices
[ x, ns ] = shiftdim( x );
[ y, ns ] = shiftdim( y );
[ n, c ] = size( x );

% temporarily shift data to mean of vertices for improved accuracy
xm = mean(x);
ym = mean(y);
x = x - xm*ones(n,1);
y = y - ym*ones(n,1);

% delta x and delta y
dx = x( [ 2:n 1 ] ) - x;
dy = y( [ 2:n 1 ] ) - y;

% summations for CW boundary integrals
A = sum( y.*dx - x.*dy )/2;
Axc = sum( 6*x.*y.*dx -3*x.*x.*dy +3*y.*dx.*dx +dx.*dx.*dy )/12;
Ayc = sum( 3*y.*y.*dx -6*x.*y.*dy -3*x.*dy.*dy -dx.*dy.*dy )/12;
Ixx = sum( 2*y.*y.*y.*dx -6*x.*y.*y.*dy -6*x.*y.*dy.*dy ...
          -2*x.*dy.*dy.*dy -2*y.*dx.*dy.*dy -dx.*dy.*dy.*dy )/12;
Iyy = sum( 6*x.*x.*y.*dx -2*x.*x.*x.*dy +6*x.*y.*dx.*dx ...
          +2*y.*dx.*dx.*dx +2*x.*dx.*dx.*dy +dx.*dx.*dx.*dy )/12;
Ixy = sum( 6*x.*y.*y.*dx -6*x.*x.*y.*dy +3*y.*y.*dx.*dx ...
          -3*x.*x.*dy.*dy +2*y.*dx.*dx.*dy -2*x.*dx.*dy.*dy )/24;
P = sum( sqrt( dx.*dx +dy.*dy ) );

% check for CCW versus CW boundary
if A < 0,
  A = -A;
  Axc = -Axc;
  Ayc = -Ayc;
  Ixx = -Ixx;
  Iyy = -Iyy;
  Ixy = -Ixy;
end

% centroidal moments
xc = Axc / A;
yc = Ayc / A;
Iuu = Ixx - A*yc*yc;
Ivv = Iyy - A*xc*xc;
Iuv = Ixy - A*xc*yc;
J = Iuu + Ivv;

% replace mean of vertices
x_cen = xc + xm;
y_cen = yc + ym;
Ixx = Iuu + A*y_cen*y_cen;
Iyy = Ivv + A*x_cen*x_cen;
Ixy = Iuv + A*x_cen*y_cen;

% principal moments and orientation
I = [ Iuu  -Iuv ;
     -Iuv   Ivv ];
[ eig_vec, eig_val ] = eig(I);
I1 = eig_val(1,1);
I2 = eig_val(2,2);
ang1 = atan2( eig_vec(2,1), eig_vec(1,1) );
ang2 = atan2( eig_vec(2,2), eig_vec(1,2) );

% return values
geom = [ A  x_cen  y_cen  P ];
iner = [ Ixx  Iyy  Ixy  Iuu  Ivv  Iuv ];
cpmo = [ I1  ang1  I2  ang2  J ];

% end of function POLYGEOM


text:
% Script to reproduce the experiments leading to the results provided in the
% Table 2 of the paper "Deep Scattering Spectrum" by J. Andén and S. Mallat.

% M=2 scattering

run_name = 'DSS_Table2_GTZAN_m2';

N=5*2^17;

src=gtzan_src('/path/to/gtzan');

fparam.filter_type = {'gabor_1d','morlet_1d'};
fparam.Q = [8 2];
fparam.J = T_to_J(8192,fparam);

options.M = 2;

Wop = wavelet_factory_1d(N, fparam, options);

feature_fun = {@(x)(format_scat(log_scat(renorm_scat(scat(x,Wop)))))};

db = prepare_database(src,feature_fun);
db.features = single(db.features);
db = svm_calc_kernel(db,'gaussian','square',1:2:size(db.features,2));

rs = RandStream.create('mt19937ar','Seed',floor(pi*1e9));
RandStream.setGlobalStream(rs);
[train_set{1}, test_set{1}] = create_partition([src.objects.class], 0.9);
for k = 2:10
	[train_set{k}, test_set{k}] = ...
		next_fold([src.objects.class], train_set{k-1}, test_set{k-1});
end

optt.kernel_type = 'gaussian';
optt.C = 2.^[0:4:8];
optt.gamma = 2.^[-16:4:-8];
optt.search_depth = 3;
optt.full_test_kernel = 1;

for k = 1:10
	[dev_err_grid,C_grid,gamma_grid] = ...
		svm_adaptive_param_search(db,train_set{k},[],optt);

	[dev_err(k),ind] = min(mean(dev_err_grid{end},2));
	C(k) = C_grid{end}(ind);
	gamma(k) = gamma_grid{end}(ind);

	optt1 = optt;
	optt1.C = C(k);
	optt1.gamma = gamma(k);

	model = svm_train(db,train_set{k},optt1);
	labels(:,k) = svm_test(db,model,test_set{k});
	err(k) = classif_err(labels(:,k),test_set{k},db.src);

	fprintf('dev err = %f, test err = %f\n',dev_err(k),err(k));

	save([run_name '.mat'],'labels','dev_err','err','C','gamma');
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function [ node_xyz, face_node, face_normal ] = stla_read ( ...
  input_file_name, node_num, face_num )

%*****************************************************************************80
%
%% STLA_READ reads graphics information from an ASCII StereoLithography file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    3D Systems, Inc,
%    Stereolithography Interface Specification,
%    October 1989.
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the name of the input file.
%
%    Input, integer NODE_NUM, the number of vertices defined.
%
%    Input, integer FACE_NUM, the number of faces defined.
%
%    Output, real NODE_XYZ(3,NODE_NUM), the coordinates of points.
%
%    Output, integer FACE_NODE(3,FACE_NUM), the nodes that make up each face.
%
%    Output, real FACE_NORMAL(3,FACE_NUM), the normal vector
%    at each face.
%
  node_xyz = [];
  face_node = [];
  face_normal = [];
%
%  If no file input, try to get one from the user.
%
  if ( nargin < 1 )
    input_file_name = input ( 'Enter the name of the ASCII STL file.' );
    if ( isempty ( input_file_name ) )
      return
    end
  end
%
%  If no size input, scan the file.
%
  if ( nargin < 3 )
    [ solid_num, node_num, face_num, text_num ] = stla_size ( ...
      input_file_name );
  end
%
%  Open the file.
%
  input_file_unit = fopen ( input_file_name, 'r' );

  if ( input_file_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'STLA_READ - Fatal error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'STLA_READ - Fatal error!' )
    return
  end
%
%  Allocate space.
%
  state = 0;
  text_num = 0;

  solid = 0;
  face = 0;
  node = 0;

  node_xyz = zeros ( 3, node_num );
  face_node = zeros ( 3, face_num );
  face_normal = zeros ( 3, face_num );
%
%  Read the next line of text.
%
  while ( 1 )

    text = fgetl ( input_file_unit );
%
%  On end of input, FGETL returns a -1 instead of a character array.
%
    if ( text == -1 )
      if ( state ~= 0 & state ~= 1 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  End-of-file, but model not finished.\n' );
        error ( 'STLA_READ - Fatal error!' );
        return
      end
      break;
    end

    text_num = text_num + 1;
    done = 1;
%
%  Read the first word in the line.
%
    [ word1, done ] = word_next_read ( text, done );

    if ( done )
      break
    end
%
%  "Doctor" the text, changing a beginning occurrence of:
%
%      END FACET to ENDFACET
%      END LOOP to ENDLOOP
%      END SOLID to ENDSOLID
%      FACET NORMAL to FACETNORMAL
%      OUTER LOOP to OUTERLOOP
%
    if ( s_eqi ( word1, 'END' ) )

      [ word2, done ] = word_next_read ( text, done );

      if ( ~s_eqi ( word2, 'FACET' ) & ...
           ~s_eqi ( word2, 'LOOP' ) & ...
           ~s_eqi ( word2, 'SOLID' ) )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  The tag END was followed by an illegal\n' );
        fprintf ( 1, '  word: "%s", when expecting\n', word );
        fprintf ( 1, '  "FACET", "LOOP", or "SOLID".\n' );
        error ( 'STLA_READ - Fatal error!' );
        return
      end

      word1 = s_cat ( word1, word2 );

    elseif ( s_eqi ( word1, 'FACET' ) )

      [ word2, done ] = word_next_read ( text, done );

      if ( ~s_eqi ( word2, 'NORMAL' ) )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  The tag FACET was followed by an illegal\n' );
        fprintf ( 1, '  word: "%s", when expecting\n', word2 );
        fprintf ( 1, '  "NORMAL".\n' );
        error ( 'STLA_READ - Fatal error!' );
        return
      end

      word1 = s_cat ( word1, word2 );

    elseif ( s_eqi ( word1, 'OUTER' ) )

      [ word2, done ] = word_next_read ( text, done );

      if ( ~s_eqi ( word2, 'LOOP' ) )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  The tag OUTER was followed by an illegal\n' );
        fprintf ( 1, '  word: "%s", when expecting\n', word2 );
        fprintf ( 1, '  "LOOP".\n' );
        error ( 'STLA_READ - Fatal error!' );
        return
      end

      word1 = s_cat ( word1, word2 );

    end
%
%  This first word tells us what to do.
%
%  SOLID - begin a new solid.
%    Valid in state 0, moves to state 1.
%  ENDSOLID - end current solid.
%    Valid in state 1, moves to state 0.
%
%  FACETNORMAL - begin a new facet.
%    Valid in state 0 or 1, moves to state 2.
%  ENDFACET - end current facet.
%    Valid in state 2, moves to state 1.
%
%  OUTERLOOP - begin a list of vertices.
%    Valid in state 2, moves to state 3, sets vertex count to 0.
%  ENDLOOP - end vertex list.
%    Valid in state 3, moves to state 2.
%
%  VERTEX - give coordinates of next vertex.
%    Valid in state 3.
%
%  End of file -
%    Valid in state 0 or 1.
%
    if ( s_eqi ( word1, 'SOLID' ) )

      if ( state ~= 0 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  A new SOLID statement was encountered, but we\n' );
        fprintf ( 1, '  have not finished processing the current solid.\n' );
        error ( 'STLA_READ - Fatal error!' );
        return
      end

      state = 1;

    elseif ( s_eqi ( word1, 'ENDSOLID' ) )

      if ( state ~= 1 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  An END SOLID statement was encountered, but\n' );
        fprintf ( 1, '  either we have not begun a solid at all, or we\n' );
        fprintf ( 1, '  are not at an appropriate point to finish the\n' );
        fprintf ( 1, '  current solid.\n' );
        error ( 'STLA_READ - Fatal error!' );
        return
      end

      state = 0;

    elseif ( s_eqi ( word1, 'FACETNORMAL' ) )

      if ( state ~= 0 & state ~= 1 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  Model not in right state for FACET.\n' );
        error ( 'STLA_READ - Fatal error!' );
        return
      end

      state = 2;
      face = face + 1;

      if ( face_num < face )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Warning!\n' );
        fprintf ( 1, '  More faces being read than expected.\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        error ( 'STLA_READ - Fatal error!' );
        return
      end

      for i = 1 : 3

        face_normal(i,face) = 0.0;

        [ word2, done ] = word_next_read ( text, done );

        if ( done )
          fprintf ( 1, '\n' );
          fprintf ( 1, 'STLA_READ - Fatal error!\n' );
          fprintf ( 1, '  File line number = %d\n', text_num );
          fprintf ( 1, '  End of information while reading a component\n' );
          fprintf ( 1, '  of the normal vector.\n' );
          error ( 'STLA_READ - Fatal error!' );
          return
        end

        [ dval, lchar, ierror ] = s_to_r8 ( word2 );

        if ( ierror ~= 0 )
          fprintf ( 1, '\n' );
          fprintf ( 1, 'STLA_READ - Fatal error!\n' );
          fprintf ( 1, '  File line number = %d\n', text_num );
          fprintf ( 1, '  Reading a component of the normal vector.\n' );
          error ( 'STLA_READ - Fatal error!' );
          return
        end

        face_normal(i,face) = dval;

      end

    elseif ( s_eqi ( word1, 'ENDFACET' ) )

      if ( state ~= 2 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  Model not in right state for ENDFACET.\n' );
        error ( 'STLA_READ - Fatal error!' );
        return
      end

      state = 1;

    elseif ( s_eqi ( word1, 'OUTERLOOP' ) )

      if ( state ~= 2 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  Model not in right state for OUTERLOOP.\n' );
        error ( 'STLA_READ - Fatal error!' );
        return
      end

      state = 3;
      vertex = 0;

    elseif ( s_eqi ( word1, 'ENDLOOP' ) )

      if ( state ~= 3 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  Model not in right state for ENDLOOP.\n' );
        error ( 'STLA_READ - Fatal error!' );
        return
      end

      state = 2;

    elseif ( s_eqi ( word1, 'VERTEX' ) )

      if ( state ~= 3 )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  Model not in right state for VERTEX.\n' );
        error ( 'STLA_READ - Fatal error!' );
        return
      end

      if ( 3 <= vertex )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Fatal error!\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        fprintf ( 1, '  More than 3 vertices specified for a face.\n' );
        error ( 'STLA_READ - Fatal error!' );
        return
      end

      for i = 1 : 3

        [ word2, done ] = word_next_read ( text, done );

        if ( done )
          fprintf ( 1, '\n' );
          fprintf ( 1, 'STLA_READ - Fatal error!\n' );
          fprintf ( 1, '  File line number = %d\n', text_num );
          fprintf ( 1, '  The value of a vertex coordinate is missing.\n' );
          error ( 'STLA_READ - Fatal error!' );
          return
        end

        [ dval, lchar, ierror ] = s_to_r8 ( word2 );

        if ( ierror ~= 0 )
          fprintf ( 1, '\n' );
          fprintf ( 1, 'STLA_READ - Fatal error!\n' );
          fprintf ( 1, '  File line number = %d\n', text_num );
          fprintf ( 1, '  The value of a vertex coordinate makes\n' );
          fprintf ( 1, '  no sense.\n' );
          error ( 'STLA_READ - Fatal error!' );
          return
        end

        temp(i) = dval;

      end

      if ( node_num <= node )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'STLA_READ - Warning!\n' );
        fprintf ( 1, '  More nodes being read than expected.\n' );
        fprintf ( 1, '  File line number = %d\n', text_num );
        error ( 'STLA_READ - Fatal error!' );
        return
      end

      node = node + 1;
      node_xyz(1:3,node) = temp(1:3);

      vertex = vertex + 1;
      face_node(vertex,face) = node;

    else

      fprintf ( 1, '\n' );
      fprintf ( 1, 'STLA_READ - Fatal error!\n' );
      fprintf ( 1, '  File line number = %d\n', text_num );
      fprintf ( 1, '  Unrecognized line in file.\n' );
      error ( 'STLA_READ - Fatal error!' );
      return

    end

  end
%
%  Close the file.
%
  fclose ( input_file_unit );

  return
end


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD 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.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

addpath(genpath('.')); init_workspace;

opt.source          = struct('camera',0,'input','_input/','bb0',[]); % camera/directory swith, directory_name, initial_bounding_box (if empty, it will be selected by the user)
opt.output          = '_output/'; mkdir(opt.output); % output directory that will contain bounding boxes + confidence

min_win             = 24; % minimal size of the object's bounding box in the scanning grid, it may significantly influence speed of TLD, set it to minimal size of the object
patchsize           = [15 15]; % size of normalized patch in the object detector, larger sizes increase discriminability, must be square
fliplr              = 0; % if set to one, the model automatically learns mirrored versions of the object
maxbbox             = 1; % fraction of evaluated bounding boxes in every frame, maxbox = 0 means detector is truned off, if you don't care about speed set it to 1
update_detector     = 1; % online learning on/off, of 0 detector is trained only in the first frame and then remains fixed
opt.plot            = struct('pex',1,'nex',1,'dt',1,'confidence',1,'target',1,'replace',0,'drawoutput',3,'draw',0,'pts',1,'help', 0,'patch_rescale',1,'save',0);

% Do-not-change -----------------------------------------------------------

opt.model           = struct('min_win',min_win,'patchsize',patchsize,'fliplr',fliplr,'ncc_thesame',0.95,'valid',0.5,'num_trees',10,'num_features',13,'thr_fern',0.5,'thr_nn',0.65,'thr_nn_valid',0.7);
opt.p_par_init      = struct('num_closest',10,'num_warps',20,'noise',5,'angle',20,'shift',0.02,'scale',0.02); % synthesis of positive examples during initialization
opt.p_par_update    = struct('num_closest',10,'num_warps',10,'noise',5,'angle',10,'shift',0.02,'scale',0.02); % synthesis of positive examples during update
opt.n_par           = struct('overlap',0.2,'num_patches',100); % negative examples initialization/update
opt.tracker         = struct('occlusion',10);
opt.control         = struct('maxbbox',maxbbox,'update_detector',update_detector,'drop_img',1,'repeat',1);

% Run TLD -----------------------------------------------------------------
%profile on;
[bb,conf] = tldExample(opt);
%profile off;
%profile viewer;

% Save results ------------------------------------------------------------
dlmwrite([opt.output '/tld.txt'],[bb; conf]');
disp('Results saved to ./_output.');

repo: OpenTLD


text:
%%*********************************************************************
%% norm_min: matrix 2-norm minimization problem
%%
%% (primal) minimize || B0 + x1*B1 + ... + xm*Bm ||_2
%%
%% (dual)   maximize    - Tr B0*Q
%%          subject to  Tr Bk*Q = 0, k=1,...,m
%%                      sum of singular values of Q is less than one.
%%
%% The matrices Bk are of size p x q.
%%
%% The problem is equivalent to the following SDP:
%%
%%    minimize    t
%%                         [0,   Bk]
%%    s.t.  sum_{k=1}^m xk*[       ]
%%                         [Bk^*, 0]
%%                               [0,     i*Bk]          [0,   B0]
%%         + sum_{k=m+1}^{2m} xk*[           ] -t*I <= -[       ]
%%                               [(i*Bk)^*, 0]          [B0^*, 0]
%%
%% Adapted from Vandenberghe and Boyd, SIAM Review 96.
%% see also igmres.m
%%--------------------------------------------------------------------
%%
%%  [blk,Avec,C,b,X0,y0,Z0,objval,x] = norm_min(B,feas,solve);
%%
%%  Input:   B = cell array such that B{k} = B_k.
%%           feas  = 1 if want feasible starting point
%%                 = 0 if otherwise.
%%           solve = 0 just to initialize
%%                 = 1 if want to solve the problem.
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%*********************************************************************

   function [blk,Avec,C,b,X0,y0,Z0,objval,x] = norm_min(B,feas,solve);

   if (nargin < 3); solve = 0; end;
   if (nargin < 2); feas = 0; end;
   if isempty(feas); feas = 0; end
%%
%% dimensions
%%
   m = size(B,2)-1;
   [p,q] = size(B{1}); isrealB = 1;
   for k = 1:m;
       [p1,q1] = size(B{k});
       if any([p q]-[p1 q1]);
          error('size of B_k is not the same for all k');
       end
       isrealB = min(isrealB, isreal(B{k}));
   end
   cmp = 1-isrealB;
%%
%% A{k} = [0 Bk; Bk' 0],  i=1,...,m+1;  A{m+1} = I
%%
   if (~cmp)
       n = p+q;
       blk{1,1} = 's'; blk{1,2} = n;
       b = [zeros(m,1); -1];
       C = -[zeros(p), B{1}; B{1}', zeros(q)];
       A = cell(1,m+1);
       for k = 1:m
           Bk = B{k+1};
           A{k} = [zeros(p), Bk; Bk', zeros(q)];
       end
       A{m+1} = -speye(n,n);
       Avec = svec(blk,A,ones(size(blk,1),1));
   elseif (cmp)
       n = p+q;
       blktmp{1,1} = 's'; blktmp{1,2} = n;
       b = [zeros(2*m,1); -1];
       Ctmp{1} = -[zeros(p),  B{1}; B{1}', zeros(q)];
       A = cell(1,2*m+1); ii = sqrt(-1);
       for k = 1:m
           Bk = B{k+1};
           A{k}   = [zeros(p), Bk; Bk', zeros(q)];
           A{k+m} = [zeros(p), ii*Bk; -ii*Bk', zeros(q)];
       end
       A{2*m+1} = -speye(n,n);
       [blk,Avec,C,b] = convertcmpsdp(blktmp,A,Ctmp,b);
   end
%%
%%
   if (feas == 1);
      if (~cmp)
         X0{1} = eye(n)/n;
      else
         X0{1} = eye(2*n)/(2*n);
      end
      y0 = 1.1*ops(C,'norm')*[zeros(length(b)-1,1); 1];
      Z0 = ops(C,'-',Atyfun(blk,Avec,[],[],y0));
   elseif (feas == 0)
      [X0,y0,Z0] = infeaspt(blk,Avec,C,b);
   end
   if (solve)
      [obj,X,y,Z] = sqlp(blk,Avec,C,b,[],X0,y0,Z0);
      objval = -mean(obj);
      if ~cmp; x = y(1:m);
      else;    x = y(1:m) +sqrt(-1)*y(m+1:2*m);
      end;
   else
      objval = []; x = [];
   end
%%*********************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function rmdirquiet(x)

% function rmdirquiet(x)
%
% <x> refers to a file or directory location
%
% Delete <x>, suppressing warnings.
% If <x> doesn't exist, we just quietly return.
%
% example:
% mkdirquiet('test');
% mkdirquiet('test');
% rmdirquiet('test');
% rmdirquiet('test');

if exist(x,'dir')
  assert(rmdir(x,'s'),sprintf('rmdir of %s failed',x));
elseif exist(x,'file')
  delete(x);
end

author: cvnlab

repo: GLMsingle

save_path: github-repos/MATLAB/cvnlab-GLMsingle


text:
function gX = ardKernGradX(kern, X, X2)

% ARDKERNGRADX Gradient of ARD kernel with respect to a point x.
% FORMAT
% DESC computes the gradient of the pre-built RBF and linear ARD
% kernel with respect to the input positions.
% ARG kern : kernel structure for which gradients are being
% computed.
% ARG x : locations against which gradients are being computed.
% RETURN g : the returned gradients. The gradients are returned in
% a matrix which is numData x numInputs x numData. Where numData is
% the number of data points and numInputs is the number of input
% dimensions in X.
%
% FORMAT
% DESC computes the gradident of the pre-built RBF and linear ARD
% kernel with respect to the input positions where both the row
% positions and column positions are provided separately.
% ARG kern : kernel structure for which gradients are being
% computed.
% ARG x1 : row locations against which gradients are being computed.
% ARG x2 : column locations against which gradients are being computed.
% RETURN g : the returned gradients. The gradients are returned in
% a matrix which is numData2 x numInputs x numData1. Where numData1 is
% the number of data points in X1, numData2 is the number of data
% points in X2 and numInputs is the number of input
% dimensions in X.
%
% SEEALSO ardKernParamInit, kernGradX, ardKernDiagGradX
%
% COPYRIGHT : Neil D. Lawrence, 2004

% KERN

gX = zeros(size(X2, 1), size(X2, 2), size(X, 1));
for i = 1:size(X, 1);
  gX(:, :, i) = ardKernGradXpoint(kern, X(i, :), X2);
end

function gX = ardKernGradXpoint(kern, x, X2)

% ARDKERNGRADXPOINT Gradient with respect to one point of x.

scales = sparse(diag(kern.inputScales));

gX = kern.linearVariance.*X2*scales;

scales = sqrt(scales);

n2 = dist2(X2*scales, x*scales);
wi2 = (.5 .* kern.inverseWidth);
rbfPart = kern.rbfVariance*exp(-n2*wi2);
for i = 1:size(x, 2)
  gX(:, i) = gX(:, i) + kern.inverseWidth*kern.inputScales(i)*(X2(:, i) - x(i)).*rbfPart;
end


text:
  function ob = block_fatrix(blocks, varargin)
%|function ob = block_fatrix(blocks, options)
%|
%| Construct block_fatrix object, a meta-Fatrix composed of Fatrix blocks,
%| such as block_diag(A_1, A_2, ..., A_M)
%| See block_fatrix_test.m for example usage.
%|
%| in
%|	blocks	{cell}	cell array of the blocks
%|
%| options
%|	'type'	char	options:
%|				'diag' for block diagonal (default)
%|				'col' for [A_1; A_2; ...; A_M]
%|				'kron' for kron(eye(Mkron), blocks{1})
%|				'row' for [A1, A2, ..., A_M]
%|				'sum' for A1 + A2 + ... + A_M
%|	'Mkron'	int	required for 'kron' type
%|	'tomo'	0|1	special support for tomography-type objects (default: 0)
%|
%| out
%|	ob [nd np]	diag,row: np = sum_m ncol(A_m)
%|
%| Copyright 05-5-12, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if nargin == 1 && streq(blocks, 'test')
	run_mfile_local block_fatrix_test
return
end

arg.blocks = blocks;

% defaults
arg.type = 'diag';
arg.chat = 0;
arg.tomo = 0;
arg.Mkron = [];

% options
arg = vararg_pair(arg, varargin);

% initialize
if ~iscell(blocks)
	fail('blocks must be cell array, not "%s"', class(blocks))
end

switch arg.type
case 'col'
	ob = block_fatrix_col(blocks, arg);
case 'diag'
	ob = block_fatrix_diag(blocks, arg);
case 'kron'
	ob = block_fatrix_kron(blocks, arg);
case 'row'
	ob = block_fatrix_row(blocks, arg);
case 'sum'
	ob = block_fatrix_sum(blocks, arg);
otherwise
	fail('unknown block type "%s"', arg.type)
end

%
% block_fatrix_col()
%
function ob = block_fatrix_col(blocks, arg)

MM = length(blocks);
arg.dims = zeros(MM, 2);
for ii=1:MM
	arg.dims(ii,:) = size(blocks{ii});
	if arg.dims(ii,2) ~= arg.dims(1,2)
		error 'all blocks must have same #cols for "col"'
	end
end
% start/end indices for selecting parts of x and y
arg.istart = cumsum([1; arg.dims(1:end-1,1)]);
arg.iend = arg.istart + arg.dims(:,1) - 1;

arg.dim = [sum(arg.dims(:,1)) arg.dims(1,2)];

if arg.tomo % prep for mat2cell later
	for ii=1:MM
		odim = arg.blocks{ii}.arg.odim;
		arg.tomo_ndim = length(odim);
		arg.tomo_dim = arg.tomo_ndim; % insist on last dimension
		if ii==1
			for id = 1:arg.tomo_ndim
				if id ~= arg.tomo_dim
					arg.mat2cell_arg{id} = odim(id);
				end
			end
		end
		arg.mat2cell_arg{arg.tomo_dim}(ii) = odim(arg.tomo_dim);
	end
end

% build Fatrix object
ob = Fatrix(arg.dim, arg, 'caller', 'block_fatrix(col)', ...
	'forw', @block_fatrix_col_forw, 'back', @block_fatrix_col_back, ...
	'free', @block_fatrix_free, 'gram', @block_fatrix_col_gram);

%
% block_fatrix_col_forw(): y = A * x
%
function y = block_fatrix_col_forw(arg, x)

MM = length(arg.blocks);
y = cell(MM,1);
for ii=1:MM
	y{ii} = arg.blocks{ii} * x;
end

if ~arg.tomo || ncol(y{1}) == 1 ...
	|| (arg.tomo && arg.dim(2) == size(x,1)) % [np (L)]
	y = cat(1, y{:});
else % handle 'tomo' case where x is not a single column
	y = cat(arg.tomo_dim, y{:});
end

%
% block_fatrix_col_back(): x = A' * y
%
function x = block_fatrix_col_back(arg, y)

MM = length(arg.blocks);
if arg.tomo
	if arg.dim(1) == size(y,1) % [nd (L)]
		arg.tomo = 0;
	else % [(odim) (L)] split into separate parts, undoing cat(?, y{:})
		tmp = size(y);
		if length(tmp) == arg.tomo_ndim
			yc = mat2cell(y, arg.mat2cell_arg{:});
		elseif length(tmp) > arg.tomo_ndim
			yc = mat2cell(y, arg.mat2cell_arg{:}, ...
				tmp((arg.tomo_ndim+1):end)); % handle (L)
		else
			error 'bug'
		end
	end
end

x = 0;
for ii=1:MM
	if ~arg.tomo % [nd (L)]
		jj = arg.istart(ii):arg.iend(ii);
		yi = y(jj,:);
	else
		yi = yc{ii};
	end
	t = arg.blocks{ii}' * yi;
	x = x + t;
end

%
% block_fatrix_col_gram()
%
function [T, reuse] = block_fatrix_col_gram(ob, W, reuse, varargin)

blocks = ob.arg.blocks;
T = cell(size(blocks));
for ii=1:length(blocks)
	A = blocks{ii};
	if isnumeric(A)
		if isempty(W)
			T{ii} = A' * A;
		else
			warn 'todo: this may not work, need piece of W!'
			T{ii} = A' * W * A;
		end
	else
		if isempty(W)
			T{ii} = build_gram(A, [], reuse, varargin{:});
		else
			if isvar('W.arg.blocks{ii}')
				T{ii} = build_gram(A, W.arg.blocks{ii}, ...
					reuse, varargin{:});
			else
				fail('block_fatrix_col_gram needs block diag W')
			end
		end
	end
end
T = fatrix_plus(T{:});

%
% block_fatrix_row()
% trick: just reuse col via transpose!
%
function ob = block_fatrix_row(blocks, arg)

for ii=1:length(blocks)
	blocks{ii} = blocks{ii}'; % trick: transpose
end

ob = block_fatrix(blocks, 'type', 'col')'; % trick: transpose

%
% block_fatrix_diag()
%
function ob = block_fatrix_diag(blocks, arg)

arg.dims = zeros(length(blocks), 2);
for ii=1:length(blocks)
	arg.dims(ii,:) = size(blocks{ii});
end
% start/end indices for selecting parts of x and y
arg.istart = cumsum([1; arg.dims(1:end-1,1)]);
arg.iend = arg.istart + arg.dims(:,1) - 1;
arg.jstart = cumsum([1; arg.dims(1:end-1,2)]);
arg.jend = arg.jstart + arg.dims(:,2) - 1;

arg.dim = sum(arg.dims, 1);

% build Fatrix object
ob = Fatrix(arg.dim, arg, 'caller', 'block_fatrix(diag)', ...
	'forw', @block_fatrix_diag_forw, 'back', @block_fatrix_diag_back, ...
	'free', @block_fatrix_free, 'gram', @block_fatrix_diag_gram, ...
	'mtimes_block', @block_fatrix_diag_mtimes_block);

%
% block_fatrix_diag_forw(): y = A * x
%
function y = block_fatrix_diag_forw(arg, x)

if nrow(x) ~= arg.dim(2)
	error('x size=%d vs dim(2)=%d', nrow(x), arg.dim(2))
end
y = [];
for ii=1:length(arg.blocks)
	t = arg.blocks{ii} * x([arg.jstart(ii):arg.jend(ii)], :);
	y = [y; t];
end

%
% block_fatrix_diag_back(): x = A' * y
% full backprojection
%
function x = block_fatrix_diag_back(arg, y)

if nrow(y) ~= arg.dim(1), error 'bad y size', end
x = [];
for ii=1:length(arg.blocks)
	t = arg.blocks{ii}' * y([arg.istart(ii):arg.iend(ii)], :);
	x = [x; t];
end

%
% block_fatrix_diag_mtimes_block()
% caution: it is quite unclear whether these are useful!
%
function y = block_fatrix_diag_mtimes_block(arg, is_transpose, x, istart, nblock)
if is_transpose
	y = block_fatrix_diag_block_back(arg, x, istart, nblock);
else
	y = block_fatrix_diag_block_forw(arg, x, istart, nblock);
end

%
% block_fatrix_diag_block_forw()
%
function y = block_fatrix_diag_block_forw(arg, x, istart, nblock)

if nrow(x) ~= arg.dim(2)
	error('x size=%d vs dim(2)=%d', nrow(x), arg.dim(2))
end
y = [];
for ii=istart:nblock:length(arg.blocks)
	t = arg.blocks{ii} * x([arg.jstart(ii):arg.jend(ii)], :);
	y = [y; t];
end

%
% block_fatrix_diag_block_back()
%
function x = block_fatrix_diag_block_back(arg, y, istart, nblock)

if nrow(y) ~= arg.dim(1), error 'bad y size', end
x = [];
for ii=istart:nblock:length(arg.blocks)
	t = arg.blocks{ii}' * y([arg.istart(ii):arg.iend(ii)], :);
	x = [x; t];
end

%
% block_fatrix_diag_gram()
%
function [T, reuse] = block_fatrix_diag_gram(ob, W, reuse, varargin)

blocks = ob.arg.blocks;
T = cell(size(blocks));
for ii=1:length(blocks)
	A = blocks{ii};
	if isnumeric(A)
		if isempty(W)
			T{ii} = A' * A;
		else
			T{ii} = A' * W * A;
		end
	else
		T{ii} = build_gram(A, W, reuse, varargin{:});
	end
end
T = block_fatrix(T, 'type', 'diag');

%
% block_fatrix_kron()
%
function ob = block_fatrix_kron(blocks, arg)

if isempty(arg.Mkron), error 'Mkron required', end
if length(blocks) ~= 1, error 'kron expects exactly one block', end

arg.dims = repmat(size(blocks{1}), [arg.Mkron 1]);
% start/end indices for selecting parts of x and y
arg.istart = cumsum([1; arg.dims(1:end-1,1)]);
arg.iend = arg.istart + arg.dims(:,1) - 1;
arg.jstart = cumsum([1; arg.dims(1:end-1,2)]);
arg.jend = arg.jstart + arg.dims(:,2) - 1;
arg.dim = sum(arg.dims,1);

% build Fatrix object
if isa(blocks{1}, 'Fatrix')
	tmp = ['F:' blocks{1}.caller];
else
	tmp = class(blocks{1});
end
tmp = sprintf('block_fatrix(kron, %s)', tmp);
ob = Fatrix(arg.dim, arg, 'caller', tmp, ...
	'forw', @block_fatrix_kron_forw, 'back', @block_fatrix_kron_back, ...
	'block_setup', @block_fatrix_kron_block_setup, ...
	'mtimes_block', @block_fatrix_kron_mtimes_block, ...
	'free', @block_fatrix_free);

%
% block_fatrix_kron_forw(): y = A * x
%
function y = block_fatrix_kron_forw(arg, x)

if nrow(x) ~= arg.dim(2)
	fail('x size=%d vs dim(2)=%d', nrow(x), arg.dim(2))
end
y = [];
for ii=1:arg.Mkron
	t = arg.blocks{1} * x([arg.jstart(ii):arg.jend(ii)], :);
	y = [y; t];
end

%
% block_fatrix_kron_back(): x = A' * y
% full backprojection
%
function x = block_fatrix_kron_back(arg, y)

if nrow(y) ~= arg.dim(1), error 'bad y size', end
x = [];
for ii=1:arg.Mkron
	t = arg.blocks{1}' * y([arg.istart(ii):arg.iend(ii)], :);
	x = [x; t];
end

%
% block_fatrix_kron_block_setup()
% apply Gblock to the base block, to prepare it for later
%
function ob = block_fatrix_kron_block_setup(ob, varargin)
ob.arg.blocks = {Gblock(ob.arg.blocks{1}, ob.nblock)};
% todo: probably need to set up some other internal variables here
% for use inside block_fatrix_kron_mtimes_block()

%
% block_fatrix_kron_mtimes_block(): y = A{ib} * x
%
function y = block_fatrix_kron_mtimes_block(arg, is_transpose, x, iblock, nblock)

bl1 = arg.blocks{1}; % base block, already put through Gblock
warn 'todo: size check not done'

if ~is_transpose % forw
%	if nrow(x) ~= arg.dim(2)
%		fail('x size=%d vs dim(2)=%d', nrow(x), arg.dim(2))
%	end
	y = [];
	for ii=1:arg.Mkron
		t = bl1{iblock} * x([arg.jstart(ii):arg.jend(ii)], :);
		y = [y; t];
	end

else % back
	y = x;
%	if nrow(y) ~= arg.dim(1), error 'bad y size', end
	x = [];
	for ii=1:arg.Mkron
		tmp = [arg.istart(ii):arg.iend(ii)]; % i list (if all data)
		% todo: we need a certain subset of that list
		fail('todo: kron subset backprojector not done')
		t = bl1{iblock}' * y(tmp, :);
		x = [x; t];
	end
	y = x;
end

%
% block_fatrix_sum()
%
function ob = block_fatrix_sum(blocks, arg)

dim = size(blocks{1});
for ii=1:length(blocks)
	if ~isequal(size(blocks{ii}), dim)
		error 'blocks must have same size for "sum"'
	end
end

% build Fatrix object
ob = Fatrix(dim, arg, 'caller', 'block_fatrix(sum)', ...
	'forw', @block_fatrix_sum_forw, 'back', @block_fatrix_sum_back, ...
	'free', @block_fatrix_free);

%
% block_fatrix_sum_forw(): y = A * x
%
function y = block_fatrix_sum_forw(arg, x)

y = 0;
for ii=1:length(arg.blocks)
	y = y + arg.blocks{ii} * x;
end

%
% block_fatrix_sum_back(): x = A' * y
% full backprojection
%
function x = block_fatrix_sum_back(arg, y)

x = 0;
for ii=1:length(arg.blocks)
	x = x + arg.blocks{ii}' * y;
end

%
% block_fatrix_free()
%
function block_fatrix_free(arg)
if arg.chat
	printm 'freeing block_fatrix object static memory'
end
for ii=1:length(arg.blocks)
	try
		free(blocks{ii})
	catch
	end
end

%
% block_fatrix_update()
%
function out = block_fatrix_update(ob, varargin)
% todo: figure out how to update, e.g., new_zmap(s), ...