text:
function [pvec, pstruct] = tapas_hgf_jget_transp(r, ptrans)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2013 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;

l = r.c_prc.n_levels;

pvec(1:l)         = ptrans(1:l);                                  % mux_0
pstruct.mux_0     = pvec(1:l);
pvec(l+1:2*l)     = exp(ptrans(l+1:2*l));                         % sax_0
pstruct.sax_0     = pvec(l+1:2*l);
pvec(2*l+1:3*l)   = ptrans(2*l+1:3*l);                            % mua_0
pstruct.mua_0     = pvec(2*l+1:3*l);
pvec(3*l+1:4*l)   = exp(ptrans(3*l+1:4*l));                       % saa_0
pstruct.saa_0     = pvec(3*l+1:4*l);
pvec(4*l+1)       = exp(ptrans(4*l+1));                           % kau
pstruct.kau       = pvec(4*l+1);
pvec(4*l+2:5*l)   = exp(ptrans(4*l+2:5*l));                       % kax
pstruct.kax       = pvec(4*l+2:5*l);
pvec(5*l+1:6*l-1) = exp(ptrans(5*l+1:6*l-1));                     % kaa
pstruct.kaa       = pvec(5*l+1:6*l-1);
pvec(6*l)         = ptrans(6*l);                                  % omu
pstruct.omu       = pvec(6*l);
pvec(6*l+1:7*l)   = ptrans(6*l+1:7*l);                            % omx
pstruct.omx       = pvec(6*l+1:7*l);
pvec(7*l+1:8*l)   = ptrans(7*l+1:8*l);                            % oma
pstruct.oma       = pvec(7*l+1:8*l);

return;

author: translationalneuromodeling

repo: tapas


text:
% Copyright 2009 - 2010 The MathWorks, Inc.
% Loading the test data
% Note that we are tracking an object whose positions are recorded
% as x and y coordinates in a Cartesian space
load position.mat
hold; grid;
for idx = 1: numPts
z = position(:,idx);
plot(z(1), z(2), 'bx');
axis([-1 1 -1 1]);
end
hold;
%% This algorithm predicts the position of a moving object based on
% its past values. It uses a Kalman filter estimator.
% Kalman filter is a recursive adaptive filter that estimates
% the state of a dynamic system from a series of
% noisy measurements. It has a broad range of application in
% signal and imge processing, control design,
% and computational finance
% Let's run the script ObjTrack.m to see how Kalman filter
% gradually estimates and tracks  a moving object.
ObjTrack01(position)
%% In this example, by assuming a zero initial condition
% ( no knowledge of the position of the object in the past)
% Kalman filter iteratively and adaptively estimates and tracks
% the object. This is depicted by the green circular traces
% as outputs  of the Kalman estimator superimposed
% on the input positions vectors depicted by blue x traces.
% Note that there are 3 sudden shifts in position and every time
% Kalman filter re-adjusts and tracks the object after few iterations
edit ObjTrack01
%% Let's look at the Kalman filter function
% Note that the Kalman estimator computes the position vector
% by computing and updating the Kalman state vector.
% We define the state vector as a 6-by-1 column vector
% that includes position (x and y) velocity (Vx Vy) and
% acceleration (Ax and Ay) measurements in a 2-dimensional
% Cartesian space. Based on the classical laws of motion:
% X = X0 + Vx*t
% Y = Y0 + Vy*t
% Vx = Vx0 + Ax*t
% Vy = Yy0 + Ay*t
% Note how these laws of motion are reflected in the Kalman
% state transion matrix A:
% x=[X Y Vx Vy Ax Ay]' %% state vector
% x=A*x;               %% recursive estimation
% dt=1;
% A=[ 1 0 dt 0 0 0;... %% state transion matrix
%     0 1 0 dt 0 0;...
%     0 0 1 0 dt 0;...
%     0 0 0 1 0 dt;...
%     0 0 0 0 1 0 ;...
%     0 0 0 0 0 1 ];
% Note how by writing only 10 lines of MATLAB code,
% we implement the Kalman estimator as if you are writing
% the mathematical formula. This is power of MATLAB!
%% Version 1 Scalar Kalman filter
% input: one sample
% output: one sample
edit kalman01.m
%% Let us now genearte C code for the kalman filter
% First we run the emlmex command to see weather there is any
% function or operator in this function that is not compliant with the
% Embedded MATLAB subset
z = position(:,1);
emlmex -eg {z} kalman01.m
%% Now that the kalman filter is compliant, Is the
% calling function (ObjTrack01) compliant with Embedded MATLAB subset?
emlmex -eg {position} ObjTrack01.m
%% No,it contains unsupported visualization function
% How do we make it complaint with the embedded MATLAB subset?
% Notice eror messages when complaining due to the use of plotting
% functions. We now showcase use of eml.extrinsic that declare
% visualization functions to be extrinsic and instruct Embedded MATLAB
% not to compile them.
edit ObjTrack02.m
%% Notice how the non-supported plotting functions now
% do not take part in compilation and thus no more eror messages
emlmex -eg {position} ObjTrack02.m
%% Since this function kalman.m is consistent with
% Embedded MATLAB subset, let us now use the emlc command to
% automatically genearte C source code for this functioon
emlc -eg {z} kalman01.m -c -T RTW -report
%% Here we note and explain the following features:
% 1. MATLAB is an untyped language, C is strictly typed.
% So we use -eg option to set the (data type, size and complexity) of the
% input variable to the function based on an example variable
% availbale in MATLAB workspace (in this case varaiable name is z).
% Based on the properties of this variable,
% Embeded MATLAB sets and infers the data type, size and complexity of
% all internal variables to the function, declares them and generate C code.
% 2. -c option indicates we just want to see the C source code
% 3. -T RTW option says generate embeddable C-code and
% 4. -report option says generate a compilation report.
% Now browse and explain the organization of Compilation report
%% Now we want to showcase how in many cases we can
% accelerate the execution speed of a function running a large data set
% by using emlmex command that compiled the function into a MEX function.
% Use of emlmex not only verifies the compliance with the Embedded MATLAB
% subset but it may potentialy accelrate the execution speed of some algorithms
% specially if the algorithm is fixed-point
%% First run the kalman algorithm with 100,000 samples without
% using emlmex to compile it to a mex function
delete *.mexw32
which kalman01
%% Use tic; toc; to record how much time it takes for processing
% 100K samples of data in this case
SIZE=100000;
Z=rand(2,SIZE);
tic;
for i=1:SIZE
    z=Z(:,i);
    y=kalman01(z);
end;
toc;
%% Second run the kalman algorithm with 100,000 samples WITH
% using emlmex to compile it to a mex function
emlmex -eg {z} kalman01.m
which kalman01
%% Use tic; toc; to record how much time it takes for processing 100K
% samples
SIZE=100000;
Z=rand(2,SIZE);
tic;
for i=1:SIZE
    z=Z(:,i);
    y=kalman01(z);
end;
toc;
%% Notice a 3 times acceleration of the speed of the execution of this
% function just by using the emlmex function
%% Version 3 Variable size Vector (Packet based)
% Kalman filter
% input: a variable-sized vector
% output: sanme size as inoput
edit ObjTrack03.m
%% Variable-size data
% Let us now genearte C code for the packet-based
% kalman filter. First we run the emlmex command to see weather there is any
% function or operator in this function that is not compliant with the
% Embedded MATLAB subset
N=100;
z = position(1:2,1:N);
eg_z=emlcoder.egs(z,[2 N]);
emlmex -eg {eg_z} kalman03.m
%% Let us run the compiled version to make sure
% the generated C code can properly handle packets of diferent sizes.
% By examining ObjTrack03.m we see that we are calling kalman03.m with
% 7 different (x-y) pairs of varying sizes ranging from 2x10 to 2x100
which kalman03
ObjTrack03
%% Since this function kalman.m is consistent with
% Embedded MATLAB subset, let us now use the emlc command to
% automatically genearte C source code for this function
% We use the maximum size of 2x100 as the upper bound for compilation to C
% and use the emlcoder.egs to indicate the maximum size of the input
% variable z.
N=100;
z = position(1:2,1:N);
eg_z=emlcoder.egs(z,[2 N]);
emlc -eg {eg_z} kalman03.m -c -T RTW -report
%% Here we note and explain the following features:
% 1. The genearted C code has just enough extra code to make it flexible
% for any size input from 2x1 to 2x100.
% 2. We showcase emlcoder.egs and explain how it works
% 3. We explain how easy it is to make an algorithm be flexible relative to
% input size:
% We show that in the function at line 36, we had o explicitly index
% the input variable z (instead of z(:,i) -> z(1:2,i) to make explicit
% which portion of variable size input us being processed at that line
%% Fixed-point data: Look at the fixed-point Kalman estimator function
% It differs from the floating-point algorithm kalman.m
% by declaring the internal variables as fixed-point
% using the fi constructor. The variable declaration (where most changes are) are separated
% from the algorithm part, which remains mostly the same relative to the
% floating point algorithm.
edit kalman04.m
%% To convert algorithm to fixed-point, you need to set numeric type
% of the fixed-point input data and specify a set of global fimath
% properties.
% As of R2009b, you can specify accumulator ad product data types
% of all intermediate variables globally, which makes the task
% of fixed-point programming easier.
%% Look at the syntax of setting the global fimath properties
% Before setting the global fimath we have the default values
resetglobalfimath
fimath
%% Now note that I am using a helper functin myglobalfimath.m
% to globaly set the fixed-point math properties of all variables in MATLAB workspace.
edit myglobalfimath.m
%%
% I want to set a global fimath properties such that SpecifyPrecision Mode is used
% for setting the Sum and Product properties, accumulators and products
% internal variables have a wordlength of 32 bits, with 20 bits for their fractional part,
% 11 bits integer part and 1 bit for sign
myglobalfimath(20);
fimath
%% Notice how it can take a long time to run fixed-point function
% without using emlmex to accelerate execution speed
% In my machine about 7 seconds for 1000 samples
delete *.mexw32
which kalman04
SIZE=1000;
Z=fi(rand(2,SIZE));
tic;
for i=1:SIZE
    z=Z(:,i);
    y=kalman04(z);
end;
toc;
%% Now you can accelerate the execution speed of fixed-point
% and in most cases even floating-point functionsby using the
% emlmex function of Fixed-Point Toolbox. It generates C code automatically
% from the MATLAB function and compiles it into MEX code.
% It will run faster. Only catch is that the kalman_fixpt.m
% must use the subset of MATLAB language allowing automatic
% MATLAB to C translation. The Embedded MATLAB subset
z = Z(:,1);
emlmex -eg {z} kalman04.m
%% Run the same function and notice the acceleration
% In my machine it run around 0.3 seconds a 20x acceleration
myglobalfimath(20);
which kalman04
SIZE=1000;
Z=fi(rand(2,SIZE));
tic;
for i=1:SIZE
    z=Z(:,i);
    y=kalman04(z);
end;
toc;
%% Now generate C source code from your fixed-point Kalman filter
% by uising the emlc command of the Real Time Workshop
% When you use the -report option it creates a link to an HTML report
% showing the automatically generated C code as kalman_fixpt.c
% Note that all variables are statically defined as native integers and
% all scaling functions are defined before the main C function
emlc -eg {z} kalman04.m -c -T RTW -report
%% Clean up
delete *.mexw32
clear mex

author: Sable

repo: mcbench-benchmarks


text:
function [Stat] = VectStat (M)
%
%   [Stat] = VectStat (M)
%
%Purpose:
%   calculates very basic stats of vectors
%
%
%Input Parameters:
%   M : MxN matrix
%      NAN or Inf values are not considered and the user is warned of their existence
%
%Output Parameters:
%
%   Stat is an Mx1 vector of structures with the following fields
%     Each Stat(i) has the following structures
%      .N : the number of elements in M(i,:)
%      .M : the mean of vector M(i,:)
%      .Md: the median of vector M(i,:)
%      .S : the (unbiased) standard deviation of M(i,:)
%      .V : the variance (S^2) of M(i,:)
%      .mx: The maximum value
%      .mn: The minimum value
%
%     if the input vector is empty , NaN values are returned in Stat
%More Info :
%
%
%
%
%     Author : Ziad Saad
%     Date : Tue Apr 20 17:13:23 CDT 1999

%Define the function name for easy referencing
FuncName = 'VectStat';

%Debug Flag
DBG = 1;
Stat.N = NaN;
Stat.M = NaN;
Stat.Md = NaN;
Stat.S = NaN;
Stat.V = NaN;

if (is_row(M) ~= -1),
	sZ = 1;
	sz2 = length(M);
	M = M(:)'; %turn M into a row vector
else
	sZ = size(M,1);
	sz2 = size(M,2);
end

if (sZ == 0 | sz2 ==0),
	ErrEval(FuncName,'Wrn_Empty M returning NaN structure for result');
	return;
end

%initialize, once
	Stat(sZ).N = 0;Stat(sZ).M = 0;Stat(sZ).Md = 0;Stat(sZ).S = 0;Stat(sZ).V = 0; Stat(sZ).mn = 0; Stat(sZ).mx = 0;

for (i=1:1:sZ),
	igood = find(isfinite(M(i,:)));
	N_igood = length(igood);
	if (N_igood < sz2),
		stmp = sprintf ('%s Warning: Inf or Nan ignored in row %g\n', FuncName, i);
	end

	Stat(i).N = N_igood;
	Stat(i).M = mean(M(i,igood));
	Stat(i).Md = median(M(i,igood));
	Stat(i).S = std(M(i,igood),0); %unbiased estimator
	Stat(i).V = Stat(i).S .^ 2;
	Stat(sZ).mn = min(M(i,igood));
	Stat(sZ).mx = max(M(i,igood));
end

return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function r = dirrand(m,n)
% DIRRAND - Uniform Dirichlet random vectors
%
%    R = DIRRAND(M) returns vector of length M chosen from
%    Dirichlet(1,...,1) distribution.
%    R = DIRRAND(M,N) returns N such vectors in MxN matrix.
%
%   References: Gelman, Carlin, Stern, Dunson, Vehtari & Rubin (2013),
%               Bayesian Data Analysis, third edition, p. 583.
%
% Copyright (c) 1998-2004 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if nargin < 2
  n=1;
end

% Generate random numbers from Gamma(1,1) distribution...
r=rand(m,n);
r=-reallog(r);
% ...and scale sum to unity
rs=sum(r);
for i1=1:m
  r(i1,:)=r(i1,:)./rs;
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function y = horzcat(varargin)
    % Horizontal concatenation for sym arrays.
    %   C = HORZCAT(A, B, ...) horizontally concatenates the symbolic
    %   variable arrays A, B, ... .  For matrices, all inputs must have the
    %   same number of rows.
    %
    %   C = HORZCAT(A,B) is called for the syntax [A B].
    %
    %   See also VERTCAT.

    % Call the superclass method
    str = horzcat@SymExpression(varargin{:});

    % create a new object with the evaluated string
    y = SymVariable(str);

end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function M = dunzip(Z)
% DUNZIP - decompress DZIP output to recover original data
%
% USAGE:
% M = dzip(Z)
%
% VARIABLES:
% Z = compressed variable to decompress
% M = decompressed output
%
% NOTES: (1) The input variable Z is created by the DZIP function and
%            is a vector of type uint8
%        (2) The decompressed output will have the same data type and
%            dimensions as the original data provided to DZIP.
%        (3) See DZIP for other notes.
%        (4) Carefully tested, but no warranty; use at your own risk.
%        (5) Michael Kleder, Nov 2005

import com.mathworks.mlwidgets.io.InterruptibleStreamCopier
a=java.io.ByteArrayInputStream(Z);
b=java.util.zip.InflaterInputStream(a);
isc = InterruptibleStreamCopier.getInterruptibleStreamCopier;
c = java.io.ByteArrayOutputStream;
isc.copyStream(b,c);
Q=typecast(c.toByteArray,'uint8');
cn = double(Q(1)); % class
nd = double(Q(2)); % # dims
s = typecast(Q(3:8*nd+2),'double')'; % size
Q=Q(8*nd+3:end);
if cn == 3
    M  = logical(Q);
elseif cn == 4
    M = char(Q);
else
    ct = {'double','single','logical','char','int8','uint8',...
        'int16','uint16','int32','uint32','int64','uint64'};
    M = typecast(Q,ct{cn});
end
M=reshape(M,s);
return

author: dsmi-lab-ntust

repo: AnomalyDetectionToolbox


text:
%sim_COF.m
%Dana Goerzen and Jamie Near, 2021.
%
% USAGE:
%  d_out=sim_COF(H,d_in,order)
%
% DESCRIPTION:
% This function nulls signal from any undesired coherences in a spin system
% Desired coherences are determined through extended phase graph analysis
% and pulse sequence design.
%
% INPUTS:
% H         = Hamiltonian operator structure.
% d_in      = input density matrix structure
% order     = desired coherence order that you wish to keep signal from
%
% OUTPUTS:
% d_out     = output density matrix structure

function d_out=sim_COF(H,d_in,order)
%initialize mask as permitting all coherences through, then iterate through coherence matrix in Hamiltonian to
%set any values that don't correspond to the desired coherence order to 0.
for n=1:length(H) %JN - Looping through the parts of the spin system:
    mask1=H(n).coherenceOrder==order;
    %zero any undesired coherences
    d_temp=mask1.*d_in{n};
    d_in{n}=d_temp;
end
d_out=d_in;
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
%%
% Calculates b based on time using Mc and plots it - called from
% cumulative window
%%

report_this_filefun(mfilename('fullpath'));

inter = 50;
incr = 10;

%%
% set step_cat as a dummy variable so that newt2 can be reassigned
% for use in mcperc_ca3.  newt2 is reset at end.
%%

step_cat = newt2;
bv2 = [];
bv3 = [] ;
me = [];
bvm = [];
ctr = 0;
i=1;
Nmin=50;
%% b=newt2;  %% initialization for bvalca3.m
day_start=0;
day_end=0;

inpr1 = 5;

nibt = 200;
win_step = 200;
for ind = 1:win_step:length(step_cat)-win_step
    newt2 = step_cat(ind:ind+nibt,:);

    %%
    % calculation based on best combination of 90% and 95% probability -- default
    %%

    %         elseif inpr1 == 5
    mcperc_ca3;
    if isnan(Mc95) == 0
        magco = Mc95;
    elseif isnan(Mc90) == 0
        magco = Mc90;
    else
        [bv magco stan av me mer me2,  pr] =  bvalca3(newt2,1,1);
    end
    l = newt2.Magnitude >= magco-0.05;
    if length(newt2(l,6)) <= Nmin
        %                disp(['%%Warning --bwithtimc--%%  less than 50 events in step ']);
    end
    [ntl,ntw] = size(newt2(l,:));
    if ntl <= 1
        disp('%%ERROR --bwithtimc--%%  Not enough data to plot');
        newt2 = step_cat;
        return
    end
    [mea bv stand,  av] = bmemag(newt2(l,:));

    days = (max(newt2.Date)-min(newt2.Date))*365.0;

    bvm = [bvm; bv step_cat(ind,3) step_cat(ind+nibt,3) magco days ind ind+nibt av];
end  %% end of for loop!!
%end

[bvml,bvmw] = size(bvm);
if bvml <= 1
    disp('%%ERROR --bwithtimc--%%  Not enough data to calculate Mc');
    newt2 = step_cat;
    return
end
if ind + win_step < length(step_cat)
    newt2 = step_cat(ind+win_step+1:length(step_cat),:);
    %%
    % calculation based on best combination of 90% and 95% probability -- default
    %%

    if inpr1 == 5
        mcperc_ca3;
        if isnan(Mc95) == 0
            magco = Mc95;
        elseif isnan(Mc90) == 0
            magco = Mc90;
        else
            [bv magco stan av me mer me2,  pr] =  bvalca3(newt2,1,1);
        end
        l = newt2.Magnitude >= magco-0.05;
        if length(newt2(l,6)) <= Nmin
            %                disp(['%%Warning --bwithtimc--%%  less than 50 events in step ']);
        end
        [ntl,ntw] = size(newt2(l,:));
        if ntl <= 1
            disp('%%ERROR --bwithtimc--%%  Not enough data to plot');
            newt2 = step_cat;
            return
        end
        [mea bv stand,  av] = bmemag(newt2(l,:));

    end
end
bvm = [bvm; bv step_cat(ind+win_step,3) step_cat(length(step_cat),3) magco days ind ind+nibt av];

%bvm(length(bvm)+1,:) = bvm(ind,:)
newt2 = step_cat;

author: CelsoReyes


text:
function [decomp] = Decomposite_eig_new(M)

% input:  X is a cell structure, each containing a PSD matrix;
%         Y is a cell structure, each containing a PSD matrix;
%         index decides the set of eigenvalues to be included;
%
dim = size(M.X_cov,1);
n_X = size(M.X_cov,3);
if(n_X == 0)
    disp('X is empty, error!');
    return;
end

% extract eigenvalues and check PSD for input;

V_x = zeros(dim,dim,n_X);
D_x = zeros(dim,dim,n_X);

for i = 1:n_X
    [V,D] = eig(M.X_cov(:,:,i));
    Sorted_D = sort(diag(D),'descend');

    V_x(:,:,i) = fliplr(V(:,:));
    D_x(:,:,i) = abs(diag(Sorted_D));
end

decomp.V = V_x;
decomp.D = D_x;

end

author: hiroyuki-kasai

repo: ClassifierToolbox


text:
% readegilocs() - look up locations for EGI EEG dataset.
%
% Usage:
%   >> EEG = readegilocs(EEG);
%   >> EEG = readegilocs(EEG, fileloc);
%
% Inputs:
%   EEG        - EEGLAB data structure
%
% Optional input:
%   fileloc    - EGI channel location file
%
% Outputs:
%   EEG        - EEGLAB data structure with channel location structure
%
% Author: Arnaud Delorme, SCCN/UCSD
%
% See also: pop_readegi(), readegi()

% Copyright (C) 12 Nov 2002 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function EEG = readegilocs(EEG, fileloc);

if nargin < 1
    help readegilocs;
    return;
end;

% importing channel locations
% ---------------------------
found = 1;
if nargin < 2
    switch EEG.nbchan
        case { 32 33 }, fileloc = 'GSN-HydroCel-32.sfp';
        case { 64 65 }, fileloc = 'GSN65v2_0.sfp';
        case { 128 129 }, fileloc = 'GSN129.sfp';
        case { 256 257 }, fileloc = 'GSN-HydroCel-257.sfp';
        otherwise, found = 0;
    end;
end;
if found
    fprintf('EGI channel location automatically detected %s ********* WARNING please check that this the proper file\n', fileloc);
    if EEG.nbchan == 64 || EEG.nbchan == 65 || EEG.nbchan == 256 || EEG.nbchan == 257
        fprintf( [ 'Warning: this function assumes you have a 64-channel system Version 2\n' ...
                   '         if this is not the case, update the channel location with the proper file' ]);
    end;
    % remove the last channel for 33 channels

    peeglab = fileparts(which('eeglab.m'));
    fileloc = fullfile(peeglab, 'sample_locs', fileloc);
    locs = readlocs(fileloc);
    locs(1).type = 'FID';
    locs(2).type = 'FID';
    locs(3).type = 'FID';
    locs(end).type = 'REF';
    if EEG.nbchan == 256 || EEG.nbchan == 257
        if EEG.nbchan == 256
            chaninfo.nodatchans = locs([end]);
            locs([end]) = [];
        end;
    elseif mod(EEG.nbchan,2) == 0,
        chaninfo.nodatchans = locs([1 2 3 end]);
        locs([1 2 3 end]) = [];
    else
        chaninfo.nodatchans = locs([1 2 3]);
        locs([1 2 3]) = [];
    end; % remove reference
    chaninfo.filename = fileloc;
    EEG.chanlocs   = locs;
    EEG.urchanlocs = locs;
    EEG.chaninfo   = chaninfo;
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function events = in_events_array(sFile, EventFile, format, EventName, isInteractive)
% IN_EVENTS_ARRAY: Read events information from a .mat or text file
%
% USAGE:  events = in_events_array(sFile, EventFile, 'times',   EventName=[ask], isInteractive=1)
%         events = in_events_array(sFile, EventFile, 'samples', EventName=[ask], isInteractive=1)
%         events = in_events_array(sFile, EventMat, ...)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, Elizabeth Bock, 2012-2023

% Parse inputs
if (nargin < 5) || isempty(isInteractive)
    isInteractive = 1;
end
if (nargin < 4) || isempty(EventName)
    EventName = [];
end

% ===== READ FILE =====
if ischar(EventFile)
    % Can read ASCII or .mat
    EventsMat = load(EventFile);
    if isstruct(EventsMat)
        fields = fieldnames(EventsMat);
        EventsMat = EventsMat.(fields{1});
        if isstruct(EventsMat)
            fields = fieldnames(EventsMat);
            EventsMat = EventsMat.(fields{1});
        end
    end
else
    EventsMat = EventFile;
end
% Force to be double
EventsMat = double(EventsMat);
% Check orientation
if (size(EventsMat,2) < size(EventsMat,1))
    EventsMat = EventsMat';
end

% ===== FORMAT EVENTS =====
% Create [samples; times] array
switch (format)
    case 'times'
        evtTimes = round(EventsMat .* sFile.prop.sfreq) ./ sFile.prop.sfreq;
    case 'samples'
        evtTimes = round(EventsMat) ./ sFile.prop.sfreq;
end

% ===== TIME OFFSET =====
if isInteractive
    % Check for offset (typical of FIF files)
    isAddOffset = 0;
    if (sFile.prop.times(1) ~= 0)
        res = java_dialog('question', ['The raw data file starts at ' num2str(sFile.prop.times(1)) ' sec.' 10 10 ...
                                      'Is this offset already added to these events?' 10 10],...
                                      'Import events', [], {'Yes', 'Add Offset','Cancel'},'Yes');
        if isempty(res) || strcmpi(res, 'Cancel')
            bst_progress('stop');
            return;
        elseif strcmpi(res, 'Add Offset')
            isAddOffset = 1;
        end
    end
    % Add a column for time
    if isAddOffset
        evtTimes = evtTimes + sFile.prop.times(1);
    end
end

% ===== CONVERT TO BRAINSTORM STRUCTURE =====
% Initialize list of events
events = db_template('event');
% Ask for a label
if ~isempty(EventName)
    events.label = EventName;
else
    res = java_dialog('input', 'Please enter a label for this event:', 'Event Label');
    if isempty(res)
        events.label = '1';
    else
        events.label = res;
    end
end
events.color      = [];
events.reactTimes = [];
events.select     = 1;
events.times      = evtTimes;
events.epochs     = ones(1, length(evtTimes));  % Epoch: set as 1 for all the occurrences
events.channels   = [];
events.notes      = [];

author: brainstorm-tools


text:
%................................................................

function outputDisplacements...
    (displacements,numberNodes,GDof)

% output of displacements in tabular form

disp('Displacements')
jj=1:numberNodes; format
f=[jj; displacements(1:2:GDof)'; displacements(2:2:GDof)'];
fprintf('Node           UX           UY\n')
fprintf('%4d   %10.4e    %10.4e\n',f)

repo: Aero-matlab


text:
function stroud_test322 ( )

%*****************************************************************************80
%
%% TEST322 tests SPHERE_CAP_AREA_3D, SPHERE_CAP_AREA_ND.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  dim_num = 3;
  ntest = 12;
  r = 1.0;
  center = [ 0.0, 0.0, 0.0 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST322\n' );
  fprintf ( 1, '  SPHERE_CAP_AREA_3D computes the volume of a\n' );
  fprintf ( 1, '    3D spherical cap, defined by a plane that cuts the\n' );
  fprintf ( 1, '    sphere to a thickness of H units.\n' );
  fprintf ( 1, '  SPHERE_CAP_AREA_ND computes the volume of an\n' );
  fprintf ( 1, '    ND spherical cap, defined by a plane that cuts the\n' );
  fprintf ( 1, '    sphere to a thickness of H units.\n' );

  area1 = sphere_area_3d ( r );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Area of the total sphere in 3D = %f\n', area1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '        R           H           Cap         Cap\n' );
  fprintf ( 1, '                                area_3d     area_nd\n' );
  fprintf ( 1, '\n' );

  for i = 0 : ntest + 1

    h = 2.0 * r * i / ntest;

    area1 = sphere_cap_area_3d ( r, h );

    area2 = sphere_cap_area_nd ( dim_num, r, h );

    fprintf ( 1, '  %12.6f  %12.6f  %12,6f  %12.6f\n', r, h, area1, area2 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = plotEarth(varargin)
%PLOTEARTH  Plots the landmasses of earth.
%   PLOTEARTH plots outlines of the landmasses of earth on a sphere using
%   black solid lines.
%
%   PLOTEARTH(LINESPEC) uses the uses the color and marker from the
%   line specification string 'LineSpec' (See PLOT for possibilities)
%
% See also spherefun/surf, spherefun/plot, spherefun/contour

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

% Land masses are stored in the data file CoastData.mat
try
    x = load('CoastData.mat', 'coast');
catch
    error('CHEBFUN:SPHEREFUN:PLOTEARTH:coastDataNotFound',...
        ['File containing the coast line data could not be found. '...
        'Try reinstalling chebfun.']);
end

if ( nargin > 0 )
    linespec = varargin{:};
else
    linespec = 'k-';
end

% Get the hold state of the current axis:
holdState = ishold;

if ( ~holdState )
    % Generate a unit sphere.
    [XX,YY,ZZ] = sphere(101);

    % Color of the sphere will be white:
    clr = [255 255 255]/255;

    % Plot the sphere, make it slightly smaller than unit so lines
    % show up more clearly.
    scl = 0.99;

    surf(scl*XX, scl*YY, scl*ZZ, 1+0*XX, 'EdgeColor', 'None', ...
        'FaceColor', clr);
    daspect([1 1 1]);
    hold on
end

h = plot3(x.coast(:,1), x.coast(:,2), x.coast(:,3), linespec);

if ( ~holdState )
    hold off
end

if ( nargout > 0 )
    varargout = { h };
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [train_msk, p] = Masks(train, train_lbl, p)
%% Matrices easying the vectorized computation in CostFunction

%% Two types of masks: one is to mask with -Inf values that are superfluous due to sentence length overhead
%% the second is to apply differentiated k-max pooling that depends on the length of each sentene; turning superfluous values to max length of first layer

%Lengths of mask independenct of tot_num_layers
p(25) = p(2)+p(4)-1; %len of masks length+pool layer 1
p(26) = indexMaskLen(p(2), 1, p) + p(6)-1; %len of masks length+pool layer 2
p(27) = indexMaskLen(p(2), 2, p) + p(36)-1; %len of masks length+pool layer 3

if p(10) == 1 %If one layer in total
    train_msk=zeros(size(train,1), 2*p(25));

    train_msk = findMasks(train, train_lbl, train_msk, p);
elseif p(10) == 2
    train_msk=zeros(size(train,1), 2*p(25)+2*p(26));

    train_msk = findMasks(train, train_lbl, train_msk, p);

elseif p(10) == 3
    train_msk=zeros(size(train,1), 2*p(25)+2*p(26)+2*p(27));
    train_msk = findMasks(train, train_lbl, train_msk, p);

end

train_msk = logical(train_msk);

end

function train_msk = findMasks(train, train_lbl, train_msk ,p)
for i=1:length(train)
    if p(60)==1
        sent_len = train_lbl(i,2); %contains length
    else
        sent_len = train_lbl{2}(i); %contains length
    end

    if p(10) == 1
        train_msk(i,sent_len+p(4)-1+1:p(25)) = 1; %layer 1 sent_lengths mask
        train_msk(i,p(25)+p(7)+1:end) = 1; %layer 1 max_pool mask
    elseif p(10) == 2
        train_msk(i,sent_len+p(4)-1+1:p(25)) = 1; %layer 1 sent_lengths mask
        train_msk(i,p(25)+indexMaskLen(sent_len,1,p)+1:2*p(25)) = 1; %layer 1 max_pool mask

        train_msk(i,2*p(25)+indexMaskLen(sent_len,1,p)+p(6)-1+1:2*p(25)+p(26)) = 1; %layer 2 sent_lengths mask
        train_msk(i,2*p(25)+p(26)+p(7)+1:end) = 1; %layer 2 max_pool mask

    elseif p(10) == 3
        train_msk(i,sent_len+p(4)-1+1:p(25)) = 1; %layer 1 sent_lengths mask
        train_msk(i,p(25)+indexMaskLen(sent_len,1,p)+1:2*p(25)) = 1; %layer 1 max_pool mask

        train_msk(i,2*p(25)+indexMaskLen(sent_len,1,p)+p(6)-1+1:2*p(25)+p(26)) = 1; %layer 2 sent_lengths mask
        train_msk(i,2*p(25)+p(26)+indexMaskLen(sent_len,2,p)+1:2*p(25)+2*p(26)) = 1; %layer 2 max_pool mask

        train_msk(i,2*p(25)+2*p(26)+indexMaskLen(sent_len,2,p)+p(36)-1+1:2*p(25)+2*p(26)+p(27)) = 1; %layer 3 sent_lengths mask
        train_msk(i,2*p(25)+2*p(26)+p(27)+p(7)+1:end) = 1; %layer 3 max_pool mask

    end
end

end

function pool_size = indexMaskLen(sent_len, num_conv_layer, p)
%Function computes pooling size for a given layer
if p(10) == 1
    pool_size = p(7); %size of final pooling layer

elseif p(10) == 2
    if num_conv_layer == 1
        pool_size = max(ceil(sent_len/2),p(7)); %half if shallow network
    end

    if num_conv_layer == 2 %not in fact possible
        pool_size = p(7);
    end
elseif p(10) == 3
    if num_conv_layer == 1
        pool_size = max(ceil((sent_len/3)*2),p(7)); %linear function
    end

    if num_conv_layer == 2
        pool_size = max(ceil(sent_len/3),p(7));
    end
else
    pool_size = -1;
end

end

author: jacoxu


text:
%% Copyright (C) 2014-2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym gamma (@var{x})
%% Symbolic gamma function.
%%
%% Example:
%% @example
%% @group
%% syms x
%% y = gamma (x)
%%   @result{} y = (sym) Γ(x)
%% @end group
%% @end example
%%
%% Note: this file is autogenerated: if you want to edit it, you might
%% want to make changes to 'generate_functions.py' instead.
%%
%% @end defmethod

function y = gamma(x)
  if (nargin ~= 1)
    print_usage ();
  end
  y = elementwise_op ('gamma', x);
end

%!error gamma (sym(1), 2)
%!assert (isequaln (gamma (sym(nan)), sym(nan)))

%!shared x, d
%! d = 1;
%! x = sym('1');

%!test
%! f1 = gamma(x);
%! f2 = gamma(d);
%! assert( abs(double(f1) - f2) < 1e-15 )

%!test
%! D = [d d; d d];
%! A = [x x; x x];
%! f1 = gamma(A);
%! f2 = gamma(D);
%! assert( all(all( abs(double(f1) - f2) < 1e-15 )))

%!test
%! % round trip
%! y = sym('y');
%! A = gamma (d);
%! f = gamma (y);
%! h = function_handle (f);
%! B = h (d);
%! assert (A, B, -eps)

repo: octsympy


text:
function [voxels,scale,translate] = read_binvox(filename)
         fid = fopen(filename,'r');
         % read header
         tline = fscanf(fid,'%s',1);
         if ~strcmp(tline,'#binvox')
             fprintf(['Error: first line reads [' tline '] instead of [#binvox]\n']);
         end
         version = fscanf(fid,'%d',1);
%          fprintf('reading binvox version %d\n',version);
         depth = -1;
         done = 0;
         while (~feof(fid)&&~done)
            tline = fscanf(fid,'%s',1);
            if strcmp(tline,'data')
                done = 1;
            elseif strcmp(tline,'dim')
                A=fscanf(fid,'%d',3);
                width = A(1);
                height= A(2);
                depth = A(3);
            elseif strcmp(tline,'translate')
                translate = fscanf(fid,'%f',3);
            elseif strcmp(tline,'scale')
                scale =fscanf(fid,'%f',1);
            end
         end

         scale = scale/width;
         totalsize = width * height * depth;
         voxels = zeros(width, height, depth);
         % read voxel data
         index = 0;end_index = 0;nr_voxels = 0;

         B = fread(fid,1,'uint8'); %read the linefeed char
         while((end_index < totalsize) && ~feof(fid))
              B = fread(fid,2,'uint8');
              value = B(1);
              count = B(2);
              end_index = index + count;

              if (end_index >=totalsize)
                  break;
              end

              for i=index:min(totalsize,end_index)
                  voxels(i+1) = value; % matlab index starts 1
              end

              if (value)
                  nr_voxels = nr_voxels+count;
              end
              index = end_index;
         end
         % matlab's strange ordering
         voxels = permute(voxels,[3,1,2]);
         fclose(fid);
end

author: jhonykaesemodel


text:
function z = r2z(R)
% Transformation of a correlation matrix to an unconstrained K(K-1)/2 vector
%
% USAGE:
%  [R] = r2z(Z)
%
% INPUTS:
%   R - A K by K correlation matrix
%
% OUTPUTS:
%   Z - K(K-1)/2 vector of values in (-inf,inf)
%
% COMMENTS:
%   See z2r for information about the transformation from Z to R.
%
% See also Z2R, R2PHI, PHI2R

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 3/27/2012

k = length(R);
C = zeros(k);
C2 = chol(R)';

for i=2:k
    rem = 1;
    for j=i-1:-1:1
        C(i,j) = C2(i,j)/sqrt(rem);
        rem = rem - C2(i,j)^2;
    end
end

C =C';
z = C(~tril(true(k)));
z = log((z+1)./(1-z));

author: bashtage

repo: mfe-toolbox


text:
classdef nncaffepool < nntest

  properties
    x
  end

  properties (TestParameter)
    type = {'avg', 'max'}
    poolx = {1 2 3}
    pooly = {1 2}
    pool = {1 2 3}
    pad = {0 1 2}
    stride = {1 2 3 4}
    stridex = {1 2 3}
    stridey = {1 2}
    padLeft = {0 1 2}
    padRight = {0 1 2}
    padTop = {0 1 2}
    padBottom = {0 1 2}
  end

  methods (TestClassSetup)
    function data(test,device)
      % make sure that all elements in x are different. in this way,
      % we can compute numerical derivatives reliably by adding a delta < .5.
      x = test.randn(15,14,3,2) ;
      x(:) = randperm(numel(x))' ;
      test.x = x ;
      test.range = 10 ;
      if strcmp(device,'gpu'), test.x = gpuArray(test.x) ; end
    end
  end

  methods (Test)
    function basic(test,poolx,pooly)
      % Test whether the avg pool output is equal to its emulation with
      % convolutional layer
      x = test.x ;
      stride = 1 ;
      pad = 0 ;
      pool = [pooly poolx] ;
      args = {'stride',stride,'pad',pad, 'method', 'avg'};
      y = vl_nncaffepool(x,pool,args{:}) ;
      y_conv = vl_nnconv(gather(x), ...
                         ones(pooly,poolx,1,size(x,3),test.currentDataType)./poolx./pooly, ...
                         zeros(1,size(x,3),test.currentDataType), ...
                         'stride', stride, ...
                         'pad', pad);
      test.eq(y, y_conv, 1e-3); % Does not pass with 1e-4
    end

    function pool_type(test, type, pool, pad, stride)
      x = test.x ;
      if pad > pool-1, return ; end
      args = {'stride',stride,'pad',pad,'method',type};
      y = vl_nncaffepool(x,pool,args{:}) ;
      dzdy = test.randn(size(y)) ;
      dzdx = vl_nncaffepool(x,pool,dzdy,args{:}) ;
      test.der(@(x) vl_nncaffepool(x,pool,args{:}), ...
               x, dzdy, dzdx, test.range * 1e-2) ;
    end

    function pool_type_and_pad(test, type, poolx, pooly, stride)
      x = test.x ;
      stride = 1 ;
      pad = 0 ;
      pool = [pooly poolx] ;
      args = {'stride',stride,'pad',pad,'method',type};
      y = vl_nncaffepool(x,pool,args{:}) ;
      dzdy = test.randn(size(y)) ;
      dzdx = vl_nncaffepool(x,pool,dzdy,args{:}) ;
      test.der(@(x) vl_nncaffepool(x,pool,args{:}), ...
               x, dzdy, dzdx, test.range * 1e-2) ;
    end

    function pool_type_and_stride(test, type, stridex, stridey)
      x = test.x ;
      pad = 0 ;
      pool = [3 2] ;
      stride = [stridey stridex] ;
      args = {'stride',stride,'pad',pad,'method',type};
      y = vl_nncaffepool(x,pool,args{:}) ;
      dzdy = test.randn(size(y)) ;
      dzdx = vl_nncaffepool(x,pool,dzdy,args{:}) ;
      test.der(@(x) vl_nncaffepool(x,pool,args{:}), ...
               x, dzdy, dzdx, test.range * 1e-2) ;
    end

    function asym_pad1(test, type, padLeft, padRight)
      x = test.x ;
      pool = [3 4] ;
      stride = [2 1] ;
      pad = [0 0 padLeft padRight] ;
      args = {'stride',stride,'pad',pad,'method',type};
      y = vl_nncaffepool(x,pool,args{:}) ;
      dzdy = test.randn(size(y)) ;
      dzdx = vl_nncaffepool(x,pool,dzdy,args{:}) ;
      test.der(@(x) vl_nncaffepool(x,pool,args{:}), ...
               x, dzdy, dzdx, test.range * 1e-2) ;
    end

    function asym_pad2(test, type, padTop, padBottom)
      x = test.x ;
      pool = [3 4] ;
      stride = [2 1] ;
      pad = [padTop padBottom 2 1] ;
      args = {'stride',stride,'pad',pad,'method',type};
      y = vl_nncaffepool(x,pool,args{:}) ;
      dzdy = test.randn(size(y)) ;
      dzdx = vl_nncaffepool(x,pool,dzdy,args{:}) ;
      test.der(@(x) vl_nncaffepool(x,pool,args{:}), ...
               x, dzdy, dzdx, test.range * 1e-2) ;
    end
  end
end

repo: MFT


text:
function [x1] = simulate_CRP(alpha,n)
% simulate labels from CRP process with concentration alpha
% IN:
%   - alpha: concentration parameter of the CRP
%   - n: size of the time series
% OUT:
%   - x1: CRP-label process (vector of size 1Xn)

% initial condition
x1 = zeros(1,n);
x1(1) = 1;

% loop over time samples
for i=2:n
   Ki = max(x1(1:i-1)); % number of classes so far
   nk = zeros(Ki,1);
   for k=1:Ki
       nk(k) = sum(x1(1:i-1)==k);
   end
   p = [nk;alpha]./(alpha+i-1);
   x1(i) = VBA_random ('Categorical', p);

end

author: MBB-team

repo: VBA-toolbox


text:
function [B,L,var,fac,E] = FA(X,s)
% FA Factor analysis (principal factoring technique)
%
% The reason to favour this method opposite to maximum likelihood
% (Statistics toolbox) is no distributional assumptions. In comparison to
% SPSS program, this provides the same results. We expect identical
% numerical solution, especially advantage of singular value decomposition.
% Initial estimates of communalities are set as squared multiple
% correlations. In iteration process (convergence criterium 0.001, limit
% for iteration 75), communalities are set to one if Heywood case occurs.
% The loadings are rotated (normalized varimax rotation, convergence
% criterium 0.00001, iteration limit 35). Estimate of rotated factors is
% based on regression approach.
%
% [B,L,var,fac,E] = FA(X,s)
%
% Inputs:
%   X is the data matrix (columns as variables, X must have more
%     than one row and more than one column).
%
%   s is the number of extracted factors (if this parameter is
%     not included, number is chosen by principal component
%     criterion with eigenvalues greater than or equal to one).
%
% Results:
%   B is the matrix of factor loadings (unrotated), last column
%     of matrix B indicates an extraction estimate of
%     communalities.
%
%   L is the varimax rotated loadings matrix.
%
%   var describes the variability proportions, last item is
%     the cumulative variance proportion.
%
%   fac is the matrix of factors.
%
%   E is the matrix of specific variances.
%
% Example: From the data containing temperature, relative humidity, wind
% speed, radiation, NO/NO2 ratio and ozone gained in boundary layer
%  [B,L,var] = FA(X,3)
%   The number of factors extracted: 3
%   Factorial number of iterations: 14
%   Rotational number of iterations: 5
%
%     B =
%        -0.962    0.145    0.088    0.954
%         0.964   -0.065   -0.230    0.986
%        -0.551    0.189    0.061    0.343
%        -0.557    0.422   -0.444    0.685
%         0.718    0.649    0.238    0.994
%        -0.975   -0.077    0.081    0.963
%
%     L =
%        -0.850   -0.308   -0.370
%         0.913    0.322    0.220
%        -0.516   -0.087   -0.262
%        -0.291   -0.083   -0.771
%         0.297    0.948    0.087
%        -0.807   -0.495   -0.258
%
%     var =
%         0.441    0.226    0.154    0.821
%
% The biggest values of rotated loadings points dependence of tempetature,
% relative humidity, wind speed and ozone in the first factor and NO/NO2
% ratio with ozone in another.
%
% Communalities show the biggest relation of NO/NO2 ratio to the others. Be
% careful, variances need not be in decreasing order.
%
% For citation ensue following format: Malec L., Skacel F., Trujillo-Ortiz
% A.(2007). FA: Factor analysis by principal factoring. A MATLAB file.
% [WWW document]. URL http://www.mathworks.com/matlabcentral/fileexchange/
% loadFile.do?objectId=14115
%
% Golub G. H., Van Loan C. F.: Matrix Computations. The Johns Hopkins
% University Press, Baltimore 1996.
% Harman H. H.: Modern Factor Analysis. The University of Chicago Press,
% Chicago 1976.
% Krzanowski W. J.: Principles of Multivariate Analysis. Oxford University
% Press, Oxford 2003.
%
% Copyright. February 26, 2007.

[m,n] = size(X);

if nargin==2 && n<s
  error 'The number of factors requested(s) is too large'
end

% Variables standardization and evaluation of correlation matrix.
e = ones(m,1);
X = (X - e * mean(X)) ./ (e * std(X));
R = cov(X);

% Maximization of variance of original variables.
a = svd(R,0);

% Number of factors according roots of principal components.
if nargin<2
   f = a>=1;
   s = sum(f);
end

fprintf ('The number of factors extracted: %.i\n', s);

% Communality estimation by coefficients of multiple correlation.
c = ones(n,1) - 1 ./ diag(R \ diag(ones(n,1)));
g = [];

% Iteration cycle which maximizes factors correlations.
for k = 1:75
    [U,D,V] = svd(R - diag(diag(R)) + diag(c),0);
    N = V * sqrt(D(:,1:s));
    p = c;
    c = sum(N.^2,2);
    g = [g find(c>1)']; % Heywood case in communality estimates.
    c(g) = 1;
    if max(abs(c - p))<0.001
       break
    end
end

fprintf ('Factorial number of iterations: %.i\n', k);

% Evaluation of factor loadings and communalities estimations.
B = [N c];

% Normalization of factor loadings.
h = sqrt(c);
N = N ./ repmat(h,1,s);

% Initial choice of eigenvalues and loadings labelling.
L = N;
z = 0;

% Iteration cycle maximizing variance of individual columns.
for l = 1:35
    [A,S,M] = svd(N' * (n * L.^3 - L * diag(sum(L.^2))),0);
    L = N * A * M';
    b = z;
    z = sum(diag(S));
    if abs(z - b)<0.00001
       break
    end
end

fprintf ('Rotational number of iterations: %.i\n', l);

% Unnormalization of factor loadings.
L = L .* repmat(h,1,s);

% Factors computation by regression and variance proportions.
t = sum(L.^2) / n;
var = [t sum(t)];
fac = X * (R \ diag(ones(n,1))) * L;

% Evaluation of given factor model variance specific matrix.
r = diag(R) - sum(L.^2,2);
E = R - L * L' - diag(r);


text:
function x = p04_x ( n )

%*****************************************************************************80
%
%% P04_X returns the least squares solution X for problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Output, real X(N,1), the least squares solution.
%
  x = [ 1.0; 2.0; 3.0; 4.0; 5.0 ];

  return
end


text:
function [P,C,vP,vC,s] = computeFluxSplits(model,mets,V, coeffSign)
% Compute relative contributions of fluxes (`V`) to the net production (`P`)
% and consumption (`C`) of a set of metabolites (`mets`)
%
% USAGE:
%
%    [P, C, vP, vC, s] = computeFluxSplits(model, mets, V, 0);
%
% INPUTS:
%    model:     a COBRA model structure with required fields
%
%               * .S:        the `m` x `n` stoichiometric matrix
%               * .mets:     an `m` x `1` cell array of metabolite identifiers
%    mets:      a list of metabolite identifiers from model.mets
%    V:         an `n` x `k` matrix of `k` flux distributions
%
% OPTIONAL INPUT:
%    coeffSign: only use the sign of the stoichiometric coefficient
%               (i.e. +1 or -1, Default = false)
%
% OUTPUTS:
%    P:     an `n` x `k` matrix of relative contributions to the production of mets
%    C:     an `n` x `k` matrix of relative contributions to the consumption of mets
%    vP:    an `n` x `k` matrix of net producing fluxes
%    vC:    an `n` x `k` matrix of net consuming fluxes
%    s:     net stoichiometry of reactions containing metabolites
%
% Example: [P,C,vP,vC,s] = computeFluxSplits(model,mets,V,1);
%
% .. Author: - Hulda S. Haraldsdottir, December 2, 2016
%            - Modified by Diana El Assal 22/7/2017

if nargin <4;
    coeffSign = false; %default
end

if coeffSign
    s = sign(sum([model.S(ismember(model.mets,mets),:); sparse(1,size(model.S,2))],1))';
else
    s = sum([model.S(ismember(model.mets,mets),:); sparse(1,size(model.S,2))],1)'; % net stoichiometry. Zero if model.mets does not contain mets.
end

W = diag(s)*V; % stoichiometrically weighted flux
vP = max(W,0); % net production
vC = -min(W,0); % net consumption
P = vP*diag(1./sum(vP)); % relative production
C = vC*diag(1./sum(vC)); % relative consumption

author: opencobra

repo: cobratoolbox


text:
function gT = rbfinfwhiteKernDiagGradX(kern, t)

% RBFINFWHITEKERNDIAGGRADX Gradient of RBF-WHITE kernel's (with integration
% limits between minus infinity and infinity) diagonal w.r.t. t.
% FORMAT
% DESC computes the gradient of the diagonal of the RBF-WHITE kernel matrix
% with respect to the elements of the input column vector given in t.
% ARG kern : the kernel structure for which gradients are being computed.
% ARG t : the input data in the form of a design matrix.
% RETURN gT : the gradients of the diagonal with respect to each element
% of t. The returned matrix has the same dimensions as t.
%
% SEEALSO : rbfinfwhiteKernParamInit, kernDiagGradX, rbfinfwhiteKernGradX
%
% COPYRIGHT : David Luengo, 2009

% KERN

if size(t, 2) > 1
  error('Input can only have one column');
end

gT = zeros(size(t));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ node_num, element_num, element_type ] = tec_zone_line_parse ( line )

%*****************************************************************************80
%
%% TEC_ZONE_LINE_PARSE parses the "ZONE" line of a TEC file.
%
%  Discussion:
%
%    The string begins with the substring "ZONE" and is followed by
%    a sequence of keywords and values joined by an equals sign.
%
%    We expect the following, but in arbitrary order, separated
%    by spaces or commas:
%
%      N = number of nodes
%      E = number of elements
%      T = optional title (we can't handle this field right now)
%      PACKING = POINT
%      ZONETYPE = FETRIANGLE or FEQUADRILATERAL or FETETRAHEDRON or FEBRICK.
%
%    Other arguments that may appear on this line will be ignore.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 February 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string LINE, a string of characters, representing the
%    "VARIABLES=" line of the file.
%
%    Output, integer NODE_NUM, the number of nodes.
%
%    Output, integer ELEMENT_NUM, the number of elements.
%
%    Output, string ELEMENT_TYPE, the element type:
%    FETRIANGLE or FEQUADRILATERAL or FETETRAHEDRON or FEBRICK..
%

%
%  Remove the initial "ZONE"
%
  line = s_behead_substring ( line, 'ZONE' );
%
%  Replace each EQUALS sign by a space.
%  Also get rid of commas and periods.
%  Do single and double quotes have to go, also?
%
  line = s_replace_ch ( line, '=', ' ' );
  line = s_replace_ch ( line, ',', ' ' );
  line = s_replace_ch ( line, '.', ' ' );
%
%  Now each pair of words represents a name and a value.
%
  node_num = -1;
  element_num = -1;
  element_type = [];

  found_num = 0;

  while ( 1 )

    [ name, line ] = s_word_extract ( line );

    if ( isempty ( name ) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'TEC_ZONE_LINE_PARSE - Fatal error!\n' );
      fprintf ( 1, '  Unexpected End of input.\n' );
      error ( 'TEC_ZONE_LINE_PARSE - Fatal error!' );
    end

    [ value, line ] = s_word_extract ( line );

    if ( isempty ( value ) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'TEC_ZONE_LINE_PARSE - Fatal error!\n' );
      fprintf ( 1, '  Unexpected End of input.\n' );
      error ( 'TEC_ZONE_LINE_PARSE - Fatal error!' );
    end

    if ( ch_eqi ( name(1), 'N' ) & node_num == -1 )

      node_num = s_to_i4 ( value );
      found_num = found_num + 1;

    elseif ( ch_eqi ( name(1), 'E' ) & element_num == -1 )

      element_num = s_to_i4 ( value );
      found_num = found_num + 1;

    elseif ( s_eqi ( name, 'DATAPACKING' ) )

      if ( ~s_eqi ( value, 'POINT' ) )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'TEC_ZONE_LINE_PARSE - Fatal error!\n' );
        fprintf ( 1, '  Value of DATAPACKING argument must be POINT.\n' );
        error ( 'TEC_ZONE_LINE_PARSE - Fatal error!' );
      end

    elseif ( s_eqi ( name, 'ZONETYPE' ) & isempty ( element_type ) )

      found_num = found_num + 1;
      element_type = value;

    else

      fprintf ( 1, '\n' );
      fprintf ( 1, '  Ignoring "%s = %s".\n', name, value );

    end

    if ( found_num == 3 )
      break
    end

  end

  return
end


text:
function [V,F] = poly2mesh(filename,minimum_angle,maximum_area)
  % POLY2MESH Triangulate interoir of polygon read from .poly  file using
  % Triangle
  %
  %
  % [V,F] = poly2mesh(filename,minimum_angle,maximum_area)
  %
  % Inputs:
  %  filename  path to .poly file
  %  minimum_angle  minimum angle parameter for Triangle
  %  maximum_area  maximum area parameter for Triangle
  %
  % Outputs:
  %   V  #vertices by 2, list of vertex positions
  %   F  #faces by 3, list of face indices
  %
  % Copyright 2011, Alec Jacobson (jacobson@inf.ethz.ch)
  %
  % See also: png2objandtga, png2poly, writePOLY, poly2VEH, triangle
  %
  warning('THIS FILE IS DEPRECATED. CALL TRIANGLE.M DIRECTLY INSTEAD');
  basename = regexprep(filename,'\.poly$','');

  minimum_angle_args = '';
  if(exist('minimum_angle'))
    minimum_angle_args = ['q' num2str(minimum_angle)];
  end
  maximum_area_args = '';
  if(exist('maximum_area'))
    maximum_area_args = ['a' num2str(maximum_area)];
  end

  % use triangle to triangulate interior
  command_line_args = ['-p' minimum_angle_args maximum_area_args];
  [V,F] = execute_triangle(command_line_args,basename);
end

author: alecjacobson

repo: gptoolbox


text:
function [ a, more ] = compnz_next ( n, k, a, more )

%*****************************************************************************80
%
%% COMPNZ_NEXT computes the compositions of the integer N into K nonzero parts.
%
%  Discussion:
%
%    A composition of the integer N into K nonzero parts is an ordered sequence
%    of K positive integers which sum to N.  The compositions (1,2,1)
%    and (1,1,2) are considered to be distinct.
%
%    The routine computes one composition on each call until there are no more.
%    For instance, one composition of 6 into 3 parts is 3+2+1, another would
%    be 4+1+1 but 5+1+0 is not allowed since it includes a zero part.
%
%    On the first call to this routine, set MORE = FALSE.  The routine
%    will compute the first element in the sequence of compositions, and
%    return it, as well as setting MORE = TRUE.  If more compositions
%    are desired, call again, and again.  Each time, the routine will
%    return with a new composition.
%
%    However, when the LAST composition in the sequence is computed
%    and returned, the routine will reset MORE to FALSE, signaling that
%    the end of the sequence has been reached.
%
%  Example:
%
%    The 10 compositions of 6 into three nonzero parts are:
%
%      4 1 1,  3 2 1,  3 1 2,  2 3 1,  2 2 2,  2 1 3,
%      1 4 1,  1 3 2,  1 2 3,  1 1 4.
%
%  Modified:
%
%    01 December 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Albert Nijenhuis and Herbert Wilf,
%    Combinatorial Algorithms,
%    Academic Press, 1978, second edition,
%    ISBN 0-12-519260-6.
%
%  Parameters:
%
%    Input, integer N, the integer whose compositions are desired.
%
%    Input, integer K, the number of parts in the composition.
%    K must be no greater than N.
%
%    Input, integer A(K), the previous composition.  On the first call,
%    with MORE = FALSE, set A = [].  Thereafter, A should be the
%    value of A output from the previous call.
%
%    Input, logical MORE.  The input value of MORE on the first
%    call should be FALSE, which tells the program to initialize.
%    On subsequent calls, MORE should be TRUE, or simply the
%    output value of MORE from the previous call.
%
%    Output, integer A(K), the next composition.
%
%    Output, logical MORE, will be TRUE unless the composition
%    that is being returned is the final one in the sequence.
%
  persistent h;
  persistent t;
%
%  We use the trick of computing ordinary compositions of (N-K)
%  into K parts, and adding 1 to each part.
%
  if ( n < k )
    more = 0;
    a(1:k) = -1;
    return
  end

  if ( ~more )

    t = n - k;
    h = 0;
    a(1) = n - k;
    a(2:k) = 0;

  else

    a(1:k) = a(1:k) - 1;

    if ( 1 < t )
      h = 0;
    end

    h = h + 1;
    t = a(h);
    a(h) = 0;
    a(1) = t - 1;
    a(h+1) = a(h+1) + 1;

  end

  more = ( a(k) ~= ( n - k ) );

  a(1:k) = a(1:k) + 1;


text:
function normals = vertexNormal(vertices, faces)
%VERTEXNORMAL Compute normals to a mesh vertices
%
%   N = vertexNormal(V, F)
%   Computes vertex normals of the mesh given by vertices V and F.
%   V is a vertex array with 3 columns, F is either a NF-by-3 or NF-by-4
%   index array, or a cell array with NF elements.
%
%   Example
%     % Draw the vertex normals of a sphere
%     s = [10 20 30 40];
%     [v f] = sphereMesh(s);
%     drawMesh(v, f);
%     view(3);axis equal; light; lighting gouraud;
%     normals = vertexNormal(v, f);
%     drawVector3d(v, normals);
%
%   See also
%     meshes3d, faceNormal, triangulateFaces
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-12-19,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

nv = size(vertices, 1);
nf = size(faces, 1);

% unit normals to the faces
faceNormals = normalizeVector3d(faceNormal(vertices, faces));

% compute normal of each vertex: sum of normals to each face
normals = zeros(nv, 3);
if isnumeric(faces)
    for i = 1:nf
        face = faces(i, :);
        for j = 1:length(face)
            v = face(j);
            normals(v, :) = normals(v,:) + faceNormals(i,:);
        end
    end
else
    for i = 1:nf
        face = faces{i};
        for j = 1:length(face)
            v = face(j);
            normals(v, :) = normals(v,:) + faceNormals(i,:);
        end
    end
end

% normalize vertex normals to unit vectors
normals = normalizeVector3d(normals);

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function lambda = gear_eigenvalues ( ii, jj, n )

%*****************************************************************************80
%
%% GEAR_EIGENVALUES returns the eigenvalues of the GEAR matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer II, JJ, define the two special entries.
%
%    Input, integer N, the order of A.
%
%    Output, real LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );
%
%  Separate the sign and value.
%
  alpha = zeros ( n, 1 );

  j = abs ( ii );
  js = i4_sign ( ii );

  k = abs ( jj );
  ks = i4_sign ( jj );

  if ( 0 < js & 0 < ks )

    w = 0;

    phi = n - floor ( ( j + k ) / 2 );
    for p = 1 : phi
      w = w + 1;
      alpha(w) = 2 * p * pi / ( 2 * n + 2 - j - k );
    end

    phi = floor ( ( j - 1 ) / 2 );
    for p = 1 : phi
      w = w + 1;
      alpha(w) = 2 * p * pi / j;
    end

    phi = floor ( ( k - 1 ) / 2 );
    for p = 1 : phi
      w = w + 1;
      alpha(w) = 2 * p * pi / k;
    end

    w = w + 1;
    alpha(w) = 0.0;

    if ( i4_is_even ( j ) & i4_is_even ( k ) )
      w = w + 1;
      alpha(w) = pi;
    end

  elseif ( 0 < js & ks < 0 )

    w = 0;

    phi = n + 1 - floor ( (j+k+1) / 2 );
    for p = 1 : phi
      w = w + 1;
      alpha(w) = ( 2 * p - 1 ) * pi / ( 2 * n + 2 - j - k );
    end

    phi = floor ( ( j - 1 ) / 2 );
    for p = 1 : phi
      w = w + 1;
      alpha(w) = 2 * p * pi / j;
    end

    phi = floor ( k / 2 );
    for p = 1 : phi
      w = w + 1;
      alpha(w) = ( 2 * p - 1 ) * pi / k;
    end

    if ( i4_is_even ( j ) & i4_is_odd ( k ) )
      w = w + 1;
      alpha(w) = pi;
    end

  elseif ( js < 0 & 0 < ks )

    w = 0;

    phi = n + 1 - floor ( ( j + k + 1 ) / 2 );
    for p = 1 : phi
      w = w + 1;
      alpha(w) = ( 2 * p - 1 ) * pi / ( 2 * n + 2 - j - k );
    end

    phi = floor ( j / 2 );
    for p = 1 : phi
      w = w + 1;
      alpha(w) = ( 2 * p - 1 ) * pi / j;
    end

    phi = floor ( ( k - 1 ) / 2 );
    for p = 1 : phi
      w = w + 1;
      alpha(w) = 2 * p * pi / k;
    end

    if ( i4_is_odd ( j ) & i4_is_even ( k ) )
      w = w + 1;
      alpha(w) = pi;
    end

  elseif ( js < 0 & ks < 0 )

    w = 0;

    phi = n - floor ( ( j + k ) / 2 );
    for p = 1 : phi
      w = w + 1;
      alpha(w) = 2 * p * pi / ( 2 * n + 2 - j - k );
    end

    phi = floor ( j / 2 );
    for p = 1 : phi
      w = w + 1;
      alpha(w) = ( 2 * p - 1 ) * pi / j;
    end

    phi = floor ( k / 2 );
    for p = 1 : phi
      w = w + 1;
      alpha(w) = ( 2 * p - 1 ) * pi / k;
    end

    if ( i4_is_odd ( j ) & i4_is_odd ( k ) )
      w = w + 1;
      alpha(w) = pi;
    end

  end

  for w = 1 : n
    lambda(w,1) = 2.0 * cos ( alpha(w) );
  end

  return
end


text:
function tck = dtiExportFibersMrtrix(fg, tck_filename)
%
%   tck = dtiExportFibersMrtrix(fg, tck_filename)
%
%   Saves a FG structure to disk in MRTRIX format (.tck)
%
% INPUTS:
%     fg -           vistasoft fiber group structure
%     tck_filename - file path to output .ick
%
% OUTPUT:
%     tck - reconstructed TCK file structure.
%
%   Based on MATLAB functions distributed with mrtrix 0.3
%
%   Note 1. The .tck HEADER is reconstructed before writing to disk from
%   information found in fg.params. The HEADER does not seem to be necessary
%   to for the final .tck file to be compatible to MRTRIX and mrview.
%
%   Note 2. The HEADER can contain information aboutThe ROIs used for tracking.
%   However, if multiple ROIs were used for tracking currently the MATLAB code
%   only stores the last ROI name used. This is a limitation existent in the
%   original MATLAB code provided with MRTRIX. We could try to solve this but
%   have not yet.
%
% Brent McPherson, Indiana University (C), 2017
% Updated by Lindsey Kitchell, Indiana University 2017
% Edits to comments by Franco Pestilli

% initialize output
tck = dtiGetFgParams(fg);

% transpose the streamline order
tck.data = fg.fibers';

for ii = 1:length(tck.data)
    % transpose the node order
    tck.data{ii} = tck.data{ii}';

end
clear ii

% save the file w/ mrtrix fxn
write_mrtrix_fibers(tck, tck_filename);

return
end

%%%%%%%%%%%%%%%%%%%%%%
% AUXILIARY FUNCTION %
%%%%%%%%%%%%%%%%%%%%%%
function write_mrtrix_fibers (fibers, filename)
%
% function: write_mrtrix_fibers (fibers, filename)
%
% writes the track data stored as a cell array in the 'data' field of the
% fibers variable to the MRtrix format track file 'filename'. All other fields
% of the fibers variable will be written as text entries in the header, and are
% expected to supplied as character arrays.
%
%
% This function was originally distributed with mrtrix 0.2/0.3

if ~isfield (fibers, 'data')
  disp ('ERROR: input fibers variable does not contain required ''data'' field');
  return;
end

if ~iscell (fibers.data)
  disp ('ERROR: input fibers.data variable should be a cell array');
  return;
end

f = fopen (filename, 'w', 'ieee-le');
if (f < 1)
  disp (['error opening ' filename ]);
  return;
end

fprintf (f, 'mrtrix fibers\ndatatype: Float32LE\ncount: %d\n', prod(size(fibers.data)));
names = fieldnames(fibers);
for i=1:size(names)
  if strcmpi (names{i}, 'data'), continue; end
  if strcmpi (names{i}, 'count'), continue; end
  if strcmpi (names{i}, 'datatype'), continue; end
  fprintf (f, '%s: %s\n', names{i}, getfield(fibers, names{i}));
end
fprintf (f, '%s: %s\n', 'timestamp: ', num2str(now)); % This allows visualizing the streamlines in mrview
data_offset = ftell (f) + 20;
fprintf (f, 'file: . %d\nEND\n', data_offset);

fwrite (f, zeros(data_offset-ftell(f),1), 'uint8');
for i = 1:prod(size(fibers.data))
  fwrite (f, fibers.data{i}', 'float32');
  fwrite (f, [ nan nan nan ], 'float32');
end

fwrite (f, [ inf inf inf ], 'float32');
fclose (f);

end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% Copyright (C) 2015, 2016, 2018, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym ipermute (@var{B}, @var{iperm})
%% Invert a permutation the indices of a symbolic array.
%%
%% Example:
%% @example
%% @group
%% A = sym([1 2 pi; 4 5 6]);
%% B = permute(A, [2 1]);
%% ipermute(B, [2 1])
%%   @result{} ans = (sym 2×3 matrix)
%%       ⎡1  2  π⎤
%%       ⎢       ⎥
%%       ⎣4  5  6⎦
%% @end group
%% @end example
%%
%% @seealso{@@sym/permute}
%% @end defmethod

function A = ipermute(B, iperm)

  if (nargin < 2)
    print_usage ();
  end

  A = permute(B, iperm);

end

%!error permute (sym(1))
%!error permute (sym(1), 2, 3)

%!test
%! syms x
%! A = [1 x];
%! perm = [2 1];
%! B = permute(A, perm);
%! C = ipermute(B, perm);
%! assert (isequal(C, A))


text:
function [names,data,data_endo,data_endo_a,data_endo_c,data_endo_c_lags,data_exo,data_exo_a,data_exo_p,data_exo_c,data_exo_c_lags,Fperiods,Fcomp,Fcperiods,Fcenddate,ar,priorexo,lambda4_2,favar]=...
    gensample(startdate,enddate,VARtype,Fstartdate,Fenddate,Fendsmpl,endo,exo,frequency,lags,F,CF,ar,lambda4,PriorExcel,priorsexogenous,pref,favar,IRFt, numendo)

% Phase 1: data loading and error checking

% if we have a FAVAR: read information data, data transformation, create indices, compute factors (PC)
if favar.FAVAR==1
    [informationstartlocation,informationendlocation,favar]=bear.favar_gensample1(startdate,enddate,favar,pref);
end

% first read the data from Excel
[data,names]=xlsread(pref.excelFile,'data');

% now, as a preliminary step: check if there is any Nan in the data; if yes, return an error since the model won't be able to run with missing data
% a simple way to test for NaN is to check for "smaller or equal to infinity": Nan is the only number for which matlab will return 'false' when asked so
%[r,c]=size(data);

% identify the date strings
datestrings=names(2:end,1);
% identify the position of the string corresponding to the start period
startlocationData=find(strcmp(datestrings,startdate));
% identify the position of the string corresponding to the end period
endlocation=find(strcmp(datestrings,enddate));
startlocation=startlocationData;
if favar.FAVAR==1 % in case we transform the data to first or second differences, we have a different startlocation
    if favar.transformation==1
        startlocation=informationstartlocation;
        endlocation=informationendlocation;
        % check is correct
    else
        startlocation=startlocationData;
    end
end

% save the whole sample temporarily
data1=data;
% adjust data to startdate and enddate
% data1=data(startlocationData:endlocation,:);
data=data(startlocationData:endlocation,:);

% identify the variable strings, endogenous and exogenous
variablestrings=names(1,2:end);

% FAVAR: augment data and variablestrings with factors
if favar.FAVAR==1
    [data,variablestrings,favar]=bear.favar_gensample2(data1,endo,variablestrings,startlocation,lags,favar);
end

% if either the start date or the date date is not recognised, return an error message
if isempty(startlocation)
    error('bear:BEARmain:UnknownStartDate', ...
        'Error: unknown start date for the sample. Please check your sample start date (remember that names are case-sensitive).');
elseif isempty(endlocation)
    error('bear:BEARmain:UnknownEndDate', ...
        'Error: unknown end date for the sample. Please check your sample end date (remember that names are case-sensitive).');
end
% also, if the start date is posterior to the end date, obviously return an error
if startlocation>=endlocation==1
    error('bear:BEARmain:InconsistentStartEndDates', ...
        'Error: inconsistency between the start and end dates. The start date must be anterior to the end date.');
end

% identify the position of the strings corresponding to the endogenous variables
% for each variable, find the corresponding string
for ii=1:numendo
    % check first that the variable ii in endo appears in the list of variable strings
    % if not, the variable is unknown: return an error
    var=endo{ii,1};
    check=find(strcmp(variablestrings,var));
    if isempty(check)==1
        message=['Error: endogenous variable ' var ' cannot be found on the excel data spreadsheet.'];
        error('BEARmain:gensample:EndoVarNotFound', message);
    end
    % if the variable is known, go on
    endolocation(ii,1)=find(strcmp(variablestrings,endo(ii,1)));
end

% identify the position of the strings corresponding to the exogenous variables
% proceed similarly to the endogenous variables, but account for the fact that exogenous may be empty
% so check first whether there are exogenous variables altogether
if isempty(exo)
    numexo=0;
else
    % if not empty, repeat what has been done with the exogenous
    numexo=size(exo,1);
    % for each variable, find the corresponding string
    for ii=1:numexo
        % check first that the variable ii in endo appears in the list of variable strings
        % if not, the variable is unknown: return an error
        var=exo{ii,1};
        check=find(strcmp(variablestrings,var));
        if isempty(check)==1
            message=['Error: exogenous variable ' var ' cannot be found on the excel data spreadsheet.'];
            error('BEARmain:gensample:ExoVarNotFound', message);
        end
        % if the variable is known, go on
        exolocation(ii,1)=find(strcmp(variablestrings,exo(ii,1)));
    end
end

% Phase 2: creation of the data matrices data_endo and data_exo

% now create the matrix of endogenous variables for the estimation sample
% it is simply the concatenation of the vectors of each endogenous variables, over the selected sample dates
data_endo=[];
% loop over endogenous variables
for ii=1:numendo
    data_endo=[data_endo data(:,endolocation(ii,1))];
end

% Similarly, create the matrix of exogenous variables for the estimation sample
data_exo=[];
for ii=1:numexo
%     data_exo=[data_exo data(startlocation:endlocation,exolocation(ii,1))];
    data_exo=[data_exo data(:,exolocation(ii,1))];
end

%% FAVAR

if favar.FAVAR==1
    if favar.transformation==1
        % correct order of transformationindex_endo
        favar.transformationindex_endo=[];
        for ii=1:size(endo,1)
            favar.transformationindex_endo=[favar.transformationindex_endo favar.transformationindex_endo_temp(:,endolocation(ii,1))];
        end
    end
    % correct order of stddev
    favar.data_exfactors_stddev=[];
    for ii=1:size(endo,1)
        favar.data_exfactors_stddev=[favar.data_exfactors_stddev favar.data_exfactors_stddev_temp(:,endolocation(ii,1))];
    end

    % compute new loadings
    %favar.L=(mvregress(data_endo,favar.X,'algorithm','cwls'))'; %appears to be not feasible for very large X
    %favar.L=(olssvd(favar.X,data_endo))';

    % IRF shock to plot
    if favar.IRFplot==1
        favar.IRF.npltXshck=size(favar.IRF.pltXshck,1);
        if IRFt==1||IRFt==2||IRFt==3
            plotXshock_indexlogical=ismember(endo,favar.IRF.pltXshck);
            favar.IRF.plotXshock_index=find(plotXshock_indexlogical==1)';
            if favar.IRF.npltXshck==0
                % error if no shock to plot is found, otherwise code crashes at a later stage
                message=['Error: Shock(' favar.IRF.npltXshck ') cannot be found.'];
                error('BEARmain:gensample:favar_IRF_npltXshck_error',message);
            end
        end
        % for IRFt 4 & 6 this step is done in loadsignres
    end

    % rotate Factors, compute new loadings for onestep, twostep estimation
    [data_endo,favar]=bear.favar_gensample3(data_endo,favar);

%%%% stationarity test
endostrings=endo(favar.variablestrings_exfactors,1);
allstrings=[favar.informationvariablestrings,endostrings'];
NSTindex=[];
NST2index=[];
catallstrings=[];
catallstrings2=[];
count=0;
count2=0;
for ii=1:size(favar.XY,2)
    adf(ii,1)=adftest(favar.XY(:,ii));
    kpss(ii,1)=kpsstest(favar.XY(:,ii));
    if adf(ii,1)==0 %|| kpss(ii,1)==1
        count=count+1;
        NSTallstrings{count,1}=allstrings(1,ii);
        NSTindex=[NSTindex;ii];
        catallstrings=strcat(catallstrings,', ',NSTallstrings{count,1});
    end
    if kpss(ii,1)==1
        count2=count2+1;
        NST2allstrings{count2,1}=allstrings(1,ii);
        NST2index=[NST2index;ii];
        catallstrings2=strcat(catallstrings2,', ',NST2allstrings{count2,1});
    end
end
if size(NSTindex,1) > 0
fprintf('%d%s%s%s\n',size(NSTindex,1),' series in X are not stationary (ADF test): ',catallstrings{1,1},'.');
end
if size(NST2index,1) > 0
fprintf('%d%s%s%s\n',size(NST2index,1),' series in X are not stationary (KPSS test): ',catallstrings2{1,1},'.');
end
%%%%

end

%% Bens check for missing values in exo and endo
[t,g]=size(data_endo);
for ii=1:t
    for jj=1:g
        temp=data_endo(ii,jj);
        if (temp<=inf)==0
            % identify the variable and the date %%%%% why is this commented
            %       NaNvariable=names{1,jj+1};
            %       NaNdate=names{ii+1,1};
            %       message=['Error: variable ' NaNvariable ' at date ' NaNdate ' (and possibly other sample entries) is identified as NaN. Please check your Excel spreadsheet: entry may be blank or non-numerical.'];
            %       msgbox(message);
            error('programme termination: An endogenous variable entry is missing');
        end
    end
end
%
%
[t,g]=size(data_exo);
for ii=1:t
    for jj=1:g
        temp=data_exo(ii,jj);
        if (temp<=inf)==0
            error('programme termination: An exogenous variable entry is missing');
        end
    end
end

% Prior values settings for the AR coefficients (either default same for all or individually in Excel in AR prior sheet)
if PriorExcel==0
    ar_default=NaN(numendo,1);
    ar_default(:,1)=ar;
    ar=ar_default;
else
    [ar]=xlsread(pref.excelFile,'AR priors');
end

if priorsexogenous==0
    exo_default=zeros(numendo,numexo+1);
    lambda4_default=zeros(numendo,numexo+1);
    for ii=1:numendo
        for jj=1:numexo+1
            priorexo(ii,jj)=0;
            lambda4_2(ii,jj)=lambda4;
        end
    end
else
    [priorexo]=xlsread(pref.excelFile,'exo mean priors');
    [lambda4]=xlsread(pref.excelFile,'exo tight priors');
    priorexo=priorexo(1:numendo,1:numexo+1);
    lambda4_2=lambda4(1:numendo,1:numexo+1);
end

% Phase 3: determination of the position of the forecast start and end periods

% if both unconditional and conditional forecasts were not selected, there is no need for all the forecast-specific matrices: simply return empty matrices
if (VARtype==1 && F==0) || ((VARtype==2 || VARtype==3 || VARtype==5 || VARtype==6 ) && (F==0 && CF==0))
    data_endo_a=[];
    data_exo_a=[];
    data_exo_p=[];
    Fperiods=[];
    Fcomp=[];
    Fcperiods=[];
    data_endo_c=[];
    data_endo_c_lags=[];
    data_exo_c=[];
    data_exo_c_lags=[];
    Fcenddate=[];

    % if forecast were selected, create all the required elements
else

    % preliminary tasks
    % first, identify the date strings, and the variable strings
    datestrings=names(2:end,1);
    if favar.FAVAR==1 % in case we transform the data to first or second differences, we have a different datestrings
        if favar.transformation==1
            datestrings=names(1+favar.informationstartlocation:end,:); %first row are labels
        end
    end

    % identify the location of the last period in the dataset
    dataendlocation=size(datestrings,1);

    % identify the position of the start period (the easy part)
    % if the start period has been selected as the first period after the sample end, identifies it directly
    if Fendsmpl==1
        Fstartlocation=find(strcmp(datestrings,enddate))+1;
        % if the start period has not been selected as the first period after the sample end, it must be within the sample: look for it
    elseif Fendsmpl==0
        Fstartlocation=find(strcmp(datestrings,Fstartdate));
        % if the start date is not recognised, return an error message
        if isempty(Fstartlocation)==1
            msgbox('Error: unknown start date for the forecasts. Select a date within a sample, or select "Start forecasts after last sample period"');
            error('unknown start date for the forecasts');
        end
    end

    % identify the position of the final forecast period (the hard part)
    % this period can be in or outside the sample, depending on the user's choice

    % if the data is yearly
    if frequency==1
        Fendlocation=str2num(Fenddate(1,1:end-1))-str2num(datestrings{1,1}(1,1:end-1))+1;

        % if the data is quarterly
    elseif frequency==2
        % first identify the year and quarter of the initial date in the whole data set (not just the sample)
        datastartyear=str2num(datestrings{1,1}(1,1:4));
        datastartquarter=str2num(datestrings{1,1}(1,6));
        % convert this date into quarters only
        datastart=datastartyear*4+datastartquarter;
        % then identify the year and quarter of the final forecast date
        forecastendyear=str2num(Fenddate(1,1:4));
        forecastendquarter=str2num(Fenddate(1,6));
        % convert this date into quarters only
        forecastend=forecastendyear*4+forecastendquarter;
        % finally, compute the number of periods that separate the two dates
        Fendlocation=forecastend-datastart+1;

        % if the data is monthly
    elseif frequency==3
        % first identify the year and month of the initial date in the whole data set (not just the sample)
        temp=datestrings{1,1};
        datastartyear=str2num(temp(1,1:4));
        temp(1,5)=' ';
        [~,datastartmonth]=strtok(temp);
        % convert this date into months only
        datastart=datastartyear*12+str2num(datastartmonth);
        % then identify the year and month of the final forecast date
        temp=Fenddate;
        forecastendyear=str2num(temp(1,1:4));
        temp(1,5)=' ';
        [~,forecastendmonth]=strtok(temp);
        % convert this date into months only
        forecastend=forecastendyear*12+str2num(forecastendmonth);
        Fendlocation=forecastend-datastart+1;

        % if the data is weekly
    elseif frequency==4
        % then identify the year and week corresponding to this final period
        temp=datestrings{end,1};
        dataendyear=str2num(temp(1,1:4));
        temp(1,5)=' ';
        [~,dataendweek]=strtok(temp);
        dataendweek=str2num(dataendweek);
        % identify the year and week corresponding to the end of the forecast period
        temp=Fenddate;
        Fendyear=str2num(temp(1,1:4));
        temp(1,5)=' ';
        [~,Fendweek]=strtok(temp);
        Fendweek=str2num(Fendweek);
        % determine whether the forecast period ends within the data set, or after the end of the data set
        if Fendyear<dataendyear
            insmpl=1;
        elseif Fendyear==dataendyear && Fendweek<=dataendweek
            insmpl=1;
        else
            insmpl=0;
        end
        % if the forecast end period lies within the dataset, simply detect its location
        if insmpl==1
            Fendlocation=find(strcmp(Fenddate,datestrings));
            % if it is outside the data set, complete until the forecast end date is reached (assuming 52 weeks per year)
        elseif insmpl==0
            % Consider two cases separately
            % first case: if the end year of the forecast is the same as the end year of the dataset
            % in this case, simply complete the missing periods
            if dataendyear==Fendyear
                complement=Fendweek-dataendweek;
                Fendlocation=dataendlocation+complement;
                % if the end year of the forecasts is posterior to the last year of the data
            elseif dataendyear<Fendyear
                % complete the first year (the one shared that ends the data set)
                complement=52-dataendweek;
                % complete the following years before the last one (if any)
                for ii=dataendyear+1:Fendyear-1
                    complement=complement+52;
                end
                % complete the final forecast year
                complement=complement+Fendweek;
                Fendlocation=dataendlocation+complement;
            end
        end

        % if the data is daily
    elseif frequency==5
        % then identify the year and day corresponding to this final period
        temp=datestrings{end,1};
        dataendyear=str2num(temp(1,1:4));
        temp(1,5)=' ';
        [~,dataendday]=strtok(temp);
        dataendday=str2num(dataendday);
        % identify the year and day corresponding to the end of the forecast period
        temp=Fenddate;
        Fendyear=str2num(temp(1,1:4));
        temp(1,5)=' ';
        [~,Fendday]=strtok(temp);
        Fendday=str2num(Fendday);
        % determine whether the forecast period ends within the data set, or after the end of the data set
        if Fendyear<dataendyear
            insmpl=1;
        elseif Fendyear==dataendyear && Fendday<=dataendday
            insmpl=1;
        else
            insmpl=0;
        end
        % if the forecast end period lies within the dataset, simply detect its location
        if insmpl==1
            Fendlocation=find(strcmp(Fenddate,datestrings));
            % if it is outside the data set, complete until the forecast end date is reached (assuming 261 opening days per year)
        elseif insmpl==0
            % Consider two cases separately
            % first case: if the end year of the forecast is the same as the end year of the dataset
            % in this case, simply complete the missing periods
            if dataendyear==Fendyear
                complement=Fendday-dataendday;
                Fendlocation=dataendlocation+complement;
                % if the end year of the forecasts is posterior to the last year of the data
            elseif dataendyear<Fendyear
                % complete the first year (the one shared that ends the data set)
                complement=261-dataendday;
                % complete the following years before the last one (if any)
                for ii=dataendyear+1:Fendyear-1
                    complement=complement+261;
                end
                % complete the final forecast year
                complement=complement+Fendday;
                Fendlocation=dataendlocation+complement;
            end
        end

        % finally, if the data is undated
    elseif frequency==6
        Fendlocation=str2num(Fenddate(1,1:end-1))-str2num(datestrings{1,1}(1,1:end-1))+1;
    end

    % from this, conclude the total number of forecast periods
    Fperiods=Fendlocation-Fstartlocation+1;

    if Fperiods<0
        msgbox('Error: The forecast start date needs to be prior to the forecast end date');
        error('invalid forecast start or end date');
    end

    % Phase 4: generation of the forecast-specific matrices
    if favar.FAVAR==0
        % load the full sample in this case
        data=data1;
    elseif favar.FAVAR==1
        data=favar.data_full;
    end
    % now create the matrix of endogenous variables for the pre-forecast period
    % it is simply the concatenation of the vectors of each endogenous variables, over the selected sample dates
    data_endo_a=[];
    % loop over endogenous variables
    for ii=1:numendo
        data_endo_a=[data_endo_a data(1:Fstartlocation-1,endolocation(ii,1))];
    end
    % also, create the matrix of exogenous variables for the pre-forecast period
    data_exo_a=[];
    for ii=1:numexo
        data_exo_a=[data_exo_a data(1:Fstartlocation-1,exolocation(ii,1))];
    end

    % create the matrix of endogenous variables for the period common to actual data and forecasts (for forecast evaluation)
    % first, check that there are such common periods: it is the case if the beginning of the forecast period is anterior to the end of the dataset
    if Fstartlocation<=dataendlocation
        % return a scalar value to indicate that forecast evaluation is possible
        Fcomp=1;
        % compute the number of common periods
        % if the forecast period ends before the end of the data set, the common periods end with the end of the forecasts
        if Fendlocation<=dataendlocation
            Fcperiods=Fperiods;
            % record the end date of the common periods
            Fcenddate=Fenddate;
            % if the forecast period ends later than the data set, the common periods end at the end of the data set
        elseif Fendlocation>dataendlocation
            Fcperiods=dataendlocation-Fstartlocation+1;
            % record the end date of the common periods
            Fcenddate=datestrings{end,1};
        end

        % create a matrix of endogenous data for the common periods
        data_endo_c=[];
        for ii=1:numendo
            data_endo_c=[data_endo_c data(Fstartlocation:min(dataendlocation,Fendlocation),endolocation(ii,1))];
        end

        % create a lagged matrix of endogenous data prior to the common periods
        % the number of values is equal to "lags"; this will be used for computation of the log predictive score
        data_endo_c_lags=[];
        for ii=1:numendo
            data_endo_c_lags=[data_endo_c_lags data(Fstartlocation-lags:Fstartlocation-1,endolocation(ii,1))];
        end

        % create a matrix of exogenous data for the common periods
        data_exo_c=[];
        for ii=1:numexo
            data_exo_c=[data_exo_c data(Fstartlocation:min(dataendlocation,Fendlocation),exolocation(ii,1))];
        end

        % create a lagged matrix of exogenous data prior to the common periods
        % the number of values is equal to "lags"; this will be used for computation of the log predictive score
        data_exo_c_lags=[];
        for ii=1:numexo
            data_exo_c_lags=[data_exo_c_lags data(Fstartlocation-lags:Fstartlocation-1,exolocation(ii,1))];
        end
        % if there are no common periods, return a scalar value to indicate that forecast evaluation is not possible
    else
        Fcomp=0;
        Fcperiods=0;
        data_exo_c=[];
        data_endo_c=[];
        Fcenddate=[];
        data_endo_c_lags=[];
        data_exo_c_lags=[];
    end

    % now create the matrix data_exo_p
    % two possible cases

    % if there are no exogenous variables, simply create an empty matrix
    if isempty(exo)
        data_exo_p=[];

        % if there are exogenous variables, load from excel
    else
        % load the data from Excel
        [num txt strngs]=xlsread(pref.excelFile,'pred exo');

        % obtain the row location of the forecast start date
        [Fstartlocation,~]=find(strcmp(strngs,Fstartdate));
        % check that the start date for the forecast appears in the sheet; if not, return an error
        if isempty(Fstartlocation)
            message=['Error: a forecast application is selected for a model that uses exogenous variables. Hence, predicted exogenous values should be supplied over the forecast periods. Yet the start date for forecasts (' Fstartdate ') cannot be found on the ''pred exo'' sheet of the Excel data file. Please verify that this sheet is properly filled, and remember that dates are case-sensitive.'];
            msgbox(message);
            error('programme termination: data error');
        end
        % obtain the row location of the forecast end date
        [Fendlocation,~]=find(strcmp(strngs,Fenddate));
        % check that the end date for the forecast appears in the sheet; if not, return an error
        if isempty(Fendlocation)
            message=['Error: a forecast application is selected for a model that uses exogenous variables. Hence, predicted exogenous values should be supplied over the forecast periods. Yet the end date for forecasts (' Fenddate ') cannot be found on the ''pred exo'' sheet of the Excel data file. Please verify that this sheet is properly filled, and remember that dates are case-sensitive.'];
            msgbox(message);
            error('programme termination: data error');
        end

        % identify the strings for the exogenous variables
        % loop over exogenous
        for ii=1:numexo
            % try to find a column match for exogenous variable ii
            [~,location]=find(strcmp(strngs,exo{ii,1}));
            % if no match is found, return an error
            if isempty(location)
                message=['Error: a forecast application is selected for a model that uses exogenous variables. Hence, predicted exogenous values should be supplied over the forecast periods. Yet the exogenous variable ''' exo{ii,1} ''' cannot be found on the ''pred exo'' sheet of the Excel data file. Please verify that this sheet is properly filled, and remember that variable names are case-sensitive.'];
                msgbox(message);
                error('programme termination: data error');
                % else, record the value
            else
                pexolocation(ii,1)=location;
            end
        end

        % if everything was fine, reconstitute the matrix data_exo_p
        % initiate
        data_exo_p=[];
        % loop over exogenous variables
        for ii=1:numexo
            % initiate the predicted values for exogenous variable ii
            predexo=[];
            % loop over forecast periods
            for jj=1:Fperiods
                temp=strngs{Fstartlocation+jj-1,pexolocation(ii,1)};
                % if this entry is empty or NaN, return an error
                if (isempty(temp) || (temp<=inf)==0)
                    message=['Error: the predicted value for exogenous variable ' exo{ii,1} ' at forecast period ' strngs{Fstartlocation+jj,1} ' (and possibly other entries) is either empty or NaN. Please verify that the ''pred exo'' sheet of the Excel data file is properly filled.'];
                    msgbox(message);
                    error('programme termination: data error');
                    % if this entry is a number, record it
                else
                    predexo=[predexo;temp];
                end
            end
            % concatenate
            data_exo_p=[data_exo_p predexo];
        end

        % also, record the exogenous values on Excel
        % replace NaN entries by blanks
        strngs(cellfun(@(x) any(isnan(x)),strngs))={[]};
        % then save on Excel
        if pref.results==1
            bear.xlswritegeneral(fullfile(pref.results_path, [pref.results_sub '.xlsx']),strngs,'pred exo','A1');
        end
    end

end

author: european-central-bank

repo: BEAR-toolbox