text:
function deposit(varargin)
%*************************************************
%DATE: 11/8/2007 (created); 	23/8/2007 (modified)
%FUNCTION:
%INPUTS:
%OUTPUTS:
%*************************************************
materialName 	= varargin{1};
thickness		= varargin{2};

shiftFrom	= 0;

%=======================material color table=========================
global Layerprofile;
global layerMatrix;

switch materialName
	case 'Si'
		materialNrm = 1;
	case 'SiO2'
		materialNrm = 2;
	case 'metal'
		materialNrm = 3;
	case 'SiNi'
		materialNrm = 4;
	case 'PoSi'
		materialNrm = 5;
	case 'PR'
		materialNrm = 9;
end

newLayerMatrixHeight = round (thickness);
newLayerMatrix = zeros(newLayerMatrixHeight, 96);

if sum(sum(layerMatrix)) == 0 % if it is empty
layerMatrix =  newLayerMatrix + materialNrm;
else
layerMatrix = vertcat(newLayerMatrix, layerMatrix); % unite
	for pointIndex = 1:96
		depositFrom = find(layerMatrix(:, pointIndex), 1,'last');
      	layerMatrix (((depositFrom+1 ) : (depositFrom+1+newLayerMatrixHeight)), pointIndex) = materialNrm;
	end
end

layerMatrix = rounddeposit(layerMatrix, materialNrm);

function layerMatrix = rounddeposit(layerMatrix, materialNrm)
%*************************************************
%DATE: 23/8/2007 (created); 	23/8/2007 (modified)
%FUNCTION:
%INPUTS:
%OUTPUTS:
%*************************************************

	for pointIndex = 1:96
		depositFromArray(pointIndex) = find(layerMatrix(:, pointIndex), 1,'last');
	end

	oldDepositFromArray = depositFromArray;
save('old', 'oldDepositFromArray');
	for pointIndex = 1:95
		diff = oldDepositFromArray(pointIndex) - oldDepositFromArray(pointIndex+1);
		if (diff>0)
		addDepositArray = calculateroundarray(diff);

		index = pointIndex + 1;
		nrm   = 1;
			while (index <= 96 & nrm <= size(addDepositArray, 2))
			depositFromArray(index) = depositFromArray(index) + addDepositArray(nrm);
			index = index + 1;
			nrm   = nrm + 1;
			end
		end
	end

	for pointIndexI = 0:94
		pointIndex = 96 - pointIndexI;
		diff = oldDepositFromArray(pointIndex) - oldDepositFromArray(pointIndex-1);
		if (diff>0)
		addDepositArray = calculateroundarray(diff);

			index = pointIndex - 1;
			nrm   = 1;
			while (index > 0 & nrm <= size(addDepositArray, 2))
			depositFromArray(index) = depositFromArray(index) + addDepositArray(nrm);
			index = index - 1;
			nrm   = nrm + 1;
			end
		end
	end

	for pointIndex = 1:96
		layerMatrix(oldDepositFromArray(pointIndex):depositFromArray(pointIndex), pointIndex) = materialNrm;
	end

function roundArray = calculateroundarray(diff)
%*************************************************
%FUNCTION:
%INPUTS:
%OUTPUTS:
%*************************************************
roundArray(1) = 0;
		a = -4/diff;
		b = diff;
		for pointDiff=1:round(diff/2)
		temp = floor (a * pointDiff^2 + b);
			if (temp < 0)
				return;
			end
			roundArray(pointDiff) = temp;
		end

repo: mcbench-benchmarks


text:
function J = rstepJacobian(f,w)

x0=w;
y0 = f(x0);
n = length(x0);
m = length(y0);
J = zeros(m,n);
dx = sqrt(eps);
for i=1:n;
    x = x0;
    x(i) = x0(i) + dx;
    y1 = f(x);

    x = x0;
    x(i) = x0(i) - dx;
    y0 = f(x);
    J(:,i) = (y1-y0)/(2*dx);
end

author: nesl


text:
function volume = cube01_volume ( )

%*****************************************************************************80
%
%% CUBE01_VOLUME returns the volume of the unit cube in 3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real VOLUME, the volume.
%
  volume = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% MAKE_B creates a weighted symmetric filter matrix encoding the
% neighborhood for isotropic label diffusion.
%
% Usage:
%
% B = make_B(1,4);
% Returns a 3 x 3 symmetric filter matrix encoding a 4-neighborhood
% corresponding to 4 equally spaced points on a 1 point radius of the
% center point.
%
% B = make_B(2,16);
% Returns a 5 x 5 symmetric filter matrix encoding 16 equally spaced points
% on a 2 point radius of the center point.
%
% Inputs:
%
% radius: distance from the center point used to determine filter matrix.
% num_points: number of neighbors used in the filter matrix.
%
% Output:
%
% B: filter matrix (approximately) encoding num_points equally spaced points
% on the speccified radius from the center point.
%
% Copyright (c) Roman Garnett, 2012--2014.

function B = make_B(radius, num_points)

  factor = 100;

  small_length = 2 * radius + 1;
  big_length = factor * small_length;

  big_B = zeros(big_length);

  offset = big_length / 2;
  side = factor / 2;

  for i = 1:num_points
    angle = (i - 1) / num_points * 2 * pi;

    x = round(offset + factor * radius * cos(angle));
    y = round(offset + factor * radius * sin(angle));

    x_lower = max(x - side, 1);
    x_upper = min(x + side, big_length);
    y_lower = max(y - side, 1);
    y_upper = min(y + side, big_length);

    big_B(x_lower:x_upper, y_lower:y_upper) = ...
        big_B(x_lower:x_upper, y_lower:y_upper) + 1;
  end

  B = zeros(small_length, small_length);

  for i = 1:small_length
    for j = 1:small_length
      x_range = (1 + (i - 1) * factor):(i * factor);
      y_range = (1 + (j - 1) * factor):(j * factor);

      B(i, j) = sum(sum(big_B(x_range, y_range)));
    end
  end

  B = B + rot90(B, 1) + rot90(B, 2) + rot90(B, 3);
  B = B / sum(B(:));

  B = round(B * 100) / 100;
end

author: muhanzhang

repo: DGCNN

save_path: github-repos/MATLAB/muhanzhang-DGCNN


text:
function test_suite = test_distancePointMesh
%TEST_DISTANCEPOINTTRIMESH  One-line description here, please.
%
%   output = test_distancePointMesh(input)
%
%   Example
%   test_distancePointMesh
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2018-03-08,    using Matlab 9.3.0.713579 (R2017b)
% Copyright 2018 INRA - Cepia Software Platform.

% corps du fichier contenant la suite de tests

% testSuite = buildFunctionHandleTestSuite(localfunctions);
test_suite = functiontests(localfunctions);

function testTetraedron(testCase) %#ok<*DEFNU>

[v, f] = createOctahedron;

testCase.assertEqual(1, distancePointMesh([2 0 0], v, f), .0001);
testCase.assertEqual(1, distancePointMesh([0 2 0], v, f), .0001);
testCase.assertEqual(1, distancePointMesh([0 0 2], v, f), .0001);
testCase.assertEqual(1, distancePointMesh([-2 0 0], v, f), .0001);
testCase.assertEqual(1, distancePointMesh([0 -2 0], v, f), .0001);
testCase.assertEqual(1, distancePointMesh([0 0 -2], v, f), .0001);

exp = sqrt(2) / 2;
testCase.assertEqual(exp, distancePointMesh([1 1 0], v, f), .0001);
testCase.assertEqual(exp, distancePointMesh([1 0 1], v, f), .0001);
testCase.assertEqual(exp, distancePointMesh([0 1 1], v, f), .0001);

function testSingleTriangle(testCase)

v = [10 10 10 ; 20 10 10 ; 10 20 10];
f = [1 2 3];

% around lower left corner
testCase.assertEqual(2, distancePointMesh([ 10  8 10 ], v, f), .0001);
testCase.assertEqual(2, distancePointMesh([  8 10 10 ], v, f), .0001);
testCase.assertEqual(2, distancePointMesh([ 12  8 10 ], v, f), .0001);
testCase.assertEqual(2, distancePointMesh([  8 12 10 ], v, f), .0001);
testCase.assertEqual(2*sqrt(2), distancePointMesh([  8  8 10 ], v, f), .0001);
% more around lower left corner
testCase.assertEqual(sqrt(5), distancePointMesh([  8  9 10 ], v, f), .0001);
testCase.assertEqual(sqrt(5), distancePointMesh([  9  8 10 ], v, f), .0001);

% around right corner
testCase.assertEqual(2, distancePointMesh([ 20  8 10 ], v, f), .0001);
testCase.assertEqual(2, distancePointMesh([ 18  8 10 ], v, f), .0001);
testCase.assertEqual(2, distancePointMesh([ 22 10 10 ], v, f), .0001);
testCase.assertEqual(2*sqrt(2), distancePointMesh([ 22  8 10 ], v, f), .0001);
testCase.assertEqual(sqrt(5), distancePointMesh([ 22  9 10 ], v, f), .0001);

% around top left corner
testCase.assertEqual(sqrt(5), distancePointMesh([  9 22 10 ], v, f), .0001);

function test_ortho_diag(testCase)
% bug found when updating distancePointMesh function
%
% mesh is a single face, query point is above the largest side of the
% triangle.
v = [0 0 0; 10 0 0; 0 10 0];
f = [1 2 3];
point111 = [11 1 0];
d = distancePointMesh(point111, v, f, 'algorithm', 'vectorized');
testCase.assertEqual(sqrt(2), d, 0.0001);

function testSingleTriangleAllSubregions(testCase)
% The algorithm described in "David Eberly, 'Distance Between Point and
% Triangle in 3D', Geometric Tools, LLC, (1999)" splits the plane of the
% triangle into region 0 to 6. The regions here are not numbered as in the
% original paper of D. Eberly but as in distancePointTrimesh_vectorized.
%
%        /\ t
%        |
%   \ R5 |
%    \   |
%     \  |
%      \ |
%       \| P3
%        *
%        |\
%        | \
%   R1   |  \   R4
%        |   \
%        | R0 \
%        |     \
%        | P1   \ P2
%  ------*-------*------> s
%        |        \
%   R3   |   R2    \   R6
%
% In the code the regions split into subregions. The regions over the
% vertices into 5, some of which only occur if the vertex has an abuse
% angle. The regions over edges split into 3 subregions, some of which
% only occur if the neighbouring vertices have acute angles.
% The following code uses triangle and one set of test points and
% switches the vertex numbers, so that every subcase is tested.

v = [36 20 10; 40 20 10; 30 30 10];

% The regions for the testPoints in the comment apply for faces = [1 2 3].
% For faces = [2 3 1] R3 becomes R5 and R4 becomes R2.
% For faces = [3 1 2] R3 becomes R6 and R4 becomes R1.
% Also the terms for the subcases change, but still all subcases are
% tested.
testPoints = [36 22 10; %R0
    42 5 10; %R3 (d<0 & a<=d)
    38 5 10; %R3 (d<0 & a>d)
    30 5 10; %R3 (d>=0 & e>=0)
    5 18 10; %R3 (d>=0 & c<=-e)
    13 13 10; %R3 (c>-e & e<0)
    30 40 10; %R4 ((c+e)-(b+d)<=0)
    50 25 10; %R4 ((c+e)-(b+d)>=a-2*b+c)
    40 30 10]; %R4 ((c+e)-(b+d)>0 & (c+e)-(b+d)<a-2*b+c)

projResults = [36 22 10;
    40 20 10;
    38 20 10;
    36 20 10;
    30 30 10;
    33 25 10;
    30 30 10;
    40 20 10;
    35 25 10];

f = [1 2 3];

[dist_s, proj_s] = distancePointMesh(testPoints, v, f, 'algorithm', 'sequential');
testCase.assertEqual(proj_s, projResults, 'AbsTol', 1e-6)
[dist_v, proj_v] = distancePointMesh(testPoints, v, f, 'algorithm', 'vectorized');
testCase.assertEqual(proj_v, projResults, 'AbsTol', 1e-6)

f = [2 3 1];

[dist_s, proj_s] = distancePointMesh(testPoints, v, f, 'algorithm', 'sequential');
testCase.assertEqual(proj_s, projResults, 'AbsTol', 1e-6)
[dist_v, proj_v] = distancePointMesh(testPoints, v, f, 'algorithm', 'vectorized');
testCase.assertEqual(proj_v, projResults, 'AbsTol', 1e-6)

f = [3 1 2];

[dist_s, proj_s] = distancePointMesh(testPoints, v, f, 'algorithm', 'sequential');
testCase.assertEqual(proj_s, projResults, 'AbsTol', 1e-6)
[dist_v, proj_v] = distancePointMesh(testPoints, v, f, 'algorithm', 'vectorized');
testCase.assertEqual(proj_v, projResults, 'AbsTol', 1e-6)

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [y,fs,wmode,fidx]=readwav(filename,mode,nmax,nskip)
%READWAV  Read a .WAV format sound file [Y,FS,WMODE,FIDX]=(FILENAME,MODE,NMAX,NSKIP)
%
% Input Parameters:
%
%	FILENAME gives the name of the file (with optional .WAV extension) or alternatively
%                 can be the FIDX output from a previous call to READWAV
%	MODE		specifies the following (*=default):
%
%    Scaling: 's'    Auto scale to make data peak = +-1
%             'r'    Raw unscaled data (integer values)
%             'q'    Scaled to make 0dBm0 be unity mean square
%             'p' *	 Scaled to make +-1 equal full scale
%             'o'    Scale to bin centre rather than bin edge (e.g. 127 rather than 127.5 for 8 bit values)
%                     (can be combined with n+p,r,s modes)
%             'n'    Scale to negative peak rather than positive peak (e.g. 128.5 rather than 127.5 for 8 bit values)
%                     (can be combined with o+p,r,s modes)
%             'g'    Scale by the gain written by the "g" option in "writewav" to restore original level
%     Offset: 'y' *	 Correct for offset in <=8 bit PCM data
%             'z'    No offset correction
%   File I/O: 'f'    Do not close file on exit
%             'd'    Look in data directory: voicebox('dir_data')
%   Display;  'h'    Print header information
%             'w'    Plot waveform
%             'W'    Plot spectrogram
%             'a'    play audio (max 10 seconds)
%             'A'    play all audio even if very long
%
%	NMAX     maximum number of samples to read (or -1 for unlimited [default])
%	NSKIP    number of samples to skip from start of file
%               (or -1 to continue from previous read when FIDX is given instead of FILENAME [default])
%
% Output Parameters:
%
%	Y        data matrix of dimension (samples,channels)
%	FS       sample frequency in Hz
%	WMODE    mode string needed for WRITEWAV to recreate the data file
%	FIDX     Information row vector containing the element listed below.
%
%           (1)  file id
%		    (2)  current position in file
%           (3)  dataoff	byte offset in file to start of data
%           (4)  nsamp	number of samples
%           (5)  nchan	number of channels
%           (6)  nbyte	bytes per data value
%           (7)  bits	number of bits of precision
%           (8)  code	Data format: 1=PCM, 2=ADPCM, 3=floating point, 6=A-law, 7=Mu-law
%           (9)  fs	    sample frequency
%           (10) mask   channel mask
%           (11) gain   gain factor in dB
%
%   If no output parameters are specified, header information will be printed.
%
%   For stereo data, y(:,1) is the left channel and y(:,2) the right
%   The mask, if specified, is a bit field giving the channels present in the following order:
%   0=FL, 1=FR, 2=FC, 3=W, 4=BL, 5=BR, 6=FLC, 7=FRC, 8=BC, 9=SL, 10=SR, 11=TC, 12=TFL, 13=TFC, 14=TFR, 15=TBL, 16=TBC, 17=TBR
%   where F=front, L=left, C=centre, W=woofer (low frequency), B=back, LC=left of centre, RC=right of centre, S=side, T=top
%
%   See also WRITEWAV.

%   *** Note on scaling ***
%   If we want to scale signal values in the range +-1 to an integer in the
%   range [-128,127] then we have four plausible choices corresponding to
%   scale factors of (a) 127, (b) 127.5, (c) 128 or (d) 128.5 but each choice
%   has disadvantages.
%   For forward scaling: (c) and (d) cause clipping on inputs of +1.
%   For reverse scaling: (a) and (b) can generate output values < -1.
%   Any of these scalings can be selected via the mode input: (a) 'o', (b) default, (c) 'on', (d) 'n'

%	   Copyright (C) Mike Brookes 1998-2011
%      Version: $Id: readwav.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Bugs/suggestions:

if nargin<1
    error('Usage: [y,fs,wmode,fidx]=READWAV(filename,mode,nmax,nskip)'); end
if nargin<2
    mode='p';
else
    mode = [mode(:).' 'p'];
end
k=find((mode>='p') & (mode<='s'));
mno=all(mode~='o');                      % scale to input limits not output limits
sc=mode(k(1));
z=128*all(mode~='z');
info=zeros(1,11);
if ischar(filename)
    if any(mode=='d')
        filename=fullfile(voicebox('dir_data'),filename);
    end
    fid=fopen(filename,'rb','l');
    if fid == -1
        fn=[filename,'.wav'];
        fid=fopen(fn,'rb','l');
        if fid ~= -1
            filename=fn;
        end
    end
    if fid == -1
        error('Can''t open %s for input',filename);
    end
    info(1)=fid;
else
    info=filename;
    fid=info(1);
end
getdat= nargout>0 || any(lower(mode)=='w') || any(lower(mode)=='a');
mh=any(mode=='h') || ~getdat;
if ~info(3)
    fseek(fid,8,-1);						% read riff chunk
    header=fread(fid,4,'*char')';
    if ~strcmp(header,'WAVE')
        fclose(fid);
        error('File does not begin with a WAVE chunck');
    end
    if mh
        fprintf('\nWAVE file: %s\n',filename);
    end
    fmtlen=-1;
    datalen=-1;
    instlen=-1;
    factlen=-1;
    riffmt='e';  % default is original wave file format
    while datalen<0						% loop until FMT and DATA chuncks both found
        header=fread(fid,4,'*char');
        len=fread(fid,1,'ulong');
        if mh
            fprintf('  %s chunk: %d bytes\n',header,len);
        end
        if strcmp(header','fmt ')					% ******* found FMT chunk *********
            fmtlen=len;                             % remember the length
            if len>16
                riffmt='x';  % might be WAVEFORMATEX format
            end
            wavfmt=fread(fid,1,'short');			% format: 1=PCM, 6=A-law, 7-Mu-law
            info(8)=wavfmt;
            info(5)=fread(fid,1,'ushort');			% number of channels
            fs=fread(fid,1,'ulong');				% sample rate in Hz
            info(9)=fs;                             % sample rate in Hz
            rate=fread(fid,1,'ulong');				% average bytes per second (ignore)
            align=fread(fid,1,'ushort');			% block alignment in bytes (container size * #channels)
            bps=fread(fid,1,'ushort');			% bits per sample
            info(7)=bps;
            %             info(6)=ceil(info(7)/8);                % round up to a byte count
            info(6)=floor(align/info(5));                       % assume block size/channels = container size
            if info(8)==-2   % wave format extensible
                cb=fread(fid,1,'ushort');			% extra bytes must be >=22
                riffmt='X';                         % WAVEFORMATEXTENSIBLE format
                wfxsamp=fread(fid,1,'ushort');			% samples union
                if wfxsamp>0
                    info(7)=wfxsamp;     % valid bits per sample
                end
                info(10)=fread(fid,1,'ulong');				% channel mask
                wfxguida=fread(fid,1,'ulong');				% GUID
                wfxguidb=fread(fid,2,'ushort');				% GUID
                wfxguidc=fread(fid,8,'uchar');				% GUID
                if wfxguida<65536
                    info(8)=wfxguida; % turn it into normal wav format
                end
                fseek(fid,len-40,0);				    % skip to end of header
            else
                if align>0 && align<(info(6)+4)*info(5)
                    info(6)=ceil(align/info(5));
                end
                fseek(fid,len-16,0);				    % skip to end of header
            end
            if mh
                fmttypes={'?' 'PCM' 'ADPCM' 'IEEE-float' '?' '?' 'A-law' 'µ-law' '?'};
                fprintf('        Format: %d = %s',info(8),fmttypes{1+max(min(info(8),8),0)});
                if wavfmt==-2
                    fprintf(' (%08x-%04x-%04x-%02x%02x-%02x%02x%02x%02x%02x%02x)\n',wfxguida,wfxguidb,wfxguidc);
                else
                    fprintf('\n');
                end
                fprintf('        %d channels at %g kHz sample rate (%d kbytes/s)\n',info(5),fs/1000,rate/1000);
                fprintf('        Mask=%x:',info(10));
                spkpos={'FL' 'FR' 'FC' 'W' 'BL' 'BR' 'FLC' 'FRC' 'BC' 'SL' 'SR' 'TC' 'TFL' 'TFC' 'TFR' 'TBL' 'TBC' 'TBR'};
                for i=1:18
                    if mod(floor(info(10)*pow2(1-i)),2)
                        fprintf([' ' spkpos{i}]);
                    end
                end
                fprintf('\n        %d valid bits of %d per sample (%d byte block size)\n',info(7),bps,align);
            end
        elseif strcmp(header','fact')				% ******* found FACT chunk *********
            factlen=len;
            if len<4
                error('FACT chunk too short');
            end
            nsamp=fread(fid,1,'ulong');				% number of samples
            fseek(fid,len-4,0);                     % skip to end of header
            if mh
                fprintf('        %d samples = %.3g seconds\n',nsamp,nsamp/fs);
            end
        elseif strcmp(header','inst')				% ******* found INST chunk *********
            instlen=len;
            if len<7
                error('INST chunk too short');
            end
            inst=fread(fid,3,'schar');
            info(11)=double(inst(3));                          % gain in dB
            if mh
                fprintf('        Gain = %d dB\n',info(11));
            end
            fseek(fid,len-3,0);                     % skip to end of header
        elseif strcmp(header','data')				% ******* found DATA chunk *********
            if fmtlen<0
                fclose(fid);
                error('File %s does not contain a FMT chunck',filename);
            end
            if factlen>3 && nsamp >0
                info(4)=nsamp;   % take data length from FACT chunk
            else
                info(4) = fix(len/(info(6)*info(5)));  % number of samples
            end
            info(3)=ftell(fid);                    % start of data
            datalen=len;
            if mh
                fprintf('        %d samples x %d channels x %d bytes/samp',info(4:6));
                if prod(info(4:6))~=len
                    fprintf(' + %d padding bytes',len-prod(info(4:6)));
                end
                fprintf(' = %g sec\n',info(4)/fs);
            end
        else							% ******* found unwanted chunk *********
            fseek(fid,len,0);
        end
    end
else
    fs=info(9);
end
if nargin<4 || nskip<0
    nskip=info(2);  % resume at current file position
end

ksamples=info(4)-nskip; % number of samples remaining
if nargin>2
    if nmax>=0
        ksamples=min(nmax,ksamples);
    end
elseif ~getdat
    ksamples=min(5,ksamples); % just read 5 samples so we can print the first few data values
end
if ksamples>0
    info(2)=nskip+ksamples;
    fseek(fid,info(3)+info(6)*info(5)*nskip,-1);
    nsamples=info(5)*ksamples;
    if any(info(8)==3)  % floating point format
        pk=1;  % peak is 1
        switch info(6)
            case 4
                y=fread(fid,nsamples,'float32');
            case 8
                y=fread(fid,nsamples,'float64');
            otherwise
                error('cannot read %d-byte floating point numbers',info(6));
        end
    else
        if ~any(info(8)==[1 6 7])
            sc='r';  % read strange formats in raw integer mode
        end
        pk=pow2(0.5,8*info(6))*(1+(mno/2-all(mode~='n'))/pow2(0.5,info(7)));  % use modes o and n to determine effective peak
        switch info(6)
            case 1
                y=fread(fid,nsamples,'uchar');
                if info(8)==1
                    y=y-z;
                elseif info(8)==6
                    y=pcma2lin(y,213,1);
                    pk=4032+mno*64;
                elseif info(8)==7
                    y=pcmu2lin(y,1);
                    pk=8031+mno*128;
                end
            case 2
                y=fread(fid,nsamples,'short');
            case 3
                y=fread(fid,3*nsamples,'uchar');
                y=reshape(y,3,nsamples);
                y=([1 256 65536]*y-pow2(fix(pow2(y(3,:),-7)),24)).';
            case 4
                y=fread(fid,nsamples,'long');
            otherwise
                error('cannot read %d-byte integers',info(6));
        end
    end
    if sc ~= 'r'
        if sc=='s'
            sf=1/max(abs(y(:)));
        elseif sc=='p'
            sf=1/pk;
        else
            if info(8)==7
                sf=2.03761563/pk;
            else
                sf=2.03033976/pk;
            end
        end
        y=sf*y;
    else                             % mode = 'r' - output raw values
        if info(8)==1
            y=y*pow2(1,info(7)-8*info(6));  % shift to get the bits correct
        end
    end
    if any(mode=='g') && info(11)~=0
        y=y*10^(info(11)/20);   % scale by the gain
    end
    if info(5)>1
        y = reshape(y,info(5),ksamples).';
    end
else
    y=[];
end
if all(mode~='f')
    fclose(fid);
end
if nargout>2  % sort out the mode input for writing this format
    wmode=char([riffmt sc 'z'-z/128]);
    if factlen>0
        wmode=[wmode 'E'];
    end
    if info(6)>1 && info(6)<5
        cszc=' cCL';
        wmode=[wmode cszc(info(6))];
    end
    switch info(8)
        case 1                                    % PCM modes
            if ~mno
                wmode=[wmode 'o'];
            end
            if any(mode=='n')
                wmode=[wmode 'n'];
            end
            wmode=[wmode num2str(info(7))];
        case 3
            if info(7)<=32
                wmode = [wmode 'v'];
            else
                wmode = [wmode 'V'];
            end
        case 6
            wmode = [wmode 'a'];
        case 7
            wmode = [wmode 'u'];
    end
    fidx=info;
end
[ns,nchan]=size(y);
if mh && ns>0
    nsh=min(ns,5);  % print first few samples
    for i=1:nsh
        fprintf('        %d:',i);
        fprintf(' %.3g',y(i,:));
        fprintf('\n');
    end
end

if ns>0.01*fs
    if any(lower(mode)=='a')
        nsh=min(ns,10*fs+ns*any(mode=='A'));
        soundsc(y(1:nsh,1:min(nchan,2)),fs);
    end
    if any(mode=='W')
        spm='pJcbf ';
        if any(mode=='w')
            spm(end)='w';
        end
        clf;
        if nchan>1
            for i=nchan:-1:1
                subplot(nchan,1,i)
                spgrambw(y(:,i),fs,spm);
            end
        else
            spgrambw(y,fs,spm);
        end
    elseif any(mode=='w')
        clf;
        if nchan>1
            for i=nchan:-1:1
                subplot(nchan,1,i)
                plot((1:ns)/fs,y(:,i));
                ylabel(['Chan ' num2str(i)]);
                if i==nchan
                    xlabel('Time (s)');
                end
            end
        else
            plot((1:ns)/fs,y);
            xlabel('Time (s)');
        end

    end
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function inspect_ft_electroderealign

% MEM 3gb
% WALLTIME 00:10:00
% DEPENDENCY ft_electroderealign ft_read_mri ft_read_sens ft_prepare_mesh ft_warp_apply

%% load mri, segmentation and electrode definition
mri = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/Subject01.mri'));
load(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/headmodel_eeg/segmentedmri.mat'));
elec = ft_read_sens(dccnpath('/home/common/matlab/fieldtrip/template/electrode/standard_1020.elc'));
temp = ft_read_sens(dccnpath('/home/common/matlab/fieldtrip/template/electrode/standard_1005.elc'));

% create a bem and a fem mesh

cfg = [];
cfg.tissue = {'brain', 'skull', 'scalp'};
cfg.numvertices = [3000 2000 1000];
bem = ft_prepare_mesh(cfg, segmentedmri);

cfg = [];
cfg.method = 'hexahedral';
cfg.tissue = {'brain', 'skull', 'scalp'};
fem = ft_prepare_mesh(cfg, segmentedmri);

%% method: fiducial

nas = mri.hdr.fiducial.mri.nas;
lpa = mri.hdr.fiducial.mri.lpa;
rpa = mri.hdr.fiducial.mri.rpa;

% they are in voxels, hence need to be transformed to head coordinates
nas = ft_warp_apply(mri.transform, nas, 'homogenous');
lpa = ft_warp_apply(mri.transform, lpa, 'homogenous');
rpa = ft_warp_apply(mri.transform, rpa, 'homogenous');

fiducials.chanpos  = [nas; lpa; rpa];
fiducials.label    = {'Nz', 'LPA', 'RPA'};
fiducials.unit     = mri.unit;
fiducials.coordsys = mri.coordsys;

cfg = [];
cfg.method = 'fiducial';
cfg.template = fiducials;
cfg.elec = elec;
cfg.channel = 'all';
cfg.fiducial = {'Nz', 'LPA', 'RPA'};
elec_realigned1 = ft_electroderealign(cfg);

figure
ft_plot_sens(elec_realigned1, 'label', 'on');
ft_plot_axes(elec_realigned1, 'fontcolor', 'k');

%% method: template

cfg = [];
cfg.method = 'template';
cfg.template = elec_realigned1;
cfg.elec = elec;
elec_realigned2 = ft_electroderealign(cfg);

figure
ft_plot_sens(elec_realigned2, 'label', 'on');
ft_plot_axes(elec_realigned2, 'fontcolor', 'k');

%% method: interactive

% rotate    [0 0 -90]
% translate [35 0 40]

cfg = [];
cfg.method = 'interactive';
cfg.headshape = bem(3);
cfg.elec = elec;
elec_realigned3 = ft_electroderealign(cfg);

figure
ft_plot_sens(elec_realigned3, 'label', 'on');
ft_plot_axes(elec_realigned3, 'fontcolor', 'k');

%% method: headshape

cfg = [];
cfg.method = 'headshape';
cfg.headshape = bem(1);
cfg.elec = elec_realigned3;
elec_realigned4 = ft_electroderealign(cfg);

figure
ft_plot_sens(elec_realigned4, 'label', 'on');
ft_plot_axes(elec_realigned4, 'fontcolor', 'k');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
clear all
close all

disp('Random Effects demo with');
disp('two-region Neural Mass Models');
disp(' ');
disp('Initial conditions are known');
disp('Flow parameters are treated as random effects');

% Backward connection ?
back=1;

% Observation noise SD
sd=0.01;

% Number of subjects
N=3;

% Number of parameters
Np=2;

m=[1,1]';
C=0.1^2*eye(2);
w_true = spm_normrnd(m,C,N);
for n=1:N,
    [M{n},U{n}] = mci_nmm_struct(back,sd,Np);
    P=w_true(:,n);
    Y{n}.y = mci_nmm_gen(M{n},U{n},P);
end

% Assign init, flow and output parameters
assign.init_par='known';
assign.flow_par='random';
assign.out_par='known';
MCI.pout0=[]; % there are no 'output' parameters
MCI.assign=assign;

MCI.M=M; MCI.U=U; MCI.Y=Y;
MCI.verbose=1;
MCI.total_its=16;
MCI.rinit=0.25;

S.prior.P=Np;
S.prior.a=Np/2;
S.prior.B=M{1}.pC;
S.prior.beta=1;
S.prior.m=zeros(Np,1);
S.N=N;
MCI.S=S;

tic;
MCI = spm_mci_mfx_dynamic (MCI);
toc

disp('Mean true params:');
mwt=mean(w_true,2)

disp('MCI second level posterior:');
MCI.sm_mean

Np=M{1}.Npflow;
for j=1:Np,
    h=figure;
    set(h,'Name',sprintf('w(%d)',j));

    min_w=min(w_true(j,:));
    max_w=max(w_true(j,:));

    plot(w_true(j,:),MCI.sw_mean(j,:),'kx','MarkerSize',10);
    hold on
    plot([min_w max_w],[min_w max_w],'k-','LineWidth',2);
    set(gca,'FontSize',18);
    xlabel('True');
    ylabel('Estimated');
    grid on
end

e=MCI.sw_mean-w_true;
sse2=trace(e'*e);
disp(sprintf('Final first level SSE=%1.2f',sse2));

figure;
plot(MCI.sm');
grid on
xlabel('RFX iteration');
ylabel('Population Level Parameters');

h=figure;
set(h,'Name',sprintf('Subject Level'));
for j=1:Np,
    subplot(Np,1,j);
    plot(squeeze(MCI.sw(j,:,:))');
    grid on
    xlabel('RFX iteration');
    ylabel(sprintf('w(%d)',j));
end

author: spm


text:
function [response, x, s, f, v, q] = tapas_huge_bold( A, B, C, D, tau, ...
    kappa, epsilon, R, u, L, E_0, r_0, V_0, vartheta_0, alpha, gamma, TR, ...
    TE, dt)
% Integrates the DCM forward equations to generate the predicted fMRI bold
% time series.
%
% INPUTS:
%   A, B, C, D - DCM connectivity matrices.
%   tau        - Venous transit time.
%   kappa      - Decay of vasodilatory signal.
%   epsilon    - Ratio of intra- and extravascular signal.
%   R          - Number of regions.
%   u          - Experimental stimuli.
%   L          - Number of experimental stimuli.
%   E_0        - Resting oxygen extraction fraction.
%   r_0        - Slope of intravascular relaxation rate.
%   V_0        - Resting venous volume.
%   vartheta_0 - Frequency offset at the outer surface of magnetized
%                vessels (Hz).
%   alpha      - Grubb's exponent.
%   gamma      - rate constant of feedback regulation.
%   TR         - Repetition time.
%   TE         - Echo time.
%   dt         - Sampling interval of inputs.
%
% OUTPUTS:
%   response - matrix of predicted response for each region
%                  (column-wise)
%   x        - time series of neuronal states
%   s        - time series of vasodilatory signal
%   f1       - time series of flow
%   v1       - time series of blood volume
%   q1       - time series of deoxyhemoglobin content.
%

%
% REFERENCE:
%   Klaas Enno Stephan, Nikolaus Weiskopf, Peter M. Drysdale, Peter A.
%   Robinson, Karl J. Friston (2007). Comparing hemodynamic models with
%   DCM. NeuroImage, 38: 387-401
%
% https://doi.org/10.1016/j.neuroimage.2007.07.040
%

% Author: Yu Yao (yao@biomed.ee.ethz.ch), Sudhir Shankar Raman
% Copyright (C) 2019 Translational Neuromodeling Unit
%                    Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of TAPAS, which is released under the terms of the GNU
% General Public Licence (GPL), version 3. For further details, see
% <https://www.gnu.org/licenses/>.
%
% This software is provided "as is", without warranty of any kind, express
% or implied, including, but not limited to the warranties of
% merchantability, fitness for a particular purpose and non-infringement.
%
% This software is intended for research only. Do not use for clinical
% purpose. Please note that this toolbox is under active development.
% Considerable changes may occur in future releases. For support please
% refer to:
% https://github.com/translationalneuromodeling/tapas/issues
%

nt = size(u, 1);
rSmp = TR/dt;
C = C'/rSmp;
if isempty(D)
    D = zeros(R, R, R);
end

tau     = 2.*exp(tau);
kappa   = .64.*exp(kappa);
epsilon = exp(epsilon);

% resting oxygen extraction fraction
E_0 = repmat(E_0, 1, R);

k1  = 4.3*vartheta_0*TE*E_0;
k2  = epsilon.*(r_0*E_0*TE);
k3  = 1 - epsilon;

% Integrate the dynamical system
[x,s,f,v,q]  = tapas_huge_int_euler(...
    A',...
    full(u*C),...
    full(u),...
    permute(B,[2 1 3]),...
    permute(D,[2 1 3]),...
    E_0,...
    1/alpha,...
    tau,...
    gamma,...
    kappa,...
    [dt,nt,R,L,0,any(B(:)),any(D(:))]);

% generate the BOLD response
response = V_0*( ...
    bsxfun(@times,k1,...
        (1 - (q(rSmp:rSmp:end,:)))) +...
    bsxfun(@times,k2,...
        (1 - (q(rSmp:rSmp:end,:)./...
        v(rSmp:rSmp:end,:)))) +...
    bsxfun(@times,k3,...
        (1-v(rSmp:rSmp:end,:))));

% demean response
response = bsxfun(@minus, response, mean(response, 1));

end

author: translationalneuromodeling

repo: tapas


text:
function view_components(Y,A,C,b,f,Cn,options)

% plot spatial components and temporal traces against filtered background
% Y:        raw data
% A:        spatial footprints
% C:        temporal components
% b:        spatial background
% f:        temporal background
% Cn:       background image (default: mean image)
% options:  options structure

% Written by:
% Eftychios A. Pnevmatikakis, Simons Foundation, 2015

defoptions = CNMFSetParms;
memmaped = isobject(Y);
if nargin < 7 || isempty(options); options = defoptions; end
if ~isfield(options,'d1') || isempty(options.d1); d1 = input('What is the total number of rows? \n'); else d1 = options.d1; end          % # of rows
if ~isfield(options,'d2') || isempty(options.d2); d2 = input('What is the total number of columns? \n'); else d2 = options.d2; end          % # of columns
% if ~isfield(options,'normalize') || isempty(options.normalize); options.normalize = ones(size(A,1),1); end
%     sn = options.normalize;
if ~isfield(options,'plot_df') || isempty(options.plot_df); options.df = defoptions.plot_df; end
    plot_df = options.plot_df;
if ~isfield(options,'make_gif') || isempty(options.make_gif); options.make_gif = defoptions.make_gif; end
    make_gif = options.make_gif;
if ~isfield(options,'save_avi') || isempty(options.save_avi); options.save_avi = defoptions.save_avi; end
    save_avi = options.save_avi;
if ~isfield(options,'sx') || isempty(options.sx); options.sx = defoptions.sx; end
    sx = min([options.sx,floor(d1/2),floor(d2/2)]);
if ~isfield(options,'pause_time') || isempty(options.pause_time); options.pause_time = defoptions.pause_time; end
    pause_time = options.pause_time;
if isfield(options,'name') && ~isempty(options.name);
    name = [options.name,'_components'];
else
    name = [defoptions.name,'_components'];
end
if ~isfield(options,'full_A') || isempty(options.full_A); full_A = defoptions.full_A; else full_A = options.full_A; end

T = size(C,2);
if ndims(Y) == 3
    Y = reshape(Y,d1*d2,T);
end
if nargin < 6 || isempty(Cn);
    Cn = reshape(mean(Y,2),d1,d2);
end

nA = sqrt(sum(A.^2))';
[K,~] = size(C);
A = A/spdiags(nA,0,K,K);    % normalize spatial components to unit energy
C = spdiags(nA,0,K,K)*C;

nr = size(A,2);     % number of ROIs
nb = size(f,1);     % number of background components
%nA = full(sum(A.^2))';  % energy of each row
%Y_r = spdiags(nA,0,nr,nr)\(A'*Y- (A'*A)*C - (A'*full(b))*f) + C;
step = 5e3;
if memmaped
    AY = zeros(K,T);
    for i = 1:step:d
        AY = AY + A(i:min(i+step-1,d),:)'*double(Y.Yr(i:min(i+step-1,d),:));
    end
else
    if issparse(A) && isa(Y,'single')
        if full_A
            AY = full(A)'*Y;
        else
            AY = A'*double(Y);
        end
    else
        AY = A'*Y;
    end
end
Y_r = (AY- (A'*A)*C - full(A'*double(b))*f) + C;

if plot_df
    %[~,Df] = extract_DF_F(Y,[A,b],[C;f],size(A,2)+1);
    [~,Df] = extract_DF_F(Y,A,C,[],options);
else
    Df = ones(size(A,2)+1,1);
end

if save_avi
    vidObj = VideoWriter([name,'.avi']);
    set(vidObj,'FrameRate',1);
    open(vidObj);
end
thr = 0.9;
fig = figure;
    set(gcf,'Position',2*[300,300,960,480]);
    set(gcf,'PaperPosition',2*[300,300,960,480]);
    int_x = zeros(nr,2*sx);
    int_y = zeros(nr,2*sx);
    cm = com(A,d1,d2);
for i = 1:nr+nb
    if i <= nr
    subplot(3,2,5);
        Atemp = reshape(A(:,i),d1,d2);
        int_x(i,:) = round(cm(i,1)) + (-(sx-1):sx);
        if int_x(i,1)<1
            int_x(i,:) = int_x(i,:) + 1 - int_x(i,1);
        end
        if int_x(i,end)>d1
            int_x(i,:) = int_x(i,:) - (int_x(i,end)-d1);
        end
        int_y(i,:) = round(cm(i,2)) + (-(sx-1):sx);
        if int_y(i,1)<1
            int_y(i,:) = int_y(i,:) + 1 - int_y(i,1);
        end
        if int_y(i,end)>d2
            int_y(i,:) = int_y(i,:) - (int_y(i,end)-d2);
        end
        Atemp = Atemp(int_x(i,:),int_y(i,:));
        imagesc(int_x(i,:),int_y(i,:),Atemp); axis square;
    end
    subplot(3,2,[1,3]);
    if i <= nr
        imagesc(2*Cn); axis equal; axis tight; axis off; hold on;
        A_temp = full(reshape(A(:,i),d1,d2));
        A_temp = medfilt2(A_temp,[3,3]);
        A_temp = A_temp(:)/norm(A_temp(:));
        [temp,ind] = sort(A_temp(:).^2,'ascend');
        temp =  cumsum(temp);
        ff = find(temp > (1-thr)*temp(end),1,'first');
        if ~isempty(ff)
            [~,ww] = contour(reshape(A_temp,d1,d2),[0,0]+A_temp(ind(ff)),'LineColor','k');
            ww.LineWidth = 2;
        end
        title(sprintf('Component %i ',i),'fontsize',16,'fontweight','bold'); drawnow; %pause;
    else
        imagesc(reshape(b(:,i-nr),d1,d2)); axis equal; axis tight;
        title('Background component','fontsize',16,'fontweight','bold'); drawnow;
    end
    subplot(3,2,[2,4,6]);
    if i <= nr
        plot(1:T,Y_r(i,:)/Df(i),'linewidth',2); hold all; plot(1:T,C(i,:)/Df(i),'linewidth',2);
        if plot_df
            title(sprintf('Component %i (calcium DF/F value)',i),'fontsize',16,'fontweight','bold');
        else
            title(sprintf('Component %i (calcium raw value)',i),'fontsize',16,'fontweight','bold');
        end
        leg = legend('Raw trace (filtered)','Inferred');
        set(leg,'FontSize',14,'FontWeight','bold');
        drawnow;
        hold off;
        if make_gif
            frame = getframe(fig); %getframe(1);
              im = frame2im(frame);
              [imind,clm] = rgb2ind(im,256);
              if i == 1;
                  imwrite(imind,clm,[name,'.gif'],'gif', 'Loopcount',inf);
              else
                  imwrite(imind,clm,[name,'.gif'],'gif','WriteMode','append');
              end
        else
            if i < nr+nb && ~save_avi
                fprintf('component %i. Press any key to continue.. \n', i);
                if pause_time == Inf;
                    pause;
                else
                    pause(pause_time);
                end
            end
        end
    else
        plot(1:T,f(i-nr,:)); title('Background activity','fontsize',16,'fontweight','bold');
        drawnow;
        if make_gif
            frame = getframe(fig); %getframe(1);
              im = frame2im(frame);
              [imind,clm] = rgb2ind(im,256);
              if i == 1;
                  imwrite(imind,clm,[name,'.gif'],'gif', 'Loopcount',inf);
              else
                  imwrite(imind,clm,[name,'.gif'],'gif','WriteMode','append');
              end
        else
            if i < nr+nb && ~save_avi
                fprintf('background component %i. Press any key to continue.. \n', i-nr);
                if pause_time == Inf;
                    pause;
                else
                    pause(pause_time);
                end
            end
        end
    end
    if save_avi
        currFrame = getframe(fig);
        writeVideo(vidObj,currFrame);
    else
        pause(0.05);
    end

end

if save_avi
    close(vidObj);
end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function [id, stop, val] = list(item, spec, tropts, fn)

% function [id, stop, val] = list(item, spec, tropts, fn)
% Find items in a cfg tree rooted at item that match a specification spec.
% By default, the filled configuration tree is searched (i.e. the
% val-branches of cfg_repeat and cfg_choice nodes).
% See MATCH for help about spec data structure.
%
% Traversal options
% struct with fields
% stopspec - match spec to stop traversal
% dflag    - traverse val or values tree
% clvl     - current level in tree
% mlvl     - maximum level to traverse - range 1 (top level only) to
%            Inf (all levels)
% cnt      - #items found so far
% mcnt    - max #items to find
% List will stop descending into subtrees if one of the conditions
% following conditions are met: item matches stopspec, clvl >= mlvl, cnt >=
% mcnt. Flag stop is true for nodes where traversal has stopped
% (i.e. items where tropts has stopped further traversal).
%
% A cell list of subsref ids to matching nodes will be returned. The id of
% this node is returned before the id of its matching children.
% If the root node of the tree matches, the first id returned will be an
% empty substruct.
% If a cell list of fieldnames is given, then the contents of these fields
% will be returned in the cell array val. If one of the fields does not
% exist, a cell with an empty entry will be returned.
% There are five pseudo-fieldnames which allow to obtain information useful
% to build e.g. a user interface for cfg trees:
% 'class' - returns the class of the current item
% 'level' - returns the level in the tree. Since data is collected
%           pre-order, children are listed after their parents. Identical
%           levels of subsequent nodes denote siblings, whereas decreasing
%           levels of subsequent nodes denote siblings of the parent node.
% 'all_set' - return all_set status of subtree rooted at item, regardless
%             whether list will descend into it or not
% 'all_set_item' - return all_set_item status of current node (i.e. whether
%                  all integrity conditions for this node are fulfilled).
%                  For in-tree nodes this can be different from all_set.
% 'showdoc' - calls showmydoc to display the help text and option hints for
%             the current item (without recursive calls for .val/.values
%             items).
%
% This function is identical for all cfg_intree classes.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: list.m 5678 2013-10-11 14:58:04Z volkmar $

rev = '$Rev: 5678 $'; %#ok

if match(item, spec)
    id = {struct('type', {}, 'subs', {})};
    tropts.cnt = tropts.cnt + 1;
    stop = false;
    if nargin > 3
        val = cell(1,numel(fn));
        for k = 1:numel(fn)
            switch fn{k}
                case 'class'
                    val{k} = {class(item)};
                case 'level'
                    val{k} = {tropts.clvl};
                case 'all_set'
                    val{k} = {all_set(item)};
                case 'all_set_item'
                    val{k} = {all_set_item(item)};
                case 'showdoc'
                    val{k} = {showmydoc(item,'')};
                case fieldnames(item)
                    val{k} = {subsref(item, substruct('.', fn{k}))};
                otherwise
                    val{k} = {{}};
            end;
        end;
    else
        val = {};
    end;
    matched = true;
else
    id = {};
    stop = [];
    if nargin > 3
        val = cell(size(fn));
        [val{:}] = deal({});
    else
        val = {};
    end;
    matched = false;
end;

if (isempty(tropts.stopspec) || ~match(item, tropts.stopspec)) && (tropts.cnt < tropts.mcnt) && (tropts.clvl < tropts.mlvl)
    tname = treepart(item, tropts.dflag);
    % subscripts to child nodes
    idsubs = substruct('.', tname, '{}', {});
    citems = subsref(item, idsubs(1));
    tropts.clvl = tropts.clvl + 1;
    for k = 1:numel(citems)
        if nargin > 3
            [id1, stop1, val1] = list(citems{k}, spec, tropts, fn);
            for l = 1:numel(fn)
                val{l} = [val{l}(:); val1{l}(:)]';
            end;
        else
            [id1, stop1] = list(citems{k}, spec, tropts);
        end;
        if ~isempty(id1)
            idsubs(2).subs = {k};
            for l = 1:numel(id1)
                id = [id(:); {[idsubs id1{l}]}]';
            end;
            stop = [stop stop1];
            tropts.cnt = tropts.cnt + numel(id1);
            if tropts.cnt >= tropts.mcnt
                if matched
                    stop(1) = true;
                end;
                break;
            end;
        end;
    end;
else
    if matched
        stop = true;
    end;
end;


text:
%kvadil 'Conditionally Translation Invariant Anti-Dilation'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vadil.pane file
%
% Parameters:
% InputFile: i 'Input Image', required: 'input image'
% InputFile: str 'Struct Element', required: 'Structuring Element'
% OutputFile: o 'Output Image ', required: 'resulting output image'
%
% Example: o = kvadil({i, str}, {'i','';'str','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vadil - Conditionally Translation Invariant Anti-Dilation
%
%  DESCRIPTION
% .I vadil
% will do the anti-dilation of an input image by a set structuring element.
%
% The output image data type will be the same as that of the highest data type
% of the input image.
%
%
%
%  EXAMPLES
% vadil -i ball.xv -str b272.str -o outimage.xv
%
% Will do the anti-dilation of image "ball.xv" by the structuring set b272.str, with the resulting image written to "outimage.xv".
%
%  "SEE ALSO"
% vcast1(1)
%
%  RESTRICTIONS
% .I vadil
% can be defined for all data types supported by Khoros, but at the moment it has been implemented just for the bit and unsigned char types.
% The structuring elements are subsets of the 3x3 matrix and the origin is always at the center of this matrix.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993-1997 Junior Barrera, Roberto Lotufo.  All rights reserved.
%

function varargout = kvadil(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kvadil(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'str', '__input';'o', '__output'};
maxval={0,0,0};
minval={0,0,0};
istoggle=[0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'vadil"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%Kalman Demo 2 (Mauna Loa Periodic CO2 Readings)
%
%  In this demo we apply the state space inference methods to the
%  well-known time series data consisting of atmospheric CO2
%  concentration readings in parts per million (ppm) by volume from
%  air samples collected at the Mauna Loa observatory, Hawaii (see
%  [1] for details and further references).
%
%  The benefit from the state space formulation is that the
%  computational complexity is linear with respect to the number
%  of data points. Due to efficient matrix solvers in Matlab
%  (favoring the traditional GP solution over sequential looping),
%  the advantages in speed start to show in datasets with thousands
%  of data points.
%
%  The take-home message from this demo is that the GP can be
%  set up exactly as any other GP regression model in GPstuff,
%  and then solved by Kalman filtering methods by specifying the
%  'type' in the GP structure to be 'KALMAN'.
%
%  The methods in this demo are based on the paper:
%
%  [1] Arno Solin and Simo Sarkka (2014). Explicit link between periodic
%      covariance functions and state space models. Accepted for
%      publication in Proceedings of the Seventeenth International
%      Conference on Artifcial Intelligence and Statistics (AISTATS 2014).
%
% Copyright (c) 2014 Arno Solin and Jukka Koskenranta

% 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.

%% Load the data

  % Load the data
  S = which('demo_kalman2.m');
  L = strrep(S,'demo_kalman2.m','demodata/maunaloa_data.txt');

  % Set training data (x = time in years, y = CO2 observations)
  data = load(L);
  y = data(:, 2:13)';  y=y(:);
  x = (data(1,1):1/12:data(end,1)+11/12)';

  % Remove data points with missing information
  x = x(y>0); y = y(y>0);

  % Make data zero mean
  ymean = mean(y); y = y-ymean;

  % Show original data
  figure(1); clf
    plot(x,y+ymean,'-k')
    xlabel('Time (year)');
    ylabel('CO_2 concentration (PPM)')
    title('Kalman Demo 2  (Mauna Loa Periodic CO2 Readings) - Data')

%% Set up GP model

  % Noise variance prior
  ps2 = prior_logunif();

  % The likelihood model
  lik = lik_gaussian('sigma2', 1, 'sigma2_prior', ps2);

  % Covariance function hyperparameter priors
  pl = prior_logunif();
  pm = prior_logunif();

  % A squared exponential covariance function
  % to deal with the smooth long term effects
  gpcf1 = gpcf_sexp('lengthScale', 100, 'magnSigma2', 5000, ...
                    'lengthScale_prior',pl,'magnSigma2_prior',pl);

  % A quasi-periodic covariance function deals with peridic
  % variation in the data. The quasi-periodic covariance function
  % is a product of a periodic covariance function and a squared
  % exponential.
  gpcf2 = gpcf_periodic('magnSigma2',1,'lengthScale',1,'period',1, ...
                        'decay',1,'lengthScale_sexp',100, ...
                        'lengthScale_prior',pl,'magnSigma2_prior',pl, ...
                        'lengthScale_sexp_prior',pl);

  % A Matern52 covariance function deals with short term
  % non-periodic effects that remain otherwise unexplained
  gpcf3 = gpcf_matern52('lengthScale', 10, 'magnSigma2', 10, ...
                        'lengthScale_prior',pl,'magnSigma2_prior',pl);

  % Finally create the GP structure
  gp = gp_set('lik', lik, 'cf', {gpcf1,gpcf2,gpcf3});

  % Set type to KALMAN
  gp = gp_set(gp,'type','KALMAN');

%% Optimize hyperparameters and predict

  % Optimization parameters
  opt=optimset('TolFun',1e-4,'TolX',1e-4,'Display','iter');

  % Find hyperparameters by optimization (BFGS)
  gp=gp_optim(gp,x,y,'opt',opt,'optimf',@fminlbfgs);

  % Set the test points
  xt = (x(end):1/12:x(end)+10)';

  % Predict values
  [Eft,Varft] =  gp_pred(gp, x, y,'xt',xt);

  % Also predict the latent components separately
  [Eft1, Varft1] = gp_pred(gp, x, y, x, 'predcf', [1 3]);
  [Eft2, Varft2] = gp_pred(gp, x, y, x, 'predcf', [2]);

%% Visualize results

  % Plot
  figure(2); clf; hold on

    % Plot the 95% confidence interval of the predictions
    p=patch([xt; flipud(xt)], ...
       [ymean + Eft + 1.96*sqrt(Varft); ...
        flipud(ymean + Eft - 1.96*sqrt(Varft))],[0.9,0.9,0.9]);
    set(p,'EdgeColor','none')

    % Plot observations
    plot(x,ymean+y,'.k','MarkerSize',5)

    % Labels and legends
    title('Kalman Demo 2 (Mauna Loa Periodic CO2 Readings)')
    xlabel('Time (years)');
    ylabel('CO_2 concentration (PPM)')
    legend('95% confidence region', ...
           'Monthly average measurements',...
           'Location', 'NorthWest');

    % Axis options
    box on; axis tight; set(gca,'Layer','top')

%% Show components separately

  figure(3); clf;
  subplot(211)

    plot(x,Eft1,'-k')

    title('Long-term trend and short-scale variation')
    xlabel('Time (years)');
    ylabel('Effect on CO_2 concentration (PPM)')

    axis tight

  subplot(212)

    plot(x,Eft2,'-k')

    title('Quasi-periodic effect')
    xlabel('Time (years)');
    ylabel('Effect on CO_2 concentration (PPM)')

    axis tight

author: gpstuff-dev

repo: gpstuff

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


text:
function pass = test_emptyObjects( pref )
% For empty chebfun2 objects, does each command deal with them
% appropriately?

% Check things work for empty chebfun2s.
f = chebfun2();
try
    f + f;
    2*f;
    f*2;
    f.^2;
    2.^f;
    f.^f;
    sqrt(f);
    sum(f);
    integral2(f);
    % max(f)
    norm(f);
    squeeze(f);
    diff(f);
    cos(f);
    sin(f);
    sinh(f);
    f.^f + f;
    diag(f);
    trace(f);
    mean(f);
    minandmax2(f);
    median(f);
    flipud(f);
    flipdim(f,1);
    pass = 1;
catch
    pass = 0;
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% Demo for Structured Edge Detector (please see readme.txt first).
addpath(genpath('../piotr_toolbox'));

%% set opts for training (see edgesTrain.m)
opts=edgesTrain();                % default options (good settings)
opts.modelDir='models/';          % model will be in models/forest
opts.modelFnm='modelBsds';        % model name
opts.nPos=5e5; opts.nNeg=5e5;     % decrease to speedup training
opts.useParfor=0;                 % parallelize if sufficient memory

%% train edge detector (~20m/8Gb per tree, proportional to nPos/nNeg)
tic, model=edgesTrain(opts); toc; % will load model if already trained

%% set detection parameters (can set after training)
model.opts.multiscale=0;          % for top accuracy set multiscale=1
model.opts.sharpen=2;             % for top speed set sharpen=0
model.opts.nTreesEval=4;          % for top speed set nTreesEval=1
model.opts.nThreads=4;            % max number threads for evaluation
model.opts.nms=0;                 % set to true to enable nms

%% evaluate edge detector on BSDS500 (see edgesEval.m)
if(0), edgesEval( model, 'show',1, 'name','' ); end

%% detect edge and visualize results
I = imread('peppers.png');
tic, E=edgesDetect(I,model); toc
figure(1); im(I); figure(2); im(1-E);

repo: interpretableCNN


text:
function [ t_start, y_start ] = p33_start ( neqn )

%*****************************************************************************80
%
%% P33_START returns the starting point for problem p33.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_START, Y_START(NEQN,1), the initial data.
%
  y_start = zeros ( neqn, 1 );

  t_start = 0.0;
  y_start = [ 2.0; 2.0 ];

  return
end


text:
function c=children(t,j)
%CHILDEREN Child nodes.
%   C=CHILDREN(T) returns an N-by-2 array C containing the numbers of the
%   child nodes for each node in the tree T, where N is the number of
%   nodes in the tree.  The children for leaf nodes (those with no children)
%   are returned as 0.
%
%   P=CHILDREN(T,J) takes an array J of node numbers and returns the children
%   for the specified nodes.
%
%   See also CLASSREGTREE, CLASSREGTREE/NUMNODES, CLASSREGTREE/PARENT.

%   Copyright 2006-2007 The MathWorks, Inc.
%   $Revision: 1.1.6.2 $  $Date: 2007/02/15 21:48:00 $

if nargin>=2 && ~validatenodes(t,j)
    error('stats:classregtree:children:InvalidNode',...
          'J must be an array of node numbers or a logical array of the proper size.');
end

if nargin<2
    c = t.children;
else
    c = t.children(j,:);
end

author: Cloud-CV

repo: object-proposals


text:
function vt = volToTalairach(volCoords, v2t)
% talairachCoords = volToTalairach(volCoords, v2t)
%
% PURPOSE:
%   Converts arbitrary coordinates (volCoords) to Talairach space, given
%   the transform info in the v2t struct. The coords must be in 'n X 3'
%   row-vector form. The v2t struct must have a 4x4 affine transform matrix
%
%       v2t.transRot
%
%   And 7 scale factors:
%
%       v2t.superiorAcScale
%       v2t.inferiorAcScale
%       v2t.rightAcScale
%       v2t.leftAcScale
%       v2t.anteriorAcScale
%       v2t.betweenAcPcScale
%       v2t.posteriorPcScale
%
%   See computeTalairach.m for an example of how to compute these.
%
% HISTORY:
%   2001.09.14 RFD (bob@white.stanford.edu) Wrote it, with help from
%              Alex Wade and Jochem Rieger.

% MrLoadRet coords are axial,coronal,sagittal, but we need coronal,axial,sagittal.
% Oh well...

% We have to make the coords homogeneous (ie. 4d)
vt = ones(size(volCoords,1),size(volCoords,2)+1);
vt(:,1:3) = [volCoords(:,2),volCoords(:,1),volCoords(:,3)];

% now, apply translation and rotation
vt = vt * v2t.transRot;

% v shoule now be on the talairach axes:
%   X = sagittal slice (right is +)
%   Y = coronal slice (anterior is +)
%   Z = axial slice (superior is +)
%
% We just need to apply the appropriate scale factor
for(ii=1:size(vt,1))
    if(vt(ii,1)>0)
        vt(ii,1) = vt(ii,1) * v2t.rightAcScale;
    else
        vt(ii,1) = vt(ii,1) * v2t.leftAcScale;
    end
    if(vt(ii,3)>0)
        vt(ii,3) = vt(ii,3) * v2t.superiorAcScale;
    else
        vt(ii,3) = vt(ii,3) * v2t.inferiorAcScale;
    end
     if(vt(ii,2)>0)
        vt(ii,2) = vt(ii,2) * v2t.anteriorAcScale;
    else
        acpc = vt(ii,2) * v2t.betweenAcPcScale;
        % we have to apply separate scale factors if it goes beyond the PC
        if(acpc<-24)
            beyondPC = vt(ii,2) + 24/v2t.betweenAcPcScale;
            vt(ii,2) = -24 + beyondPC*v2t.posteriorPcScale;
        else
            vt(ii,2) = acpc;
        end
    end
end
vt = vt(:,1:3);
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [A_hat, E_hat, iter] = exact_alm_rpca(D, lambda, tol, maxIter)

% Oct 2009
% This matlab code implements the augmented Lagrange multiplier method for
% Robust PCA.
%
% D - m x n matrix of observations/data (required input)
%
% lambda - weight on sparse error term in the cost function
%
% tol - tolerance for stopping criterion.
%     - DEFAULT 1e-7 if omitted or -1.
%
% maxIter - maximum number of iterations
%         - DEFAULT 1000, if omitted or -1.
%
% Initialize A,E,Y,u
% while ~converged
%   minimize
%     L(A,E,Y,u) = |A|_* + lambda * |E|_1 + <Y,D-A-E> + mu/2 * |D-A-E|_F^2;
%   Y = Y + \mu * (D - A - E);
%   \mu = \rho * \mu;
% end
%
% Minming Chen, October 2009. Questions? v-minmch@microsoft.com ;
% Arvind Ganesh (abalasu2@illinois.edu)
%
% Copyright: Perception and Decision Laboratory, University of Illinois, Urbana-Champaign
%            Microsoft Research Asia, Beijing

[m,n] = size(D);

if(nargin < 2) lambda = 1 / sqrt(m); end
if(nargin < 3) tol = 1e-7; elseif(tol == -1) tol = 1e-7; end
if(nargin < 4) maxIter = 1000; elseif(maxIter == -1) maxIter = 1000; end

% initialize
Y = sign(D);
norm_two = norm(Y, 2); %norm_two = lansvd(Y, 1, 'L');
norm_inf = norm( Y(:), inf) / lambda;
dual_norm = max(norm_two, norm_inf);
Y = Y / dual_norm;

A_hat = zeros( m, n);
E_hat = zeros( m, n);
dnorm = norm(D, 'fro');
tolProj = 1e-6 * dnorm;
total_svd = 0;
mu = .5/norm_two; % this one can be tuned
rho = 6;          % this one can be tuned

iter = 0;
converged = false;
stopCriterion = 1;
sv = 5;
svp = sv;
while ~converged
    iter = iter + 1;
    % solve the primal problem by alternative projection
    primal_converged = false;
    primal_iter = 0;
    sv = sv + round(n * 0.1);
    while primal_converged == false
        temp_T = D - A_hat + (1/mu)*Y;
        temp_E = max( temp_T - lambda/mu,0) + min( temp_T + lambda/mu,0);

        %if choosvd(n, sv) == 1
        %    [U S V] = lansvd(D - temp_E + (1/mu)*Y, sv, 'L');
        %else
            %[U,S,V] = svd(D - temp_E + (1/mu)*Y, 'econ');
            [U,S,V] = svdecon(D - temp_E + (1/mu)*Y); % fastest
        %end

        diagS = diag(S);
        svp = length(find(diagS > 1/mu));
        if svp < sv
            sv = min(svp + 1, n);
        else
            sv = min(svp + round(0.05*n), n);
        end
        temp_A = U(:,1:svp)*diag(diagS(1:svp)-1/mu)*V(:,1:svp)';

        if norm(A_hat - temp_A, 'fro') < tolProj && norm(E_hat - temp_E, 'fro') < tolProj
            primal_converged = true;
        end
        A_hat = temp_A;
        E_hat = temp_E;
        primal_iter = primal_iter + 1;
        total_svd = total_svd + 1;
    end

    Z = D - A_hat - E_hat;
    Y = Y + mu*Z;
    mu = rho * mu;

    %% stop Criterion
    stopCriterion = norm(Z, 'fro') / dnorm;
    if stopCriterion < tol
        converged = true;
    end

    disp(['Iteration' num2str(iter) ' #svd ' num2str(total_svd) ' r(A) ' num2str(svp)...
        ' |E|_0 ' num2str(length(find(abs(E_hat)>0)))...
        ' stopCriterion ' num2str(stopCriterion)]);

    if ~converged && iter >= maxIter
        disp('Maximum iterations reached') ;
        converged = 1 ;
    end
end

if nargin == 5
    fclose(fid);
end

author: andrewssobral

repo: lrslibrary


text:
function [all_theta] = oneVsAll(X, y, num_labels, lambda)
%ONEVSALL trains multiple logistic regression classifiers and returns all
%the classifiers in a matrix all_theta, where the i-th row of all_theta
%corresponds to the classifier for label i
%   [all_theta] = ONEVSALL(X, y, num_labels, lambda) trains num_labels
%   logisitc regression classifiers and returns each of these classifiers
%   in a matrix all_theta, where the i-th row of all_theta corresponds
%   to the classifier for label i

% Some useful variables
m = size(X, 1);
n = size(X, 2);

% You need to return the following variables correctly
all_theta = zeros(num_labels, n + 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: You should complete the following code to train num_labels
%               logistic regression classifiers with regularization
%               parameter lambda.
%
% Hint: theta(:) will return a column vector.
%
% Hint: You can use y == c to obtain a vector of 1's and 0's that tell use
%       whether the ground truth is true/false for this class.
%
% Note: For this assignment, we recommend using fmincg to optimize the cost
%       function. It is okay to use a for-loop (for c = 1:num_labels) to
%       loop over the different classes.
%
%       fmincg works similarly to fminunc, but is more efficient when we
%       are dealing with large number of parameters.
%
% Example Code for fmincg:
%
%     % Set Initial theta
%     initial_theta = zeros(n + 1, 1);
%
%     % Set options for fminunc
%     options = optimset('GradObj', 'on', 'MaxIter', 50);
%
%     % Run fmincg to obtain the optimal theta
%     % This function will return theta and the cost
%     [theta] = ...
%         fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), ...
%                 initial_theta, options);
%

	for c = 1 : num_labels,
		initial_theta = zeros(n + 1 , 1);
		options = optimset('GradObj' , 'on' , 'MaxIter' , 50);
		[theta] = ...
			fmincg(@(t)(lrCostFunction(t , X , (y == c) , lambda)), ...
					initial_theta , options);
		all_theta(c,:) = theta';

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

end

author: zhouxc

repo: Stanford-Machine-Learning-Course


text:
function zx = cvx_expand_dim( z, dim, nx )

%CVX_EXPAND_DIM   Expands an N-D array along a specified dimension.
%   CVX_EXPAND_DIM( X, DIM, NX ) stacks NX copies of the matrix X along the
%   dimension NX. It is equivalent to CAT( DIM, X, X, ..., X ), where X is
%   repeated NX times.
%
%   This is an internal CVX function, and as such no checking is performed to
%   insure that the arguments are valid.

zdims = cell( 1, max( ndims(z), dim ) );
[ zdims{:} ] = deal( ':' );
zdims{dim} = ones( 1, nx );
zx = z( zdims{:} );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function element_node = grid_q16_element ( nelemx, nelemy )

%*****************************************************************************80
%
%% GRID_Q16_ELEMENT produces a grid of 16 node quadrilaterals.
%
%  Example:
%
%    Input:
%
%      NELEMX = 2, NELEMY = 2
%
%    Output:
%
%      ELEMENT_NODE =
%         1,  2,  3,  4,  8,  9, 10, 11, 15, 16, 17, 18, 22, 23, 24, 25;
%         4,  5,  6,  7, 11, 12, 13, 14, 18, 19, 20, 21, 25, 26, 27, 28;
%        22, 23, 24, 25, 29, 30, 31, 32, 36, 37, 38, 39, 43, 44, 45, 46;
%        25, 26, 27, 28, 32, 33, 34, 35, 39, 40, 41, 42, 46, 47, 48, 49.
%
%  Grid:
%
%   43-44-45-46-47-48-49
%    |        |        |
%    |        |        |
%   36 37 38 39 40 41 42
%    |        |        |
%    |        |        |
%   29 30 31 32 33 34 35
%    |        |        |
%    | 3      | 4      |
%   22-23-24-25-26-27-28
%    |        |        |
%    |        |        |
%   15 16 17 18 19 20 21
%    |        |        |
%    |        |        |
%    8  9 10 11 12 13 14
%    |        |        |
%    | 1      | 2      |
%    1--2--3--4--5--6--7
%
%  Reference Element Q16:
%
%    |
%    1 13--14--15--16
%    |  |   :   :   |
%    |  |   :   :   |
%    |  9..10..11..12
%    S  |   :   :   |
%    |  |   :   :   |
%    |  5...6...7...8
%    |  |   :   :   |
%    |  |   :   :   |
%    0  1---2---3---4
%    |
%    +--0-----R-----1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NELEMX, NELEMY, the number of elements along the
%    X and Y directions.  The number of elements generated will be
%    NELEMX * NELEMY.
%
%    Output, integer ELEMENT_NODE(16,NELEMX*NELEMY), the nodes that form
%    each element.
%
  element = 0;

  for j = 1 : nelemy
    for i = 1 : nelemx

      base = ( j - 1 ) * 3 * ( 3 * nelemx + 1 ) + 3 * i - 2;

      element = element + 1;

      element_node( 1,element) = base;
      element_node( 2,element) = base                          + 1;
      element_node( 3,element) = base                          + 2;
      element_node( 4,element) = base                          + 3;
      element_node( 5,element) = base +     ( 3 * nelemx + 1 );
      element_node( 6,element) = base +     ( 3 * nelemx + 1 ) + 1;
      element_node( 7,element) = base +     ( 3 * nelemx + 1 ) + 2;
      element_node( 8,element) = base +     ( 3 * nelemx + 1 ) + 3;
      element_node( 9,element) = base + 2 * ( 3 * nelemx + 1 );
      element_node(10,element) = base + 2 * ( 3 * nelemx + 1 ) + 1;
      element_node(11,element) = base + 2 * ( 3 * nelemx + 1 ) + 2;
      element_node(12,element) = base + 2 * ( 3 * nelemx + 1 ) + 3;
      element_node(13,element) = base + 3 * ( 3 * nelemx + 1 );
      element_node(14,element) = base + 3 * ( 3 * nelemx + 1 ) + 1;
      element_node(15,element) = base + 3 * ( 3 * nelemx + 1 ) + 2;
      element_node(16,element) = base + 3 * ( 3 * nelemx + 1 ) + 3;

    end
  end

  return
end


text:
function y=sinhSlope3(x)
%%SINHSLOPE3 Evaluate the function 6*(sinh(x)-x)/x^3 avoiding
%            problems at the point x=0, where the limit equals 1. This
%            function arises in range enclosure algorithms for interval
%            arithmetic.
%
%INPUTS: x A vector of matrix of real values at which the function should
%          be evaluated.
%
%OUTPUTS: y The value(s) of the function 6*(sinh(x)-x)/x^3 at the
%           point(s) in x.
%
%For values of x having magnitude over 1e-2 or for complex values, the
%function is evaluated as written. For smaller values, a Taylor series
%expansion is used to avoid the singularity near the origin.
%
%This is two times one of the slope function in Table 10.6 of Section
%10.6.2 of [1]. As the standard says the value at 0 is 1, but the value of
%the function in the standard at 0 is 1/2, it is assumed that there is a
%typo in the standard. Hence this function is twice that listed in [1].
%This file implements the function for non-intervals. See the Interval
%class for an implementation for Interval arithmetic.
%
%REFERENCES:
%[1] IEEE Standard for Interval Arithmetic," in IEEE Std 1788-2015,
%    pp.1-97, June 30 2015.
%
%November 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%Allocate space
y=zeros(size(x));

%For large enough values, use the standard function.
sel=abs(x)>1e-2|~isreal(x);
x1=x(sel);
y(sel)=6*(sinh(x1)-x1)./x1.^3;

%For smaller magnitude values, use an sixth-order Taylor series expansion.
%The relative error at 10^-2 is on the order of 1e-23, which is less than
%the precision of a double floating point number at that value. The series
%is given in Horner form.
x1=x(~sel);
y(~sel)=1+x1.^2.*(1/20+x1.^2.*(1/840+x1.^2/60480));

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% some statistical tests

path(path, '../../images/');

if not(exist('name'))
name = 'line_vertical_blurred';
name = 'sawtooth';
name = 'line_vertical';
name = 'flow2';
name = 'flow2';
name = 'noise';
name = 'turbulence';
name = 'corral';
name = 'lena';
name = 'barb';
name = 'reptilskin';
name = 'hair';
name = 'disk';
name = 'square';
end
if not(exist('wavtype'))
wavtype = 'redun';
wavtype = 'ortho';
end

rep = ['results/stats/'];

save_images = 1;

if exist(rep)~=7
    mkdir(rep);
end

n = 512;
n0 = [];
if strcmp(name, 'disk') || strcmp(name, 'square')
    n0 = n;
end

M = load_image(name, n0); M = sum(M,3);
M = rescale(crop(M,n));
n = size(M,1);

options.wavelet_type = 'biorthogonal';
options.wavelet_vm = 4;
Jmin = 4;
Jmax = log2(n)-1;

switch wavtype
    case 'ortho'
        MW = perform_wavelet_transform(M, Jmin, +1, options);
        %% scale selection
        j = Jmax;
        [selx,sely] = compute_quadrant_selection(j,2);
        MW1 = MW(selx,sely);
        % father
        [selx,sely] = compute_quadrant_selection(j-1,2);
        MWf = MW(selx,sely);
        sel = floor( 1:0.5:nj/2+0.5 );
        MWf = MWf(sel,sel);
        % brother
        [selx,sely] = compute_quadrant_selection(j,1);
        MWb = MW(selx,sely);
    case 'redun'
        MW = perform_atrou_transform(M,Jmin,options);
        MW1 = MW{1};
        MWf = MW{4};
        MWb = MW{2};
end
nj = size(MW1,1);

if save_images
    warning off;
%    imwrite(rescale(M),[rep name '-original.png'], 'png');
    warning off;
end

%% compute 1D histogram
nb_bins = 501;
[N,X] = histo(MW1(:),nb_bins);
N = N/prod( size(MW1) );
clf;
plot( X, log2(N) );
t = max(abs(MW1(:)))*0.8;
axis( [-t t -16 max(log2(N))] );

if save_images
%    saveas(gcf, [rep name '-histo.eps'], 'epsc');
%    saveas(gcf, [rep name '-histo.png'], 'png');
end

% parameters for display
pdisp = 800;
ms = 15;
nbins = 50;

sibling_list = { 'same','same','same','same','same','father','brother' };
delta_list = { [1 0] [0 1] [2 0] [0 2] [4 0] [0 0] [0 0] };
ntests = length(sibling_list);
nrows = 3;
ncols = ceil( (ntests+2)/nrows );

clf;
subplot(nrows, ncols, 1);
imageplot( M, 'Original' );
subplot(nrows, ncols, 2);
imageplot( MW1, 'Wavelets coefs' );

for itest = 1:ntests
    sibling = sibling_list{itest};
    d = delta_list{itest};

    switch sibling
        case 'same'
            MW2 = MW1;
        case 'father'
            MW2 = MWf;
        case 'brother'
            MW2 = MWb;
    end

    %% select neighboring coefficients
    [Y,X] = meshgrid(1:nj,1:nj);
    I = find( X+d(1)>0 & X+d(1)<=nj & Y+d(2)>0 & Y+d(2)<=nj );
    v1 = MW1(I);
    J = X(I)+d(1) + (Y(I)+d(2)-1)*nj;
    v2 = MW2(J);

    %* compute covariance matrix
    w1 = v1-mean(v1);
    w2 = v2-mean(v2);
    kappa1 = 2*mean(w1.*w2)/(mean(w1.^2)+mean(w2.^2));
    w1 = w1/sqrt(sum(w1.^2));
    w2 = w2/sqrt(sum(w2.^2));
    kappa = sum(w1.*w2);

    % bounds
    t1 = max(abs(v1))*0.7;
    t2 = max(abs(v2))*0.7;
    t1 = std(v1)*3;
    t2 = std(v2)*3;

    %% display normalized pairwise
    x1 = linspace(-t1,t1,nbins)';
    x2 = linspace(-t2,t2,nbins)';
    h = compute_pairwise_histogram( cat(2,v1,v2), x1, x2, options);
    H = h./repmat( max(h,[], 1), [nbins 1] );
    H(isnan(H)) = min(H(:));

    str = sibling;
    if d(1)~=0 || d(2)~=0
        str = [str ',\delta=(' num2str(d(1)) ',' num2str(d(2)) ')'];
    end

    subplot(nrows, ncols, itest+2);
    hold on;
    imageplot(-H);
    sel = randperm(length(v1)); sel = sel(1:pdisp);
    % aa = plot( v1(sel),v2(sel), '.' );  axis([-t1 t1 -t2 t2]);
    % set(aa, 'MarkerSize', ms);
    title([str ',c=' num2str(kappa*100, 2) '%']);
    hold off;

end
saveas(gcf, [rep name '-' wavtype '-histo-pairw.png'], 'png');

author: gpeyre

repo: matlab-toolboxes


text:
function ihist = imghist(img)
% Author: Javier Montoya (jmontoyaz@gmail.com).
%         http://www.lis.ic.unicamp.br/~jmontoya
%
% IMGHIST calculates the histogram of a given image.
% Input parameters:
%    img: image I (passed as a bidimensional matrix).
% Ouput parameters:
%    ihist: histogram.
%
% Usage:
%    I     = imread('tire.tif');
%    ihist = imghist(I);
%    figure; stem(ihist); title('Image Histogram');

   if exist('img', 'var') == 0
      error('Error: Specify an input image.');
   end

   ihist       = [];
   [rows,cols] = size(img);
   maxgval     = 255;
   ihist       = zeros(1,maxgval);

   for i=0:maxgval
      ihist(i+1) = sum(img(:)==i);
   end
end

author: Sable

repo: mcbench-benchmarks


text:
function pyObj = bst_mat2py(matArray, isForced)
% BST_MAT2PY: Converts Matlab matrices to Python arrays
%
% INPUTS:
%    - matArray: Matlab array to convert to a Python array
%    - isForced: Boolean, if 1: force a manual conversion of the array.
%                Useful for preventing Matlab from transposing the row vectors into column vectors.
%
% @=============================================================================
% 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: Martin Cousineau, 2019
%          Francois Tadel, 2020

% Parse inputs
if (nargin < 2) || isempty(isForced)
    isForced = 0;
end

% Get Matlab version
MatlabVersion = bst_get('MatlabVersion');
% Empty objects
if isempty(matArray)
    pyObj = py.None;
% Matlab >= 2020b: Can convert everything
elseif (MatlabVersion >= 909) && ~isForced
    pyObj = matArray;
% Matlab <= 2020a: Must do some manual conversions
else
    switch class(matArray)
        case {'double', 'single', 'logical', 'int8', 'int16', 'int32', 'int64', 'uint8', 'uint16', 'uint32', 'uint64'}
            % Convert numerical matrices to Numpy arrays
            shape = int64(size(matArray));
            pyObj = py.numpy.array(reshape(matArray', 1, []));
            pyObj = pyObj.reshape(shape);

        case 'cell'
            % Convert cell arrays to Python lists (which support mixed types)
            pyObj = py.list(matArray(:)');

        otherwise
            error(['Unsupported class: ' class(matArray)]);
    end
end

author: brainstorm-tools


text:
function [sigma_draw, Ft_draw, F_draw, Fcomp_draw] = VARdrawpost(VAR)
% =======================================================================
% Draw from the posterior distribution of a VAR model
% =======================================================================
% [sigma_draw, Ft_draw] = VARdrawpost(VAR)
% -----------------------------------------------------------------------
% INPUT
%   - VAR: structure, result of VARmodel function
% -----------------------------------------------------------------------
% OUPUT
%   - sigma_draw: draw from the posterior of VCV matrix of VAR residuals
%   - Ft_draw: draw from the posterior of Ft of VAR
% -----------------------------------------------------------------------
% EXAMPLE
%   - See VARToolbox_Code.m in "../Primer/"
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2012. Updated November 2020
% -----------------------------------------------------------------------

%% Get relevant parameters from VAR structure
%===============================================
nobs = VAR.nobs;
k    = [];
X    = VAR.X;
const= VAR.const;
nvar = VAR.nvar;
nlag = VAR.nlag;

%% OLS estimates
%============================================
Ft_hat = VAR.Ft;
sigma_hat = VAR.sigma;
inv_sigma_hat = inv(sigma_hat);

%% Draw the VCV matrix (sigma)
%============================================
inv_sigma_draw = wishrnd(inv_sigma_hat/nobs,nobs);
sigma_draw     = inv(inv_sigma_draw);

%% Draw the coeffiecient matrix (Ft)
%============================================
% Compute VCV of coefficient matrix
aux1 = (X'*X)\eye(size(X,2));
aux2 = kron(sigma_draw, aux1);
aux2 = (aux2 + aux2')/2; % Force symmetry (might fail bc of numerical error)
Fthat_vec = Ft_hat(:);
Ftdraw = mvnrnd(Fthat_vec,aux2);
Ft_draw = reshape(Ftdraw,k,nvar);
F_draw = Ft_draw';
Fcomp_draw = [F_draw(:,1+const:nvar*nlag+const); eye(nvar*(nlag-1)) zeros(nvar*(nlag-1),nvar)];

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function S = querylpv(lpv, p)
%QUERYLPV Query an LPV model at given parameter point
%	S = QUERYLPV(lpv, p)
%
%	lpv   - LPV model
%	p     - given parameter
%
%	S     - system matrix at the given point

% TODO: convert to lti sys

S = zeros(size(lpv));
for i = 1:size(lpv,1)
	for j = 1:size(lpv,2)
		S(i,j) = lpv{i,j}(p);
	end
end


text:
function showStereoMatches(I1, I2, matches, current_time_flag, style, lineWidth,  ptSyle, ptWidth)
%SHOWSTEREOMATCHES Plot matches detected in two images using MATLAB
%
% INPUT:
%   - I2_l: left image at time t
%   - I2_r: right image at time t
%   - matches(1, N): array of structure conatining keypoints
%   - current_time_flag: flag for matches in previous time (=0) or current time (=1)
%   - style: defines color and type of plotted line
%   - lineWidth: defines thickness of plotted line
%   - ptSyle: defines type of points; eg: 'o'
%   - ptWidth: defines thickness of plotted points

% false color
imgOverlay = imfuse(I2, I1);
% create red-cyan image instead of the imfuse default
imgOverlay(:,:,1) = imgOverlay(:,:,2);
imgOverlay(:,:,2) = imgOverlay(:,:,3);

% show image
imshow(imgOverlay);
hold on

% for all matches
x_from = zeros(1, size(matches, 2));
x_to = zeros(1, size(matches, 2));
y_from = zeros(1, size(matches, 2));
y_to = zeros(1, size(matches, 2));
num = 0;
for match = matches
    num = num + 1;
    if(current_time_flag)
        x_from(num) = match.pt2_r.location(1);
        x_to(num)   = match.pt2_l.location(1);
        y_from(num) = match.pt2_r.location(2);
        y_to(num)   = match.pt2_l.location(2);
    else
        x_from(num) = match.pt1_r.location(1);
        x_to(num)   = match.pt1_l.location(1);
        y_from(num) = match.pt1_r.location(2);
        y_to(num)   = match.pt1_l.location(2);
    end
end

% plot line
plot([y_from; y_to], [x_from; x_to], style, 'Linewidth', lineWidth);
% plot end points
scatter(y_from, x_from, ptSyle, 'r', 'Linewidth', ptWidth);
scatter(y_to, x_to, ptSyle, 'g', 'Linewidth', ptWidth);
end

repo: Stereo-Odometry-SOFT


text:
%data   :3*length increment type equally spaced sensor output.
%alg    :Type of the coning algorithm to be used
function [inc corr]=coning_minor(data, alg)
inlen=size(data,2);
switch (alg)
    case(0)
        %ignagni(1990):Algorithm A (no correction - constant approximation)
        inc=data;
        corr=zeros(size(data));
    case(1)
        %ignagni(1990):Algorithm D (quadratic approximation to 3 points)
        outlen=floor((inlen-3)/3)+1;
        inc=zeros(3,outlen);
        corr=zeros(3,outlen);

        ind=1;
        for i=3:3:inlen
            inc(:,ind)=sum(data(:,i-2:i),2);
            corr(:,ind)=(33/80)*cross(data(:,i-2), data(:,i))+(57/80)*cross(data(:,i-1),data(:,i)-data(:,i-2));
            ind=ind+1;
        end
    case(2)
        %ignagni(1990):Algorithm E
        outlen=floor((inlen-3)/3)+1;
        inc=zeros(3,outlen);
        corr=zeros(3,outlen);

        ind=1;
        for i=3:3:inlen
            inc(:,ind)=sum(data(:,i-2:i),2);
            corr(:,ind)=(9/20)*cross(data(:,i-2), data(:,i))+(27/40)*cross(data(:,i-1),data(:,i)-data(:,i-2));
            ind=ind+1;
        end
    case(3)
        %ignagni(1996):Algorithm 1
        outlen=floor((inlen-4)/2)+1;
        inc=zeros(3,outlen);
        corr=zeros(3,outlen);

        ind=1;
        for i=4:2:inlen
            vr_a=sum(data(:,i-1:i),2);
            vr_b=sum(data(:,i-3:i-2),2);
            inc(:,ind)=var_a;
            corr(:,ind)=(32/45)*cross(data(:,i-1), data(:,i))+(-1/180)*cross(vr_b,vr_a);
            ind=ind+1;
        end
    case(4)
        %ignagni(1996):Algorithm 2
        outlen=floor((inlen-3)/2)+1;
        inc=zeros(3,outlen);
        corr=zeros(3,outlen);

        ind=1;
        for i=3:2:inlen
            inc(:,ind)=sum(data(:,i-1:i),2);
            corr(:,ind)=cross(((-1/30)*data(:,i-2)+(11/15)*data(:,i-1)), data(:,i));
            ind=ind+1;
        end
    case (5)
        %ignagni(1996):Algorithm 3 - Ignagni (1990):Algorithm F
        outlen=floor((inlen-3)/3)+1;
        inc=zeros(3,outlen);
        corr=zeros(3,outlen);

        ind=1;
        for i=3:3:inlen
            inc(:,ind)=sum(data(:,i-2:i),2);
            corr(:,ind)=cross(((9/20)*data(:,i-2)+(27/20)*data(:,i-1)), data(:,i));
            ind=ind+1;
        end
    case (6)
        %ignagni(1996):Algorithm 4
        outlen=floor((inlen-6)/3)+1;
        inc=zeros(3,outlen);
        corr=zeros(3,outlen);

        ind=1;
        for i=6:3:inlen
            vr_a=sum(data(:,i-2:i),2);
            vr_b=sum(data(:,i-5:i-3),2);
            inc(:,ind)=vr_a;
            corr(:,ind)=cross(((243/560)*data(:,i-2)+(1539/1120)*data(:,i-1)), data(:,i))+(1/3360)*cross(vr_b,vr_a);
            ind=ind+1;
        end
    case (7)
        %ignagni(1996):Algorithm 5
        outlen=floor((inlen-4)/3)+1;
        inc=zeros(3,outlen);
        corr=zeros(3,outlen);

        ind=1;
        for i=4:3:inlen
            inc(:,ind)=sum(data(:,i-2:i),2);
            corr(:,ind)=cross(((3/280)*data(:,i-3)+(57/140)*data(:,i-2)+(393/280)*data(:,i-1)), data(:,i));
            ind=ind+1;
        end
    case (8)
        %ignagni(1996):Algorithm 6
        outlen=floor((inlen-5)/3)+1;
        inc=zeros(3,outlen);
        corr=zeros(3,outlen);

        ind=1;
        for i=5:3:inlen
            inc(:,ind)=sum(data(:,i-2:i),2);
            corr(:,ind)=cross(((-1/420)*data(:,i-4)+(1/40)*data(:,i-3)+(157/420)*data(:,i-2)+(1207/840)*data(:,i-1)), data(:,i));
            ind=ind+1;
        end
    case (9)
        %ignagni(1996):Algorithm 7
        outlen=floor((inlen-4)/4)+1;
        inc=zeros(3,outlen);
        corr=zeros(3,outlen);

        ind=1;
        for i=4:4:inlen
            inc(:,ind)=sum(data(:,i-3:i),2);
            corr(:,ind)=cross(((54/105)*data(:,i-3)+(92/105)*data(:,i-2)+(214/105)*data(:,i-1)), data(:,i));
            ind=ind+1;
        end
    otherwise
        disp('Undefined Algorithm (Perhaps, someone else adds it later)');
end

repo: ekfmonoslam


text:
%% housekeeping
close all
home()
%% Choose a model type: see cell "create the structural VAR model" below
model_type=5;

%% Create dataset
do_plot=true;
scale=100;

[db,varlist]=create_dataset(scale,do_plot);

%% Create the structural VAR model
close()

% first we create a template structure
% ------------------------------------
tpl=svar.template();

% we update the fields of the structure
% --------------------------------------
tpl.endogenous=varlist;
tpl.nlags=2;

% create restrictions on parameters as well as markov chains
%------------------------------------------------------------
switch model_type
    case 0
        % constant-parameter model
        [restrictions,tpl]=create_restrictions_and_markov_chains0(tpl);
    case 1
        % Coefficients are switching regimes across all equations
        % (synchronized case)
        [restrictions,tpl]=create_restrictions_and_markov_chains1(tpl);
    case 2
        % Coefficients and variances have different chains, different
        % regimes, and different durations
        [restrictions,tpl]=create_restrictions_and_markov_chains2(tpl);
    case 3
        % Only coefficients in monetary policy equation are changing
        [restrictions,tpl]=create_restrictions_and_markov_chains3(tpl);
    case 4
        % Only variance in monetary policy equation is changing
        [restrictions,tpl]=create_restrictions_and_markov_chains4(tpl);
    case 5
        % Both coefficients and variances in monetary policy equation
        % change with two independent Markov processes
        [restrictions,tpl]=create_restrictions_and_markov_chains5(tpl);
    case 6
        % Only variances in ALL three equations switch
        [restrictions,tpl]=create_restrictions_and_markov_chains6(tpl);
    otherwise
        error('the coded model types are 0, 1, 2, 3, 4 and 6')
end

% finally we create a svar object by pushing the structure into svar
%--------------------------------------------------------------------
m=svar(tpl,'data',db,'estim_linear_restrictions',restrictions);

%% Find posterior mode

mest=estimate(m,'estim_start_date','1960Q1');

%% Printing out A0hat, A1hat, and A2hat in a form compatible with Zha''s original Matlab code

output_in_original_zha_matlab_code(mest,scale)

% uncomment this if you want to capture output
% [A0hat,A1hat,A2hat]=output_in_original_zha_matlab_code(mest,scale)

%% Markov chain Monte Carlo
% Note that because of the linear restrictions, not all parameters are
% estimated. Hence, the effective number of estimated parameters is smaller
% than the number of parameters declared by the user. This is reflected in
% the dimensions of lb,ub,x0 and SIG below. The user does not have to be
% concerned about those.
[objective,lb,ub,x0,SIG]=pull_objective(mest);

ndraws_mcmc         = 1500;  % number of parameter draws through MCMC.
ndraws_burnin       = floor(0.1*ndraws_mcmc); % number of parameter draws to be burned
mcmc_options=struct('burnin',ndraws_burnin,'N',ndraws_mcmc,'thin',1);

Results=mh_sampler(objective,lb,ub,mcmc_options,x0,SIG);

%% Marginal data density
% pick yours: 'bridge','mhm','mueller','swz','is','ris','cj'
tic
log_mdd = mcmc_mdd(Results.pop,lb,ub,...
    struct('log_post_kern',objective,... % function to MINIMIZE !!!
    'algorithm','swz',... % MDD algorithm
    'L',2000 ... % Number of i.i.d. draws of the proposal density function
));
minutes_MDD_Took = toc/60;

laplace_approx=mest.estimation.posterior_maximization.log_marginal_data_density_laplace;
disp(['log(MDD): Laplace Approximation: ',num2str(laplace_approx)])
disp(['log(MDD): Sims-Waggoner-Zha Approximation: ',num2str(log_mdd)])
disp([' Minutes Sims-Waggoner-Zha Approximation took: ',num2str(minutes_MDD_Took)])

%% Impulse responses
myirfs=irf(mest);

do_plot_irfs(myirfs,mest);

%% do smoothed probabilities

do_plot_smoothed_probabilities(mest)

%% Out-of sample forecasts at the mode

mycast=forecast(mest);

do_plot_unconditional_forecasts(mycast,mest)

%% Conditional forecast on ygap
% conditional information
%-------------------------
ygap=scale*(-0.025:0.005:8*0.005).';
cond_db=struct('ygap',ts('2015Q3',ygap));

% options for the exercise
%--------------------------
myoptions=struct('cbands',[10,20,50,80,90],'do_plot',true,'nsteps',20,...
    'param_uncertainty',false,'shock_uncertainty',true,'ndraws',200);

% do it
%-------
tic
[fkst,bands,hdl]=do_conditional_forecasts(mest,db,cond_db,Results.pop,myoptions);
fprintf('\n\n Computing conditional forecasts took %0.4f minutes\n\n',toc/60);
%% Median conditional forecasts
% we take advantage of the fact that in the bands above we specified 50 in
% the bands above
%--------------------------------------------------------------------------
figure('name','Median Conditional Forecasts');
nvars=mest.endogenous.number;
for ivar=1:nvars
    thisname=mest.endogenous.name{ivar};
    subplot(nvars,1,ivar)
    plot(bands.(thisname)('ci_50'),'linewidth',2)
    title(mest.endogenous.tex_name{ivar})
end

author: jmaih


text:
function y = test_h(x,index)
%TEST_H       Some test functions collected from Coconut
%
%   y = test_h(x,index)
%
%Input index specifies test function. Copyright see below.
%

% written  04/04/04     S.M. Rump
%

switch index

  case 1     % source see bottom of file
    y = x(3)-1 + x(1).^2 + x(2).^2 + (x(3)+x(4)).^2 + sin(x(3)).^2 + x(1).^2*x(2).^2 + ...
        x(4)-3 + sin(x(3)).^2 + (x(4)-1).^2 + x(2).^4 + x(3).^4 + (x(4)+x(1)).^2 + ...
        (x(1)-4 + sin(x(4)).^2 + x(2).^2*x(3).^2).^2 + sin(x(4)).^4;

  case 2     % source see bottom of file
    N = length(x);      % model problem: N = 1000, initial x=ones(N,1);
    I = 1:N-4;
    y = sum( (-4*x(I)+3.0).^2 ) + sum( ( x(I).^2 + 2*x(I+1).^2 + ...
              3*x(I+2).^2 + 4*x(I+3).^2 + 5*x(N).^2 ).^2 );

end

% function 1 taken from   http://www.sor.princeton.edu/~rvdb/ampl/nlmodels/cute/allinitu.mod
% # AMPL Model by Hande Y. Benson
% #
% # Copyright (C) 2001 Princeton University
% # All Rights Reserved
% #
% # Permission to use, copy, modify, and distribute this software and
% # its documentation for any purpose and without fee is hereby
% # granted, provided that the above copyright notice appear in all
% # copies and that the copyright notice and this
% # permission notice appear in all supporting documentation.
%
% #   Source:
% #   N. Gould, private communication.
%
% #   SIF input: Nick Gould, June 1990.
%
% #   classification OUR2-AY-4-0
%
% var x{1..4};
%
% minimize f:
% x[3]-1 +
% x[1]^2+
% x[2]^2 + (x[3]+x[4])^2 +
% sin(x[3])^2 + x[1]^2*x[2]^2 + x[4]-3 +
% sin(x[3])^2 +
% (x[4]-1)^2 +
% (x[2]^2)^2+
% (x[3]^2 + (x[4]+x[1])^2)^2 +
% (x[1]-4 + sin(x[4])^2 + x[2]^2*x[3]^2)^2 +
% sin(x[4])^4;
%
% solve;
% display f;
% display x;

% function 2 taken from   http://www.sor.princeton.edu/~rvdb/ampl/nlmodels/cute/bdqrtic.mod
% # AMPL Model by Hande Y. Benson
% #
% # Copyright (C) 2001 Princeton University
% # All Rights Reserved
% #
% # Permission to use, copy, modify, and distribute this software and
% # its documentation for any purpose and without fee is hereby
% # granted, provided that the above copyright notice appear in all
% # copies and that the copyright notice and this
% # permission notice appear in all supporting documentation.
%
% #   Source: Problem 61 in
% #   A.R. Conn, N.I.M. Gould, M. Lescrenier and Ph.L. Toint,
% #   "Performance of a multifrontal scheme for partially separable
% #   optimization",
% #   Report 88/4, Dept of Mathematics, FUNDP (Namur, B), 1988.
%
% #   SIF input: Ph. Toint, Dec 1989.
%
% #   classification SUR2-AN-V-0
%
% param N:=1000;
% var x{1..N} := 1.0;
%
% minimize f:
% sum {i in 1..N-4} (-4*x[i]+3.0)^2 + sum {i in 1..N-4} (x[i]^2+2*x[i+1]^2+3*x[i+2]^2+4*x[i+3]^2+5*x[N]^2)^2;
%
% solve;
% display f;
% display x;

repo: OpenMAS