text:
function varargout = slinvoke(invoke_descr, varargin)
%SLINVOKE Invokes a function
%
% $ Syntax $
%   - [y1, y2, ..., yn] = slinvoke(invoke_descr, x1, x2, ..., xm)
%
% $ Arguments $
%   - invoke_descr:       the invoking descriptor
%   - x1, x2, ..., xm:    the input arguments
%   - y1, y2, ..., yn:    the output arguments
%
% $ Description $
%   - [y1, y2, ..., yn] = slinvoke(invoke_descr, placepos, x1, x2, ..., xm)
%     invokes some function with a specified manner. The invoking descriptor
%     can be in following form:
%     1. a string representing the function name, a function handle or an
%        inline object. In this case, no fixed argument is binded. The x1,
%        x2, ..., xm are directly fed to the function.
%     2. a cell array like {func, var_pos, fa1, fa2, ... fak}. Here
%        func is the function name, function handle or inline object that
%        is to be invoked. var_pos is the positions of variable arguments.
%        fa1, fa2, ... are the fixed(binded) arguments. If var_pos is [],
%        the variable arguments will be put first.
%
% $ Examples $
%   - Invoke a function without fixed arguments
%     \{
%           slinvoke('plot', x, y);
%           y = slinvoke({'sin', []}, x);
%     \}
%
%   - Invoke a function with fixed arguments binded, the following pairs
%     of statements are equivalent.
%     \{
%           K = slinvoke({'strfind', 1, 'pattern'}, s);
%           K = strfind(s, 'pattern');
%
%           slinvoke({'plot', [1 2], 'r-'}, x, y);
%           plot(x, y, 'r-');
%
%           slinvoke({'plot', 2, x, 'b+'}, y);
%           plot(x, y, 'b+');
%
%           b = slinvoke({'isequal', [], 0}, x);
%           b = isequal(x, 0);
%     \}
%
% $ History $
%   - Created by Dahua Lin on Dec 28th, 2005
%

%% parse input arguments
if ~iscell(invoke_descr)
    func = invoke_descr;
    var_pos =[];
    fix_args = {};
else
    len_invoke_descr = length(invoke_descr);
    func = invoke_descr{1};
    if len_invoke_descr >= 2
        var_pos = invoke_descr{2};
    else
        var_pos = [];
    end
    if len_invoke_descr >= 3
        fix_args = invoke_descr(3:end);
    else
        fix_args = {};
    end
end

if isempty(var_pos) && ~isempty(varargin)
    var_pos = 1:length(varargin);
end
n_fix_args = length(fix_args);
n_var_args = length(var_pos);
if length(varargin) ~= n_var_args
    error('sltoolbox:argmismatch', ...
        'The variable input arguments do not match that described in invoking descriptor');
end

%% organize input arguments
if n_fix_args == 0 && n_var_args == 0
    if nargout == 0
        feval(func);
    else
        [varargout{1:nargout}] = feval(func);
    end
elseif n_var_args == 0
    if nargout == 0
        feval(func, fix_args{:});
    else
        [varargout{1:nargout}] = feval(func, fix_args{:});
    end
else
    n_args = max(n_fix_args + n_var_args, max(var_pos));
    input_args = cell(1, n_args);
    indicator_var = false(1, n_args);
    indicator_var(var_pos) = true;
    fix_pos = find(~indicator_var);

    for i = 1 : n_var_args
        input_args{var_pos(i)} = varargin{i};
    end

    for i = 1 : n_fix_args
        input_args{fix_pos(i)} = fix_args{i};
    end

    if nargout == 0
        feval(func, input_args{:});
    else
        [varargout{1:nargout}] = feval(func, input_args{:});
    end
end

author: lmthang


text:
function [pvec, pstruct] = tapas_condhalluc_obs2_transp(r, ptrans)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2016 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

pvec    = NaN(1,length(ptrans));
pstruct = struct;

pvec(1)    = exp(ptrans(1));         % be
pstruct.be = pvec(1);

pvec(2)    = exp(ptrans(2));         % nu
pstruct.nu = pvec(2);

return;

author: translationalneuromodeling

repo: tapas


text:
filename = 'CantileverSquare';%'CantileverSquareSYmmetricMesh';%'CantileverSquareNewFine';%'LshapeTriFine';%'CantileverSquareNew';'ArchTriFine';%'CantileverSquare';%'ArchTriFine';%'CantileverSquare';%'Lshape';%'LshapeTriFine';%'CantileverSquareSmall';%'';%'Lshape';'CantileverSquare';'Lshape';%'CantileverSquare';%'LshapeFine';%'Bridge_quad_coarse';%'Arch_quad_coarse';%'BridgeCool_Quadrilateral_Bilinear_Structured_Coarse';%'Bridge';%'CantileverSquareSmall';
ptype = 'MACRO';
initial_case = 'given';
m1 = 0.0101;
m2 = 0.0101;
%cost = {'compliance'};
cost = {'stressNorm'};
%cost = {'stressNorm','compliance'};
%weights = [0.55,0.45];
weights = 1;
constraint = {'volumeConstraint'};
filterType = 'PDE';
constraint_case = 'EQUALITY';

Vfrac_initial = 0.3;
optimality_initial = 1e-5;
constr_initial = 1e-5;

Vfrac_final = 0.3;
optimality_final = 1e-5;
constr_final = 1e-5;

stressNormExponent_initial = 2;
stressNormExponent_final = 64;

optimizer = 'DualNestedInPrimal';
optimizerUnconstrained = 'PROJECTED GRADIENT';

designVariable = 'MicroParams';
ub = 0.989;
lb = 0.011;
homegenizedVariablesComputer = 'ByVademecum';
%
vademecumFileName = 'SuperEllipseQMax';
%vademecumFileName = 'SuperEllipseQ2';
%vademecumFileName = 'SuperEllipseQOptAnalytic';
%
% designVariable = 'Density';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';

kfrac = 2;
nsteps = 20;

plotting = false;
printing = true;
monitoring = true;
monitoring_interval = 1;
maxiter = 500;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function varargout = erfc(varargin)
%ERFC (overloaded)

switch class(varargin{1})

    case 'double'
        error('Overloaded SDPVAR/ERFC CALLED WITH DOUBLE. Report error')

    case 'sdpvar'
        varargout{1} = InstantiateElementWise(mfilename,varargin{:});

    case 'char'

        operator = struct('convexity','none','monotonicity','decreasing','definiteness','positive','model','callback');
        operator.bounds = @bounds;
        operator.range = [-1 1];
        operator.derivative =@(x)-exp(-x.^2)*2/sqrt(pi);
        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error('SDPVAR/ERF called with CHAR argument?');
end

function [L,U] = bounds(xL,xU)
L = erfc(xL);
U = erfc(xU);

author: shamilmamedov


text:
function [XS] = dataclean_XS(XS)

for i=1:size(XS,1)
  string = XS(i,3);
  string = string{1};
  if(strcmp(string,'Utilities ')==1),
    XS(i,3) = {'Utilities'};
  elseif(strcmp(string,'Services ')==1),
    XS(i,3) = {'Services'};
  elseif(strcmp(string,'#Technology')==1),
    XS(i,3) = {'Technology'};
  end
end

end

author: gudbrandtandberg


text:
function [ primal_obj ] = compute_primal(X, Y, W, Omega, lambda)
% Inputs
%   Xtrain: input training data (m-length cell)
%   Ytrain: output training data (m-length cell)
%   W: current models (d x m)
%   Omega: precision matrix (m x m)
%   lambda: regularization parameter
% Output
%   primal objective

% compute primal
total_loss = 0;
for t=1:length(X)
    preds = Y{t}.*(X{t}*W(:, t));
    total_loss = total_loss + mean(max(0.0, 1.0 - preds));
