text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% SPHARM-MAT is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function userdataHandle = initOptions(userdataHandle)

% Start initialization
% Not implemented: 'ParamQuad', 'ParamConf'
userdataHandle.keyword.list = {'ParamCALD', 'ParamPDM', 'ParamQuad', 'ParamConf', 'ExpLSF', 'ExpIRF', ...
    'AligSHREC', 'AligFOE', 'AligLandmark', 'ExpAlig', 't_map', 'PCA', 't_map+PCA', ...
    'PDM_stat', 'Import', 'DisplayRes', 'bim2gipl', 'fix2gipl', 'gipl2bim', 'smo2surf_para_meta', ...
    'surf_para_meta2smo', 'des2meta_coef', 'meta_coef2des', 'reg2procalign_meta', 'InHouse_Fix', ...
    'PDM_Fix', 'Config'};
userdataHandle.keyword.selected = '';

userdataHandle.mesh={'orig';'quad32';'quad64';'quad128';'quad256';'quad512';'icosa1';'icosa2'; ...
    'icosa3';'icosa4';'icosa5';'icosa6'};
userdataHandle.space={'object';'param';'both'};
userdataHandle.shade={'solid';'mesh';'both'};
userdataHandle.overlay={'none';'adc_paramap'};
userdataHandle.colormap={'jet';'hot';'summer';'cool';'autumn';'winter';'spring'};
userdataHandle.export={'screen';'png';'both'};

% CALD Parameterization
userdataHandle.ParamCALD.vars = {'MeshGridSize','MaxSPHARMDegree','Tolerance','Smoothing','Iteration', ...
        'LocalIteration','t_major','SelectDiagonal','OutDirectory'};
userdataHandle.ParamCALD.args = [1 1 1 1 1 1 10 10 200];  % Encode the number and type of each variables
userdataHandle.ParamCALD.inFilter = {...
    '*_bim.mat;*_fix.mat;*_obj.mat', 'Binary and Mesh Objects (*_bim.mat;*_fix.mat;*_obj.mat)'; ...
    '*_obj.mat','Mesh Objects (*_obj.mat)'; ...
    '*_bim.mat;*_fix.mat', 'Bianry Objects (*_bim.mat, *_fix.mat)'; ...
    '*.mat', 'All Objects (*.mat)' ...
    };
userdataHandle.ParamCALD.default={'50', '6', '2', '2', '100', '10', '', '','./'};

for i=1:length(userdataHandle.ParamCALD.vars)
    varName = sprintf('userdataHandle.ParamCALD.%s', userdataHandle.ParamCALD.vars{i});
    if userdataHandle.ParamCALD.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.ParamCALD.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.ParamCALD.default{%d});', varName, i);
        eval(strCC);
    end
end

% PDM Parameterization (GenParamMesh)
userdataHandle.ParamPDM.vars = {'iter', 'label', 'OutDirectory', 'others'};
userdataHandle.ParamPDM.args = [1 1 200 20];
userdataHandle.ParamPDM.inFilter = {'*.gipl;*.hdr', 'Bianry Objects (*.gipl, *hdr)'};
userdataHandle.ParamPDM.default = {'500','1','./',''};
userdataHandle.ParamPDM.path = '';
userdataHandle.ParamPDM.command = 'GenParaMesh';
userdataHandle.ParamPDM.exist = 0;

for i=1:length(userdataHandle.ParamPDM.vars)
    varName = sprintf('userdataHandle.ParamPDM.%s', userdataHandle.ParamPDM.vars{i});
    if userdataHandle.ParamPDM.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.ParamPDM.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.ParamPDM.default{%d});', varName, i);
        eval(strCC);
    end
end

% Parameterization using a similar method to Brechbuhler's methd
userdataHandle.ParamQuad.vars = {'NRS_step','OutDirectory'};
userdataHandle.ParamQuad.args = [1 200];  % Encode the number and type of each variables
userdataHandle.ParamQuad.inFilter = {'*_bim.mat;*_fix.mat', 'Bianry Objects (*_bim.mat, *_fix.mat)'; ...
   '*_obj.mat','Mesh Objects (*_obj.mat)'; '*.mat', 'All Objects (*.mat)'};
userdataHandle.ParamQuad.default={'2','./'};

for i=1:length(userdataHandle.ParamQuad.vars)
    varName = sprintf('userdataHandle.ParamQuad.%s', userdataHandle.ParamQuad.vars{i});
    if userdataHandle.ParamQuad.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.ParamQuad.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.ParamQuad.default{%d});', varName, i);
        eval(strCC);
    end
end

% Expansion in the sense of least mean square
userdataHandle.ExpLSF.vars = {'MaxSPHARMDegree','OutDirectory'};
userdataHandle.ExpLSF.args = [1 200];
userdataHandle.ExpLSF.inFilter = {'*_ini.mat;*_smo.mat', 'Parameterized Objects (*_ini.mat, *_smo.mat)'; '*.mat', 'All Objects (*.mat)'};
userdataHandle.ExpLSF.default = {'15', './'};

for i=1:length(userdataHandle.ExpLSF.vars)
    varName = sprintf('userdataHandle.ExpLSF.%s', userdataHandle.ExpLSF.vars{i});
    if userdataHandle.ExpLSF.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.ExpLSF.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.ExpLSF.default{%d});', varName, i);
        eval(strCC);
    end
end

% SHREC Alignment
userdataHandle.AligSHREC.vars = {'Template', 'MaxSPHARMDegree', 'GroupAlpha', 'NormalizeSize', ...
    'BaseRes', 'HierarchyStep', 'HierarchyDepth', 'Top_K', 'GammaRes', 'OutDirectory'};
userdataHandle.AligSHREC.args = [100 1 1 10 1 1 1 1 1 200];
userdataHandle.AligSHREC.inFilter = {'*_des.mat;*_prm.mat', 'SPHARM Descriptors (*_des.mat,*_prm.mat)'; '*.mat', 'All Objects (*.mat)'};
userdataHandle.AligSHREC.default = {'','12','100','','1','1','3','1','2','./'};

for i=1:length(userdataHandle.AligSHREC.vars)
    varName = sprintf('userdataHandle.AligSHREC.%s', userdataHandle.AligSHREC.vars{i});
    if userdataHandle.AligSHREC.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.AligSHREC.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.AligSHREC.default{%d});', varName, i);
        eval(strCC);
    end
end

% First order ellipsoid alignment
userdataHandle.AligFOE.vars = {'CPoint', 'NPole', 'MaxSPHARMDegree', 'OutDirectory'};
userdataHandle.AligFOE.args = [10 10 1 200];
userdataHandle.AligFOE.inFilter = {'*_des.mat', 'SPHARM Descriptors (*_des.mat)'; '*.mat', 'All Objects (*.mat)'};
userdataHandle.AligFOE.default = {'', '', '12', './'};

for i=1:length(userdataHandle.AligFOE.vars)
    varName = sprintf('userdataHandle.AligFOE.%s', userdataHandle.AligFOE.vars{i});
    if userdataHandle.AligFOE.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.AligFOE.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.AligFOE.default{%d});', varName, i);
        eval(strCC);
    end
end

% PDM expansion and alignment (ParaToSPHARMMesh)
userdataHandle.ExpAlig.vars = {'flipTemplate', 'subdivLevel', 'spharmDegree', ...
        'regTemplate', 'FinalFlip', 'OutDirectory', 'others'};
userdataHandle.ExpAlig.args = [100 1 1 100 1 200 20];
userdataHandle.ExpAlig.inFilter = {'*_surf.meta', 'Parameterized Objects (*_surf.meta)'};
userdataHandle.ExpAlig.default = {'','20','12','','0','./',''};
userdataHandle.ExpAlig.path = '';
userdataHandle.ExpAlig.command = 'ParaToSPHARMMesh';
userdataHandle.ExpAlig.exist = 0;

for i=1:length(userdataHandle.ExpAlig.vars)
    varName = sprintf('userdataHandle.ExpAlig.%s', userdataHandle.ExpAlig.vars{i});
    if userdataHandle.ExpAlig.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.ExpAlig.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.ExpAlig.default{%d});', varName, i);
        eval(strCC);
    end
end

% t-statistical analysis
userdataHandle.t_map.vars = {'Atlas', 'Smoothing_FWHM', 'EqualVariance', 'Signal', 'SampleMesh', 'OutputNamePrefix','OutDirectory','GroupIDs'};
userdataHandle.t_map.args = [100 1 10 10 10 10 200 10];
userdataHandle.t_map.inFilter = {'*.mat', 'Registered Objects (*.mat)'; '*.mat', 'All Objects (*.mat)'};
userdataHandle.t_map.default = {'./atlas.mat','5','','','','t_map','./','Ctrl,PT'};

for i=1:length(userdataHandle.t_map.vars)
    varName = sprintf('userdataHandle.t_map.%s', userdataHandle.t_map.vars{i});
    if userdataHandle.t_map.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.t_map.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.t_map.default{%d});', varName, i);
        eval(strCC);
    end
end

% PCA analysis
userdataHandle.PCA.vars = {'GroupID','OutputName', 'OutDirectory'};
userdataHandle.PCA.args = [10 10 200];
userdataHandle.PCA.inFilter = {'*.mat', 'Registered Objects (*.mat)'; '*.mat', 'All Objects (*.mat)'};
userdataHandle.PCA.default = {'All','PCA_stat.mat','./'};

for i=1:length(userdataHandle.PCA.vars)
    varName = sprintf('userdataHandle.PCA.%s', userdataHandle.PCA.vars{i});
    if userdataHandle.PCA.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.PCA.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.PCA.default{%d});', varName, i);
        eval(strCC);
    end
end

% Display T-stat maps
userdataHandle.res_t_map.vars = {'Threshold_p_value', 'Overlay','Colormap'};
userdataHandle.res_t_map.args = [1 10 10];
userdataHandle.res_t_map.inFilter = {'t_map_*.mat',  'T-Stat Results (t_map_*.mat)'; ...
    '*.mat', 'All Objects (*.mat)'};
userdataHandle.res_t_map.default = {'0.05', '', ''};

for i=1:length(userdataHandle.res_t_map.vars)
    varName = sprintf('userdataHandle.res_t_map.%s', userdataHandle.res_t_map.vars{i});
    if userdataHandle.res_t_map.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.res_t_map.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.res_t_map.default{%d});', varName, i);
        eval(strCC);
    end
end

% Display PCA results
userdataHandle.res_PCA.vars = {'Level','Sigma','Mesh','MaxSPHARMDegree'};
userdataHandle.res_PCA.args = [1 1 10 1];
userdataHandle.res_PCA.inFilter = {'PCA_*.mat',  'PCA Results (PCA_*.mat)'; ...
    '*.mat', 'All Objects (*.mat)'};
userdataHandle.res_PCA.default = {'4','3','','15'};

for i=1:length(userdataHandle.res_PCA.vars)
    varName = sprintf('userdataHandle.res_PCA.%s', userdataHandle.res_PCA.vars{i});
    if userdataHandle.res_PCA.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.res_PCA.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.res_PCA.default{%d});', varName, i);
        eval(strCC);
    end
end

% PDM topology fix (SegPostProcess)
userdataHandle.PDM_Fix.vars = {'space','OutDirectory', 'others'};
userdataHandle.PDM_Fix.args = [3 200 20];
userdataHandle.PDM_Fix.inFilter = {'*.gipl;*.hdr', 'Bianry Objects (*.gipl, *hdr)'};
userdataHandle.PDM_Fix.default = {'0.75,0.75,0.75','./',''};
userdataHandle.PDM_Fix.path = '';
userdataHandle.PDM_Fix.command = 'SegPostProcess';
userdataHandle.PDM_Fix.exist = 0;

for i=1:length(userdataHandle.PDM_Fix.vars)
    varName = sprintf('userdataHandle.PDM_Fix.%s', userdataHandle.PDM_Fix.vars{i});
    if userdataHandle.PDM_Fix.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.PDM_Fix.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.PDM_Fix.default{%d});', varName, i);
        eval(strCC);
    end
end

% InHouse topology fixing method
userdataHandle.InHouse_Fix.vars = {'Connectivity','Epsilon','OutDirectory'};
userdataHandle.InHouse_Fix.args = [10 1 200];
userdataHandle.InHouse_Fix.inFilter = {'*_bim.mat;*_fix.mat', 'Bianry Objects (*_bim.mat, *_fix.mat)'; ...
    '*.mat', 'All Objects (*.mat)'};
userdataHandle.InHouse_Fix.default = {'','1.5','./'};

for i=1:length(userdataHandle.InHouse_Fix.vars)
    varName = sprintf('userdataHandle.InHouse_Fix.%s', userdataHandle.InHouse_Fix.vars{i});
    if userdataHandle.InHouse_Fix.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.InHouse_Fix.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.InHouse_Fix.default{%d});', varName, i);
        eval(strCC);
    end
end

% System configuration
userdataHandle.Config.vars = {'PDM_path_pc','PDM_path_unix'};
userdataHandle.Config.args = [200,200];
userdataHandle.Config.default = {'./','./'};
userdataHandle.Config.filename = 'sysConfig.mat';

for i=1:length(userdataHandle.Config.vars)
    varName = sprintf('userdataHandle.Config.%s', userdataHandle.Config.vars{i});
    if userdataHandle.Config.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.Config.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.Config.default{%d});', varName, i);
        eval(strCC);
    end
end

userdataHandle.Config.loaded = 0;

if exist(fullfile(pwd, userdataHandle.Config.filename), 'file')
    load(fullfile(pwd, userdataHandle.Config.filename));
    if ispc
        PDM_path = PDM_path_pc;
    elseif isunix
        PDM_path = PDM_path_unix;
    end

    userdataHandle.Config.PDM_path_pc = PDM_path_pc;
    userdataHandle.Config.PDM_path_unix = PDM_path_unix;

    if exist('PDM_path','var')
        userdataHandle.Config.loaded = 1;
        userdataHandle.Config.PDM_path = PDM_path;
        userdataHandle.ParamPDM.path = PDM_path;
        userdataHandle.PDM_Fix.path = PDM_path;
        userdataHandle.ExpAlig.path = PDM_path;
%         userdataHandle.PDM_stat.path = PDM_path;
    end

    com1 = [PDM_path '/' userdataHandle.ParamPDM.command];
    com2 = [PDM_path '/' userdataHandle.ExpAlig.command];
%     com3 = [PDM_path '/' userdataHandle.PDM_stat.command];
    com4 = [PDM_path '/' userdataHandle.PDM_Fix.command];

    if ispc
        com1 = [PDM_path '/' userdataHandle.ParamPDM.command '.exe'];
        com2 = [PDM_path '/' userdataHandle.ExpAlig.command '.exe'];
        com4 = [PDM_path '/' userdataHandle.PDM_Fix.command '.exe'];
    end
    if exist(com1, 'file')
        userdataHandle.ParamPDM.exist =1;
    else
        userdataHandle.ParamPDM.exist =0;
    end
    if exist(com2, 'file')
        userdataHandle.ExpAlig.exist =1;
    else
        userdataHandle.ExpAlig.exist =0;
    end
    if exist(com4, 'file')
        userdataHandle.PDM_Fix.exist =1;
    else
        userdataHandle.PDM_Fix.exist =0;
    end
end

% Display objects
userdataHandle.DisplayObjs.vars = {'Space','Mesh','Shade','Overlay','Export','Degree','Template'};
userdataHandle.DisplayObjs.args = [10 10 10 10 10 1 100];
userdataHandle.DisplayObjs.inFilter = {...
    '*.mat', 'All Objects (*.mat)'; ...
    '*_bim.mat; *_fix.mat', 'Bianry Objects (*_bim.mat, *_fix.mat)'; ...
    '*_obj.mat', 'Mesh Objects (*_obj.mat)'; ...
    '*_ini.mat; *_smo.mat', 'Parameterized Objects (*_ini.mat, *_smo.mat)'; ...
    '*_des.mat', 'SPHARM Descriptors (*_des.mat)'; ...
    '*_reg.mat; *_nor.mat',  'Registered Objects (*_reg.mat, *_nor.mat)' ...
    };
userdataHandle.DisplayObjs.default = {'','','','','','',''};

for i=1:length(userdataHandle.DisplayObjs.vars)
    varName = sprintf('userdataHandle.DisplayObjs.%s', userdataHandle.DisplayObjs.vars{i});
    if userdataHandle.DisplayObjs.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.DisplayObjs.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.DisplayObjs.default{%d});', varName, i);
        eval(strCC);
    end
end

%  Scale objects
userdataHandle.ScaleObjs.vars = {'ScalingFactor','OutDirectory'};
userdataHandle.ScaleObjs.args = [100 200];
userdataHandle.ScaleObjs.inFilter = {'*_des.mat; *_prm.mat; *_reg.mat', 'SPHARM Descriptors (*_des.mat, *_prm.mat, *_reg.mat)'; ...
    '*_obj.mat; *_ini.mat; *_smo.mat',  'Meshed Objects (*_obj.mat, *_ini.mat, *_smo.mat)'; ...
    '*.mat', 'All Objects (*.mat)'};
userdataHandle.ScaleObjs.default = {'','./'};

for i=1:length(userdataHandle.ScaleObjs.vars)
    varName = sprintf('userdataHandle.ScaleObjs.%s', userdataHandle.ScaleObjs.vars{i});
    if userdataHandle.ScaleObjs.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.ScaleObjs.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.ScaleObjs.default{%d});', varName, i);
        eval(strCC);
    end
end

% Create average object
userdataHandle.AverageObjs.vars = {'OutputName','OutDirectory'};
userdataHandle.AverageObjs.args = [10 200];
userdataHandle.AverageObjs.inFilter = {'*_des.mat; *_prm.mat; *_reg.mat', 'SPHARM Descriptors (*_des.mat, *_prm.mat, *_reg.mat)'; ...
    '*.mat', 'All Objects (*.mat)'};
userdataHandle.AverageObjs.default = {'atlas.mat','./'};

for i=1:length(userdataHandle.AverageObjs.vars)
    varName = sprintf('userdataHandle.AverageObjs.%s', userdataHandle.AverageObjs.vars{i});
    if userdataHandle.AverageObjs.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.AverageObjs.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.AverageObjs.default{%d});', varName, i);
        eval(strCC);
    end
end

% Import objects
userdataHandle.Import.vars = {'ResampleFactor','OutDirectory'};
userdataHandle.Import.args = [1 200];
userdataHandle.Import.inFilter = {'*.stl', 'Mesh Objects (*.stl)'; ...
    '*.m', 'Surface-mesh Objects (*.m)'; ...
    '*.hdr; *.nii', 'Binary Analyze/NIfTI Objects (*.hdr, *.nii)'};% ...
%    '*.stl; *.m; *.hdr; *.nii', 'All Eligible Files (*.stl, *.m, *.hdr, *.nii)'};
userdataHandle.Import.default = {'1','./'};

for i=1:length(userdataHandle.Import.vars)
    varName = sprintf('userdataHandle.Import.%s', userdataHandle.Import.vars{i});
    if userdataHandle.Import.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.Import.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.Import.default{%d});', varName, i);
        eval(strCC);
    end
end

% Convert file formats
userdataHandle.bim2gipl.inFilter = {'*_bim.mat', 'Bianry Objects (*_bim.mat)'; ...
    '*.mat', 'All Objects (*.mat)'};
userdataHandle.bim2gipl.args=[200];
userdataHandle.bim2gipl.vars = {'OutDirectory'};
userdataHandle.bim2gipl.OutDirectory = './';

userdataHandle.fix2gipl.inFilter = {'*_fix.mat', 'Bianry Objects (*_fix.mat)'; ...
    '*.mat', 'All Objects (*.mat)'};
userdataHandle.fix2gipl.args=[200];
userdataHandle.fix2gipl.vars = {'OutDirectory'};
userdataHandle.fix2gipl.OutDirectory = './';

userdataHandle.gipl2bim.inFilter = {'*.gipl', 'Bianry Objects (*.gipl)'};
userdataHandle.gipl2bim.args=[200];
userdataHandle.gipl2bim.vars = {'OutDirectory'};
userdataHandle.gipl2bim.OutDirectory = './';

userdataHandle.smo2surf_para_meta.inFilter = {'*_smo.mat', 'Parameterized Objects (*_smo.mat)'};
userdataHandle.smo2surf_para_meta.args=[200];
userdataHandle.smo2surf_para_meta.vars = {'OutDirectory'};
userdataHandle.smo2surf_para_meta.OutDirectory = './';

userdataHandle.surf_para_meta2smo.inFilter = {'*_surf.meta;*_para.meta', 'Parameterized Objects (*_surf.meta, *_para.meta)'};
userdataHandle.surf_para_meta2smo.args=[200];
userdataHandle.surf_para_meta2smo.vars = {'OutDirectory'};
userdataHandle.surf_para_meta2smo.OutDirectory = './';

userdataHandle.des2meta_coef.inFilter = {'*_des.mat', 'SPHARM Descriptors (*_des.mat)'};
userdataHandle.des2meta_coef.args=[200];
userdataHandle.des2meta_coef.vars = {'OutDirectory'};
userdataHandle.des2meta_coef.OutDirectory = './';

userdataHandle.meta_coef2des.inFilter = {'*_surfSPHARM.meta;*_surfSPHARM.coef', 'SPHARM Descriptors (*_surfSPHARM.meta, *_surfSPHARM.coef)'};
userdataHandle.meta_coef2des.args=[200];
userdataHandle.meta_coef2des.vars = {'OutDirectory'};
userdataHandle.meta_coef2des.OutDirectory = './';

userdataHandle.ellalign_meta_coef2reg.inFilter = {'*_ellalign.meta;*_ellalign.coef', 'Registered SPHARM Descriptors (*_ellalign.meta, *_ellalign.coef)'};
userdataHandle.ellalign_meta_coef2reg.args=[200];
userdataHandle.ellalign_meta_coef2reg.vars = {'OutDirectory'};
userdataHandle.ellalign_meta_coef2reg.OutDirectory = './';

% reg2procalign_meta option (mesh selection)
userdataHandle.reg2procalign_meta.vars = {'SampleMesh','OutDirectory'};
userdataHandle.reg2procalign_meta.args = [10,200];
userdataHandle.reg2procalign_meta.inFilter = {'*_reg.mat;*_nor.mat', 'Registered Objects (*_reg.mat, *_nor.mat)'};
userdataHandle.reg2procalign_meta.default = {'','./'};

for i=1:length(userdataHandle.reg2procalign_meta.vars)
    varName = sprintf('userdataHandle.reg2procalign_meta.%s', userdataHandle.reg2procalign_meta.vars{i});
    if userdataHandle.reg2procalign_meta.args(i) < 10
        strCC = sprintf('%s=str2num(userdataHandle.reg2procalign_meta.default{%d});', varName, i);
        eval(strCC);
    else
        strCC = sprintf('%s=(userdataHandle.reg2procalign_meta.default{%d});', varName, i);
        eval(strCC);
    end
end

userdataHandle.others.inFilter = {'*.txt', 'List of Objects (*.txt)'};  %'t_map' 'PCA' 't_map+PCA'

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%WATERSHED  Performs a marker-based image segmentation using the watershed algorithm
%
%     marker = cv.watershed(image, marker)
%
% ## Input
% * __image__ Input 8-bit 3-channel image.
% * __marker__ Input 32-bit single-channel image (map) of markers. It should
%   have the same size as `image`.
%
% ## Output
% * __marker__ Output 32-bit single-channel image (map) of markers. It has the
%   same size as `image`.
%
% The function implements one of the variants of watershed, non-parametric
% marker-based segmentation algorithm, described in [Meyer92].
%
% Before passing the image to the function, you have to roughly outline the
% desired regions in the image markers with positive (>0) indices. So, every
% region is represented as one or more connected components with the pixel
% values 1, 2, 3, and so on. Such markers can be retrieved from a binary mask
% using cv.findContours and cv.drawContours. The markers are "seeds" of the
% future image regions. All the other pixels in markers, whose relation to the
% outlined regions is not known and should be defined by the algorithm, should
% be set to 0's. In the function output, each pixel in markers is set to a
% value of the "seed" components or to -1 at boundaries between the regions.
%
% ### Note
% Any two neighbor connected components are not necessarily separated by
% a watershed boundary (-1 pixels); for example, they can touch each other in
% the initial marker image passed to the function.
%
% ## References
% [Meyer92]:
% > Fernand Meyer. "Color image segmentation". In Image Processing and its
% > Applications, 1992, International Conference on, pages 303-306. IET, 1992.
%
% See also: cv.grabCut, cv.findContours, watershed, boundarymask
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function res = imLUT(img, lut)
%IMLUT apply a lut to a gray-scale image.
%
%   IM2 = imLUT(IMG, LUT).
%   IMG is a gray-scale image, 1, 8 or 16 bits image, LUT is a double array
%   with 2**Nbits elements.
%   Each element x in IMG will by replaced by the value of the (x+1)-th
%   element in the LUT.
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 25/10/2004.
%

%   HISTORY

res = zeros(size(img));
for i=0:length(lut)-1
    res(img==i) = lut(i+1);
end

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
% DEMBRENDAN1 Use the GP-LVM to model the Frey face data with back constraints.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'brendan';
experimentNo = 1;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('fitc');
options.optimiser = 'conjgrad';
options.back = 'kbr';
options.backOptions = kbrOptions(Y);
options.backOptions.kern = kernCreate(Y, 'rbf');
options.backOptions.kern.inverseWidth = 0.00001;

latentDim = 2;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
capName = dataSetName;;
capName(1) = upper(capName(1));
save(['dem' capName num2str(experimentNo) '.mat'], 'model');

if exist('printDiagram') & printDiagram
  fgplvmPrintPlot(model, lbls, capName, experimentNo);
end

% Load the results and display dynamically.
fgplvmResultsDynamic(dataSetName, experimentNo, 'image', [20 28], 1, 0, 1)

% compute the nearest neighbours errors in latent space.
errors = fgplvmNearestNeighbour(model, lbls);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function val = simAssetPrice_2(TimeSeries, Period, nYears, nTrials)
% SIMASSETPRICE_2 simulates the asset prices series in a Monte Carlo
% simulation.  Is similar to simAssetPrice_orig, but removes dependence on
% the GBM object (at some expense to readibility).

[expReturn, expCov] = calcExpMoments(TimeSeries, Period);
nAssets = size(expCov,1);
StartAssetPrice = TimeSeries(end, :)';

T = cholcov(expCov);
randNums = randn(nTrials*nYears, nAssets);
randNums = permute(reshape(randNums * T, nTrials, nYears, nAssets), [2,3,1]);
returns = bsxfun(@plus, randNums, expReturn);
StartAssetPrice = repmat(StartAssetPrice', [1,1,nTrials]);
val = cumprod([StartAssetPrice; 1 + returns]);

val(val<0) = 0;

end

function [expReturn, expCov] = calcExpMoments(TimeSeries, Period)
% Estimates and annualizes the expected returns and covariances for the
% given time series
returns     = price2ret(TimeSeries);
expReturn   = mean(returns);
expCov      = cov(returns);

% Annualize monthly, weekly and daily returns and correlations
if Period == 'm'
    [expReturn, expCov] = arith2geom(expReturn, expCov, 12);
elseif Period == 'w'
    % Not implemented
elseif Period == 'd'
    % Not implemented
else
    % Do nothing
end

end

author: Sable

repo: mcbench-benchmarks


text:
addpath(genpath('D:\GitHub\KiloSort2')) % path to kilosort folder
addpath('D:\GitHub\npy-matlab')

pathToYourConfigFile = 'D:\GitHub\KiloSort2\configFiles'; % take from Github folder and put it somewhere else (together with the master_file)
run(fullfile(pathToYourConfigFile, 'configFileBench384.m'))

% common options for every probe
ops.chanMap     = 'D:\GitHub\KiloSort2\configFiles\neuropixPhase3A_kilosortChanMap.mat';
ops.trange      = [0 Inf]; % TIME RANGE IN SECONDS TO PROCESS

 % these settings overwrite any settings from config
ops.Th       = [10 10];     % threshold on projections (like in Kilosort1)
ops.lam      = 40^2;   % weighting on the amplitude penalty (like in Kilosort1, but it has to be much larger)

ops.ThS      = [8 8];  % lower bound on acceptable single spike quality
ops.momentum = [20 400]; % number of samples to average over
ops.minFR    = 1/50; % minimum spike rate (Hz)

ops.sigmaMask  = 30;

ops.Nfilt       = 1024; % max number of clusters
ops.nfullpasses = 3; % how many forward backward passes to do
ops.nPCs        = 3; % how many PCs to project the spikes into

ops.NchanTOT = 385;
ops.useRAM = 0;
ops.ccsplit = .99;

% rootZ = 'D:\DATA\ALLEN\mouse366119\probeC_2018-03-02_15-18-32_SN619041624\experiment1\recording1\continuous\Neuropix-120.0\';
% rootZ = 'H:\DATA\Spikes\Robbins\ZNP1';
rootZ = 'H:\DATA\Spikes\WillAllen';
rootH = 'H:\DATA\Spikes\temp\';

fs = dir(fullfile(rootZ, '*.bin'));

fname = fs(1).name;

ops.fbinary     = fullfile(rootZ,  fname);
ops.fproc       = fullfile(rootH, 'temp_wh2.dat'); % residual from RAM of preprocessed data

% preprocess data
rez = preprocessDataSub(ops);

rez = clusterSingleBatches(rez);

% figure(191);
% imagesc(rez.ccb(rez.iorig, rez.iorig), [20 100])
%%
learnAndSolve8
% rez = learnAndSolve8(rez);

rez2    = splitAllClusters(rez);

rezToPhy(rez2, rootZ);

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
classdef nnbnorm < nntest
  properties (TestParameter)
    rows = {2 8 13}
    cols = {2 8 17}
    numDims = {1 3 4}
    batchSize = {2 7}
  end
  methods (Test)
    function basic(test, rows, cols, numDims, batchSize)
      r = rows ;
      c = cols ;
      nd = numDims ;
      bs = batchSize ;
      x = test.randn(r, c, nd, bs, 'single') ;
      g = test.randn(1, 1, nd, 1, 'single');
      b = test.randn(1, 1, nd, 1, 'single');
      g = test.randn(nd, 1, 'single');
      b = test.randn(nd, 1, 'single');

      y = vl_nnbnorm(x,g,b) ;
      dzdy = test.randn(size(y), 'single') ;
      [dzdx,dzdg,dzdb] = vl_nnbnorm(x,g,b,dzdy) ;

      test.der(@(x) vl_nnbnorm(x,g,b), x, dzdy, dzdx, test.range * 1e-2, -1e-3) ;
      test.der(@(g) vl_nnbnorm(x,g,b), g, dzdy, dzdg, test.range * 1e-2, -1e-3) ;
      test.der(@(b) vl_nnbnorm(x,g,b), b, dzdy, dzdb, test.range * 1e-2, -1e-3) ;
    end
  end
end

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function q = cubic_integrate ( a, b )

%*****************************************************************************80
%
%% CUBIC_INTEGRATE integrates the cubic from A to B.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 February 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the integration interval.
%
%    Output, real Q, the integral from A to B.
%
  q = cubic_antiderivative ( b ) - cubic_antiderivative ( a );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ppt, op] = pptOpen(pptFile);
% For Windows machines only: Create an Active X Server object used
% for talking with PowerPoint, writing to the specified powerpoint
% presentation. Creates it if needed, otherwise opens and gets ready
% to append to it.
%
% Usage:
%  [ppt, op] = pptOpen(pptFile);
%
% pptFile: path to the powerpoint file to open.
% ppt: ActiveX Server for PowerPoint.
% op: the operation for the specific ppt file.
%
% ras, 10/2005.
if ~ispc
    warning('Sorry, pptOpen is Windows-only right now.');
    return
end

% Start an ActiveX session with PowerPoint:
ppt = actxserver('PowerPoint.Application');
ppt.Visible = 1;

if ~exist(pptFile,'file');
  % Create new presentation:
  op = invoke(ppt.Presentations,'Add');
else
  % Open existing presentation:
  op = invoke(ppt.Presentations,'Open',pptFile);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% [B1_] = tensor_fibers_column(A)
%
function [B1_] = tensor_fibers_column(A)
  [n1, n2, n3] = size(A);
  for j = 1:n2
    for k = 1:n3
      B1_{j,k} = A(:,j,k);
    end
  end
end

author: andrewssobral

repo: lrslibrary


text:
function h=ref_pfilt(f,g,a)
%REF_PFILT  Reference pfilt handling structs
%   Usage:  h=pfilt(f,g);
%           h=pfilt(f,g,a,dim);
%
%   `pfilt(f,g)` applies the filter *g* to the input *f*. If *f* is a
%   matrix, the filter is applied along each column.
%
%   `pfilt(f,g,a)` does the same, but downsamples the output keeping only
%   every a'th sample (starting with the first one).
%
%   `pfilt(f,g,a,dim)` filters along dimension dim. The default value of
%   [] means to filter along the first non-singleton dimension.
%
%   The filter *g* can be a vector, in which case the vector is treated
%   as a zero-delay FIR filter.
%
%   The filter *g* can be a cell array. The following options are
%   possible:
%
%     * If the first element of the cell array is the name of one of the
%       windows from |firwin|, the whole cell array is passed onto
%       |firfilter|.
%
%     * If the first element of the cell array is `'bl'`, the rest of the
%     cell array is passed onto |blfilter|.
%
%     * If the first element of the cell array is `'pgauss'`, `'psech'`,
%       the rest of the parameters is passed onto the respective
%       function. Note that you do not need to specify the length *L*.
%
%   The coefficients obtained from filtering a signal *f* by a filter *g* are
%   defined by
%
%   ..          L-1
%      c(n+1) = sum f(l+1) * g(an-l+1)
%               l=0
%
%   .. math:: c\left(n+1\right)=\sum_{l=0}^{L-1}f\left(l+1\right)g\left(an-l+1\right)
%
%   where $an-l$ is computed modulo $L$.
%
%   See also: pconv

if nargin<3
    a=1;
end;

[L,W]=size(f);

l=(0:L-1).'/L;
if isstruct(g)
    if isfield(g,'h')
        g_time=circshift(postpad(g.h,L),g.offset).*exp(2*pi*1i*(round(g.fc*L/2))*l);
        G=fft(g_time);
    elseif isfield(g,'H')
        G=circshift(postpad(g.H(L),L),g.foff(L)).*exp(-2*pi*1i*round(g.delay)*l);
    else
       error('%s: Unknown filter definition.',upper(mfilename));
    end;

else
    G=fft(fir2long(g,L));
end;

if numel(a) > 1
    % This is possibly a fractional subsampling

    afrac = a(1)/a(2);
    if a(1) ~= L
        error('%s: The length in a(1) is not equal to L.',upper(mfilename));
    end
    N = L/afrac;
    if abs(N-round(N))>1e-5
        error('%s: Output length is not integer.',upper(mfilename));
    else
        N = round(N);
    end

    foff = g.foff(L);

    h=zeros(N,W,assert_classname(f));
    for w=1:W
        h(:,w) = blfilt(f(:,w),G,N,foff,afrac);

        if isstruct(g) && isfield(g,'realonly') && g.realonly
            G2 = involute(G);
            supp = numel(g.H(L));
            foff2 = -L+mod(L-foff-supp,L)+1;
            h(:,w) = (h(:,w) + blfilt(f(:,w),G2,N,foff2,afrac))/2;
        end;

    end;

else
    % This is regular subsampling case
    if isstruct(g) && isfield(g,'realonly') && g.realonly
        G=(G+involute(G))/2;
    end;

    N=L/a;
    h=zeros(N,W,assert_classname(f));
    for w=1:W
        F=fft(f(:,w));
        h(:,w)=ifft(sum(reshape(F.*G,N,a),2))/a;
    end;
end

function h = blfilt(f,G,N,foff,afrac)
L = size(f,1);

align = foff - floor(foff/N)*N;

F = circshift(fft(f).*G,-foff);

F = postpad(F,ceil(L/N)*N);

F = circshift(F,align);

F = sum(reshape(F,N,numel(F)/N),2);

h = ifft(F)/afrac;

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function bg = ball_grid ( n, r, c, ng )

%*****************************************************************************80
%
%% BALL_GRID computes grid points inside a ball.
%
%  Discussion:
%
%    The grid is defined by specifying the radius and center of the ball,
%    and the number of subintervals N into which the horizontal radius
%    should be divided.  Thus, a value of N = 2 will result in 5 points
%    along that horizontal line.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of subintervals.
%
%    Input, real R, the radius of the ball.
%
%    Input, real C(3), the coordinates of the center of the ball.
%
%    Input, integer NG, the number of grid points, as determined by
%    BALL_GRID_COUNT.
%
%    Output, real BG(3,NG), the grid points inside the ball.
%
  bg = zeros ( 3, ng );

  p = 0;

  for i = 0 : n
    x = c(1) + r * 2 * i / ( 2 * n + 1 );

    for j = 0 : n
      y = c(2) + r * 2 * j / ( 2 * n + 1 );

      for k = 0 : n

        z = c(3) + r * 2 * k / ( 2 * n + 1 );

        if ( r * r < ( x - c(1) ).^2 ...
                   + ( y - c(2) ).^2 ...
                   + ( z - c(3) ).^2 )
          break
        end
        p = p + 1;
        bg(1:3,p) = [ x, y, z ]';
        if ( 0 < i )
          p = p + 1;
          bg(1:3,p) = [ 2 * c(1) - x, y, z ]';
        end
        if ( 0 < j )
          p = p + 1;
          bg(1:3,p) = [ x, 2 * c(2) - y, z ]';
        end
        if ( 0 < k )
          p = p + 1;
          bg(1:3,p) = [ x, y, 2 * c(3) - z ]';
        end
        if ( 0 < i && 0 < j )
          p = p + 1;
          bg(1:3,p) = [ 2 * c(1) - x, 2 * c(2) - y, z ]';
        end
        if ( 0 < i && 0 < k )
          p = p + 1;
          bg(1:3,p) = [ 2 * c(1) - x, y, 2 * c(3) - z ]';
        end
        if ( 0 < j && 0 < k )
          p = p + 1;
          bg(1:3,p) = [ x, 2 * c(2) - y, 2 * c(3) - z ]';
        end
        if ( 0 < i && 0 < j && 0 < k )
          p = p + 1;
          bg(1:3,p) = [ 2 * c(1) - x, 2 * c(2) - y, 2 * c(3) - z ]';
        end
      end
    end
  end

  return
end


text:
%% Machine Learning  Online Class === Exercise 7:K-Means Clustering and Principle Component Analysis
close all;clc

%% Optional (ungraded) Exercise: PCA for Visualization--1
A = double(imread('lena.png'));
A = A / 255;
[m,n,l] = size(A);
X = reshape(A,m*n, l);
K = 16;
max_iters = 10;
init_centroids = kMeansInitCentroids(X, K);
[centroids, idx] = runkMeans(X, init_centroids, max_iters);

palette = hsv(K);
colors = palette(idx, :);

figure(1);
scatter3(X(:, 1), X(:, 2), X(:,3), 10, colors);
title('Pixel dataset plotted in 3D. Color shows centroid memberships');
hold off;

%% Optional (ungraded) Exercise: PCA for Visualization--2
[X_norm, mu, sigma] = featureNormalize(X);

[U, S] = pca(X_norm);
Z = projectData(X_norm, U, 2);

figure(2);
plotDataPoints(Z,idx, K);
title('Pixel dataset plotted in 2D, using PCA for dimensionality reduction');

author: ShiMengjie

repo: Machine-Learning-Andrew-Ng


text:
function uv=spherAng2Uv(azEl,systemType,includeW,Ms,Muv)
%%SPHERANG2UV Convert azimuth and elevation in spherical coordinates into a
%             [u;v] or [u;v;w] direction cosine vector in 3D. Direction
%             cosines u and v are the first two elements of a unit vector
%             in 3D.
%
%INPUTS: azEl A 2XN set of N points in [azimuth;elevation] in radians to
%           convert to direction cosines.
%  systemType An optional parameter specifying the axis from which the
%           angles are measured in radians. Possible values are
%           0 (The default if omitted) Azimuth is measured
%             counterclockwise from the x-axis in the x-y plane.
%             Elevation is measured up from the x-y plane (towards the
%             z-axis). This is consistent with common spherical
%             coordinate systems for specifying longitude (azimuth) and
%             geocentric latitude (elevation).
%           1 Azimuth is measured counterclockwise from the z-axis in the
%             z-x plane. Elevation is measured up from the z-x plane
%             (towards the y-axis). This is consistent with some spherical
%             coordinate systems that use the z-axis as the boresight
%             direction of the radar.
%           2 This is the same as 0 except instead of being given
%             elevation, one is given the angle away from the z-axis, which
%             is (pi/2-elevation).
%           3 This is the same as 0 except azimuth is measured clockwise
%             from the y-axis in the x-y plane instead of counterclockwise
%             from the x-axis. This coordinate system often arises when
%             given "bearings" in a local East-North-Up coordinate system,
%             where the bearing directions are measured East of North.
%  includeW An optional boolean value indicating whether a third direction
%           cosine component should be included. The u and v direction
%           cosines are two parts of a 3D unit vector. Generally, one might
%           assume that the target is in front of the sensor, so the third
%           component would be positive and is not needed. However, the
%           third component can be included if ambiguity exists. The
%           default if this parameter is omitted or an empty matrix is
%           passed is false.
%    Ms,Muv If either the spherical coordinate system or the u-v coordinate
%           system is rotated compared to the global Cartesian coordinate
%           system, these optional 3X3 matrices provide the rotations. Ms
%           is a 3X3 matrix to go from the alignment of a global
%           Cartesian coordinate system to that in which the spherical
%           coordinates are computed. Similarly, Muv is a rotation matrix
%           to go from the alignment of a global Cartesian cordinate system
%           to that in which the u-v(-w) coordinates are computed. If
%           either of these in omitted or an empty matrix is passed, then
%           the missing one is replaced with the identity matrix.
%
%OUTPUTS: uv A 2XN (without w) or 3XN (with w) set of direction cosines
%            values corresponding to the specified angles.
%
%Direction cosines and spherical coordinate systems are discussed in [1].
%
%EXAMPLE:
%In this example, a spherical value is converted to uv coordiantes and then
%back, demonostrating the consistency of the functions spherAng2Uv and
%uv2SpherAng. The relative error should be about zero.
% zSpher=[0.4;0.7];
% systemType=3;
% includeW=true;
% Ms=Euler2Ang2RotMat(0.1,0.2,'xy');
% Muv=Euler1Ang2RotMat(0.25,'z');
% convBack=uv2SpherAng(spherAng2Uv(zSpher,systemType,includeW,Ms,Muv),systemType,Ms,Muv);
% relErr=max(abs((convBack-zSpher)./zSpher))
%
%REFERENCES:
%[1] D. F. Crouse, "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%
%June 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<5||isempty(Muv))
    Muv=eye(3,3);
end

if(nargin<4||isempty(Ms))
    Ms=eye(3,3);
end

if(nargin<3||isempty(includeW))
    includeW=false;
end

if(nargin<2||isempty(systemType))
    systemType=0;
end

azimuth=azEl(1,:);
elevation=azEl(2,:);

if(systemType==2)
    elevation=pi/2-elevation;
    systemType=0;
elseif(systemType==3)
    azimuth=pi/2-azimuth;
    systemType=0;
end

switch(systemType)
    case 0
        uv=Muv*Ms'*[cos(azimuth).*cos(elevation);
                    sin(azimuth).*cos(elevation);
                    sin(elevation)];
    case 1
        uv=Muv*Ms'*[sin(azimuth).*cos(elevation);
                    sin(elevation);
                    cos(azimuth).*cos(elevation)];
    otherwise
        error('Invalid system type specified.')
end

if(includeW==false)
    uv=uv(1:2,:);
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function cur_shift = WardGetExpShift(img1, img2, ward_percentile, shift_bits)
%
%
%       cur_shift = WardGetExpShift(img1, img2, ward_percentile, shift_bits)
%
%       This function computes the Ward's MTB.
%
%       Input:
%           -img1: the target image
%           -img2: the image that needs to be aligned to img1
%
%       Output:
%           -cur_shift: shifting vector for aligning img2 into img1.
%
%     Copyright (C) 2012-15  Francesco Banterle
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%

if(~exist('ward_percentile', 'var'))
    ward_percentile = 0.5;
end

if(~exist('shift_bits', 'var'))
    shift_bits = 6;
end

cur_shift = zeros(2, 1);
shift_ret = zeros(2, 1);

while(shift_bits > 0)
    %computing MTB
    sml_img1 = imresize(img1, 2^(-shift_bits), 'bilinear');
    sml_img2 = imresize(img2, 2^(-shift_bits), 'bilinear');

    [tb1, eb1] = WardComputeThreshold(sml_img1, ward_percentile);
    [tb2, eb2] = WardComputeThreshold(sml_img2, ward_percentile);

    min_err = size(sml_img1, 1) * size(sml_img1, 2);

    tb1 = logical(tb1);
    eb1 = logical(eb1);

    for i=-1:1
        for j=-1:1
            xs = cur_shift(1) + i;
            ys = cur_shift(2) + j;

            shifted_tb2 = logical(imShift(tb2, [xs, ys]));
            shifted_eb2 = logical(imShift(eb2, [xs, ys]));

            diff_b = bitxor(tb1, shifted_tb2);
            diff_b = diff_b & eb1;
            diff_b = diff_b & shifted_eb2;

            err = sum(diff_b(:));

            if (err < min_err)
                shift_ret = [xs; ys];
                min_err = err;
            end
        end
    end

    shift_bits = shift_bits - 1;
    cur_shift  = shift_ret * 2;
end

end

author: banterle


text:
function cardan_poly_coef_test ( )

%*****************************************************************************80
%
%% CARDAN_POLY_COEF_TEST tests CARDAN_POLY_COEF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  n_max = 10;

  s = 1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CARDAN_POLY_COEF_TEST\n' );
  fprintf ( 1, '  CARDAN_POLY_COEF returns the coefficients of a\n' );
  fprintf ( 1, '  Cardan polynomial.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  We use the parameter S = %f\n', s );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Table of polynomial coefficients:\n' );
  fprintf ( 1, '\n' );

  for n = 0 : n_max
    c = cardan_poly_coef ( n, s );
    fprintf ( 1, '  %2d  ', n );
    for i = 0 : n
      fprintf ( 1, '  %9f', c(i+1) );
    end
    fprintf ( 1, '\n' );
  end

  return
end


text:
function varargout = svmAverageGUI(varargin)
% svmAverageGUI
%   Given an svm structure, opens a GUI displaying the contents with the
%   ability to select, modify, and average across fields describing the
%   data.
%
% Usage:
%   svmAverageGUI(svm)
%
% Proper comments and ability to return the structure (even more
% important!) forthcoming.
%
% [renobowen@gmail.com 2010]
%

    % Begin initialization code - DO NOT EDIT
    gui_Singleton = 1;
    gui_State = struct('gui_Name',       mfilename, ...
                       'gui_Singleton',  gui_Singleton, ...
                       'gui_OpeningFcn', @svmAverageGUI_OpeningFcn, ...
                       'gui_OutputFcn',  @svmAverageGUI_OutputFcn, ...
                       'gui_LayoutFcn',  [] , ...
                       'gui_Callback',   []);
    if nargin && ischar(varargin{1})
        gui_State.gui_Callback = str2func(varargin{1});
    end

    if nargout
        [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
    else
        gui_mainfcn(gui_State, varargin{:});
    end
    % End initialization code - DO NOT EDIT
end

function svmAverageGUI_OpeningFcn(hObject, eventdata, handles, svm, varargin)
    if ~(exist('svm','var'))
        svm = [];
    end

    % Choose default command line output for svmAverageGUI
    handles.output      = hObject;
    handles.svm         = svm;
    handles.saveToVar   = [];

    for i = 1:2:length(varargin)
        switch (lower(varargin{i}))
            case {'savetovar'}
                handles.saveToVar = varargin{i + 1};
            otherwise
                fprintf(1, 'Unrecognized option ''%s''\n', varargin{i});
        end
    end

    handles.changed     = 0;

    % Update handles structure
    guidata(hObject, handles);

    if (~isempty(svm))
        svmAverageGUI_InitFcn(hObject, handles);
    end
end

function svmAverageGUI_InitFcn(hObject, handles)
    handles = DataList_InitFcn(hObject, handles);
    handles = RunList_InitFcn(hObject, handles);
    handles = GroupList_InitFcn(hObject, handles);
    handles = TrialList_InitFcn(hObject, handles);
    guidata(hObject, handles);
end

function varargout = svmAverageGUI_OutputFcn(hObject, eventdata, handles)

    if (isempty(handles.svm))
        fprintf(1, 'No svm structure specified - closing GUI.\n');
        close(handles.AverageGUI);
    end

end

function DataList_UpdateSelections(hObject, handles)
    runValues       = unique(handles.svm.run);
    groupValues     = unique(handles.svm.group);
    trialValues     = unique(handles.svm.trial);

    selectedRuns    = runValues(get(handles.RunList, 'Value'));
    selectedGroups  = groupValues(get(handles.GroupList, 'Value'));
    selectedTrials  = trialValues(get(handles.TrialList, 'Value'));

    mask =  ismember(handles.svm.run, selectedRuns) & ...
            ismember(handles.svm.group, selectedGroups) & ...
            ismember(handles.svm.trial, selectedTrials);

    set(handles.DataList, 'Value', find(mask));
    guidata(hObject, handles);
end

function handles = DataList_InitFcn(hObject, handles)
    strings = num2str([handles.svm.run handles.svm.group handles.svm.trial]);
    set(handles.DataList, 'String', strings);
    set(handles.DataList, 'Max', size(strings,1));
    set(handles.DataList, 'Value', 1);
    set(handles.DataList, 'ListboxTop', 1);

end

function handles = RunList_InitFcn(hObject, handles)
    strings = num2str(unique(handles.svm.run));
    set(handles.RunList, 'String', strings);
    set(handles.RunList, 'Max', size(strings,1));
    set(handles.RunList, 'Value', 1);
    set(handles.RunList, 'ListboxTop', 1);
end

function handles = GroupList_InitFcn(hObject, handles)
    strings = handles.svm.grouplabel;
    set(handles.GroupList, 'String', strings);
    set(handles.GroupList, 'Max', length(strings));
    set(handles.GroupList, 'Value', 1);
    set(handles.GroupList, 'ListboxTop', 1);
end

function handles = TrialList_InitFcn(hObject, handles)
    strings = num2str(unique(handles.svm.trial));
    set(handles.TrialList, 'String', strings);
    set(handles.TrialList, 'Max', size(strings,1));
    set(handles.TrialList, 'Value', 1);
    set(handles.TrialList, 'ListboxTop', 1);

end

function ResetDetailSelectors(hObject, handles)
    set(handles.RunList, 'Value', 1);
    set(handles.RunList, 'ListboxTop', 1);
    set(handles.GroupList, 'Value', 1);
    set(handles.GroupList, 'ListboxTop', 1);
    set(handles.TrialList, 'Value', 1);
    set(handles.TrialList, 'ListboxTop', 1);
    guidata(hObject, handles);
end

function UpdateLists(hObject, handles)
    svmAverageGUI_InitFcn(hObject, handles);
end

function handles = RepairIndices(hObject, handles)
    nLabels = length(handles.svm.grouplabel);
    markedForDeletion = zeros(1,nLabels);
    subFromInds = 0;
    for i = 1:length(handles.svm.grouplabel)
        if (sum(ismember(handles.svm.group, i)) == 0)
            markedForDeletion(i) = 1;
            subFromInds = subFromInds + 1;
        else
            inds = ismember(handles.svm.group, i);
            handles.svm.group(inds) = i - subFromInds;
        end
    end

    handles.svm.grouplabel = handles.svm.grouplabel(~markedForDeletion);
end

%% List callbacks
function DataList_Callback(hObject, eventdata, handles)
% DataList_Callback
%   Selecting within this window resets the selected run/group/trial in the
%   other lists (we're now manually editing the selection and the remaining
%   lists will no longer contain info representative of what we've
%   selected.

    ResetDetailSelectors(hObject, handles);
end

function RunList_Callback(hObject, eventdata, handles)
    DataList_UpdateSelections(hObject, handles);
end

function GroupList_Callback(hObject, eventdata, handles)
    DataList_UpdateSelections(hObject, handles);
end

function TrialList_Callback(hObject, eventdata, handles)
    DataList_UpdateSelections(hObject, handles);
end

%% Action button callbacks - used to modify data
function RelabelButton_Callback(hObject, eventdata, handles)
    % hObject    handle to AverageButton (see GCBO)
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)

    run         = str2num(get(handles.RunInput, 'String'));
    trial       = str2num(get(handles.TrialInput, 'String'));
    groupLabel  = get(handles.GroupInput, 'String');

    if (isempty(run) && isempty(trial) && isempty(groupLabel))
        fprintf(1, 'Empty assignment fields!\n');
        return;
    end

    mask        = false(get(handles.DataList, 'Max'),1);
    inds        = get(handles.DataList, 'Value');
    mask(inds)  = true;

    if (~isempty(run))
        handles.svm.run(mask) = run;
    end

    if (~isempty(trial))
        handles.svm.trial(mask) = trial;
    end

    if (~isempty(groupLabel))
        group = cellfind(handles.svm.grouplabel, groupLabel);
        if (isempty(group))
            group = max(handles.svm.group) + 1;
            handles.svm.grouplabel{group} = groupLabel;
        end
        handles.svm.group(mask) = group;
        handles = RepairIndices(hObject, handles);
    end

    handles.changed = 1;
    UpdateLists(hObject, handles);
end

function AverageButton_Callback(hObject, eventdata, handles)

    run         = str2num(get(handles.RunInput, 'String'));
    trial       = str2num(get(handles.TrialInput, 'String'));
    groupLabel  = get(handles.GroupInput, 'String');

    if (isempty(run) || isempty(trial) || isempty(groupLabel))
        fprintf(1, 'Empty assignment fields!\n');
        return;
    end

    group = cellfind(handles.svm.grouplabel, groupLabel);
    if (isempty(group))
        group = max(handles.svm.group) + 1;
        handles.svm.grouplabel{group} = groupLabel;
    end

    mask        = false(get(handles.DataList, 'Max'),1);
    inds        = get(handles.DataList, 'Value');
    mask(inds)  = true;

    newsvm          = handles.svm;
    newsvm.data     = newsvm.data(~mask,:);
    newsvm.run      = newsvm.run(~mask,:);
    newsvm.group    = newsvm.group(~mask,:);
    newsvm.trial    = newsvm.trial(~mask,:);

    newsvm.data     = [newsvm.data; mean(handles.svm.data(mask,:), 1)];
    newsvm.run      = [newsvm.run; run];
    newsvm.group    = [newsvm.group; group];
    newsvm.trial    = [newsvm.trial; trial];

    handles.svm     = newsvm;
    handles         = RepairIndices(hObject, handles);
    handles.changed = 1;
    UpdateLists(hObject, handles);
end

function SaveAsButton_Callback(hObject, eventdata, handles)
    if (handles.changed)
        default = {'svm'};
        if (~isempty(handles.saveToVar))
            default = handles.saveToVar;
        end
        varName = inputdlg('Assign to variable: ', 'svmAverageGUI', 1, default);
        if (isempty(varName)) return; end

        assignin('base', varName{1}, handles.svm);

        fprintf(1, 'Saved data to variable ''%s''.\n', varName{1});

        handles.changed = 0;
        guidata(hObject, handles);
    else
        fprintf(1, 'No changes to save.\n');
    end
end

function SaveButton_Callback(hObject, eventdata, handles)
    if (handles.changed)
        if (isempty(handles.saveToVar))
            handles.saveToVar = inputdlg('Assign to variable: ', 'svmAverageGUI', 1, {'svm'});
            if (isempty(handles.saveToVar)) return; end
        end
        assignin('base', handles.saveToVar{1}, handles.svm);

        fprintf(1, 'Saved data to variable ''%s''.\n', handles.saveToVar{1});

        handles.changed = 0;
        guidata(hObject, handles);
    else
        fprintf(1, 'No changes to save.\n');
    end
end

function AverageGUI_DeleteFcn(hObject, eventdata, handles)
    if (handles.changed)
        buttonChosen = questdlg('Save changes?', 'svmAverageGUI', 'Yes', 'No', 'Yes');
        if (strcmp(buttonChosen,'Yes'))
            SaveButton_Callback(hObject, eventdata, handles);
        end
    end
end

function HelpButton_Callback(hObject, eventdata, handles)
    help svmAverageGUI;
end

%% Generated by GUIDE, necessary to run
% Create functions
function DataList_CreateFcn(hObject, eventdata, handles), end
function RunList_CreateFcn(hObject, eventdata, handles), end
function GroupList_CreateFcn(hObject, eventdata, handles), end
function TrialList_CreateFcn(hObject, eventdata, handles), end
function RunInput_CreateFcn(hObject, eventdata, handles), end
function GroupInput_CreateFcn(hObject, eventdata, handles), end
function TrialInput_CreateFcn(hObject, eventdata, handles), end

% Callbacks I didn't need
function RunInput_Callback(hObject, eventdata, handles), end
function TrialInput_Callback(hObject, eventdata, handles), end
function GroupInput_Callback(hObject, eventdata, handles), end


text:
function keep = newtonpolytope(exponent_m,exponent_p)
%NEWTONPOLYTOPE  Internal function to remove monimials in SOS programs using Newton polytope

% WARNING : THIS CODE SUCKS AND IS ONLY USED AS A BACK UP PLAN
% IF EVERYTHING ELSE FAILS (CRASHING LP SOLVERS ETC)

% *************************************
% TRY TO CALCULATE CONVEX HULL
% *************************************
try
    cnvhull = convhulln(full(exponent_p));
catch
    keep = 1:size(exponent_m,1);
    return
end

% ***************************************
% GET THE UNIQUE POINTS OF IN CONVEX HULL
% ***************************************
unique_points = unique(cnvhull);

% ***************************************
% CALCULATE A POINT IN INTERIOR
% ***************************************
p_c = sum(exponent_p(unique_points,:),1)'/length(unique_points);

% ***************************************
% CALCULATE HYPER-PLANES Ai^T(x-bi)=0
% ***************************************
j = 1;
A = [];
for i = 1:size(cnvhull,1)
    X = exponent_p(cnvhull(i,:),:)';
    y = X(:,1);
    dX = X(:,2:end)-repmat(X(:,1),1,size(X,2)-1);
    Atemp = null(full(dX'));
    % Is this a full-dimensional facet
    if size(Atemp,2)==1
        direction = (p_c-y)'*Atemp;
        if direction > 0
           A{j}=-Atemp;
       else
           A{j}=Atemp;
       end
        b{j}=y;
        j=j+1;
    end
end

% ***************************************
% CHECK IF MONOMIAL IS IN NEWTON POLYTOPE
% ***************************************
if ~isempty(A)
    keep = [];
    for j = 1:size(exponent_m,1)
        inside = 1;
        y = exponent_m(j,:)';
        if isempty(findrows(exponent_p,y)) % Numerically safe on border
            i = 1;
            while inside & (i<=length(A))
                inside = inside & ((A{i}'*(y-b{i}))<=1e-9);
                i = i+1;
            end
        end
        if inside
            keep = [keep j];
        end
    end
else
    keep = 1:size(exponent_m,1);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [err,time,solver,eqn] = femStokes(node,elem,pde,bdFlag,option,varargin)
%% FEMSTOKES solve the Stokes equation by various finite element methods
%
%   FEMSTOKES computes approximations to the Stokes equation on a
%   sequence of meshes obtained by uniform refinement of a input mesh.
%
% Created by Ming Wang at Nov., 2012.
%
% See also femPoisson femrateStokes
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

%% Default setting of mesh and pde data
if ~exist('node','var') || ~exist('elem','var')
    [node,elem] = squaremesh([0,1,0,1],0.125);  % default mesh is a square
end
if ~exist('option','var'), option = []; end
if ~exist('pde','var')
    pde = Stokesdata1;                          % default data
end
if ~exist('bdFlag','var')
    bdFlag = setboundary(node,elem,'Dirichlet');
end

%% Parameters
option = femoption(option);
maxIt = option.maxIt;   maxN = option.maxN; L0 = option.L0;
option = femStokesoption(option);
elemType = option.elemType; refType = option.refType;

%% Generate an initial mesh
for k = 1:L0
    if strcmp(option.refType,'red')
        [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
    elseif strcmp(option.refType,'bisect')
        [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
    end
end

%% Initialize err
erruH1 = zeros(maxIt,1); errpL2 = zeros(maxIt,1);
errTime = zeros(maxIt,1); solverTime = zeros(maxIt,1);
assembleTime = zeros(maxIt,1); meshTime = zeros(maxIt,1);
itStep = zeros(maxIt,1);  stopErr = zeros(maxIt,1); flag = zeros(maxIt,1);
N = zeros(maxIt,1);

%% Finite Element Method
for k = 1:maxIt
    % solve the equation
    switch upper(elemType)
        case 'P2P1'
            [u,p,edge,A,eqn,info] = StokesP2P1(node,elem,pde,bdFlag,option);
        case 'P2P0'
            [u,p,edge,A,eqn,info] = StokesP2P0(node,elem,pde,bdFlag,option);
        case 'ISOP2P1'
            [u,p,edge,A,eqn,info] = StokesisoP2P1(node,elem,pde,bdFlag,option);
        case 'ISOP2P0'
            [u,p,edge,A,eqn,info] = StokesisoP2P0(node,elem,pde,bdFlag,option);
        case 'CRP0'
            [u,p,edge,A,eqn,info] = StokesCRP0(node,elem,pde,bdFlag,option);
        case 'CRP1'
            [u,p,edge,A,eqn,info] = StokesCRP1(node,elem,pde,bdFlag,option);
        case 'MINI'
            [u,p,edge,A,eqn,info] = StokesMini(node,elem,pde,bdFlag,option);
        case 'P1BP1'
            [u,p,edge,A,eqn,info] = StokesP1bP1(node,elem,pde,bdFlag,option);
    end
    % compute error
    tic;
%     if strcmp(elemType,'P1BP1')
%         uI = Lagrangeinterpolate(pde.exactu,node,elem);
%     elseif strcmp(elemType(1:2),'CR')
%         uI = Lagrangeinterpolate(pde.exactu,node,elem,'CR',edge);
%     end
    if strcmp(elemType,'P2P0') || strcmp(elemType,'P2P1') || ...
       strcmp(elemType,'isoP2P0') || strcmp(elemType,'isoP2P1')
        uI = pde.exactu([node; (node(edge(:,1),:)+node(edge(:,2),:))/2]);
    elseif strcmp(elemType,'CRP0') || strcmp(elemType,'CRP1')
        uI = pde.exactu((node(edge(:,1),:)+node(edge(:,2),:))/2);
    elseif strcmp(elemType,'Mini')
        uI = pde.exactu(node);
    elseif strcmp(elemType,'P1bP1')
        % bubble part won't be taken in the error computation
        Nv = size(node,1); NT = size(elem,1);
        u0 = pde.exactu(node);
        uI = u;
        uI([(1:Nv)'; NT+Nv+(1:Nv)']) = u0(:);
    end
    erruH1(k) = sqrt((u-uI(:))'*A*(u-uI(:)));
    errpL2(k) = getL2error(node,elem,pde.exactp,p);
    errTime(k) = toc;
    % record time
    solverTime(k) = info.solverTime;
    assembleTime(k) = info.assembleTime;
    if option.printlevel>1
        fprintf('Time to compute the error %4.2g s \n H1 err %4.2g    L2err %4.2g \n',...
            errTime(k),erruH1(k), errpL2(k));
    end
    % record solver information
    itStep(k) = info.itStep;
%    stopErr(k) = info.stopErr;
%    flag(k) = info.flag;
    % plot
    N(k) = size(node,1);
    if option.plotflag && N(k) < 2e3 % show mesh and solution for small size
        if length(p) == size(elem,1) % piecewise constant function
            p = recoverP02P1(node,elem,p);
        end
        figure(1);  showresult(node,elem,p);
    end
    if N(k) > maxN
        break;
    end
    % refine mesh
    tic;
    if strcmp(refType,'red')
        [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
    elseif strcmp(refType,'bisect')
        [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
    end
    meshTime(k) = toc;
end
h = 1./sqrt(N(1:k));

%% Plot convergence rates
if option.rateflag
    figure;
    set(gcf,'Units','normal');
    set(gcf,'Position',[0.25,0.25,0.55,0.4]);
    subplot(1,2,1)
    showrateh(h,erruH1(1:k),1,'k-+','|u_I-u_h|_1');
    subplot(1,2,2)
    showrateh(h,errpL2(1:k),1,'m-+','||p-p_h||');
end

%% Output
err = struct('N',N,'uH1',erruH1(1:k),'pL2',errpL2(1:k));
time = struct('N',N,'err',errTime(1:k),'solver',solverTime(1:k), ...
              'assmble',assembleTime(1:k),'mesh',meshTime(1:k));
solver = struct('N',N(1:k),'itStep',itStep(1:k),'time',solverTime(1:k),...
                'stopErr',stopErr(1:k),'flag',flag(1:k));

%% Display error
ts = zeros(k,6); ts = char(ts);
display('#nodes     |u_I-u_h|_1      ||p-p_h||  ');
% format shorte
display([num2str(err.N) ts num2str(err.uH1,'%0.5e') ts num2str(err.pL2,'%0.5e')]);

repo: Aero-matlab


text:
function [env] = env_secant(x_data, y_data, view, side)
% Function call: env_secant(x_data, y_data, view, side)
% Calculates the top envelope of data <y_data> over <x_data>.
% Method used: 'secant-method'
% env_secant() observates the max. slope of about <view> points,
% and joints them to the resulting envelope.
% An interpolation over original x-values is done finally.
% <side> ('top' or 'bottom') defines which side to evolve.
% Author: Andreas Martin, Volkswagen AG, Germany

side = strcmpi( {'top','bottom'}, side ) * [ 1 ; -1 ];

assert( view > 1, ...
       'Parameter <view> too small!' );
assert( ndims (x_data) == 2, ...
       'Parameter <x_data> has to be vector type!' );
assert( size (x_data, 1) == 1 || size (x_data, 2) == 1, ...
       'Parameter <x_data> has to be vector type (Nx1)!' );
assert( ndims (y_data) == 2, ...
       'Parameter <y_data> has to be vector type (Nx1)!' );
assert( size (y_data, 1) == 1 || size (y_data, 2) == 1, ...
       'Parameter <y_data> has to be vector type (Nx1)!' );
assert( length (x_data) == length (y_data), ...
       'Parameters <x_data> and <y_data> must have same length!' );
assert( side ~= 0, ...
       'Parameter <side> must be ''top'' or ''bottom''' );

y_data = y_data(:);
data_len = length( y_data );
x_new = [];
y_new = [];

i = 1;
while i < data_len;
    ii = i+1:min( i + view, data_len );
    [ m, idx ] = max( ( y_data(ii) - y_data(i) ) ./ (ii-i)' .* side );

    % Equidistant x_data assumed! Use next row instead, if not:
    %[ m, idx ] = max( ( y_data(ii) - y_data(i) ) ./ ( x_data(ii) - x_data(i) ) * side );

    % New max. slope: store new "observation point"
    i = i + idx;
    x_new = [ x_new x_data(i) ];
    y_new = [ y_new y_data(i) ];
end;

env = interp1( x_new, y_new, x_data, 'linear', 'extrap' );


text:
function Y = prtRvUtilLaplaceCdf(X,mu,theta)
% Y = laplacecdf(X,mu,theta)

if numel(X) > length(X)
    error('X must have 1 singleton dimension')
end
Y = zeros(size(X));

Y(X>=mu) = 1-1/2*exp(-(X(X>=mu)-mu)./theta);
Y(X<mu) = 1/2*exp((X(X<mu)-mu)./theta);

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
classdef OptimizerToDesignVariableTranslator <  handle

    methods (Access = public, Static)

        function n = translate(optimizer)
            switch optimizer
                case {'SLERP', 'PROJECTED SLERP', 'HAMILTON-JACOBI'}
                    n = 'LevelSet';
                case {'PROJECTED GRADIENT', 'MMA', 'IPOPT'}
                    n = 'Density';
            end

        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function cadj=ref_spreadadj_4(coef)
%REF_SPREADADJ_4  Symbol of adjoint spreading function.
%   Usage: cadj=ref_spreadadj_4(c,number);
%
%   Development version by FJ for comparison of different implementations
%   cadj=SPREADADJ(c) will compute the symbol cadj of the spreading
%   operator that is the adjoint of the spreading operator with symbol c.
%
%   Improved implementation of the direct formula with higher memory
%   needs (but better speed)
%
%   This implementation uses the same improvements as case 3, but
%   also avoid the use of loop by using matrix pointwise
%   multiplication, which improves the computation time, but also
%   requires more memory due to the contruction of the (L-1)x(L-1)
%   matrix temp
%

L=size(coef,1);

cadj=zeros(L);

% Proceesing for ii==0 or jj==0
cadj(1,1)=conj(coef(1,1));
cadj(2:end,1)=conj(coef(end:-1:2,1));
cadj(1,2:end,1)=conj(coef(1,end:-1:2));

% Processing for ii~=0 and jj~=0

% Precomputation for exponential term

% Optimization note : As said in note of case 3 ,we are computing
% the Lth root of unity which have special properties and symetries
% that could be exploited to highly reduce this computation
temp2=exp((-i*2*pi/L)*(0:L-1));

% Optimization note : Here we are computing indexes for all the
% exponential terms, which leads to a highly structured matrix
% which strcture can be formalized (notably it is symetric) and
% used to reduce the computation cost
temp=mod((1:L-1)'*(1:L-1),L)+1;

% Optimization note : Finaly we construct the matrix containing all
% the needed exponential terms.
% This matrix is known as the DFT matrix and appears in the matrix
% formulation of the dft, but in our case it is used for pointwise
% multiplication instead of matrix mulitplication.
% There might be optimized algorithms for the computation of this
% matrix described in the context of fft.
temp=temp2(temp);

cadj(2:L,2:L)=conj(coef(L:-1:2,L:-1:2)).*temp;


text:
function versionPhysio = tapas_physio_version()
% returns version number of this copy of PhysIO
%
%   versionPhysio = tapas_physio_version()
%
% IN
%
% OUT
%
% EXAMPLE
%   tapas_physio_version
%
%   See also

% Author: Lars Kasper
% Created: 2018-10-19
% Copyright (C) 2018 TNU, Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS PhysIO Toolbox, which is released under
% the terms of the GNU General Public License (GPL), version 3. You can
% redistribute it and/or modify it under the terms of the GPL (either
% version 3 or, at your option, any later version). For further details,
% see the file COPYING or <http://www.gnu.org/licenses/>.
%
versionPhysio = 'R2022a-v8.1.0';

author: translationalneuromodeling

repo: tapas


text:
% AUTHOR:         Aaron Nicolson
% AFFILIATION:    Signal Processing Laboratory, Griffith University
%
% This Source Code Form is subject to the terms of the Mozilla Public
% License, v. 2.0. If a copy of the MPL was not distributed with this
% file, You can obtain one at http://mozilla.org/MPL/2.0/.

clear all; close all; clc

ver = {''};

for i = 1:length(ver)
    T = readtable(['./iter/', ver{i}, '.csv']);
    iter = 1:height(T);
    plot(iter, T.loss); hold on;
    xlabel('Iteration');
    ylabel('Loss');
end
legend(ver);
hold off;

author: anicolson

repo: DeepXi

save_path: github-repos/MATLAB/anicolson-DeepXi


text:
function[varargout]=vshift(varargin)
% VSHIFT  Cycles the elements of an array along a specified dimension.
%
%   Y=VSHIFT(X,N,DIM) cycles the elements of X N places along dimension DIM.
%
%   Example: x=[1 2 3 4 5];
%            vshift(x,+1,2)=[2 3 4 5 1]
%            vshift(x,-1,2)=[5 1 2 3 4]
%
%   Note shifting by N and then by -N recovers the original array.
%
%   [Y1,Y2,...YN]=VSHIFT(X1,X2,...XN,N,DIM) also works.
%
%   VSHIFT(X1,X2,...XN,N,DIM); with no arguments overwrite the original
%   input variables.
%
%   Note that VSHIFT is similar to Matlab's CIRCSHIFT, but has the opposite
%   convection for positive and negative shifts.
%
%   ------------------------------------------------------------------
%   Y=VSHIFT(X,N,DIM,INDEX,DIM2) applies this shift selectively, only to
%   that subset of X obtained by indexing X with INDEX along DIM2, i.e.
%
%		    1 2      DIM2     DIMS(X)
%		    | |        |         |
%		  X(:,:, ... INDEX, ..., :)
%
%   is cycled N places along dimension DIM, but the remainder of X is not.
%   DIM and DIM2 cannot be the same.  The above extensions to multiple
%   output varibles work in this case as well.
%   ------------------------------------------------------------------
%
%   See also: VINDEX, CIRCSHIFT.
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2001--2016 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--t')
  vshift_test,return
end

%/********************************************************
%Sort out input arguments
nax=2;
if nargin>4
  if  length(varargin{end-3}(:))==1
     nax=4;
     dim2=(varargin{end});
     jj=varargin{end-1}(:);
     dim=varargin{end-2};
     n=varargin{end-3};
     if dim==dim2
       error('DIM and DIM2 cannot be the same.')
     end
     %n, jj, dim,dim2
  end
end

if nax==2
  dim=varargin{end};
  n=varargin{end-1};
end
%\********************************************************

for i=1:length(varargin)-nax
  if nax==2
    varargout{i}=vshift1(varargin{i},n,dim);
  else
    varargout{i}=vshift1(varargin{i},n,dim,jj,dim2);
  end
end
eval(to_overwrite(nargin-nax))

function[y]=vshift1(x,n,ndim,jj,ndim2)

N=size(x,ndim);
if n>0
    ii=[(n+1:N) (1:n)];
elseif n<0
    n=-n;
    ii=[(N-(n-1):N) (1:N-n)];
elseif n==0
    ii=(1:N);
end

if nargin==3
%   Same as using circshift... these are the same speed
    %array=zeros(max(ndim,numel(size(x))),1);
    %array(ndim)=-n;
    %ndim
 %   y=circshift(x,-n,ndim);
    y=vindex(x,ii,ndim);
else
    y1=vindex(x,jj,ndim2);
    y1=vindex(y1,ii,ndim);
    %vsize(x,y1,jj,ndim2);
    %x,y1,ii,jj,ndim2
    y=vindexinto(x,y1,jj,ndim2);
end

function[]=vshift_test
x=(1:10);
ans1=[(2:10) 1];
reporttest('VSHIFT col case', aresame(vshift(x,1,2),ans1))

x=(1:10)';
ans1=[10 (1:9)]';
reporttest('VSHIFT row case', aresame(vshift(x,-1,1),ans1))

clear x ans1
x(:,:,1)=[1 2; 3 4];
x(:,:,2)=2*[1 2; 3 4];
ans1(:,:,2)=x(:,:,1);
ans1(:,:,1)=x(:,:,2);
reporttest('VSHIFT mat case', aresame(vshift(x,1,3),ans1))

clear x ans1
x(:,:,1)=[1 2; 3 4];
x(:,:,2)=2*[1 2; 3 4];
ans1(:,:,1)=[3 2;1 4];
ans1(:,:,2)=2*[3 2;1 4];
reporttest('VSHIFT mat selective case one', aresame(vshift(x,1,1,1,2),ans1))

clear x ans1
x(:,:,1)=[1 2; 3 4];
x(:,:,2)=2*[1 2; 3 4];
ans1(:,:,1)=[3 4;1 2];
ans1(:,:,2)=2*[3 4;1 2];
reporttest('VSHIFT mat selective case two', aresame(vshift(x,1,1,1:2,2),ans1))

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function n=grid2n(i,j,Nj)
n=i+(j-1)*Nj;


text:
classdef surfPrinter < figurePrinter

     properties (Access = private)

    end

    methods (Access = public)

        function obj = surfPrinter(figureID)
            obj.init(figureID);
            obj.setFontSizes();
            obj.setPositionAndSize();
        end

    end

    methods (Access = private)

        function init(obj,figureID)
            obj.figID = figureID;
        end

    end

end


text:
function C = long(d)
%LONG         Long class constructor
%
%  C = long(d)
%
%  C is long scalar or column vector, array input d forced to be column vector
%
%Long representation
%
%  C = C.sign * sum( C.mantissa(i)*beta^(-i) ) * beta^C.exponent +/- C.error
%
%where C.error = C.error.mant * beta^C.error.exp and with
%  summation from 1 to precision (=size(C.mantissa,2)) with
%
%  1 <= precision <= INTLAB_LONG_PRECISION
%
%and integers
%
%  C.sign       in {-1,1}
%  C.mantissa   in 0 .. beta-1
%  C.exponent   representable integer ( -2^52+1 .. 2^52-1 )
%  C.error      nonnegative double, stored by C.error.mant and C.error.exp
%
%Computations can be executed with or w/o error term, see longinit. For
%  computational speed, comparison, min/max and absolute value refer to
%  midpoint.
%  To compare, for example, intervals (for computation with error term)
%  use inf(A)>sup(B) instead of A>B, and so forth.
%
%For control of working precision, see help longprecision. Base beta is a
%  power of 2 so that double precision floating point numbers are stored
%  without error.
%An example of long arithmetic with big cancellation is
%
%  x = long(-20);
%  y = long(1);  t = long(1);  i = 0;
%  while abs(t)>1e-20
%    i = i+1;
%    t = t*x/i;
%    y = y+t;
%  end
%  format long
%  Y = long2intval(y)
%
%producing
%
%  intval Y =
%    1.0e-008 *
%     0.20612_________
%
%The poor accuracy improves with more precision. After
%  longprecision(40);
%and the same statements as above we obtain
%
%  intval Y =
%    1.0e-008 *
%     0.20611536224378
%
%For more information try demolong .
%

% written  12/30/98     S.M. Rump
% modified 02/09/01     S.M. Rump  performance improvement
% modified 09/29/02     S.M. Rump  care for NaN components
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 08/26/12     S.M. Rump  global variables removed
%

  INTLAB_LONG_LOGBETA = getappdata(0,'INTLAB_LONG_LOGBETA');
  INTLAB_LONG_BETA = getappdata(0,'INTLAB_LONG_BETA');
  INTLAB_LONG_ERROR = getappdata(0,'INTLAB_LONG_ERROR');

  if nargin==0
    C.sign = [];
    C.exponent = [];
    C.mantissa = [];
    C.error.mant = [];
    C.error.exp = [];
    C = class(C,'long');
    return
  end

  if isa(d,'long')
    C = d;
  else
    if isempty(d)
      C.sign = [];
      C.exponent = [];
      C.mantissa = [];
      C.error.mant = [];
      C.error.exp = [];
      C = class(C,'long');
      return
    end
    sized = size(d);
    n = prod(sized);
    if n~=sized(1)
      warning('input array for long forced to be column vector')
      d = d(:);
    end
    indexNaN = isnan(d);
    d(indexNaN) = [];
    [s e m] = splitdble(d);

    % get sign
    C.sign = s;

    % get exponent
    q = ceil(e/INTLAB_LONG_LOGBETA);
    C.exponent = q;

    % get mantissa digits
    C.mantissa = zeros(size(d,1),ceil(53/INTLAB_LONG_LOGBETA)+1);
    p = 0;
    while any(m)
      p = p+1;
      m = m*INTLAB_LONG_BETA;
      C.mantissa(:,p) = floor(m);
      m = m - C.mantissa(:,p);
    end

    % treat zero components
    index = ( d==0 );
    if any(index)
      C.sign(index) = 1;
      C.exponent(index) = -inf;
    end

    % adjust mantissa digits to exponent
    r = INTLAB_LONG_LOGBETA*q - e;
    index = ( r~=0 );
    if any(index)
      C.mantissa(index,:) = shiftright(C.mantissa(index,:),r(index));
    end

    % omit trailing zeros (improves performance)
    [m index] = max(C.mantissa(:,end:-1:1)~=0,[],2);
    index(all(C.mantissa'==0)) = size(C.mantissa,2);
    if min(index)~=1
      C.mantissa = C.mantissa(:,1:end-min(index)+1);
    end

    if any(indexNaN)
      Csign = C.sign;
      Cexponent = C.exponent;
      Cmantissa = C.mantissa;
      C.sign = zeros(n,1);
      C.exponent = zeros(n,1);
      C.mantissa = zeros(n,size(Cmantissa,2));
      C.sign(indexNaN) = NaN;
      C.exponent(indexNaN) = NaN;
      C.mantissa(indexNaN) = NaN;
      C.sign(~indexNaN) = Csign;
      C.exponent(~indexNaN) = Cexponent;
      C.mantissa(~indexNaN) = Cmantissa;
    end

    % set error
    C.error.mant = zeros(n,1);
    C.error.exp = zeros(n,1);

    C = class(C,'long');

  end

repo: OpenMAS


text:
function utmp = IEC_MANN(Nx,Ny,Nz,delta,UHub,HubHt,IECturbC)
%%INITIALISATION

% definition of constants
twopi=2*pi;
%constants and derived parameters from IEC
gamma = 3.9; %IEC, (B.12)
alpha = 0.2; %IEC, sect. 6.3.1.2

%set delta1 according to guidelines (chap.6)
if HubHt<=60,
    delta1=0.7*HubHt;
else
    delta1=42;
end;

%IEC, Table 1, p.22
if IECturbC == 'A',
    Iref=0.16;
elseif IECturbC == 'B',
    Iref=0.14;
elseif IECturbC == 'C',
    Iref=0.12;
else
    error('IECturbC can be equal to A,B or C;adjust the input value')
end;

%IEC, sect. 6.3.1.3
b=5.6;
sigma1=Iref*(0.75*UHub+b);
sigma2=0.7*sigma1;
sigma3=0.5*sigma1;
%derived constants
l=0.8*delta1; %IEC, (B.12)
sigmaiso=0.55*sigma1; %IEC, (B.12)

Cij2=zeros(3,3,Nx,Ny,Nz);
ikx=cat(2,-Nx/2:-1,1:Nx/2-1);
[x y z]=ndgrid(ikx,1:Ny,1:Nz);
k1=(x)*l/(Nx*delta)*twopi;
k2=(y-Ny/2)*l/(Ny*delta)*twopi;
k3=(z-Nz/2)*l/(Nz*delta)*twopi;
kabs=sqrt(k1.^2+k2.^2+k3.^2);
beta= gamma./(kabs.^(2/3).*real(pfq([1/3 17/6],4/3,-kabs.^(-2))));
k03=k3+beta.*k1;
k0abs=sqrt(k1.^2+k2.^2+k03.^2);
Ek0=1.453*k0abs.^4./(1+k0abs.^2).^(17/6);
C1=beta.*k1.^2.*( k0abs.^2 - 2*k03.^2 + beta.*k1.*k03 )./( kabs.^2.*( k1.^2 + k2.^2 ));
C2=k2.*k0abs.^2./ (exp( (3/2).*log( k1.^2 + k2.^2 ) )) .* atan2( beta.*k1.* sqrt( k1.^2 + k2.^2 ) ,( k0abs.^2 - k03.*k1.*beta));
xhsi1=C1 - k2.*C2./k1;
xhsi2=k2.*C1./k1 + C2;
CC=sigmaiso*sqrt(twopi*pi*l^3.*Ek0./(Nx*Ny*Nz*delta^3.*k0abs.^4));
Cij2(1,1,1:size(k1,1),:,:)= CC.*( k2.*xhsi1);
Cij2(1,2,1:size(k1,1),:,:)= CC.*( k3 - k1.*xhsi1 + beta.*k1);
Cij2(1,3,1:size(k1,1),:,:)= CC.*( -k2);
Cij2(2,1,1:size(k1,1),:,:)= CC.*( k2.*xhsi2 - k3 - beta.*k1);
Cij2(2,2,1:size(k1,1),:,:)= CC.*( -k1.*xhsi2);
Cij2(2,3,1:size(k1,1),:,:)= CC.*( k1);
Cij2(3,1,1:size(k1,1),:,:)= CC.*( k0abs.^2.*k2 ./ (kabs.^2));
Cij2(3,2,1:size(k1,1),:,:)= CC.*( -k0abs.^2.*k1 ./ (kabs.^2));
%Cij2(isnan(Cij2))=0;
%% FOURIER COEFFICIENTS
n=complex(normrnd(0,1,[3 1 Nx Ny Nz]),normrnd(0,1,[3 1 Nx Ny Nz]));
H = reshape(mtimesx(Cij2,n,'speed'),[],Nx,Ny,Nz);
%% TIME SERIES
u(1,:,:,:)=real(fftn(fftshift(H(1,:,:,:))));
u(2,:,:,:)=real(fftn(fftshift(H(2,:,:,:))));
u(3,:,:,:)=real(fftn(fftshift(H(3,:,:,:))));

stdDv=zeros(Ny*Nz,3);
stdDv(:,1)=reshape(std(u(1,:,:,:)),Ny*Nz,1); % Standard deviation for u-component
stdDv(:,2)=reshape(std(u(2,:,:,:)),Ny*Nz,1); % Standard deviation for v-component
stdDv(:,3)=reshape(std(u(3,:,:,:)),Ny*Nz,1); % Standard deviation for w-component

% Introduction of a Corrective Factor SF in order to adjust the statistics
% as advised by the OC4 project DTUs report

sf=zeros(Ny*Nz,3);
sf(:,1)=sqrt(sigma1^2./(stdDv(:,1).^2));
sf(:,2)=sqrt((sigma2)^2./(stdDv(:,2).^2));
sf(:,3)=sqrt((sigma3)^2./(stdDv(:,3).^2));

utmp=zeros(3,Nx,Ny,Nz);
for i=1:length(sf),
    for j=1:Ny,
        for k=1:Nz,
            utmp(1,:,j,k)=sf(i,1).*u(1,:,j,k);
            utmp(2,:,j,k)=sf(i,2).*u(2,:,j,k);
            utmp(3,:,j,k)=sf(i,3).*u(3,:,j,k);
        end;
    end;
end;
for i=1:Nz,
    thisZ=HubHt-(Nz*delta)/2+1.*i*(Nz*delta)/(Nz-1);
    profile=UHub*(thisZ/HubHt)^alpha;
    utmp(1,:,:,i)=(utmp(1,:,:,i)+profile);
end;

repo: mcbench-benchmarks


text:
%compute speed in the direction perpendicular to tail-central direction
function [data,units]=compute_velpertc(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
velpertc=cell(1,numlarvae);
%tailcentralangperp=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    %tailcentralangperp{1,i}=trx(larva).tailcentralang-pi/2;
    velpertc{1,i}=trx(larva).velmag_ctr.*(cos(trx(larva).velang).*cos(trx(larva).tailcentralang(1,1:end-1)+pi/2)+sin(trx(larva).velang).*sin(trx(larva).tailcentralang(1,1:end-1)+pi/2));
end

units=parseunits('mm/s');
data=velpertc;

author: kristinbranson

repo: JAABA


text:
function [params, names] = vargpTimeDynamicsExtractParam(model)

% VARGPTIMEDYNAMICSEXTRACTPARAM Extract parameters from the GP time dynamics model.
%
%	Description:
%
%	PARAMS = VARGPTIMEDYNAMICSEXTRACTPARAM(MODEL) extracts the model
%	parameters from a structure containing the information about a
%	Gaussian process dynamics model.
%	 Returns:
%	  PARAMS - a vector of parameters from the model.
%	 Arguments:
%	  MODEL - the model structure containing the information about the
%	   model.
%	DESC does the same as above, but also returns parameter names.
%	ARG model : the model structure containing the information about
%	the model.
%	RETURN params : a vector of parameters from the model.
%	RETURN names : cell array of parameter names.
%
%
%
%	See also
%	GPEXTRACTPARAM, GPTIMEDYNAMICSCREATE, GPTIMEDYNAMICSEXPANDPARAM, MODELEXTRACTPARAM

%	Based on GPTIMEDYNAMICSEXTRACTPARAM
% VARGPLVM

if nargout > 1
  returnNames = true; % Also return parameter names in an array
else
  returnNames = false;
end

% Parameters are extracted in the following order: (notation: % parameter{size})
% [dynamicsVardistParams{dynamics.vardist.nParams}        % mu_bar, lambda
%  dynamics.kernParams{dynamics.kern.nParams}]            % sigmaRBF, lt,
%  sigmaWhite

% Variational parameters (reparametrized means and covariances)
if returnNames
  [varParams, varNames] = modelExtractParam(model.vardist);
  names = varNames;
else
  %varParams = vardistExtractParam(model.vardist);
  varParams = modelExtractParam(model.vardist);
end
params = varParams;

% Kernel parameters
if returnNames
  [kernParams, kernParamNames] = kernExtractParam(model.kern);
  for i = 1:length(kernParamNames)
    kernParamNames{i} = ['Kernel, ' kernParamNames{i}];
  end
  names = {names{:}, kernParamNames{:}};
else
  kernParams = kernExtractParam(model.kern);
end

params = [params kernParams];

author: SheffieldML

repo: vargplvm


text:
% Chapter 4 - Electromagnetic Waves and Optical Resonators.
% Program_4b - Bifurcation Diagram for a Nonlinear Optical Resonator.
% Copyright Birkhauser 2013. Stephen Lynch.

% Bifurcation diagram for a simple fiber resonator (Figures 4.13 & 4.16(a)).
clear
halfN=1999;N=2*halfN+1;N1=1+halfN;
format long;E1=zeros(1,N);E2=zeros(1,N);
Esqr=zeros(1,N);Esqr1=zeros(1,N);ptsup=zeros(1,N);
C=0.345913;
E1(1)=0;kappa=0.0225;Pmax=120;phi=0;

% Ramp the power up
for n=1:halfN
    E2(n+1)=E1(n)*exp(1i*(abs(C*E1(n))^2-phi));
    E1(n+1)=1i*sqrt(1-kappa)*sqrt(n*Pmax/N1)+sqrt(kappa)*E2(n+1);
    Esqr(n+1)=abs(E1(n+1))^2;
end

% Ramp the power down
for n=N1:N
    E2(n+1)=E1(n)*exp(1i*(abs(C*E1(n))^2-phi));
    E1(n+1)=1i*sqrt(1-kappa)*sqrt(2*Pmax-n*Pmax/N1)+sqrt(kappa)*E2(n+1);
    Esqr(n+1)=abs(E1(n+1))^2;
end

for n=1:halfN
    Esqr1(n)=Esqr(N+1-n);
    ptsup(n)=n*Pmax/N1;
end

% Plot the bifurcation diagrams
fsize=15;
subplot(2,1,1)
plot(Esqr(1:N),'.','MarkerSize',1)
xlabel('Number of Ring Passes','FontSize',fsize);
ylabel('Output','FontSize',fsize);

subplot(2,1,2)
hold on
plot(ptsup(1:halfN),Esqr(1:halfN),'.','MarkerSize',1);
plot(ptsup(1:halfN),Esqr1(1:halfN),'.','MarkerSize',1);
xlabel('Input Power','FontSize',fsize);
ylabel('Output Power','FontSize',fsize);
hold off

% End of Program_4b.

author: Sable


text:
classdef HolonomicConstraint < handle
    % HolonomicConstraint: represents a group of holonomic constraints of a
    % dynamical systems.
    %
    % A holonomic constraint is defined as a function of state variables
    % that equals to constants, i.e.,
    %   h(x) == hd
    %
    % @author ayonga @date 2017-04-20
    %
    % Copyright (c) 2016-2017, AMBER Lab
    % All right reserved.
    %
    % Redistribution and use in source and binary forms, with or without
    % modification, are permitted only in compliance with the BSD 3-Clause
    % license, see
    % http://www.opensource.org/licenses/bsd-license.php

    % properties determined internally
    properties (SetAccess=protected, GetAccess=public)
       % The dimension of the virtual constraints
        %
        % @type integer
        Dimension

        % The symbolic representation of constant parameter value 'hd' that
        % the holonomic constraints associated wtih.
        %
        % @type SymVariable
        Param

        % The symbolic representation of the input variables associated
        % with the holonomic constraints
        %
        % @type SymVariable
        Input
    end

    % properties must be determined by the users
    properties (SetAccess=protected, GetAccess=public)
        % The name of the virtual constraints
        %
        % @type char
        Name

        % The label of the holonomic constraint
        %
        % @type char
        ConstrLabel

        % The highest order of derivatives to enforce the holonomic
        % constraints. This is the same as the relative degree of a virtual
        % constraint.
        %
        % @type integer
        DerivativeOrder

    end

    properties (Dependent)
        % The holonomic constraint expression
        %
        % @type SymFunction
        ConstrExpr

        % The Jacobian matrix of the holonomic constraitns
        %
        % @type SymFunction
        ConstrJac

         % The Jacobian matrix of the holonomic constraitns
        %
        % @type SymFunction
        ConstrJacDot

        % The name of the parameter variable 'hd' associated with the
        % holonomic constraints
        %
        % @type char
        ParamName

        % The name othe input variables (external forces, etc.) associated
        % with the holonomic constraints
        %
        % @type char
        InputName
    end

    %% GET methods
    methods
        function cstr = get.ConstrExpr(obj)
            cstr = obj.h_;
        end
        function jac = get.ConstrJac(obj)
            jac = obj.Jh_;
        end
        function jacdot = get.ConstrJacDot(obj)
            jacdot = obj.dJh_;
        end
        function name = get.ParamName(obj)
            name = ['p' obj.Name];
        end
        function name = get.InputName(obj)
            name = ['f' obj.Name];
        end
        function name = get.h_name(obj)
            name = ['h_' obj.Name '_' obj.Model.Name];
        end

        function name = get.Jh_name(obj)
            name = ['Jh_' obj.Name '_' obj.Model.Name];
        end

        function name = get.dJh_name(obj)
            name = ['dJh_' obj.Name '_' obj.Model.Name];
        end

        function name = get.dh_name(obj)
            name = ['dh_' obj.Name '_' obj.Model.Name];
        end

        function name = get.ddh_name(obj)
            name = ['ddh_' obj.Name '_' obj.Model.Name];
        end
    end

    properties (Dependent)

        h_name

        Jh_name

        dh_name

        dJh_name

        ddh_name
    end

    properties (Access = protected)
        % The dynamical system model
        %
        % @type DynamicalSystem
        Model

        % The holonomic constraint expression
        %
        % @type SymFunction
        h_

        % The Jacobian matrix of the holonomic constraint expression
        %
        % @type SymFunction
        Jh_

        % The first order derivative of the holonomic constraint expression
        %
        % @type SymFunction
        dh_
        % The first order derivative of the Jacobian matrix of the
        % holonomic constraint expression
        %
        %
        % @type SymFunction
        dJh_

        % The second order derivatives of the holonomic constraints
        % expression
        %
        % @type SymFunction
        ddh_
    end

    methods

        function obj = HolonomicConstraint(model, h, name, varargin)
            % The class constructor function
            %
            % Parameters:
            % model: the dynamical system model in which the virtual
            % constraints are defined @type DynamicalSystem
            % h: the symbolic expression of the constraints
            % name: the name of the virtual constraints @type char
            % @type SymExpression
            % varargin: optional parameters. In details
            %  ConstrLabel: labels for constraints @type cellstr
            %  Jacobian: the custom Jacobian matrix @type SymExpression
            %  DerivativeOrder: the degree of holonomic constraints
            %  @type integer

            if nargin == 0
                return;
            end

            % validate (model) argument
            validateattributes(model, {'ContinuousDynamics'},...
                {'scalar'},...
                'HolonomicConstraint','model');
            obj.Model = model;
            x = model.States.x;

            % validates (name) argument
            validateName(obj, name);
            obj.Name = name;

            % parse the input options
            args = struct(varargin{:});
            assert(isscalar(args),...
                'The values of optional properties are must be scalar data.');

            if ~isempty(h)
                % validate (ya) argument
                validateattributes(h,{'SymExpression'},...
                    {'nonempty','vector'},...
                    'HolonomicConstraint','h');
                if isrow(h) % convert to column vector if it is a row vector
                    h = vertcat(h(:));
                end
            elseif isfield(args, 'LoadPath') && ~isempty(args.LoadPath)
                h = SymExpression([]);
                h = load(h, args.LoadPath, obj.h_name);
            else
                error(['Unable to create the HolonomicConstraint object. ',...
                    'Either the expression is empty or the load path is not specified.'],...
                    'HolonomicConstraint');
            end

            if isfield(args, 'LoadPath') && ~isempty(args.LoadPath)
                is_loaded = true;
            else
                is_loaded = false;
            end

            dim = length(h);
            obj.Dimension = dim;

            hd = SymVariable(obj.ParamName, [dim,1]);
            obj.Param = hd;
            obj.Input = SymVariable(obj.InputName,[dim,1]);

            if is_loaded % the loaded expression is h:= h_a - h_d
                obj.h_ = SymFunction(obj.h_name, h, {x, hd});
            else
                obj.h_ = SymFunction(obj.h_name, h-hd, {x, hd});
            end

            % validate and assign the desired outputs
            if isfield(args, 'ConstrLabel')
                obj.setConstrLabel(args.ConstrLabel);
            end

            if is_loaded
                Jh = SymFunction(obj.Jh_name, [], {model.States.x});
                obj.Jh_ = load(Jh,args.LoadPath);
            else
                if isfield(args, 'Jacobian')
                    obj.setJacobian(args.Jacobian);
                end
            end

            if isfield(args, 'DerivativeOrder')
                obj.setDerivativeOrder(args.DerivativeOrder);
            end

            if is_loaded
                obj.configure(args.LoadPath);
            else
                obj.configure();
            end

        end

    end

    methods

        % configure and compile the holonomic constraints
        obj = configure(obj, load_path);

        % enforce as NLP constraints
        nlp = imposeNLPConstraint(obj, nlp);

        % calculate the Jacobian matrix of the holonomic constraints
        [Jh, dJh] = calcJacobian(obj, x, dx);

        h = calcConstraint(obj, x);

        % save the symbolic expressions
        saveExpression(obj, export_path, varargin);

    end

    % set functions
    methods

        function export(obj, export_path, varargin)
            % export the symbolic expressions of the constraints matrices and
            % vectors and compile as MEX files.
            %
            % Parameters:
            %  export_path: the path to export the file @type char
            %  varargin: variable input parameters @type varargin
            %   Vars: a list of symbolic variables @type SymVariable
            %   File: the (full) file name of exported file @type char
            %   ForceExport: force the export @type logical
            %   BuildMex: flag whether to MEX the exported file @type logical
            %   Namespace: the namespace of the function @type char

            export(obj.h_,export_path, varargin{:});
            export(obj.Jh_,export_path, varargin{:});

            if ~isempty(obj.dJh_)
                export(obj.dJh_,export_path, varargin{:});
            end

        end

        % Jacobian matrix
        function obj = setJacobian(obj, jac)
            % sets the Jacobian matrix of the holonomic constraints if it
            % is provided directly by the users
            %
            % Parameters:
            % jac: the jacobian matrix @type SymExpression

            model = obj.Model;
            if nargin > 1
                validateattributes(jac,{'SymExpression'},...
                    {'2d','size',[obj.Dimension,model.numState]},...
                    'HolonomicConstraint','Jacobian');
            else
                jac = jacobian(obj.h_, model.States.x);
            end
            obj.Jh_ = SymFunction(obj.Jh_name, jac, {model.States.x});
        end

        % RelativeDegree
        function obj = setDerivativeOrder(obj, degree)
            % sets the highest derivative order of holonomic constraints in
            % order to enforce as bilateral constraints of the system
            %
            % Parameters:
            % degree: derivative order @type integer

            if degree == 1
                error('Currently we do not support J(x)dx = 0 type of holonomic constraitns.')
            end

            validateattributes(degree, {'double'},...
                {'nonempty','scalar','positive','integer','>=',1,'<=',2},...
                'HolonomicConstraint','DerivativeOrder');
            obj.DerivativeOrder = degree;

        end

        % OutputLabel
        function obj = setConstrLabel(obj, label)
            % sets the naming labels of outputs
            %
            % Parameters:
            % label: the cell array of labels @type cellstr

            validateattributes(label,{'cell'},...
                {'nonempty','numel',obj.Dimension,'row'},...
                'HolonomicConstraint','ConstrLabel');
            cellfun( @(x) validateattributes(...
                x, {'char'},{}), label);

            obj.ConstrLabel = label;
        end

        % Name
        function name = validateName(~, name)
            validateattributes(name, {'char'},...
                {'nonempty','scalartext'},...
                'HolonomicConstraint','Name');

            assert(isempty(regexp(name, '\W', 'once')) || ~isempty(regexp(name, '\$', 'once')),...
                'HolonomicConstraint:invalidSymbol', ...
                'Invalid symbol string, can NOT contain special characters.');

            assert(isempty(regexp(name, '_', 'once')),...
                'HolonomicConstraint:invalidSymbol', ...
                'Invalid symbol string, can NOT contain ''_''.');

        end
    end
end

author: ayonga

repo: frost-dev

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


text:
%--- help for mdd/mhm ---
%
%mdd/mhm is a function.
%    log_mhm = mhm(obj, mhm_tau, opts)
%

author: jmaih


text:
function [gal] = in32gal(in3)
% Convert volume from cubic inches to US liquid gallons.
% Chad Greene 2012
gal = in3*0.004329004329 ;


text:
function C = blkdiageye(X,k)
% construct block-diagonal matrix with k copies of X on diagonal
% Equivalent to C = kron(eye(k),X) but much faster
%
% Author: Tim Mullen, 2011 (C) SCCN/INC/UCSD

ss = repmat({X},1,k);
C = blkdiag(ss{:});

% alternate form
% ss = repmat('X,',1,k); ss(end)=[];
% eval(sprintf('C = blkdiag(%s);',ss));

repo: REST