end
primal_obj = total_loss + lambda / 2 * trace(W * Omega * W');

end

author: gingsmith

repo: fmtl

save_path: github-repos/MATLAB/gingsmith-fmtl


text:
function [D,D1ols,q,C]=irfiv_ols_for_bootstrap_GK(txt,names,IVrotate,betadraw,n,Xdraw,Ydraw,k,p,enddate,startdate,cut1,cut2,cut3,cut4)
%% Copyright Ben Schumann
% function [D, gamma]=bear.irfiv_ols(names, betahat,sigmahat, n,X,Y,k,p,enddate,startdate)
% instrumental variable identification in an OLS setting
% inputs:  - matrix 'betahat': vec(OLS estimates of the reduced form)
%          - matrix 'sigmahat': vec(OLS estimates of sigma)
%          - matrix 'X': Independend Variable
%          - matrix 'Y': Dependend Variable
%          - matrix 'IV': Rotated IV from wild bootstrap
%          - integer 'IRFperiods': number of periods for IRFs
%          - integer 'n': number of endogenous variables in the VAR model (defined p 7 of technical guide)
%          - integer 'm': number of exogenous variables in the VAR model (defined p 7 of technical guide)
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - integer 'k': number of coefficients to estimate for each equation in the BVAR model (defined p 7 of technical guide)
%          - integer 'T': number of observations
%          - string  'stardate': VAR startdate
%          - string  'enddate': VAR enddate
% outputs: - matrix 'D': record of the accepted draws for the structural matrix D
%          - matrix 'gamma': record of the draws for the structural disturbances variance-covariance matrix gamma
IV = IVrotate;
%% Preparation for first stage regression
%get reduced form residuals
beta = betadraw;
B    = reshape(beta,k,n);
EPS  = Ydraw-Xdraw*B;

[EPSIV,IVcut] = bear.cut_EPS_IV_GK_new(txt, names, EPS, IV, cut1, cut2, cut3, cut4, startdate, enddate, p);

%% Imposing the covariance restrictions
%E_1 = EPSIV'*IVcut/length(IVcut);
%E11 = E_1(1,:);
%E21 = E_1(2:end,:);
%Mu = E21*E11^(-1); %relative impulse vector
%% normalize to a one standard deviation shock
sigmahatIV=(1/(length(EPSIV)-k))*(EPSIV'*EPSIV);

%get the gamma vector
%partition the reduced form VCV
%Sigma11 = sigmahatIV(1,1);
%Sigma12 = sigmahatIV(1,2:end);
%Sigma21 = sigmahatIV(2:end,1);
%Sigma22 = sigmahatIV(2:end,2:end);

%Gamma = Sigma22 + Mu*Sigma11*Mu' - Sigma21*Mu' - Mu*Sigma21'; %%%%%Gamma ouput is not used after this
%get b12 as in Michelle Piffers notes
%b12b12t = (Sigma21-Mu*Sigma11)'*Gamma^(-1)*(Sigma21-Mu*Sigma11);
%b11b11t = Sigma11 - b12b12t;
% b11 = chol(bear.nspd(b11b11t)); %%this is the scaling vector

%% first stage regression (this results in the same vector as Mu)

%step 2: Regress the first reduced form shock on the instrument
Shock = EPSIV(:,1);
[nobs , ~] = size(IVcut);
XX = [ones(nobs,1) IVcut];
[~, nvar] = size(XX);
%get OLS estimate
XpXi = (XX'*XX)\eye(nvar);
betaIV=XpXi*(XX'*Shock);
%get predicted value
IVpred = XX*betaIV;

%% second stage regression
%step 3: Regress the other reduced form shocks on the predicted value
ImpactIRFIV = zeros(n,1);
ImpactIRFIV(1,1) = 1;

for hh=2:n
Shock = EPSIV(:,hh);
[nobs,~]= size(IVpred);
IVpredtemp = [ones(nobs,1) IVpred];
[~,nvar] = size(IVpredtemp);
IVpIVi = (IVpredtemp'*IVpredtemp)\eye(nvar);
betaIV2=IVpIVi*(IVpredtemp'*Shock);
ImpactIRFIV(hh,1) = betaIV2(2,1); %should be equal to Mu from 2:end
end

% step 5: Create the structural matrix and only fill the first column as
% this is the only one identified
D=zeros(n,n);
% Step 6: Replace the first Column in the Cholesky Decomposition by
%the structural impact matrix computed above
%%another way to retrieve b11 (the scalar that scales the IRF to be a 1sdt Shock) is simply
C=chol(bear.nspd(sigmahatIV),'lower');
b=ImpactIRFIV;
%%Recover the vector q that maps the first column of C into b such that Cq=b;
q = C\b;
%%b11 is the euclidian length of q
b11q = 1/norm(q);

%
D(1:end,1) = ImpactIRFIV*b11q;
D1ols = ImpactIRFIV*b11q; % for the IRFt6 TakeOLS option

author: european-central-bank

repo: BEAR-toolbox


text:
    function [y1,y2,y3,y4] = organize(x1,x2,x3,x4,x5,x6)
                     dis1 = (x1-x3)^2+(x2-x4)^2;
                     dis2 = (x1-x5)^2+(x2-x6)^2;
                       if dis1 <= dis2
                          y1 = x3;
                          y2 = x4;
                          y3 = x5;
                          y4 = x6;
                      else
                          y1 = x5;
                          y2 = x6;
                          y3 = x3;
                          y4 = x4;
                      end;

author: Sable

repo: mcbench-benchmarks


text:
function img = imvoronoi2d(varargin)
%IMVORONOI2D Generate a 2D voronoi image from a set of points
%
%   IMG = imvoronoi2d(DIM, POINTS);
%   where DIM is a 1*2 array containing the size of the image, and POINTS
%   is a N*2 array of double containing germs coordinate, build the voronoi
%   image corresponding to given germs.
%
%   IMG = imvoronoi2d(LX, LY, POINTS);
%   Specifes the coordinates of the pixels with the two row vectors LX and
%   LY.
%
%   imvoronoi2d(DIM, POINTS, BORDER)
%   Also specifies edge condition. BORDER can be one of:
%   'free'      does not modify result image
%   'periodic'  consider germs are repeated in each direction
%   'remove'    remove cells touching borders of image
%
%   IMG = imvoronoi2d(DIM, N);
%   Creates an image of size DIM (1*2 array), using N germs with uniform
%   distribution in the image space.
%
%   imvoronoi2d
%   without argument return a 100*100 image with 10 random points as germs.
%
%
%   Example:
%   img = imvoronoi2d([100 100], [20 50; 80 30;70 90]);
%   imshow(img);
%
%   See also
%     imShapes, imVoronoi3d, imPowerDiagram, dilatedVoronoi
%

%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 28/04/2004.
%

%   HISTORY
%   use faster algorithm for distance function (loop on points and not on
%   pixels).
%   2008/10/10 update doc, code clean up
%   2009/03/04 switch coordinate to comply with (x,y) of germs
%   29/05/2009 use more possibilities for specifying grid

%% extract input arguments

% default values
points      = [];       % points array
N           = 20;       % default number of germs
edgecond    = 'free';   % edge condition
conn        = 4;        % 4-connectivity by default in dimension 2

% compute coordinate of image voxels
[lx, ly, varargin] = parseGridArgs(varargin{:});

% extraction of points, or number of points
if ~isempty(varargin)
    var = varargin{1};

    if length(var)>1    % use an array of germs
        points = var;
    else                % get number of germs
        N = var;
    end
end

% extraction of edge condition
if length(varargin)>1
    edgecond = varargin{2};
end

%% initialisations

% size in each direction
Ny = length(ly);
Nx = length(lx);

% create array of points, if it does not exist
if isempty(points)
    % generation of germs array
    points = zeros(N, 2);
    for i=1:N
        points(i,1) = rand*lx(end) + lx(1);
        points(i,2) = rand*ly(end) + ly(1);
    end
end

% create periodic conditions if needed:
% add the same point in each direction, including diagonals
if strcmp(edgecond, 'periodic')
    points = duplicateGerms(lx, ly, points);
end

% total number of points to consider
N = size(points, 1);

%% Main algorithm
% * first create distance function: each pixel get the distance to the
%   closest point
% * then perform watershed

% fill distance with a default value
dist = Nx * Ny * ones([Ny Nx]);

% update distance map with distance to closest point
for p = 1:N
    dx = repmat(reshape(lx-points(p, 1), [1 Nx]), [Ny 1]);
    dy = repmat(reshape(ly-points(p, 2), [Ny 1]), [1 Nx]);
    dist = min(dist, hypot(dx, dy));
end

% use a watershed to generate cells, with 4-connectivity
img = watershed(dist, conn) ~= 0;

%% post processing
% ---------------------------------------------

% if edge condition is 'remove', kill all cells touching borders
if strcmp(edgecond, 'remove')
    img = removeBorderRegions(img);
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function [frq,cr] = v_cent2frq(c)
%V_FRQ2ERB  Convert Hertz to Cents frequency scale [C,CR]=(FRQ)
%	frq = v_frq2mel(c) converts a vector of frequencies in cents
%	to the corresponding values in Hertz.
%   100 cents corresponds to one semitone and 440Hz corresponds to 5700
%   cents.
%   The optional cr output gives the gradient in Hz/cent.
%
%	The relationship between cents and frq is given by:
%
%	c = 1200 * log2(f/(440*(2^((3/12)-5)))
%
%	Reference:
%
%     [1] Ellis, A.
%         On the Musical Scales of Various Nations
%         Journal of the Society of Arts, 1885, 485-527

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: v_cent2frq.m 10865 2018-09-21 17:22:45Z 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent p q
if isempty(p)
    p=1200/log(2);
    q=5700-p*log(440);
end
% c = 1200*sign(frq).*log2(frq/(440*2^((3/12)-5)));
af=(exp((abs(c)-q)/p));
frq=sign(c).*af;
cr=af/p;
if ~nargout
    plot(c,frq,'-x');
    ylabel(['Frequency (' v_yticksi 'Hz)']);
    xlabel(['Frequency (' v_xticksi 'Cents)']);
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function value = gpl_dimension ( filename )

%*****************************************************************************80
%
%% GPL_DIMENSION determines the "dimension" of the data in a GPL file.
%
%  Discussion:
%
%    A satisfactory GPL file contains only three kinds of records:
%    * comment lines begin with a "#" and are ignored;
%    * blank lines;
%    * lines with exactly N numeric values;
%
%    This function attempts to determine the value of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILENAME, the name of the file.
%
%    Output, integer VALUE, the "dimension" of the data.
%    -2, inconsistent data was found.
%    -1, no data was found at all.
%     0, only empty lines were read.
%     n, every nonempy line contains N values.
%
  file = fopen ( filename, 'rt' );

  d = -1;

  while ( 1 )

    s = fgetl ( file );

    if ( ~ ischar ( s ) )
      break
    end

    if ( length ( s ) < 1 )
      continue
    end

    if ( s(1) == '#' )
      continue;
    end

    d2 = s_word_count ( s );

    if ( d == -1 )
      d = d2;
    elseif ( d == d2 )

    elseif ( d == 0 )
      d = d2;
    else
      d = -2;
      break
    end

  end

  fclose ( file );

  value = d;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function trans = createTranslation3d(varargin)
%CREATETRANSLATION3D Create the 4x4 matrix of a 3D translation.
%
%   usage:
%   TRANS = createTranslation3d(DX, DY, DZ);
%   return the translation corresponding to DX and DY.
%   The returned matrix has the form :
%   [1 0 0 DX]
%   [0 1 0 DY]
%   [0 0 1 DZ]
%   [0 0 0  1]
%
%   TRANS = createTranslation3d(VECT);
%   return the translation corresponding to the given vector [x y z].
%
%
%   See also:
%   transforms3d, transformPoint3d, transformVector3d,
%   createRotationOx, createRotationOy, createRotationOz, createScaling3d
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 06/04/2004.
%

%   HISTORY
%   22/04/2009 rename as createTranslation3d

if isempty(varargin)
    % assert translation with null vector
    dx = 0;
    dy = 0;
    dz = 0;
elseif length(varargin)==1
    % translation vector given in a single argument
    var = varargin{1};
    dx = var(1);
    dy = var(2);
    dz = var(3);
else
    % translation vector given in 3 arguments
    dx = varargin{1};
    dy = varargin{2};
    dz = varargin{3};
end

% create the translation matrix
trans = [1 0 0 dx ; 0 1 0 dy ; 0 0 1 dz; 0 0 0 1];

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [hp] = kW2hp(kW)
% Convert power from kilowatts to mechanical horsepower.
% Chad A. Greene 2012
hp = kW*1.34102209;


text:
function varargout=ceil(varargin)
%FLOOR (overloaded)

switch class(varargin{1})

    case 'sdpvar' % Overloaded operator for SDPVAR objects. Pass on args and save them.

        x = varargin{1};
        [n,m] = size(x);
        x = reshape(x,n*m,1);
        y = [];
        for i = 1:n*m
            y = [y;yalmip('addextendedvariable',mfilename,extsubsref(x,i))];
        end
        y = reshape(y,n,m);
        varargout{1} = y;

    case 'char' % YALMIP send 'graph' when it wants the epigraph or hypograph
        switch varargin{1}
            case {'milp','graph'}
                % Description using epigraphs
                t = varargin{2};
                X = varargin{3};

                c = intvar(1,1);
                F = (x-1 <= c <= x);

                varargout{1} = F;
                varargout{2} = struct('convexity','milp','monotonicity','milp','definiteness','milp');
                varargout{3} = X;

            otherwise
                error('SDPVAR/SORT called with CHAR argument?');
        end
    otherwise
        error('Strange type on first argument in SDPVAR/SORT');
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function FPtable = loadFPtable(fn)
% load the lookup table from CCF to FP slices

[~, fnBase] = fileparts(fn);

fid = fopen(fn, 'r');

titles = textscan(fid, '%s%s%s%s%s%s%s%s%s%s%s', 1, 'delimiter', ',');
titles = cellfun(@(x)x{1}, titles, 'uni', false);
data = textscan(fid, '%s%s%d%d%d%d%d%d%s%s%s', 'delimiter', ',');

FPtable = table(data{:}, 'VariableNames', titles);

fclose(fid);

author: cortex-lab

repo: allenCCF

save_path: github-repos/MATLAB/cortex-lab-allenCCF


text:
% DEMHSVARGPLVMREGRESSION A script to run deep GP regression.
%
% This is a generic demo. You can replace the data used here with your own
% data and run it (ie Ytr{1} has to be the observed data and inpX has to be
% your observed labels).
%
% To configure the deepGP used here, do the following:
%   1. If any of the fields you want to change appear in this demo in the
%   "configuring the deep gp" section, they change it directly there.
%   2. If the field you're looking for is not there, then check the available
%    configuration options in hsvargplvm_init.m. The way this works, is that
%    you just need to overwrite the corresponding workspace variables.
%    e.g. if in hsvargplvm_init.m you see a field "fixInducing", you can
%    overwrite this field by just calling this demo as:
%    >> fixInducing = true; demHsvargplvmClassification
%   3. If the field you're looking for is not in hsvargplvm_init.m, then
%   check (in this order) svargplvm_init.m and vargplvm_init.m, again
%   overwritting the configuration by specifying the variable to exist, as
%   above.
%
% SEEALSO: demHsvargplvmClassification.m
%
% COPYRIGHT: Andreas Damianou, 2014
% DEEPGP

%% ------ CONFIGURING THE DEEP GP
%--- Mandatory configurations
if ~exist('Ytr', 'var'), error('You need to specify your outputs in Ytr{1}=...'); end
if ~exist('inpX', 'var'), error('You need to specify your inputs in inpX=...'); end

%--- Optional configurations: Whatever configuration variable is not already set (ie does not exist
% as a variable in the workspace) is set to a default value.
if ~exist('experimentNo'), experimentNo = 404; end
if ~exist('K'), K = 30; end
if ~exist('Q'), Q = 6; end
if ~exist('baseKern'), baseKern = 'rbfardjit'; end % {'rbfard2','white','bias'}; end
% This is called "dynamics" and "time" for historical reasons.. It actually refers to a coupling GP in the uppermost level
if ~exist('dynamicsConstrainType'), dynamicsConstrainType = {'time'}; end
stackedOpt = [];
if exist('stackedInitVardistIters', 'var'), stackedOpt.stackedInitVardistIters=stackedInitVardistIters; end
if exist('stackedInitIters', 'var'), stackedOpt.stackedInitIters=stackedInitIters; end
if exist('stackedInitSNR', 'var'), stackedOpt.stackedInitSNR=stackedInitSNR; end
if exist('stackedInitK', 'var'), stackedOpt.stackedInitK=stackedInitK; end
if ~exist('initXOptions', 'var'), initXOptions = []; end

% Initialise script based on the above variables. This returns a struct
% "globalOpt" which contains all configuration options
hsvargplvm_init;

% Automatically calibrate initial variational covariances - better to not change that
globalOpt.vardistCovarsMult = [];

[options, optionsDyn] = hsvargplvmOptions(globalOpt, inpX);

%% ------------- Initialisation and model creation
% Initialise latent spaces, unless the user already did that
if ~(iscell(options.initX) && prod(size(options.initX{1})) > 1)
    [globalOpt, options, optionsDyn, initXOptions] = hsvargplvmRegressionInitX(globalOpt, options, optionsDyn, inpX, Ytr, stackedOpt);
end

% Create the deep GP based on the model options, global options
% (configuration) and options for initialising the latent spaces X
model = hsvargplvmModelCreate(Ytr, options, globalOpt, initXOptions);

% Since we do regression, we need to add a GP on the parent node. This GP
% couples the inputs and is parametrised by options in a struct "optionsDyn".
model = hsvargplvmAddParentPrior(model, globalOpt, optionsDyn);

%-- We have the option to not learn the inducing points and/or fix them to
% the given inputs.
% Learn inducing points? (that's different to fixInducing, ie tie them
% to X's, if learnInducing is false they will stay in their original
% values, ie they won't constitute parameters of the model).
if exist('learnInducing') && ~learnInducing
    model = hsvargplvmPropagateField(model, 'learnInducing', false);
end
%--

if globalOpt.fixInducing && globalOpt.fixInducing
    model = hsvargplvmPropagateField(model, 'fixInducing', true);
    for m=1:model.layer{end}.M % Not implemented yet for parent node
        model.layer{end}.comp{m}.fixInducing = false;
    end
end

%!!!!!!!!!!!!!!!!!!!!!!!!-----------------------
if exist('DEBUG_entropy','var') && DEBUG_entropy
    model.DEBUG_entropy = true;for itmp=1:model.H, model.layer{itmp}.DEBUG_entropy = true; end
end

params = hsvargplvmExtractParam(model);
model = hsvargplvmExpandParam(model, params);
model.globalOpt = globalOpt;
% Computations can be made in parallel, if option is activated
model.parallel = globalOpt.enableParallelism;

fprintf('# Scales after init. latent space:\n')
hsvargplvmShowScales(model,false);
%% OPTIMISATION
[model,modelPruned, modelInitVardist] = hsvargplvmOptimiseModel(model, true, true);

% If you decide to train for further iterations...
% modelOld = model; [model,modelPruned, ~] = hsvargplvmOptimiseModel(model, true, true, [], {0, [100]});

author: SheffieldML

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


text:
function drawSurfPatch(u, v, z, varargin)
%DRAWSURFPATCH Draw a 3D surface patch, with 2 parametrized surfaces.
%
%   usage:
%   drawSurfPatch(u, v, zuv)
%   where u, v, and zuv are three matrices the same size, u and
%   corresponding to each parameter, and zuv being equal to a function of u
%   and v.
%
%   drawSurfPatch(u, v, zuv, p0)
%   If p0 is specified, two lines with u(p0(1)) and v(p0(2)) are drawn on
%   the surface, and corresponding tangent are also shown.
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 24/05/2005.
%

%   HISTORY
%   2005-06-08 add doc.
%   2007-01-04 remove unused variables and change function name
%   2010-03-08 code cleanup, use drawPolyline3d

% prepare figure
hold on;

% draw the surface interior
surf(u, v, z, 'FaceColor', 'g', 'EdgeColor', 'none');

% draw the surface boundaries
drawPolyline3d(u(1,:), v(1,:), z(1,:))
drawPolyline3d(u(end,:), v(end,:), z(end,:))
drawPolyline3d(u(:,end), v(:,end), z(:,end))
drawPolyline3d(u(:,1), v(:,1), z(:,1))

% eventually draw two perpendicular lines on the surface
if ~isempty(varargin)
    pos = varargin{1};
    drawPolyline3d(u(pos(1),:), v(pos(1),:), z(pos(1),:));
    drawPolyline3d(u(:,pos(2)), v(:,pos(2)), z(:,pos(2)));
end


text:
function theta = GFM_rand_init( tree, M, min_prob, min_std )

% Initial values for GFM EM algorithm
% Use the 'small' EM algorithm of Fan & Xia to initialize the 'big' EM
% algorithm of Crouse et al/Durand et al.
%
% Syntax:
%   theta = GFM_rand_init( tree, M, min_prob, min_std )
%
% Input:
%   tree  : Wavelet transform from DWT2_TO_TREE with L levels and D
%           dimensions
%
%   M     : The number of mixtures.
%
%	min_p : Mininimum value of the probabilities
%
%	min_s : Mininimum value of the standard devs
%
%
% Output:
%   theta : L-by-3-by-D cell with parameters. For the third dimension row l
%           are the parameters for level l:
%           state probs, transition probs, standard deviations
%
%
% See also: GFM_INIT, GFM_EM_WRAPPER, GFM_EM

L = length(tree) - 1;
D = length(tree{2});

% Initialize output
theta = cell( L, 3, D );

for l = 1:L
    for d = 1:D
		% State probs
        p = max( min_prob, rand(1,2) );
        theta{l, 1, d} = p / sum(p);

		if l > 1
			% Trans probs
    		p = max( min_prob, rand(2) );
    		theta{l, 2, d} = p ./ repmat( sum(p), [M 1] );

			% Update state probs to be consistent with the trans probs
			theta{l,1,d} = theta{l-1,1,d} * theta{l,2,d}';
		end

		% Standard devs
        theta{l, 3, d} = max( min_std, rand(1,M)*std(tree{l+1}{d}(:)) );
    end
end

end


text:
function [p,q] = bandwidth(A)
%BANDWIDTH    Upper and lower bandwidth of matrix A
%
%   A   = 0  for  i-j > p  or j-i > q
%    ij
%
%    [p,q] = bandwidth(A)
%

% written  10/16/98     S.M. Rump
% modified 09/21/02     S.M. Rump  one output argument corrected
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  if A.complex
    [p1,q1] = bandwidth(A.mid);
    [p2,q2] = bandwidth(A.rad);
  else
    [p1,q1] = bandwidth(A.inf);
    [p2,q2] = bandwidth(A.sup);
  end

  p = max(p1,p2);
  q = max(q1,q2);

  if nargout<=1
    p = max(abs(p),abs(q));
  end

repo: OpenMAS


text:
% THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION IS RELEASED "AS IS."  THE U.S. GOVERNMENT MAKES NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, CONCERNING THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION, INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT WILL THE U.S. GOVERNMENT BE LIABLE FOR ANY DAMAGES, INCLUDING ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE, OR INABILITY TO USE, THIS SOFTWARE OR ANY ACCOMPANYING DOCUMENTATION, EVEN IF INFORMED IN ADVANCE OF THE POSSIBILITY OF SUCH DAMAGES.
%
% change log:
%   040204 smw stub file
%
% file: get_r_ci_na_1s.m
%
function varargout = get_r_ci_na_1s(varargin)
  error('get_r_ci_na_1s is not part of current release.');
return;


text:
function [qp]=zn3tr(input)
% [qp]=zn3tr(input)
% Ziegler-Nichols PID controller for processes of 3rd order.
% This function computes parameters of the controller (q0, q1, q2, p1, p2).
% Controller is based on trapezoidal method of discretization.
% Transfer function of the controller is as follows:
%
%            q0 + q1*z^-1 + q2*z^-2     q0 + q1*z^-1 + q2*z^-2
% G(z^-1) = ------------------------ = ------------------------
%                  1 - z^-1              1 + p1*z^-1 + p2*z^-2
%
% where p1=-1 and p2=0.
%
% Transfer function of the controlled system is:
%
%               b1*z^-1 + b2*z^-2 + b3*z^-3
% Gs(z^-1) = ---------------------------------
%             1 + a1*z^-1 + a2*z^-2 + a3*z^-3
%
% Input: input ... input parameters
%                  input(1) ... a1
%                  input(2) ... b1
%                  input(3) ... a2
%                  input(4) ... b2
%                  input(5) ... a3
%                  input(6) ... b3
%                  input(7) ... sample time T0
% Output: qp ... controller parameters
%                qp(1) ... q0
%                qp(2) ... q1
%                qp(3) ... q2
%                qp(4) ... p1 (-1)
%                qp(5) ... p2 (0)

a1 = input(1);
b1 = input(2);
a2 = input(3);
b2 = input(4);
a3 = input(5);
b3 = input(6);
T0 = input(7);

% compute ultimate gain and frequency
[Kpu, Tu] =  ultim([b1 b2 b3],[a1 a2 a3],T0);

Kp = 0.6*Kpu;
Ti = Tu/2;
Td = Tu/8;

q0 = Kp*(1 + T0/(2*Ti) + Td/T0);
q1 = -Kp*(1 - T0/(2*Ti) + 2*Td/T0);
q2 = Kp*(Td/T0);
p1 = -1;
p2 = 0;

qp=[q0; q1; q2; p1; p2];


text:
function [ FIFF ] = fiff_define_constants()
%
%    [ FIFF ] = fiff_define_constants()
%
%    Defines structure containing the FIFF constants
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%    Revision 1.23  2009/03/30 11:37:37  msh
%    Added copying of measurement info blocks from the original like in mne_browse_raw
%
%    Revision 1.22  2009/02/16 04:38:21  msh
%    Fixed comments on coordinate transforms
%
%
%    Revision 1.20  2008/11/16 21:31:23  msh
%    Added mne_transform_coordinates and new coordinate frame definitions
%
%    Revision 1.19  2008/09/30 17:22:39  msh
%    Added possibility to load unprocessed MaxShield data
%
%    Revision 1.18  2008/06/07 21:22:10  msh
%    Added FIFF_FORWARD_SOLUTION_GRAD constant
%
%    Revision 1.17  2008/05/26 10:49:26  msh
%    Update to incorporate the already weighted lead field basis
%
%    Revision 1.16  2008/05/06 20:40:56  msh
%    Fixed ordering of output for compatibility with maxfilter averager
%
%    Revision 1.15  2008/04/16 22:24:57  msh
%    Added megacq parameters to the measurement info
%
%    Revision 1.14  2008/04/16 17:29:28  msh
%    Added acquisition parameter tags.
%
%    Revision 1.13  2008/03/13 19:18:06  msh
%    Read and write FIFF_MEAS_DATE from/to FIFFB_MEAS_INFO as appropriate
%
%    Revision 1.12  2007/01/29 21:21:21  msh
%    Added reading of the additional source prior covariances.
%
%    Revision 1.11  2006/09/25 19:48:16  msh
%    Added projection item kinds to fiff_define_constants
%    Changed some fields to int32 in surface structures
%
%    Revision 1.10  2006/09/24 18:52:43  msh
%    Added FIFFV_REF_MEG_CH to fiff_define_constants.
%    Added coord_frame field to dig point structure.
%
%    Revision 1.9  2006/06/22 21:22:46  msh
%    Take into account the possibility of calibrated compensation matrices
%
%    Revision 1.8  2006/05/05 03:50:40  msh
%    Added routines to compute L2-norm inverse solutions.
%    Added mne_write_inverse_sol_stc to write them in stc files
%    Several bug fixes in other files
%
%    Revision 1.7  2006/04/28 18:01:27  msh
%    Added mne_read_bem_surfaces
%    Fixed errors in mne_read_source_spaces and improved triangle data computation.
%
%    Revision 1.6  2006/04/26 19:50:58  msh
%    Added fiff_read_mri
%
%    Revision 1.5  2006/04/23 15:29:39  msh
%    Added MGH to the copyright
%
%    Revision 1.4  2006/04/15 12:21:00  msh
%    Several small improvements
%
%    Revision 1.3  2006/04/10 23:26:53  msh
%    Added fiff reading routines
%
%    Revision 1.2  2005/12/05 20:23:21  msh
%    Added fiff_save_evoked. Improved error handling.
%
%    Revision 1.1  2005/12/05 16:01:04  msh
%    Added an initial set of fiff writing routines.
%
%

%
% Blocks
%
FIFF.FIFFB_MEAS               = 100;
FIFF.FIFFB_MEAS_INFO          = 101;
FIFF.FIFFB_RAW_DATA           = 102;
FIFF.FIFFB_PROCESSED_DATA     = 103;
FIFF.FIFFB_CONTINUOUS_DATA    = 112;
FIFF.FIFFB_EVOKED             = 104;
FIFF.FIFFB_ASPECT             = 105;
FIFF.FIFFB_SUBJECT            = 106;
FIFF.FIFFB_ISOTRAK            = 107;
FIFF.FIFFB_HPI_MEAS           = 108;
FIFF.FIFFB_HPI_RESULT         = 109;
FIFF.FIFFB_DACQ_PARS          = 117;
FIFF.FIFFB_REF                = 118;
FIFF.FIFFB_SMSH_RAW_DATA      = 119;
FIFF.FIFFB_SMSH_ASPECT        = 120;
FIFF.FIFFB_HPI_SUBSYSTEM      = 121;

FIFF.FIFFB_PROJ               = 313;
FIFF.FIFFB_PROJ_ITEM          = 314;
FIFF.FIFFB_MRI                = 200;
FIFF.FIFFB_MRI_SET            = 201;
FIFF.FIFFB_MRI_SLICE          = 202;
FIFF.FIFFB_PROCESSING_HISTORY = 900;
FIFF.FIFFB_SSS_INFO           = 502;
FIFF.FIFFB_SSS_CAL_ADJUST     = 503;
FIFF.FIFFB_SSS_ST_INFO        = 504;
FIFF.FIFFB_SSS_BASES          = 505;
%
% Of general interest
%
FIFF.FIFF_FILE_ID         = 100;
FIFF.FIFF_DIR_POINTER     = 101;
FIFF.FIFF_BLOCK_ID        = 103;
FIFF.FIFF_BLOCK_START     = 104;
FIFF.FIFF_BLOCK_END       = 105;
FIFF.FIFF_FREE_LIST       = 106;
FIFF.FIFF_FREE_BLOCK      = 107;
FIFF.FIFF_NOP             = 108;
FIFF.FIFF_PARENT_FILE_ID  = 109;
FIFF.FIFF_PARENT_BLOCK_ID = 110;
%
%  Megacq saves the parameters in these tags
%
FIFF.FIFF_DACQ_PARS      = 150;
FIFF.FIFF_DACQ_STIM      = 151;

FIFF.FIFF_SFREQ       = 201;
FIFF.FIFF_NCHAN       = 200;
FIFF.FIFF_DATA_PACK   = 202;
FIFF.FIFF_CH_INFO     = 203;
FIFF.FIFF_MEAS_DATE   = 204;
FIFF.FIFF_SUBJECT     = 205;
FIFF.FIFF_COMMENT     = 206;
FIFF.FIFF_NAVE        = 207;
FIFF.FIFF_DIG_POINT   = 213;
FIFF.FIFF_LOWPASS     = 219;
FIFF.FIFF_COORD_TRANS = 222;
FIFF.FIFF_HIGHPASS    = 223;
FIFF.FIFF_NAME        = 233;
FIFF.FIFF_DESCRIPTION = FIFF.FIFF_COMMENT;
%
% HPI result
%
FIFF.FIFF_HPI_COIL_MOMENTS       = 240;
FIFF.FIFF_HPI_FIT_GOODNESS       = 241;
FIFF.FIFF_HPI_FIT_ACCEPT         = 242;
FIFF.FIFF_HPI_FIT_GOOD_LIMIT     = 243;
FIFF.FIFF_HPI_FIT_DIST_LIMIT     = 244;
FIFF.FIFF_HPI_COIL_NO            = 245;
FIFF.FIFF_HPI_COILS_USED         = 246;
FIFF.FIFF_HPI_DIGITIZATION_ORDER = 247;
%
% Pointers
%
FIFF.FIFFV_NEXT_SEQ    = 0;
FIFF.FIFFV_NEXT_NONE   = -1;
%
% Channel types
%
FIFF.FIFFV_MEG_CH     =   1;
FIFF.FIFFV_REF_MEG_CH = 301;
FIFF.FIFFV_EEG_CH     =   2;
FIFF.FIFFV_MCG_CH     = 201;
FIFF.FIFFV_STIM_CH    =   3;
FIFF.FIFFV_EOG_CH     = 202;
FIFF.FIFFV_EMG_CH     = 302;
FIFF.FIFFV_ECG_CH     = 402;
FIFF.FIFFV_MISC_CH    = 502;
FIFF.FIFFV_RESP_CH    = 602;             % Respiration monitoring
%
% Quaternion channels for head position monitoring
%
FIFF.FIFFV_QUAT_0   = 700;                 % Quaternion parameter q0; obsolete for unit quaternion
FIFF.FIFFV_QUAT_1   = 701;                 % Quaternion parameter q1; rotation
FIFF.FIFFV_QUAT_2   = 702;                 % Quaternion parameter q2; rotation
FIFF.FIFFV_QUAT_3   = 703;                 % Quaternion parameter q3; rotation
FIFF.FIFFV_QUAT_4   = 704;                 % Quaternion parameter q4; translation
FIFF.FIFFV_QUAT_5   = 705;                 % Quaternion parameter q5; translation
FIFF.FIFFV_QUAT_6   = 706;                 % Quaternion parameter q6; translation
FIFF.FIFFV_HPI_G    = 707;                 % Goodness-of-fit in continuous hpi
FIFF.FIFFV_HPI_ERR  = 708;                 % Estimation error in continuous hpi
FIFF.FIFFV_HPI_MOV  = 709;                 % Estimated head movement speed in continuous hpi
%
% Coordinate frames
%
FIFF.FIFFV_COORD_UNKNOWN        = 0;
FIFF.FIFFV_COORD_DEVICE         = 1;
FIFF.FIFFV_COORD_ISOTRAK        = 2;
FIFF.FIFFV_COORD_HPI            = 3;
FIFF.FIFFV_COORD_HEAD           = 4;
FIFF.FIFFV_COORD_MRI            = 5;
FIFF.FIFFV_COORD_MRI_SLICE      = 6;
FIFF.FIFFV_COORD_MRI_DISPLAY    = 7;
FIFF.FIFFV_COORD_DICOM_DEVICE   = 8;
FIFF.FIFFV_COORD_IMAGING_DEVICE = 9;
%
% Needed for raw and evoked-response data
%
FIFF.FIFF_FIRST_SAMPLE       = 208;
FIFF.FIFF_LAST_SAMPLE        = 209;
FIFF.FIFF_ASPECT_KIND        = 210;
FIFF.FIFF_DATA_BUFFER        = 300;      % Buffer containing measurement data
FIFF.FIFF_DATA_SKIP          = 301;      % Data skip in buffers
FIFF.FIFF_EPOCH              = 302;      % Buffer containing one epoch and channel
FIFF.FIFF_DATA_SKIP_SAMP     = 303;      % Data skip in samples
FIFF.FIFF_MNE_BASELINE_MIN   = 304;      % Time of baseline beginning
FIFF.FIFF_MNE_BASELINE_MAX   = 305;      % Time of baseline end
%
% Different aspects of data
%
FIFF.FIFFV_ASPECT_AVERAGE       = 100;      % Normal average of epochs
FIFF.FIFFV_ASPECT_STD_ERR       = 101;      % Std. error of mean
FIFF.FIFFV_ASPECT_SINGLE        = 102;      % Single epoch cut out from the continuous data
FIFF.FIFFV_ASPECT_SUBAVERAGE    = 103;
FIFF.FIFFV_ASPECT_ALTAVERAGE    = 104;      % Alternating subaverage
FIFF.FIFFV_ASPECT_SAMPLE        = 105;      % A sample cut out by graph
FIFF.FIFFV_ASPECT_POWER_DENSITY = 106;    % Power density spectrum
FIFF.FIFFV_ASPECT_DIPOLE_WAVE   = 200;    % Dipole amplitude curve
%
% BEM surface IDs
%
FIFF.FIFFV_BEM_SURF_ID_UNKNOWN    = -1;
FIFF.FIFFV_BEM_SURF_ID_BRAIN      = 1;
FIFF.FIFFV_BEM_SURF_ID_SKULL      = 3;
FIFF.FIFFV_BEM_SURF_ID_HEAD       = 4;
%
% More of those defined in MNE
%
FIFF.FIFFV_MNE_SURF_UNKNOWN       = -1;
FIFF.FIFFV_MNE_SURF_LEFT_HEMI     = 101;
FIFF.FIFFV_MNE_SURF_RIGHT_HEMI    = 102;
%
%   These relate to the Isotrak data
%
FIFF.FIFFV_POINT_CARDINAL = 1;
FIFF.FIFFV_POINT_HPI      = 2;
FIFF.FIFFV_POINT_EEG      = 3;
FIFF.FIFFV_POINT_ECG      = FIFF.FIFFV_POINT_EEG;
FIFF.FIFFV_POINT_EXTRA    = 4;

FIFF.FIFFV_POINT_LPA = 1;
FIFF.FIFFV_POINT_NASION = 2;
FIFF.FIFFV_POINT_RPA = 3;
%
%   SSP
%
FIFF.FIFF_PROJ_ITEM_KIND         = 3411;
FIFF.FIFF_PROJ_ITEM_TIME         = 3412;
FIFF.FIFF_PROJ_ITEM_NVEC         = 3414;
FIFF.FIFF_PROJ_ITEM_VECTORS      = 3415;
FIFF.FIFF_PROJ_ITEM_CH_NAME_LIST = 3417;
%
%   MRIs
%
FIFF.FIFF_MRI_SOURCE_PATH       = 1101;
FIFF.FIFF_MRI_SOURCE_FORMAT     = 2002;
FIFF.FIFF_MRI_PIXEL_ENCODING    = 2003;
FIFF.FIFF_MRI_PIXEL_DATA_OFFSET = 2004;
FIFF.FIFF_MRI_PIXEL_SCALE       = 2005;
FIFF.FIFF_MRI_PIXEL_DATA        = 2006;
FIFF.FIFF_MRI_WIDTH             = 2010;
FIFF.FIFF_MRI_WIDTH_M           = 2011;
FIFF.FIFF_MRI_HEIGHT            = 2012;
FIFF.FIFF_MRI_HEIGHT_M          = 2013;
%
FIFF.FIFFV_MRI_PIXEL_BYTE       = 1;
FIFF.FIFFV_MRI_PIXEL_WORD       = 2;
FIFF.FIFFV_MRI_PIXEL_SWAP_WORD  = 3;
FIFF.FIFFV_MRI_PIXEL_FLOAT      = 4;
%
%   These are the MNE fiff definitions
%
FIFF.FIFFB_MNE                    = 350;
FIFF.FIFFB_MNE_SOURCE_SPACE       = 351;
FIFF.FIFFB_MNE_FORWARD_SOLUTION   = 352;
FIFF.FIFFB_MNE_PARENT_MRI_FILE    = 353;
FIFF.FIFFB_MNE_PARENT_MEAS_FILE   = 354;
FIFF.FIFFB_MNE_COV                = 355;
FIFF.FIFFB_MNE_INVERSE_SOLUTION   = 356;
FIFF.FIFFB_MNE_NAMED_MATRIX       = 357;
FIFF.FIFFB_MNE_ENV                = 358;
FIFF.FIFFB_MNE_BAD_CHANNELS       = 359;
FIFF.FIFFB_MNE_VERTEX_MAP         = 360;
FIFF.FIFFB_MNE_EVENTS             = 361;
FIFF.FIFFB_MNE_MORPH_MAP          = 362;
%
% CTF compensation data
%
FIFF.FIFFB_MNE_CTF_COMP           = 370;
FIFF.FIFFB_MNE_CTF_COMP_DATA      = 371;
FIFF.FIFFB_MNE_DERIVATIONS        = 372;

FIFF.FIFFB_MNE_EPOCHS             = 373;
FIFF.FIFFB_MNE_ICA                = 374;
%
% Fiff tags associated with MNE computations (3500...)
%
%
% 3500... Bookkeeping
%
FIFF.FIFF_MNE_ROW_NAMES              = 3502;
FIFF.FIFF_MNE_COL_NAMES              = 3503;
FIFF.FIFF_MNE_NROW                   = 3504;
FIFF.FIFF_MNE_NCOL                   = 3505;
FIFF.FIFF_MNE_COORD_FRAME            = 3506;      % Coordinate frame employed. Defaults:
%  FIFFB_MNE_SOURCE_SPACE       FIFFV_COORD_MRI
%  FIFFB_MNE_FORWARD_SOLUTION   FIFFV_COORD_HEAD
%  FIFFB_MNE_INVERSE_SOLUTION   FIFFV_COORD_HEAD
FIFF.FIFF_MNE_CH_NAME_LIST           = 3507;
FIFF.FIFF_MNE_FILE_NAME              = 3508;      % This removes the collision with fiff_file.h (used to be 3501)
%
% 3510... 3590... Source space or surface
%
FIFF.FIFF_MNE_SOURCE_SPACE_POINTS        = 3510;    % The vertices
FIFF.FIFF_MNE_SOURCE_SPACE_NORMALS       = 3511;    % The vertex normals
FIFF.FIFF_MNE_SOURCE_SPACE_NPOINTS       = 3512;    % How many vertices
FIFF.FIFF_MNE_SOURCE_SPACE_SELECTION     = 3513;    % Which are selected to the source space
FIFF.FIFF_MNE_SOURCE_SPACE_NUSE          = 3514;    % How many are in use
FIFF.FIFF_MNE_SOURCE_SPACE_NEAREST       = 3515;    % Nearest source space vertex for all vertices
FIFF.FIFF_MNE_SOURCE_SPACE_NEAREST_DIST  = 3516;    % Distance to the Nearest source space vertex for all vertices
FIFF.FIFF_MNE_SOURCE_SPACE_ID            = 3517;    % Identifier
FIFF.FIFF_MNE_SOURCE_SPACE_TYPE          = 3518;    % Surface or volume

FIFF.FIFF_MNE_SOURCE_SPACE_NTRI          = 3590;    % Number of triangles
FIFF.FIFF_MNE_SOURCE_SPACE_TRIANGLES     = 3591;    % The triangulation
FIFF.FIFF_MNE_SOURCE_SPACE_NUSE_TRI      = 3592;    % Number of triangles corresponding to the number of vertices in use
FIFF.FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES = 3593;    % The triangulation of the used vertices in the source space
FIFF.FIFF_MNE_SOURCE_SPACE_NNEIGHBORS    = 3594;     % Number of neighbors for each source space point (used for volume source spaces)
FIFF.FIFF_MNE_SOURCE_SPACE_NEIGHBORS     = 3595;    % Neighbors for each source space point (used for volume source spaces)
FIFF.FIFF_MNE_SOURCE_SPACE_VOXEL_DIMS    = 3596;    % Voxel space dimensions in a volume source space
FIFF.FIFF_MNE_SOURCE_SPACE_INTERPOLATOR  = 3597;    % Matrix to interpolate a volume source space into a mri volume
FIFF.FIFF_MNE_SOURCE_SPACE_MRI_FILE      = 3598;    % MRI file used in the interpolation

FIFF.FIFF_MNE_SOURCE_SPACE_DIST          = 3599;    % Distances between vertices in use (along the surface)
FIFF.FIFF_MNE_SOURCE_SPACE_DIST_LIMIT    = 3600;    % If distance is above this limit (in the volume) it has not been calculated
%
% 3520... Forward solution
%
FIFF.FIFF_MNE_FORWARD_SOLUTION       = 3520;
FIFF.FIFF_MNE_SOURCE_ORIENTATION     = 3521;    % Fixed or free
FIFF.FIFF_MNE_INCLUDED_METHODS       = 3522;
FIFF.FIFF_MNE_FORWARD_SOLUTION_GRAD  = 3523;
%
% 3530... Covariance matrix
%
FIFF.FIFF_MNE_COV_KIND               = 3530;      % What kind of a covariance matrix
FIFF.FIFF_MNE_COV_DIM                = 3531;      % Matrix dimension
FIFF.FIFF_MNE_COV                    = 3532;      % Full matrix in packed representation (lower triangle)
FIFF.FIFF_MNE_COV_DIAG               = 3533;      % Diagonal matrix
FIFF.FIFF_MNE_COV_EIGENVALUES        = 3534;       % Eigenvalues and eigenvectors of the above
FIFF.FIFF_MNE_COV_EIGENVECTORS       = 3535;
FIFF.FIFF_MNE_COV_NFREE              = 3536;       % Number of degrees of freedom
%
% 3540... Inverse operator
%
% We store the inverse operator as the eigenleads, eigenfields,
% and weights
%
FIFF.FIFF_MNE_INVERSE_LEADS              = 3540;   % The eigenleads
FIFF.FIFF_MNE_INVERSE_LEADS_WEIGHTED     = 3546;   % The eigenleads (already weighted with R^0.5)
FIFF.FIFF_MNE_INVERSE_FIELDS             = 3541;   % The eigenfields
FIFF.FIFF_MNE_INVERSE_SING               = 3542;   % The singular values
FIFF.FIFF_MNE_PRIORS_USED                = 3543;   % Which kind of priors have been used for the source covariance matrix
FIFF.FIFF_MNE_INVERSE_FULL               = 3544;   % Inverse operator as one matrix
% This matrix includes the whitening operator as well
% The regularization is applied
FIFF.FIFF_MNE_INVERSE_SOURCE_ORIENTATIONS = 3545;  % Contains the orientation of one source per row
% The source orientations must be expressed in the coordinate system
% given by FIFF_MNE_COORD_FRAME
%
% 3550... Saved environment info
%
FIFF.FIFF_MNE_ENV_WORKING_DIR        = 3550;     % Working directory where the file was created
FIFF.FIFF_MNE_ENV_COMMAND_LINE       = 3551;     % The command used to create the file
%
% 3560... Miscellaneous
%
FIFF.FIFF_MNE_PROJ_ITEM_ACTIVE       = 3560;     % Is this projection item active?
FIFF.FIFF_MNE_EVENT_LIST             = 3561;     % An event list (for STI 014)
FIFF.FIFF_MNE_HEMI                   = 3562;     % Hemisphere association for general purposes
%
% 3570... Morphing maps
%
FIFF.FIFF_MNE_MORPH_MAP              = 3570;     % Mapping of closest vertices on the sphere
FIFF.FIFF_MNE_MORPH_MAP_FROM         = 3571;     % Which subject is this map from
FIFF.FIFF_MNE_MORPH_MAP_TO           = 3572;     % Which subject is this map to
%
% 3580... CTF compensation data
%
FIFF.FIFF_MNE_CTF_COMP_KIND         = 3580;     % What kind of compensation
FIFF.FIFF_MNE_CTF_COMP_DATA         = 3581;     % The compensation data itself
FIFF.FIFF_MNE_CTF_COMP_CALIBRATED   = 3582;     % Are the coefficients calibrated?

%
% 3700... Real-Time Communication
%
FIFF.FIFF_MNE_RT_COMMAND            = 3700;    	% Fiff Real-Time Command
FIFF.FIFF_MNE_RT_CLIENT_ID          = 3701;     % Fiff Real-Time mne_rt_server client id

%
% 3710... Real-Time Blocks
%
FIFF.FIFFB_MNE_RT_MEAS_INFO         = 3710;     % Fiff Real-Time Measurement Info

%
% Fiff values associated with MNE computations
%
FIFF.FIFFV_MNE_FIXED_ORI            = 1;
FIFF.FIFFV_MNE_FREE_ORI             = 2;

FIFF.FIFFV_MNE_MEG                  = 1;
FIFF.FIFFV_MNE_EEG                  = 2;
FIFF.FIFFV_MNE_MEG_EEG              = 3;

FIFF.FIFFV_MNE_UNKNOWN_COV          = 0;
FIFF.FIFFV_MNE_SENSOR_COV           = 1;
FIFF.FIFFV_MNE_NOISE_COV            = 1;         % This is what it should have been called
FIFF.FIFFV_MNE_SOURCE_COV           = 2;
FIFF.FIFFV_MNE_FMRI_PRIOR_COV       = 3;
FIFF.FIFFV_MNE_SIGNAL_COV           = 4;         % This will be potentially employed in beamformers
FIFF.FIFFV_MNE_DEPTH_PRIOR_COV      = 5;         % The depth weighting prior
FIFF.FIFFV_MNE_ORIENT_PRIOR_COV     = 6;     % The orientation prior
%
% Source space types (values of FIFF_MNE_SOURCE_SPACE_TYPE)
%
FIFF.FIFFV_MNE_SPACE_UNKNOWN  = -1;
FIFF.FIFFV_MNE_SPACE_SURFACE  = 1;
FIFF.FIFFV_MNE_SPACE_VOLUME   = 2;
FIFF.FIFFV_MNE_SPACE_DISCRETE = 3;
%
% Covariance matrix channel classification
%
FIFF.FIFFV_MNE_COV_CH_UNKNOWN  = -1;  % No idea
FIFF.FIFFV_MNE_COV_CH_MEG_MAG  =  0;  % Axial gradiometer or magnetometer [T]
FIFF.FIFFV_MNE_COV_CH_MEG_GRAD =  1;  % Planar gradiometer [T/m]
FIFF.FIFFV_MNE_COV_CH_EEG      =  2;  % EEG [V]
%
% Projection item kinds
%
FIFF.FIFFV_PROJ_ITEM_NONE           = 0;
FIFF.FIFFV_PROJ_ITEM_FIELD          = 1;
FIFF.FIFFV_PROJ_ITEM_DIP_FIX        = 2;
FIFF.FIFFV_PROJ_ITEM_DIP_ROT        = 3;
FIFF.FIFFV_PROJ_ITEM_HOMOG_GRAD     = 4;
FIFF.FIFFV_PROJ_ITEM_HOMOG_FIELD    = 5;
FIFF.FIFFV_MNE_PROJ_ITEM_EEG_AVREF  = 10;
%
% Additional coordinate frames
%
FIFF.FIFFV_MNE_COORD_TUFTS_EEG   =  300;         % For Tufts EEG data
FIFF.FIFFV_MNE_COORD_CTF_DEVICE  = 1001;     % CTF device coordinates
FIFF.FIFFV_MNE_COORD_CTF_HEAD    = 1004;         % CTF head coordinates
FIFF.FIFFV_MNE_COORD_MRI_VOXEL   = 2001;         % The MRI voxel coordinates
FIFF.FIFFV_MNE_COORD_RAS         = 2002;         % Surface RAS coordinates with non-zero origin
FIFF.FIFFV_MNE_COORD_MNI_TAL     = 2003;         % MNI Talairach coordinates
FIFF.FIFFV_MNE_COORD_FS_TAL_GTZ  = 2004;         % FreeSurfer Talairach coordinates (MNI z > 0)
FIFF.FIFFV_MNE_COORD_FS_TAL_LTZ  = 2005;         % FreeSurfer Talairach coordinates (MNI z < 0)
FIFF.FIFFV_MNE_COORD_FS_TAL      = 2006;         % FreeSurfer Talairach coordinates
%
% CTF coil and channel types
%
FIFF.FIFFV_REF_MEG_CH             = 301;
%
%   Data types
%
FIFF.FIFFT_VOID                  = 0;
FIFF.FIFFT_BYTE                  = 1;
FIFF.FIFFT_SHORT                 = 2;
FIFF.FIFFT_INT                   = 3;
FIFF.FIFFT_FLOAT                 = 4;
FIFF.FIFFT_DOUBLE                = 5;
FIFF.FIFFT_JULIAN                = 6;
FIFF.FIFFT_USHORT                = 7;
FIFF.FIFFT_UINT                  = 8;
FIFF.FIFFT_ULONG                 = 9;
FIFF.FIFFT_STRING                = 10;
FIFF.FIFFT_LONG                  = 11;
FIFF.FIFFT_DAU_PACK13            = 13;
FIFF.FIFFT_DAU_PACK14            = 14;
FIFF.FIFFT_DAU_PACK16            = 16;
FIFF.FIFFT_COMPLEX_FLOAT         = 20;
FIFF.FIFFT_COMPLEX_DOUBLE        = 21;
FIFF.FIFFT_OLD_PACK              = 23;
FIFF.FIFFT_CH_INFO_STRUCT        = 30;
FIFF.FIFFT_ID_STRUCT             = 31;
FIFF.FIFFT_DIR_ENTRY_STRUCT      = 32;
FIFF.FIFFT_DIG_POINT_STRUCT      = 33;
FIFF.FIFFT_CH_POS_STRUCT         = 34;
FIFF.FIFFT_COORD_TRANS_STRUCT    = 35;
FIFF.FIFFT_DIG_STRING_STRUCT     = 36;
FIFF.FIFFT_STREAM_SEGMENT_STRUCT = 37;
%
% Units of measurement
%
FIFF.FIFF_UNIT_NONE = -1;
%
% SI base units
%
FIFF.FIFF_UNIT_M   = 1;
FIFF.FIFF_UNIT_KG  = 2;
FIFF.FIFF_UNIT_SEC = 3;
FIFF.FIFF_UNIT_A   = 4;
FIFF.FIFF_UNIT_K   = 5;
FIFF.FIFF_UNIT_MOL = 6;
%
% SI Supplementary units
%
FIFF.FIFF_UNIT_RAD = 7;
FIFF.FIFF_UNIT_SR  = 8;
%
% SI base candela
%
FIFF.FIFF_UNIT_CD  = 9;
%
% SI derived units
%
FIFF.FIFF_UNIT_HZ  = 101;
FIFF.FIFF_UNIT_N   = 102;
FIFF.FIFF_UNIT_PA  = 103;
FIFF.FIFF_UNIT_J   = 104;
FIFF.FIFF_UNIT_W   = 105;
FIFF.FIFF_UNIT_C   = 106;
FIFF.FIFF_UNIT_V   = 107;
FIFF.FIFF_UNIT_F   = 108;
FIFF.FIFF_UNIT_OHM = 109;
FIFF.FIFF_UNIT_MHO = 110;
FIFF.FIFF_UNIT_WB  = 111;
FIFF.FIFF_UNIT_T   = 112;
FIFF.FIFF_UNIT_H   = 113;
FIFF.FIFF_UNIT_CEL = 114;
FIFF.FIFF_UNIT_LM  = 115;
FIFF.FIFF_UNIT_LX  = 116;
%
% Others we need
%
FIFF.FIFF_UNIT_T_M   = 201; % T/m
FIFF.FIFF_UNIT_AM    = 202; % Am
FIFF.FIFF_UNIT_AM_M2 = 203; % Am/m^2
FIFF.FIFF_UNIT_AM_M3 = 204; % Am/m^3
%
% Multipliers
%
FIFF.FIFF_UNITM_E    = 18;
FIFF.FIFF_UNITM_PET  = 15;
FIFF.FIFF_UNITM_T    = 12;
FIFF.FIFF_UNITM_MEG  = 6;
FIFF.FIFF_UNITM_K    = 3;
FIFF.FIFF_UNITM_H    = 2;
FIFF.FIFF_UNITM_DA   = 1;
FIFF.FIFF_UNITM_NONE = 0;
FIFF.FIFF_UNITM_D    = -1;
FIFF.FIFF_UNITM_C    = -2;
FIFF.FIFF_UNITM_M    = -3;
FIFF.FIFF_UNITM_MU   = -6;
FIFF.FIFF_UNITM_N    = -9;
FIFF.FIFF_UNITM_P    = -12;
FIFF.FIFF_UNITM_F    = -15;
FIFF.FIFF_UNITM_A    = -18;

%
% MNE epochs bookkeeping
%
FIFF.FIFFB_MNE_EPOCHS_SELECTION    = 3800;  % the epochs selection
FIFF.FIFFB_MNE_EPOCHS_DROP_LOG     = 3801;  % the drop log

return

end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Sigma_Vector=alpha_generator(l,M,nk)

if l==0 || l>M-1 || nk>M
    display('wrong input parameter')
    Sigma_Vector=[];
    return
end

flag_one=0;

for k=1:nk-1

    temp(k)=exp(j*2*pi*k/nk);
    if temp(k)==exp(j*2*pi*l/M)
        flag_one=1;
    end

end

if flag_one==0
    Sigma_Vector=0;
else
    Sigma_Vector=1;
end


text:
function b = NormalizeArray(a)
% Normalizes array a. This means that the minimum value will become 0 and
% the maximum value 1.
%
% a:            Input array.
%
% b:            Normalized output array
%
%     Jasper Uijlings - 2013

minVal = min(a(:));
maxVal = max(a(:));

diffVal = maxVal - minVal;

b = a - minVal;
if diffVal ~= 0
    b = b ./ diffVal;
end

author: zouchuhang

repo: LayoutNet


text:
function [lin_restr,nonlin_restr,markov_chains]=create_restrictions_and_markov_chains5(markov_chains)
% create_restrictions_and_markov_chains5 -- creates restrictions and
% markov chains for the SVAR model in which both coefficients and variance
% for the monetary policy equation are changing with two independent Markov
% processes
%
% ::
%
%
%   [lin_restr,nonlin_restr,tpl]=create_restrictions_and_markov_chains5(tpl)
%
% Args:
%
%    - **markov_chains** [empty|struct]: structure of previously defined
%    markov chains
%
% Returns:
%    :
%
%    - **lin_restr** [cell]: cell array of restrictions (see below).
%
%    - **nonlin_restr** [cell]: cell array of inequality restrictions
%
%    - **markov_chains** [struct]: modified markov chains
%
% Note:
%
%    - The syntax to construct a restriction
%      --> ai(eqtn)
%      --> ai(eqtn,vbl)
%      --> ai(eqtn,vbl,chain_name,state)
%      --> a(eqtn)
%      --> a(eqtn,vbl)
%      --> a(eqtn,vbl,chain_name,state)
%      - **eqtn** [integer]: integer
%      - **vbl** [integer|char]: integer or variable name
%      - **i** [integer]: lag
%      - **chain_name** [char]: name of the markov chain
%      - **state** [integer]: state number
%
%    - The lag coefficients are labelled a0, a1, a2,...,ak, for a model with k
%    lags. Obviously, a0 denotes the contemporaneous coefficients.
%
%    - The constant terms labelled c_1_1, c_2_2,...,c_n_n, for a model with n
%    endogenous variables.
%
%    - The standard deviations labelled s_1_1, s_2_2,...,s_n_n, for a
%    model with n endogenous variables.
%
% Example:
%
%    See also:

if nargin==0||isempty(markov_chains)

    markov_chains=struct('name',{},...
    'states_expected_duration',{},...
    'controlled_parameters',{});

end

% The parameter restrictions are identical to those in the model
% with regime switching in the policy coefficients only. Hence the mp_coef
% markov chain will also be common to those two models.
[lin_restr,nonlin_restr,markov_chains]=create_restrictions_and_markov_chains3(markov_chains);

% We add the volatility Markov chain from the model in which only the
% volatility of the monetary policy equation changes. N.B: In the process,
% we want to make sure we do not over-write the restrictions above!
%--------------------------------------------------------------------------
[~,~,markov_chains]=create_restrictions_and_markov_chains4(markov_chains);

end

author: jmaih


text:
function y = tanhminus1(x)
%TANHMINUS1 Hyperbolic tangent minus one.
%
%   TANHMINUS1(X) is TANH(X)-1 calculated in a way that is numerically better
%   when X is large and positive.
%
%   This example illustrates the difference
%
%      x = 17.5:0.01:20.5;
%      plot(x, tanh(x)-1, 'y-', x, tanhminus1(x), 'r-');

%   Author:      Peter J. Acklam
%   Time-stamp:  2003-10-13 15:04:13 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   % check number of input arguments
   error(nargchk(1, 1, nargin));

   y = -2 ./ (1 + exp(2 * x));

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function do_plot_unconditional_forecasts(mycast,m)

vnames=m.endogenous;

nvars=numel(vnames);

figure('name','unconditional forecasts');

for ivar=1:nvars

    subplot(nvars,1,ivar);

    plot(mycast.(vnames{ivar}),'linewidth',2);

    ylabel(vnames{ivar});

end

tmp=set([mycast.ygap mycast.pi mycast.FFR],'varnames',{'ygap','pi','FFR'});

disp('********** Unconditional forecasts ***********');
display(tmp)

end


text:
% DEMO  --  Chebyshev Polynomial Interpolation and Differentiation
% UPDATED  --  October 28, 2013
% Written by Matthew Kelly, Cornell University
%
% This script demonstrates how chebyshev interpolants can be used to get
% accurate derivate information.
%

%% General Settings
clear; clc;

%What order should the approximation be?
order = 34;  %Low = 25, %Med = 1000, %High = 100000

%How many points should be used for error calculations and plotting?
nTime = 1000;

%What domain should we be looking at?
d = [0,1.3];

%Time for use in plots
time = linspace(d(1),d(2),nTime);

% %The following version of time will yield more accurate results and a
% %faster runtime because the points in time do not line up exactly with
% the chebyshev grid points.
% time = linspace(d(1)-1e-4, d(2)+2e-4,nTime);

%% DEMO  --  Fit to analytic function

%Set up analytic function
IO.domain = d;
IO.userFunc = @testFunction;

%Get the values at each of the chebyshev nodes
f = chebyshevFit(IO, order);         %Chebyshev Values
tic
[y,Dy,DDy,DDDy] = chebyshevInterpolate(f,time,d);   %Approximation
toc

%get the exact (analytic) solution for each derivative
[g,Dg,DDg,DDDg] = testFunction(time);

%Show results
figure(401); clf;
subplot(4,2,1);
    plot(time,y,'b-','LineWidth',2)
    title(['function approximation - order ' num2str(order)]);
subplot(4,2,3)
    plot(time,Dy,'b-','LineWidth',2)
    title(['derivative approximation - order ' num2str(order)]);
subplot(4,2,5)
    plot(time,DDy,'b-','LineWidth',2)
    title(['second derivative approximation - order ' num2str(order)]);
subplot(4,2,7)
    plot(time,DDDy,'b-','LineWidth',2)
    title(['third derivative approximation - order ' num2str(order)]);
subplot(4,2,2);
    semilogy(time,abs(g-y))
    title('error in function');
subplot(4,2,4)
    semilogy(time,abs(Dg-Dy))
    title('error in derivative');
subplot(4,2,6)
    semilogy(time,abs(DDg-DDy))
    title('error in second derivative');
subplot(4,2,8)
    semilogy(time,abs(DDDg-DDDy))
    title('error in third derivative');

author: MatthewPeterKelly

repo: dscTutorials


text:
function [sampled_graphs, accept_ratio, num_edges] = learn_struct_mcmc(data, ns, varargin)
% LEARN_STRUCT_MCMC  Monte Carlo Markov Chain search over DAGs assuming fully observed data
% [sampled_graphs, accept_ratio, num_edges] = learn_struct_mcmc(data, ns, ...)
%
% data(i,m) is the value of node i in case m.
% ns(i) is the number of discrete values node i can take on.
%
% sampled_graphs{m} is the m'th sampled graph.
% accept_ratio(t) = acceptance ratio at iteration t
% num_edges(t) = number of edges in model at iteration t
%
% The following optional arguments can be specified in the form of name/value pairs:
% [default value in brackets]
%
% scoring_fn - 'bayesian' or 'bic' [ 'bayesian' ]
%              Currently, only networks with all tabular nodes support Bayesian scoring.
% type       - type{i} is the type of CPD to use for node i, where the type is a string
%              of the form 'tabular', 'noisy_or', 'gaussian', etc. [ all cells contain 'tabular' ]
% params     - params{i} contains optional arguments passed to the CPD constructor for node i,
%              or [] if none.  [ all cells contain {'prior', 1}, meaning use uniform Dirichlet priors ]
% discrete   - the list of discrete nodes [ 1:N ]
% clamped    - clamped(i,m) = 1 if node i is clamped in case m [ zeros(N, ncases) ]
% nsamples   - number of samples to draw from the chain after burn-in [ 100*N ]
% burnin     - number of steps to take before drawing samples [ 5*N ]
% init_dag   - starting point for the search [ zeros(N,N) ]
%
% e.g., samples = my_learn_struct_mcmc(data, ns, 'nsamples', 1000);
%
%
% Modified by Mingyi Wang  (mingyiwang@hotmail.com) Sep 18, 2006 (based on Sonia Leach (SML)'s version ( 2/4/02, 9/5/03))
%
% Some bugs in update_ancestor_matrix() were fixed. This function can call mk_nbrs_of_digraph properly
%

[n ncases] = size(data);

% set default params
type = cell(1,n);
params = cell(1,n);
for i=1:n
 type{i} = 'tabular';
 %params{i} = { 'prior', 1};
 params{i} = { 'prior_type', 'dirichlet', 'dirichlet_weight', 1 };
end
scoring_fn = 'bayesian';
discrete = 1:n;
clamped = zeros(n, ncases);
nsamples = 100*n;
burnin = 5*n;
dag = zeros(n);

args = varargin;
nargs = length(args);
for i=1:2:nargs
 switch args{i},
  case 'nsamples',   nsamples = args{i+1};
  case 'burnin',     burnin = args{i+1};
  case 'init_dag',   dag = args{i+1};
  case 'scoring_fn', scoring_fn = args{i+1};
  case 'type',       type = args{i+1};
  case 'discrete',   discrete = args{i+1};
  case 'clamped',    clamped = args{i+1};
  case 'gconstraint', gconstraint=args{i+1};  %Added by mingyi
  case 'params',     if isempty(args{i+1}), params = cell(1,n); else params = args{i+1};  end

 end
end

% We implement the fast acyclicity check described by P. Giudici and R. Castelo,
% "Improving MCMC model search for data mining", submitted to J. Machine Learning, 2001.

% SML: also keep descendant matrix C
use_giudici = 1;
%use_giudici = 0; %Revised by MIngyi
if use_giudici
 [nbrs, ops, nodes, A] = mk_nbrs_of_digraph(dag);
else
 [nbrs, ops, nodes] = mk_nbrs_of_dag(dag);
 A = [];
end

num_accepts = 1;
num_rejects = 1;
T = burnin + nsamples;
accept_ratio = zeros(1, T);
num_edges = zeros(1, T);
sampled_graphs = cell(1, nsamples);
%sampled_bitv = zeros(nsamples, n^2);

for t=1:T
 [dag, nbrs, ops, nodes, A, accept] = take_step(dag, nbrs, ops, ...
                    nodes, ns, data, clamped, A, ...
                      scoring_fn, discrete, type, params);
 num_edges(t) = sum(dag(:));
 num_accepts = num_accepts + accept;
 num_rejects = num_rejects + (1-accept);
 accept_ratio(t) =  num_accepts/num_rejects;
 if t > burnin
   sampled_graphs{t-burnin} = dag;
   %sampled_bitv(t-burnin, :) = dag(:)';
 end
 fprintf('MCMC: %d/%d\n',t,T);
end

%%%%%%%%%

function [new_dag, new_nbrs, new_ops, new_nodes, A,  accept] = ...
   take_step(dag, nbrs, ops, nodes, ns, data, clamped, A,  ...
     scoring_fn, discrete, type, params, prior_w)

global gconstraint;   %Added by Mingyi
use_giudici = ~isempty(A);
if use_giudici
 [new_dag, op, i, j, new_A] =  pick_digraph_nbr(dag, nbrs, ops, nodes,A); % updates A
 [new_nbrs, new_ops, new_nodes] =  mk_nbrs_of_digraph(new_dag,new_A);
else
 d = sample_discrete(normalise(ones(1, length(nbrs))));
 new_dag = nbrs{d};
 op = ops{d};
 i = nodes(d, 1); j = nodes(d, 2);
 [new_nbrs, new_ops, new_nodes] = mk_nbrs_of_dag1(new_dag);
end
%For debug
% fprintf('op:%s,i:%d,j:%d\n',op,i,j);
% if ~acyclic(new_dag)
%     error('new dag must be acyclic!')
% end
% if size(find(diag(new_A)),1)>0
%   A=A
%   new_A=new_A
%   error('new A must be acyclic!')
%  end
%debug ends

bf =  bayes_factor(dag, new_dag, op, i, j, ns, data, clamped, scoring_fn, discrete, type, params);

%R = bf * (new_prior / prior) * (length(nbrs) / length(new_nbrs));
R = bf * (length(nbrs) / length(new_nbrs));
u = rand(1,1);
if u > min(1,R) % reject the move
 accept = 0;
 new_dag = dag;
 new_nbrs = nbrs;
 new_ops = ops;
 new_nodes = nodes;
else
 accept = 1;
 if use_giudici
    A = new_A; % new_A already updated in pick_digraph_nbr
 end
end

%%%%%%%%%

function bfactor = bayes_factor(old_dag, new_dag, op, i, j, ns, data, clamped, scoring_fn, discrete, type, params)

u = find(clamped(j,:)==0);
LLnew = score_family(j, parents(new_dag, j), type{j}, scoring_fn, ns, discrete, data(:,u), params{j});
LLold = score_family(j, parents(old_dag, j), type{j}, scoring_fn, ns, discrete, data(:,u), params{j});
bf1 = exp(LLnew - LLold);

if strcmp(op, 'rev')  % must also multiply in the changes to i's family
 u = find(clamped(i,:)==0);
 LLnew = score_family(i, parents(new_dag, i), type{i}, scoring_fn, ns, discrete, data(:,u), params{i});
 LLold = score_family(i, parents(old_dag, i), type{i}, scoring_fn, ns, discrete, data(:,u), params{i});
 bf2 = exp(LLnew - LLold);
else
 bf2 = 1;
end
bfactor = bf1 * bf2;

%%%%%%%% Giudici stuff follows %%%%%%%%%%

% SML: This now updates A as it goes from digraph it choses
function [new_dag, op, i, j, new_A] = pick_digraph_nbr(dag, digraph_nbrs, ops, nodes, A)

d = sample_discrete(normalise(ones(1, length(digraph_nbrs))));
%d = myunidrnd(length(digraph_nbrs),1,1);
i = nodes(d, 1); j = nodes(d, 2);
new_dag = digraph_nbrs(:,:,d);

op = ops{d};
new_A = update_ancestor_matrix(A, op, i, j, dag);
%for debug
% if op=='add'
%     if ~(dag(i,j)==0 & new_dag(i,j)==1)
%         fprintf('error add\n');
%     end
% end
% if op=='del'
%     if ~(dag(i,j)==1 & new_dag(i,j)==0)
%       fprintf('new dag del calculation is error!\n')
%     end
% end
% if op=='rev'
%     if ~(dag(i,j)==1 & dag(j,i)==0 & new_dag(i,j)==0 & new_dag(j,i)==1)
%       fprintf('new dag rev calculation is error!\n')
%     end
% end
% new_AA = reachability_graph(new_dag');
% if find(diag(new_AA)==1)
%     fprintf('cyclic\n');
% end
% if ~isequal(new_A,new_AA)
%    fprintf('new A calculation is error!\n')
% end
%debug ends

%%%%%%%%%%%%%%

function A = update_ancestor_matrix(A,  op, i, j, dag)

switch op
case 'add',
 A = do_addition(A,  op, i, j, dag);
case 'del',
 A = do_removal(A,  op, i, j, dag);
case 'rev',
 A = do_removal(A,  op, i, j, dag);
 A = do_addition(A,  op, j, i, dag);
end

%%%%%%%%%%%%

function A = do_addition(A, op, i, j, dag)

A(j,i) = 1;     % i is an ancestor of j
anci = find(A(i,:));
if ~isempty(anci)
 A(j,anci) = 1;   % all of i's ancestors are added to Anc(j)
end

descj = find(A(:,j));  %all the descendants of j are selected
if ~isempty(descj)
 for k=descj(:)'
   A(k,i) = 1;        % i is the ancestor of descj
   if ~isempty(anci)  % all of i's ancestors are also the ancestor of each descendant of j
       A(k,anci)=1;
   end
 end
end

%%%%%%%%%%%

function A = do_removal(A, op, i, j, dag)
descj = find(A(:,j));
A = update_row(A,i, j, dag);   % compute the A(j,:) row for dag i->j removal

if ~isempty(descj)
  order = topological_sort(dag);  %all the parent nodes are before to the children nodes
  [junk, perm] = sort(order);     %node i is perm(i)-TH in order
  descj_topnum = perm(descj);     %descj(i) is descj_topnum(i)-th in order

% SML: now re-sort descj by rank in descj_topnum
  [junk, perm] = sort(descj_topnum);
  descj = descj(perm);
  for k = descj(:)'
    A = old_update_row(A, k, dag);
  end
end

%%%%%%%%%

function A = update_row(A, i,j, dag)
% We compute row j of A
A(j, :) = 0;
ps = parents(dag, j);
ps=setdiff(ps,i);  % All the parents except i
if ~isempty(ps)
 A(j, ps) = 1;
end
for k=ps(:)'
 anck = find(A(k,:));
 if ~isempty(anck)
   A(j, anck) = 1;
 end
end

%%%%%%%%%

function A = old_update_row(A, j, dag)

% We compute row j of A
A(j, :) = 0;
ps = parents(dag, j);
if ~isempty(ps)
 A(j, ps) = 1;
end
for k=ps(:)'
 anck = find(A(k,:));
 if ~isempty(anck)
   A(j, anck) = 1;
 end
end

%%%%%%%%

function A = init_ancestor_matrix(dag)

order = topological_sort(dag);
A = zeros(length(dag));
for j=order(:)'
 A = update_row(A, j, dag);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
classdef StressNormShapeFuncCreator < handle

    properties (Access = private)
        outFile
        stressShape
        stressShapeDB
    end

    methods (Access = public)

        function obj = StressNormShapeFuncCreator(d)
            obj.init(d)
            obj.createStressShapeDataBase();
            obj.createStressShape();
        end

        function s = getStressNormShape(obj)
            s = obj.stressShape;
        end

    end

    methods (Access = private)

        function init(obj,d)
            obj.outFile       = d.outFile;
        end

        function createStressShape(obj)
            dB = obj.stressShapeDB();
            sF = ShFunc_StressNorm(dB);
%             sF.filter.preProcess();
            obj.stressShape = sF;
        end

        function createStressShapeDataBase(obj)
            dB.filename    = obj.outFile;
            dB.TOL         = obj.createMaterialProperties();
            dB.material    = 'ISOTROPIC';
            dB.method      = 'SIMPALL';
            dB.pdim        = '2D';
            dB.stressHomog = [1,0,0]';
            dB.filter      = 'P1';
            dB.optimizer   = 'SLERP';
            dB.pdim        = '2D';
            dB.scale       = 'MICRO';
            obj.stressShapeDB = dB;
        end
    end

    methods (Access = private, Static)

        function matProp = createMaterialProperties()
            matProp.rho_plus  =  1;
            matProp.rho_minus =  0;
            matProp.E_plus    =  1;
            matProp.E_minus   =  1.0000e-03;
            matProp.nu_plus   =  0.3333;
            matProp.nu_minus  =  0.3333;
        end

    end

end

author: SwanLab

repo: Swan


text:
classdef libad9361 < handle

    properties (Hidden)
        libad9361_win = 'https://github.com/analogdevicesinc/libad9361-iio/releases/download/v0.2/libad9361-0.2-win64.zip';
        libad9361_osx = 'https://github.com/analogdevicesinc/libad9361-iio/releases/download/v0.2/ad9361-0.2-Darwin.tar.gz';
        libad9361_linux= 'http://swdownloads.analog.com/cse/travis_builds/master_latest_libad9361-iio-trusty.tar.gz';
        os
        headerFileNames = {'ad9361.h','ad9361-wrapper.h'};
        libraryFileNames = {'ad9361','libad9361.so','libad9361.dll'};
        tmpdir = 'tmp_deps_install';
    end

    methods (Access=protected, Hidden)

        function loc = getDependencyTargetLocation(~)
            rootDir = fileparts(strtok(mfilename('fullpath'), '+'));
            loc = fullfile(rootDir,'deps','ad9361');
        end

        function [files,dlfn,fulltmp] = downloadFiles(obj)
            if ismac
                dlfn = websave('dl.tar.gz',obj.libad9361_osx);
                fulltmp = fullfile(pwd,obj.tmpdir);
                files = untar(dlfn,fulltmp);
            elseif ispc
                dlfn = websave('dl.zip',obj.libad9361_win);
                fulltmp = fullfile(pwd,obj.tmpdir);
                files = unzip(dlfn,fulltmp);
            else
                dlfn = websave('dl.tar.gz',obj.libad9361_linux);
                fulltmp = fullfile(pwd,obj.tmpdir);
                files = untar(dlfn,fulltmp);
            end
        end

    end

    methods

        function [check,lib,head] = checkForDependencies(obj)
            loc = getDependencyTargetLocation(obj);
            lib = false;
            for l = 1:length(obj.libraryFileNames)
                lib = lib || exist(fullfile(loc,obj.libraryFileNames{l}),'file');
            end
            head = true;
            for h = 1:length(obj.headerFileNames)
                head = head && exist(fullfile(loc,'include',obj.headerFileNames{h}),'file');
            end
            check = lib && head;
        end

        function download_libad9361(obj)
            if obj.checkForDependencies()
                fprintf('Dependencies already installed, not reinstalling\n');
                return
            end
            fprintf('Downloading libad9361 library\n');
            [files,dlfn,fulltmp] = downloadFiles(obj);
            target = obj.getDependencyTargetLocation();
            if ~exist(target,'dir')
                mkdir(target);
            end
            if ~exist(fullfile(target,'include'),'dir')
                mkdir(fullfile(target,'include'));
            end
            % Includes
            for f = 1:length(files)
                if contains(files{f},obj.headerFileNames)
                    copyfile(files{f},fullfile(target,'include'));
                end
            end
            % Libraries
            for f = 1:length(files)
                if exist(files{f},'dir')
                    continue;
                end
                fp = split(files{f}, '/'); fp = fp{end};
                s = dir(files{f});
                if contains(fp,obj.libraryFileNames) && s.bytes>0
                    if ismac && ~contains(fp,'.')
                        copyfile(files{f},fullfile(target,'libad9361.dylib'));
                    end
                    if ispc && contains(files{f},'.dll')
                        copyfile(files{f},target);
                    end
                    if isunix
                        copyfile(files{f},fullfile(target,'libad9361.so'));
                    end
                end
            end
            % Cleanup
            rmdir(fulltmp,'s');
            delete(dlfn);
            fprintf('Installing libad9361 complete\n');
            % Path stuff
            fprintf('Adding libad9361 to path\n');
            addpath(genpath(target));
            try
                warning('error','MATLAB:SavePath:PathNotSaved');
                savepath();
            catch
               warning(['savepath failed. libad9361 is currently on path, '...
                   'but is likely to not be across MATLAB restarts.',...
                   'Please add "',target,'" to path']);
            end
        end
    end
end

author: analogdevicesinc


text:
function [D] = spm_eeg_invert_classic_mix(D,val,Qpriors,surfind,ugainfiles)
%function [D] = spm_eeg_invert_classic_mix(D,val,Qpriors,surfind,ugainfiles)
%%
% ReML inversion of multiple posterior current variances from previous
% iterations spm_eeg_invert_classic or spm_eeg_invert
% ReML estimation of regularisation hyperparameters using the
% spatiotemporal hierarchy implicit in EEG/MEG data
%
% D contains the data and inversion parameters (see
% spm_eeg_invert.m/spm_eeg_invert_classic.m)
% val the inversion index
% Qpriors is N solutions of rows by Nd variance estimates
% surfind is N solutions long and contains indices into ugainfiles to these priors with different lead field structures
%% ugainfiles are the SPMgain matrices for the different surfaces
%
% Output D will have a solution which is optimal REML mixture of Qpriors
%
% Created by:
%               Gareth Barnes - g.barnes@ucl.ac.uk
%
% This version is for single subject single modality analysis and therefore
% contains none of the associated scaling factors.
%
% $Id: spm_eeg_invert_classic_mix.m 6077 2014-06-30 16:55:03Z spm $

Nl = length(D);

if Nl>1,
    error('function only defined for a single subject');
end;

% D - SPM data structure
%==========================================================================
if nargin > 1
    D.val = val;
elseif ~isfield(D, 'val')
    D.val = 1;
end

val=D.val;

inverse   = D.inv{val}.inverse;

% forward model
%--------------------------------------------------------------------------

% defaults
%--------------------------------------------------------------------------

try, type = inverse.type;   catch, type = 'GS';     end
try, s    = inverse.smooth; catch, s    = 0.6;      end
try, Np   = inverse.Np;     catch, Np   = 256;      end
try, Nr   = inverse.Nr;     catch, Nr   = 16;       end %% requested number of temporal modes, could be changed depending on svd
try, xyz  = inverse.xyz;    catch, xyz  = [0 0 0];  end
try, rad  = inverse.rad;    catch, rad  = 128;      end
try, hpf  = inverse.hpf;    catch, hpf  = 48;       end %% need to one day put these the correct way round
try, lpf  = inverse.lpf;    catch, lpf  = 0;        end
try, sdv  = inverse.sdv;    catch, sdv  = 4;        end
try, Han  = inverse.Han;    catch, Han  = 1;        end
try, woi  = inverse.woi;    catch, woi  = [];       end
try, Nm   = inverse.Nm;     catch, Nm   = [];       end
try, Nt   = inverse.Nt;     catch, Nt   = [];       end %% fixed number of temporal modes
try, Ip   = inverse.Ip;     catch, Ip   = [];       end
try, SHUFFLELEADS=inverse.SHUFFLELEADS;catch, SHUFFLELEADS=[];end

% defaults
%--------------------------------------------------------------------------
%type = inverse.type;    % Type of inversion scheme

% get specified modalities to invert (default to all)

%--------------------------------------------------------------------------
modalities = D.inv{val}.forward.modality;       % MEG in this case

Nmax  = 16;         % max number of temporal modes

% check lead fields and get number of dipoles (Nd) and channels (Nc)
%==========================================================================

if inverse.Nm~=length(inverse.Ic{1}),
    disp('Using reduced spatial modes');
    if length(unique(surfind))>1,
        error('Cannot merge different surface files with reduced spatial modes at the moment');
    end;
    U=inverse.U{1};
else
    disp('Using all spatial modes');
    U=eye(inverse.Nm);
end;

fprintf('Checking leadfields')

usurfind=unique(surfind);
Nsurf=length(usurfind);
if Nsurf~=size(ugainfiles,1); %% number of different surfaces
    error('need a lead field file for each surface');
end;

L=[];
g_ind=zeros(Nsurf,2); %% index start and end of each set of lead fields
for f=1:Nsurf; %% load in  all leadfield

    fname = deblank(ugainfiles(f,:));
    G = load(fullfile(fname)); % Relative path

    label = G.label;
    g_ind(f,1)=size(L,2)+1; %% number of lead fields per surface
    L     =[L G.G];
    g_ind(f,2)=size(L,2); %% number of lead fields per surface
    clear G;

end;

if size(modalities,1)>1,
    error('not defined for multiple modalities');
end;
Ic  = setdiff(D.indchantype(modalities), badchannels(D));
Nd    = size(L,2);      % Number of dipoles

fprintf(' - done\n')

% check for (e.g., empty-room) sensor components (in Qe)
%==========================================================================
QE = 1;                     % No empty room noise measurement

%==========================================================================
% Spatial projectors - keep these to map directly to channels for now
%==========================================================================

A=U;
UL=U*L;

clear L

Nm    = size(UL,1);         % Number of spatial projectors

% Report
%----------------------------------------------------------------------
fprintf('Using %d spatial modes',Nm)

% No dipole is eliminated
%--------------------------------------------------------------------------
Is    = 1:Nd;               % Accepted dipoles
%Ns    = length(Is);         % Ns = Nd in this case

%==========================================================================
% Temporal projector
%==========================================================================
AY    = {};                                      % pooled response for MVB
AYYA  = 0;                                       % pooled response for ReML

% Time-window of interest
%----------------------------------------------------------------------

if isempty(woi)
    w      = 1000*[min(D.time) max(D.time)];
else
    w=woi; %% in milliseconds
end;

It     = (w/1000 - D.timeonset)*D.fsample + 1;
It     = max(1,It(1)):min(It(end), length(D.time));
It     = fix(It);

% Peristimulus time
%----------------------------------------------------------------------
pst    = 1000*D.time;                   % peristimulus time (ms)
pst    = pst(It);                       % windowed time (ms)
dur    = (pst(end) - pst(1))/1000;      % duration (s)
dct    = (It - It(1))/2/dur;            % DCT frequencies (Hz)
Nb     = length(It);                    % number of time bins

% Serial correlations
%----------------------------------------------------------------------
K      = exp(-(pst - pst(1)).^2/(2*sdv^2));
K      = toeplitz(K);
qV     = sparse(K*K');

% Confounds and temporal subspace
%----------------------------------------------------------------------

T      = spm_dctmtx(Nb,Nb);         % use plot(T) here!

j      = find( (dct >= lpf) & (dct <= hpf) ); %% THis is the wrong way round but leave for nowfor compatibility with spm_eeg_invert
T      = T(:,j);                    % Apply the filter to discrete cosines
dct    = dct(j);                    % Frequencies accepted

%% Hanning window
%----------------------------------------------------------------------

if Han
    W  = sparse(1:Nb,1:Nb,spm_hanning(Nb)); %% use hanning unless specified
else
    W=1;
end;

% get trials or conditions
%----------------------------------------------------------------------
try
    trial = D.inv{D.val}.inverse.trials;
catch
    trial = D.condlist;
end
Ntrialtypes=length(trial);
% get temporal covariance (Y'*Y) to find temporal modes
%======================================================================
%MY    = cell(Nmod,1);                        % mean response
YTY   = sparse(0);                           % accumulator

% get (spatially aligned) data
%------------------------------------------------------------------

YY    = 0;
%    MY{m} = 0;
N=0;
badtrialind=D.badtrials;
for j = 1:Ntrialtypes,                          % pool over conditions
    c     = D.indtrial(trial{j});     % and trials
    c=setxor(c,badtrialind);
    length(c)
    Nk    = length(c);
    for k = 1:Nk
        Y     = A*D(Ic,It,c(k));

        YY    = YY + Y'*Y;
        N     = N + 1;
    end
end
YY=YY./N;

% Apply any Hanning and filtering
%------------------------------------------------------------------
YY         = W'*YY*W;     % Hanning
YTY         = T'*YY*T;     % Filter

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

if isempty(Nt),

    [U E]  = spm_svd(YTY,exp(-8));          % get temporal modes
    if isempty(U),
        warning('nothing found using spm svd, using svd');
        [U E]  = svd(YTY);          % get temporal modes
    end;
    E      = diag(E)/trace(YTY);            % normalise variance
    Nr     = min(length(E),Nmax);           % number of temporal modes
    Nr=max(Nr,1); %% use at least one mode
else
    [U E]  = svd(YTY);          % get temporal modes

    E      = diag(E)/trace(YTY);            % normalise variance
    disp('Fixed number of temporal modes');
    Nr=Nt;
end;

V      = U(:,1:Nr);                     % temporal modes
VE     = sum(E(1:Nr));                  % variance explained

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

% projection and whitening
%----------------------------------------------------------------------
S      = T*V;                           % temporal projector
Vq     = S*pinv(S'*qV*S)*S';            % temporal precision

% get spatial covariance (Y*Y') for Gaussian process model
%======================================================================

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

% loop over Ntrialtypes trial types
%----------------------------------------------------------------------
UYYU = 0;
AYYA=0;
Nn    =0;                             % number of samples
AY={};
Ntrials=0;
for j = 1:Ntrialtypes,

    UY{j} = sparse(0);
    c       = D.indtrial(trial{j});
    Nk      = length(c);

    % loop over epochs
    %------------------------------------------------------------------
    for k = 1:Nk

        % stack (scaled aligned data) over modalities
        %--------------------------------------------------------------

        Y       = D(Ic,It,c(k))*S;
        Y=A*Y; %% CHANGE

        % accumulate first & second-order responses
        %--------------------------------------------------------------
        Nn       = Nn + Nr;         % number of samples
        %Y           = spm_cat(MY);           % contribution to ERP, Y and MY are the same for 1 modality
        YY          = Y*Y';                  % and covariance
        Ntrials=Ntrials+1;
        %YYep{k}=YY;%% one for each trial and condition

        % accumulate statistics (subject-specific)
        %--------------------------------------------------------------
        UY{j}     = UY{j} + Y;           % condition-specific ERP
        UYYU     = UYYU + YY;          % subject-specific covariance

        % and pool for optimisation of spatial priors over subjects
        %--------------------------------------------------------------
        AY{end + 1} = Y;                     % pooled response for MVB
        AYYA        = AYYA    + YY;          % pooled response for ReML

    end
end

AYYA=AYYA./Nn;
AY=spm_cat(AY);

% assuming equal noise over subjects (Qe) and modalities AQ
%--------------------------------------------------------------------------
AQeA   = A*QE*A';           % Note that here it is A*A'
Qe{1}  = AQeA/(trace(AQeA)); % it means IID noise in virtual sensor space

% Inverse solution
%==========================================================================
QP     = {};
LQP    = {};
LQPL   = {};

for j=1:size(Qpriors,1),
    % Accumulate empirical priors (New set of patches for the second inversion)
    %------------------------------------------------------------------

    disp(sprintf('Adding prior set %d of %d',j,size(Qpriors,1)));
    dum=sparse(zeros(1,Nd));

    dum(g_ind(surfind(j),1):g_ind(surfind(j),2))=Qpriors(j,:);

    qp=sparse(diag(dum)); %% qp is Nd*Nd
    QP{j}=sparse(diag(qp));
    LQP{j}=UL*qp;
    LQPL{j} = LQP{j}*UL';
    clear qp
end;

%==========================================================================
% Step 2: Re-estimate for each subject separately (fusing all modalities)
%==========================================================================

fprintf('Inverting subject 1\n')

% generate sensor component (Qe) per modality
%----------------------------------------------------------------------
AQeA  = A*QE*A';                % Again it is A*A'
AQ    = AQeA/(trace(AQeA));

% using spatial priors
%----------------------------------------------------------------------
Np    = length(LQPL);       % Final number of priors
Ne    = length(Qe);         % Sensor noise prior
Q     = [Qe LQPL];

% re-do ReML (with informative hyperpriors)
% Here is performed the second inversion
%======================================================================

Q0          = exp(-2)*trace(AYYA)*AQ/trace(AQ);
[Cy,h,Ph,F] = spm_reml_sc(AYYA,[],Q,1,-4,16,Q0);

% Data ID
%----------------------------------------------------------------------
% When should I use the data ID?
% For comparison purposes it is necessary to guarantee that the IDs have
% the same value.
%
% When using the same dataset but different lead field matrices is a good
% example of fail, because the spatial projector will generate different
% virtual sensors, and therefore different data for the inversion.
ID    = spm_data_id(YY);

% Covariance: sensor space - Ce and source space - L*Cp
%----------------------------------------------------------------------
Cp    = sparse(0);
LCp   = sparse(0);
hp    = h(Ne + (1:Np));
for j = 1:Np
    Cp  =  Cp + hp(j)*QP{j};
    LCp = LCp + hp(j)*LQP{j};
end

% MAP estimates of instantaneous sources
%======================================================================
% This is equivalent to M = Cp*UL'*inv(Qe + UL*Cp*UL'))
% with Cp the posterior source covariance (with optimal h values)
M     = LCp'/Cy;

% conditional variance (leading diagonal)
% Cq    = Cp - Cp*L'*iC*L*Cp;
%----------------------------------------------------------------------
Cq    = Cp - sum(LCp.*M')';

% evaluate conditional expectation
%----------------------------------------------------------------------
% evaluate conditional expectation (of the sum over trials)
%----------------------------------------------------------------------
SSR   = 0;
SST   = 0;
J     = {};

for j = 1:Ntrialtypes

    % trial-type specific source reconstruction
    %------------------------------------------------------------------
    J{j} = M*UY{j};

    % sum of squares
    %------------------------------------------------------------------
    SSR  = SSR + sum(var((UY{j} - UL*J{j}))); %% changed variance calculation
    SST  = SST + sum(var( UY{j}));

end

% J = M*Y;
%
% % sum of squares
% %------------------------------------------------------------------
% SSR  = sum(var((Y - UL*J),0,2));
% SST  = sum(var(Y,0,2));

% accuracy; signal to noise (over sources)
%======================================================================
R2   = 100*(SST - SSR)/SST;
fprintf('Percent variance explained %.2f (%.2f)\n',full(R2),full(R2*VE));

% Save results
% DEMO: WARNING! These results are not coincident in format with
%                those generated in the SPM8
%======================================================================
inverse.type   = type;                 % inverse model
inverse.M      = M;                    % MAP projector (reduced)
inverse.J   = J;                    % Conditional expectation
inverse.Y      = Y;                    % ERP data (reduced)
inverse.L      = UL;                   % Lead-field (reduced)
inverse.qC     = Cq;                   % spatial covariance
inverse.qV     = Vq;                   % temporal correlations
inverse.T      = S;                    % temporal projector
inverse.U      = {A};                    % spatial projector
inverse.Is     = Is;                   % Indices of active dipoles
inverse.It     = It;                   % Indices of time bins
try
    inverse.Ic{1}     = Ic;                   % Indices of good channels
catch
    inverse.Ic    = Ic;                   % Indices of good channels
end;
inverse.Nd     = Nd;                   % number of dipoles
inverse.pst    = pst;                  % peristimulus time
inverse.dct    = dct;                  % frequency range
inverse.F      = F;                    % log-evidence
inverse.ID     = ID;                   % data ID
inverse.R2     = R2;                   % variance explained (reduced)
inverse.VE     = R2*VE;                % variance explained
inverse.woi    = w;                    % time-window inverted

inverse.modality = modalities;         % modalities inverted

% save in struct
%----------------------------------------------------------------------
D.inv{val}.inverse = inverse;
D.inv{val}.method  = 'Imaging';

% display
%======================================================================
spm_eeg_invert_display(D);
drawnow

return

author: spm


text:
function fx = p38_fun ( n, x )

%*****************************************************************************80
%
%% P38_FUN evaluates the integrand for problem 38.
%
%  Discussion:
%
%    The problem has a parameter ALPHA that can be set by calling
%    P38_PARAM_SET.
%
%    The integrand oscillates more strongly as ALPHA is increased.
%
%    The suggested range for ALPHA is 0 to 10.
%
%  Interval:
%
%    0 <= x <= pi
%
%  Integrand:
%
%    cos ( 2^ALPHA * sin ( x ) )
%
%  Exact Integral:
%
%    pi * J0 ( 2^ALPHA )
%
%    where J0 ( x ) is the J Bessel function of order 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Piessens, Elise de Doncker-Kapenga,
%    Christian Ueberhuber, David Kahaner,
%    QUADPACK: A Subroutine Package for Automatic Integration,
%    Springer, 1983, page 83.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  alpha = p38_param_get ( );

  fx = cos ( 2.0^alpha * sin ( x ) );

  return
end


text:
function l4vec_uniform_test ( )

%*****************************************************************************80
%
%% L4VEC_UNIFORM_TEST tests L4VEC_UNIFORM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 December 2014
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'L4VEC_UNIFORM_TEST\n' );
  fprintf ( 1, '  L4VEC_UNIFORM computes a vector of\n' );
  fprintf ( 1, '  pseudorandom logical values.\n' );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The initial seed is %d\n', seed );

  [ l, seed ] = l4vec_uniform ( n, seed );

  l4vec_print ( n, l, '  Uniform L4VEC:' );

  return
end


text:
function [s, S_l, S_t] = pluckerSegment(l,t)

% PLUCKERSEGMENT Segment from Plucker line and endpoint abscissas.
%   PLUCKERSEGMENT(L,T) returns a segment [p1;p2] of two endpoints given a
%   Plucker line L and its two endpoint abscissas in T.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

[e1,e2,E1_l,E2_l,E1_t1,E2_t2] = pluckerEndpoints(l,t(1),t(2));

s = [e1;e2];

if nargout > 1
    S_l = [E1_l;E2_l];
    if nargout > 2
        S_t = [E1_t1 E1_t2
            E2_t1;E2_t2];
    end
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb