text:
function unit = seconds2unit(seconds,seconds0,seconds1)
% Convert seconds past midnight to unit times
%
% USAGE:
%   [UNIT] = seconds2unit(SECONDS,SECONDS0,SECONDS1)
%
% INPUTS:
%   SECONDS - m by 1 column vector of times expressed as seconds past midnight
%               (e.g. 1:00:00 is 3600, 12:00:15 is 43215) where m>=2
%   SECONDS0   - Base time, maps to 0 un the unit interval
%   SECONDS1   - End time, maps to 1 un the unit interval
%
% OUTPUTS:
%   UNIT    - m by 1 column vector of times measured as fraction of time
%               between min(SECONDS) and max(SECONDS)
%
% COMMENTS:
%   This is a helper function for REALIZED_KERNEL
%
%  See also WALL2UNIT, WALL2SECONDS, SECONDS2WALL, REALIZED_KERNEL, REALIZED_VARIANCE

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 5/1/2008

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin~=3
    error('Three inputs required.')
end
if any(seconds<0)
    error('SECONDS must be non-negative');
end
if seconds0<0
    error('SECONDS0 must be non-negative');
end
if seconds1<0
    error('SECONDS1 must be non-negative');
end

if (seconds1-seconds0)<eps
    error('SECONDS1 must be larger than SECONDS0.');
end

if size(seconds,2)>size(seconds,1)
    seconds=seconds';
end

seconds = double(seconds);
seconds0 = double(seconds0);
seconds1 = double(seconds1);

if size(seconds,2)>1
    error('SECONDS must be an m by 1 column vector');
end
if length(seconds)<2
    error('SECONDS must contain at least 2 elements');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
unit=(seconds-seconds0)/(seconds1 - seconds0);

author: bashtage

repo: mfe-toolbox


text:
function output = callecos(yalmipmodel)
originalModel = yalmipmodel;
% Fix ECOS options
options = yalmipmodel.options;
options.ecos.verbose = options.verbose~=0;
if ~isempty(yalmipmodel.binary_variables)
    options.ecos.bool_vars_idx = yalmipmodel.binary_variables;
end
if ~isempty(yalmipmodel.integer_variables)
    options.ecos.int_vars_idx = yalmipmodel.integer_variables;
end
model.opts = options.ecos;

if ~isempty(yalmipmodel.ub)
    [yalmipmodel.F_struc,yalmipmodel.K] = addStructureBounds(yalmipmodel.F_struc,yalmipmodel.K,yalmipmodel.ub,yalmipmodel.lb);
end

% Write nonlinear functions using exponential cone operators, if possible
[yalmipmodel,output] = normalizeExponentialCone(yalmipmodel);
if output.problem
    return
end

% Map to ECOS syntax
cones = [];
cones.f = yalmipmodel.K.f;
cones.l = yalmipmodel.K.l;
cones.q = yalmipmodel.K.q;
cones.s = yalmipmodel.K.s;
cones.ep = yalmipmodel.K.e;
model.c = full(yalmipmodel.c);
model.A = -yalmipmodel.F_struc(1:cones.f,2:end);if isempty(model.A);model.A = [];end
model.b = full(yalmipmodel.F_struc(1:cones.f,1));if isempty(model.b);model.b = [];end
model.G = -yalmipmodel.F_struc(1+cones.f:end,2:end);if isempty(model.G);model.G = [];end
model.h = full(yalmipmodel.F_struc(1+cones.f:end,1));if isempty(model.h);model.h = [];end
if nnz(cones.l) > 0
    model.dims.l = cones.l;
end
if nnz(cones.q) > 0
    model.dims.q = cones.q;
end
if nnz(cones.ep) > 0
    model.dims.e = cones.ep;
    tempG = model.G(1:model.dims.l+sum(cones.q),:);
    temph = model.h(1:model.dims.l+sum(cones.q));
    top = 1+model.dims.l+sum(cones.q);
    for i = 1:model.dims.e
        tempG = [tempG;model.G(top + [0;2;1],:)];
        temph = [temph;model.h(top + [0;2;1],:)];
        top = top  + 3;
    end
    model.G = tempG;
    model.h = temph;
end

if options.savedebug
    save ecosdebug model
end

if options.showprogress;showprogress(['Calling ' yalmipmodel.solver.tag],options.showprogress);end
if isempty(model.A)
    solvertime = tic;
    [x,y,info,s,z] = ecos(model.c,model.G,model.h,model.dims,model.opts);
    solvertime = toc(solvertime);
else
    solvertime = tic;
    [x,y,info,s,z] = ecos(model.c,model.G,model.h,model.dims,model.A,model.b,model.opts);
    solvertime = toc(solvertime);
end

% Internal format, only keep original variablesop
if ~isempty(x)
    Primal = x(1:length(originalModel.c));
else
    Primal = nan(length(originalModel.c),1);
end
if ~isempty(yalmipmodel.evalMap)
    % No support for duals when exponential cones yet
    Dual = [];
else
    Dual   = [y;z];
end

switch info.exitflag
    case 0
        problem = 0;
    case 1
        problem = 1;
    case 2
        problem = 2;
    case -1
        problem = 3;
    case {-2,-3}
        problem = 4;
    case -7
        problem = 9;
    case {10,11}
        problem = 4;
    otherwise
        problem = 9;
end

% Save ALL data sent to solver
if options.savesolverinput
    solverinput.model = model;
else
    solverinput = [];
end

% Save ALL data from the solution?
if options.savesolveroutput
    solveroutput.x = x;
    solveroutput.y = y;
    solveroutput.info = info;
    solveroutput.s = s;
    solveroutput.z = z;
else
    solveroutput = [];
end

% Standard interface
output = createOutputStructure(Primal,Dual,[],problem,yalmipmodel.solver.tag,solverinput,solveroutput,solvertime);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function results = spm_cfg_results
% SPM Configuration file for Results Report
%__________________________________________________________________________
% Copyright (C) 2005-2018 Wellcome Trust Centre for Neuroimaging

% $Id: spm_cfg_results.m 7268 2018-02-27 16:48:13Z guillaume $

%--------------------------------------------------------------------------
% spmmat Select SPM.mat
%--------------------------------------------------------------------------
spmmat         = cfg_files;
spmmat.tag     = 'spmmat';
spmmat.name    = 'Select SPM.mat';
spmmat.help    = {'Select the SPM.mat file that contains the design specification.'};
spmmat.filter  = 'mat';
spmmat.ufilter = '^SPM\.mat$';
spmmat.num     = [1 1];

%--------------------------------------------------------------------------
% titlestr Results Title
%--------------------------------------------------------------------------
titlestr         = cfg_entry;
titlestr.tag     = 'titlestr';
titlestr.name    = 'Results Title';
titlestr.help    = {'Heading on results page - determined automatically if left empty'};
titlestr.val     = {''};
titlestr.strtype = 's';
titlestr.num     = [0 Inf];
titlestr.hidden  = true;

%--------------------------------------------------------------------------
% contrasts Contrast(s)
%--------------------------------------------------------------------------
contrasts         = cfg_entry;
contrasts.tag     = 'contrasts';
contrasts.name    = 'Contrast(s)';
contrasts.help    = {
                     'Index of contrast(s). If more than one number is entered, analyse a conjunction hypothesis.'
                     ''
                     'If only one number is entered, and this number is "Inf", then results are printed for all contrasts found in the SPM.mat file.'
}';
contrasts.strtype = 'n';
contrasts.num     = [1 Inf];

%--------------------------------------------------------------------------
% threshdesc Threshold type
%--------------------------------------------------------------------------
threshdesc        = cfg_menu;
threshdesc.tag    = 'threshdesc';
threshdesc.name   = 'Threshold type';
threshdesc.help   = {''};
threshdesc.labels = {'FWE' 'none' 'FDR'};
threshdesc.values = {'FWE' 'none' 'FDR'};
threshdesc.val    = {'FWE'};

%--------------------------------------------------------------------------
% thresh Threshold
%--------------------------------------------------------------------------
thresh         = cfg_entry;
thresh.tag     = 'thresh';
thresh.name    = 'Threshold';
thresh.help    = {''};
thresh.strtype = 'r';
thresh.num     = [1 1];
thresh.val     = {0.05};

%--------------------------------------------------------------------------
% extent Extent (voxels)
%--------------------------------------------------------------------------
extent         = cfg_entry;
extent.tag     = 'extent';
extent.name    = 'Extent (voxels)';
extent.help    = {''};
extent.strtype = 'w';
extent.num     = [1 1];
extent.val     = {0};

%--------------------------------------------------------------------------
% conjunction Conjunction Number
%--------------------------------------------------------------------------
conjunction         = cfg_entry;
conjunction.tag     = 'conjunction';
conjunction.name    = 'Conjunction number';
conjunction.help    = {'Conjunction number. Unused if a simple contrast is entered.'
    'For Conjunction Null, enter 1.'
    'For Global Null, enter the number of selected contrasts.'
    'For Intermediate, enter the number of selected contrasts minus the number of effects under the Null.'}';
conjunction.strtype = 'n';
conjunction.num     = [1 1];
conjunction.val     = {1};

%--------------------------------------------------------------------------
% contrasts Contrast(s)
%--------------------------------------------------------------------------
contrasts1         = cfg_entry;
contrasts1.tag     = 'contrasts';
contrasts1.name    = 'Contrast(s)';
contrasts1.help    = {'Index of contrast(s) for masking.'};
contrasts1.strtype = 'n';
contrasts1.num     = [1 Inf];

%--------------------------------------------------------------------------
% thresh Mask threshold
%--------------------------------------------------------------------------
thresh1         = cfg_entry;
thresh1.tag     = 'thresh';
thresh1.name    = 'Mask threshold';
thresh1.help    = {''};
thresh1.strtype = 'r';
thresh1.num     = [1 1];
thresh1.val     = {0.05};

%--------------------------------------------------------------------------
% mtype Nature of mask
%--------------------------------------------------------------------------
mtype        = cfg_menu;
mtype.tag    = 'mtype';
mtype.name   = 'Nature of mask';
mtype.help   = {''};
mtype.labels = {'Inclusive' 'Exclusive'};
mtype.values = {0 1};

%--------------------------------------------------------------------------
% mask Mask using contrast
%--------------------------------------------------------------------------
contrast      = cfg_branch;
contrast.tag  = 'contrast';
contrast.name = 'Contrast';
contrast.val  = {contrasts1 thresh1 mtype};
contrast.help = {'Masking using contrast.'};

%--------------------------------------------------------------------------
% name Mask image
%--------------------------------------------------------------------------
name         = cfg_files;
name.tag     = 'name';
name.name    = 'Mask image(s)';
name.help    = {''};
name.filter  = {'image'};
name.ufilter = '.*';
name.num     = [1 Inf];

%--------------------------------------------------------------------------
% image Mask using image
%--------------------------------------------------------------------------
image      = cfg_branch;
image.tag  = 'image';
image.name = 'Image';
image.val  = {name mtype};
image.help = {'Masking using image(s).'};

%--------------------------------------------------------------------------
% none No Masking
%--------------------------------------------------------------------------
none      = cfg_const;
none.tag  = 'none';
none.name = 'None';
none.val  = { 1 };
none.help = {'No masking.'};

%--------------------------------------------------------------------------
% mask Masking
%--------------------------------------------------------------------------
mask        = cfg_choice;
mask.tag    = 'mask';
mask.name   = 'Masking';
mask.help   = {''};
mask.values = {none contrast image};
mask.val    = {none};

%--------------------------------------------------------------------------
% conspec Contrast query
%--------------------------------------------------------------------------
conspec      = cfg_branch;
conspec.tag  = 'conspec';
conspec.name = 'Contrast query';
conspec.val  = {titlestr contrasts threshdesc thresh extent conjunction mask};
conspec.help = {''};

%--------------------------------------------------------------------------
% generic Contrasts
%--------------------------------------------------------------------------
generic        = cfg_repeat;
generic.tag    = 'generic';
generic.name   = 'Contrasts';
generic.help   = {''};
generic.values = {conspec};
generic.num    = [1 Inf];

%--------------------------------------------------------------------------
% units Units
%--------------------------------------------------------------------------
units        = cfg_menu;
units.tag    = 'units';
units.name   = 'Data type';
units.help   = {['Data type. This option is only meaningful for M/EEG '...
    'data. Keep the default ''Volumetric'' for any other kind of data.']};
units.labels = {'Volumetric (2D/3D)',...
                'Scalp-Time',...
                'Scalp-Frequency',...
                'Time-Frequency',...
                'Frequency-Frequency'}';
units.values = { 1 2 3 4 5 };
units.val    = { 1 };

%--------------------------------------------------------------------------
% basename Basename
%--------------------------------------------------------------------------
basename         = cfg_entry;
basename.tag     = 'basename';
basename.name    = 'Basename';
basename.help    = {'Enter basename of output files ''spm?_????_<basename>.ext''.'};
basename.strtype = 's';
basename.num     = [1 Inf];

%--------------------------------------------------------------------------
% background Background image
%--------------------------------------------------------------------------
background         = cfg_files;
background.tag     = 'background';
background.name    = 'Background image';
background.help    = {'Background image.'};
background.filter  = {'image'};
background.ufilter = '.*';
background.num     = [1 1];

%--------------------------------------------------------------------------
% orientation Image orientation
%--------------------------------------------------------------------------
orientation        = cfg_menu;
orientation.tag    = 'orientation';
orientation.name   = 'Image orientation';
orientation.help   = {'Image orientation.'};
orientation.labels = {'Axial','Coronal','Sagittal'}';
orientation.values = {'axial','coronal','sagittal'};

%--------------------------------------------------------------------------
% Slices Slices
%--------------------------------------------------------------------------
slices         = cfg_entry;
slices.tag     = 'slices';
slices.name    = 'Slices';
slices.help    = {'Slices to display (mm).'};
slices.strtype = 'r';
slices.num     = [1 Inf];

%--------------------------------------------------------------------------
% nsubj Number of subjects
%--------------------------------------------------------------------------
nsubj         = cfg_entry;
nsubj.tag     = 'nsubj';
nsubj.name    = 'Number of subjects';
nsubj.help    = {'Number of subjects.'};
nsubj.strtype = 'r';
nsubj.num     = [1 1];

%--------------------------------------------------------------------------
% label Label
%--------------------------------------------------------------------------
grplabel         = cfg_entry;
grplabel.tag     = 'label';
grplabel.name    = 'Label';
grplabel.help    = {'Group label.'};
grplabel.strtype = 's';
grplabel.num     = [0 Inf];

%--------------------------------------------------------------------------
% group
%--------------------------------------------------------------------------
group      = cfg_branch;
group.tag  = 'group';
group.name = 'Group';
group.val  = {nsubj grplabel};
group.help = {['Number of subjects and labels per group. ', ...
    'For a single subject analysis, enter "1" and "single subject".']};

%--------------------------------------------------------------------------
% groups
%--------------------------------------------------------------------------
groups        = cfg_repeat;
groups.tag    = 'groups';
groups.name   = 'Groups';
groups.help   = {['Number of groups. ', ...
    'For a single subject analysis, specify one group.']};
groups.values = {group};
groups.num    = [1 Inf];

%--------------------------------------------------------------------------
% modality Modality
%--------------------------------------------------------------------------
modality        = cfg_menu;
modality.tag    = 'modality';
modality.name   = 'Modality';
modality.help   = {'Modality.'};
modality.labels = {'Anatomical MRI',...
                   'Functional MRI',...
                   'Diffusion MRI',...
                   'PET',...
                   'SPECT',...
                   'EEG',...
                   'MEG'
}';
modality.values = {'AMRI','FMRI','DMRI','PET','SPECT','EEG','MEG'};

%--------------------------------------------------------------------------
% refspace Reference space
%--------------------------------------------------------------------------
refspace        = cfg_menu;
refspace.tag    = 'refspace';
refspace.name   = 'Reference space';
refspace.help   = {['Reference space. For an experiment completed only ',...
    'within SPM, choose one of the first four options.']};
refspace.labels = {'Subject space (no normalisation)',...
                   'Normalised space (using segment)',...
                   'Normalised space (using old segment)',...
                   'Customised space',...
                   'Other normalised MNI space',...
                   'Other normalised Talairach space',...
}';
refspace.values = {'subject','ixi','icbm','custom','mni','talairach'};

%--------------------------------------------------------------------------
% exports
%--------------------------------------------------------------------------
exports{1}      = cfg_branch;
exports{1}.tag  = 'tspm';
exports{1}.name = 'Thresholded SPM';
exports{1}.val  = { basename };
exports{1}.help = {'Save filtered SPM{.} as an image.'};

exports{end+1}    = cfg_branch;
exports{end}.tag  = 'binary';
exports{end}.name = 'All clusters (binary)';
exports{end}.val  = { basename };
exports{end}.help = {'Save filtered SPM{.} as a binary image.'};

exports{end+1}    = cfg_branch;
exports{end}.tag  = 'nary';
exports{end}.name = 'All clusters (n-ary)';
exports{end}.val  = { basename };
exports{end}.help = {'Save filtered SPM{.} as an n-ary image.'};

pf = spm_print('format');
for i=1:numel(pf)
    exports{end+1}    = cfg_const;
    exports{end}.tag  = pf(i).label{1};
    exports{end}.name = pf(i).name;
    exports{end}.val  = { true };
    exports{end}.help = {pf(i).name};
end
exports{end+1}    = cfg_const;
exports{end}.tag  = 'csv';
exports{end}.name = 'CSV file';
exports{end}.val  = { true };
exports{end}.help = {exports{end}.name};
if ispc
    exports{end+1}    = cfg_const;
    exports{end}.tag  = 'xls';
    exports{end}.name = 'Excel spreadsheet file';
    exports{end}.val  = { true };
    exports{end}.help = {exports{end}.name};
end

exports{end+1}    = cfg_branch;
exports{end}.tag  = 'montage';
exports{end}.name = 'Montage';
exports{end}.val  = { background orientation slices };
exports{end}.help = {'Display montage.'};

exports{end+1}    = cfg_branch;
exports{end}.tag  = 'nidm';
exports{end}.name = 'NIDM (Neuroimaging Data Model)';
exports{end}.val  = {modality refspace groups};
exports{end}.help = {exports{end}.name};

%--------------------------------------------------------------------------
% export Export results
%--------------------------------------------------------------------------
export        = cfg_repeat;
export.tag    = 'export';
export.name   = 'Export results';
export.help   = {['Select the export format you want. PostScript (PS) is '...
               'the only format that allows to append figures to the same ' ...
               'file.']};
export.values = exports;
uiprintdef    = spm_get_defaults('ui.print');
for i=1:numel(exports)
    if strcmp(uiprintdef,exports{i}.tag)
        export.val = exports(i);
        break;
    end
end

%--------------------------------------------------------------------------
% results Results Report
%--------------------------------------------------------------------------
results      = cfg_exbranch;
results.tag  = 'results';
results.name = 'Results Report';
results.val  = {spmmat generic units export};
results.help = {'Statistical inference.'};
results.prog = @spm_run_results;
results.vout = @vout_results;

%==========================================================================
function dep = vout_results(job)

dep(1)            = cfg_dep;
dep(1).sname      = 'xSPM Variable';
dep(1).src_output = substruct('.','xSPMvar');
dep(1).tgt_spec   = cfg_findspec({{'strtype','e'}});

dep(2)            = cfg_dep;
dep(2).sname      = 'TabDat Variable';
dep(2).src_output = substruct('.','TabDatvar');
dep(2).tgt_spec   = cfg_findspec({{'strtype','e'}});

for i=1:numel(job.export)
    fn = char(fieldnames(job.export{i}));
    if ismember(fn,{'tspm','binary','nary'})
        dep(end+1)          = cfg_dep;
        dep(end).sname      = 'Filtered image';
        dep(end).src_output = substruct('.','filtered');
        dep(end).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
        break; % so far, a single dependency pointing to the last exported image
    end
end

author: spm


text:
function T = build_toep( c, k, n )
%
%  T = build_toep( c, k, n );
%
%  Given:
%    c - the nonzero part of a central column of a banded Toeplitz
%        matrix
%    k - index of c containing the diagonal element of T
%    n - dimension of the banded Toeplitz matrix
%
%  The banded Toeplitz matrix is constructed explicitly.
%

%  J. Nagy  2/11/02

m = length( c );

col = zeros(n,1);
row = col';
col(1:m-k+1,1) = c(k:m);
row(1,1:k) = c(k:-1:1)';
T = toeplitz( col, row );

repo: IRtools


text:
function [ J ] = vec2jacl( vec )

tolerance = 1e-12;

if size(vec,1) == 3

    phi = vec;

    ph = norm(phi);
    if ph < tolerance
        % If the angle is small, fall back on the series representation
        J = vec2jaclSeries(phi,10);
    else
        axis = phi/ph;

        cph = (1 - cos(ph))/ph;
        sph = sin(ph)/ph;

        J = sph * eye(3) + (1 - sph) * axis * axis' + cph * hat(axis);
    end

elseif size(vec,1) == 6

    phi = vec(1:3);
    rho = vec(4:6);

    ph = norm(phi);
    if ph < tolerance
        % If the angle is small, fall back on the series representation
        J = vec2jaclSeries(phi,10);
    else
        Jsmall = vec2jacl( phi );
        Q = vec2Ql( vec );
        J = [ Jsmall Q; zeros(3) Jsmall ];
    end
end
end

author: mbrossar


text:
function [cost grad] = svmCost(w, X, y, lambda)
% cost = HingeLoss^2 + lambda*||w||^2
% 1 2 3 4 5 step
yp = X*w;
idx = find(yp.*y<1);
err = yp(idx)-y(idx);
cost = err'*err + lambda*w'*w;
grad = 2*X(idx,:)'*err + 2*lambda*w;
end

author: Grootzz

repo: GLCM-SVM

save_path: github-repos/MATLAB/Grootzz-GLCM-SVM


text:
function bessel_y1_values_test ( )

%*****************************************************************************80
%
%% BESSEL_Y1_VALUES_TEST demonstrates the use of BESSEL_Y1_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BESSEL_Y1_VALUES_TEST:\n' );
  fprintf ( 1, '  BESSEL_Y1_VALUES stores values of\n' );
  fprintf ( 1, '  the Bessel Y1 function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X            Y1(X)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = bessel_y1_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [l, sh] = TellAfni_Commands (f, fpart )
%
%   [l, sh] = TellAfni_Commands (f, field)
%
%Purpose:
%   Extract, all the commands from README.driver
%   No guarantee that you'll get all the commands extracted,
%   Use this as a guide but read the entire README.driver for the full picture
%
%Input Parameters:
%   f: Full path and name of README.driver file
%      pass no parameters if you want the function to try and
%      locate the file automatically.
%      Use '' if you want to look for the file interactively.
%   field: Particular command of interest.
%          If specified, the help section from README.driver
%          is returned in sh.
%
%Output Parameters:
%   l a cell string of the possible parameters
%   sh: help string corresponding to field
%
%
%Key Terms:
%
%More Info :
%    TellAfni
%    TellAfni_Commands
%    Test_TellAfni
%    AFNI's README.driver
%    NewCs
%
%     Author : Ziad Saad
%     Date : Wed Dec 7 15:34:57 EST 2005
%     SSCC/NIMH/ National Institutes of Health, Bethesda Maryland

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

%Debug Flag
DBG = 1;

%initailize return variables
err = 1;

if (nargin == 0), %try to find README.driver
   f = '';
   fpart = '';
elseif (nargin == 1),
   fpart = '';
end

if (isempty(f)),
   [s,w] = unix('locate README.driver');
   if (~isempty(w)), %got a hit
      i = find(isspace(w));
      f = w(1:i(1)-1);
   end
end

if (~isempty(f)),
   if (exist(f,'file') == 2),
      fprintf(1,'Using file %s for info.\n', f);
   else
      fprintf(2,'Failed to find README file %s\n', f);
      return;
   end
else
   f = '';
end

if (isempty(f)),
   [f, p] = uigetfile('README.driver', 'Pick a README.driver');
   if (f == 0),
      return;
   end
   f = sprintf('%s%c%s',p,filesep,f);
end

fid = fopen(f, 'r');
if (fid < 0),
   fprintf(2,'Failed to read %s\n', f);
   return;
end

s = fscanf(fid,'%c');
ns = length(s);
fclose(fid);

[err,inl] = FindChar(s,'NewLine');
ninl = length(inl);

ncm = 50;
cnt = 1;
S = char(zeros(ninl, ncm));
sl = [];
if (ninl),
   for (i=1:1:ninl-1),
      [err,w] = GetWord(s(inl(i):min(inl(i)+ncm, ns)),1);
      w = zdeblank(w);
      nw = length(w);
      if (nw > 1),
         if ( (sum(abs(double(upper(w))-double(w))) == 0) & sum(isletter(w)) > nw/2),
            % fprintf(2,'>%s<\t',w); pause
            S(cnt,1:nw) = w;
            %store location in s
            sl(cnt) = inl(i);
            cnt = cnt + 1;
         end
      end
   end
else
   fprintf(2,'Failed to understand file.\n');
   return;
end

if (cnt < 2),
   fprintf(2,'Failed to understand file, found nothing.\n');
   return;
end

[l, ilu] = unique(cellstr (S(1:cnt-1,:)));
ils = sl(ilu);

%do they want particular info ?
sh = '';
if (~isempty(fpart)),
   fpart = upper(zdeblank(fpart));
   cnt = 1;
   fnd = 0;
   while (cnt <= length(l)),
      if (~isempty(findstr(char(l(cnt)),fpart))), %(strncmp(char(l(cnt)),fpart,length(fpart))),
         fnd = cnt;
         %find location of next command
         dils = (ils - ils(cnt)); dils(find(dils<=0)) = max(dils + 1); [jj,imin] = min(dils);
         %here is the string
         sh = [sh s(ils(cnt)+1:min(ils(imin), length(s)))];
      end
      cnt = cnt + 1;
   end
   if (~fnd),
      fprintf(2,'Failed to find help for %s.\n', fpart);
      return;
   end
end

return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function params = init_default_params(params)

% Initialize default parameters
default_params.use_gpu = false;
default_params.gpu_id = [];
default_params.clamp_position = false;

def_param_names = fieldnames(default_params);
for k = 1:numel(def_param_names)
    param_name = def_param_names{k};
    if ~isfield(params, param_name)
        params.(param_name) = default_params.(param_name);
    end
end

params.visualization = params.visualization;

repo: STRCF


text:
function x = mono_upto_next_grlex ( m, n, x )

%*****************************************************************************80
%
%% MONO_UPTO_NEXT_GRLEX: grlex next monomial with total degree up to N.
%
%  Discussion:
%
%    We consider all monomials in an M dimensional space, with total
%    degree up to N.
%
%    For example:
%
%    M = 3
%    N = 3
%
%    #  X(1)  X(2)  X(3)  Degree
%      +------------------------
%    1 |  0     0     0        0
%      |
%    2 |  0     0     1        1
%    3 |  0     1     0        1
%    4 |  1     0     0        1
%      |
%    5 |  0     0     2        2
%    6 |  0     1     1        2
%    7 |  0     2     0        2
%    8 |  1     0     1        2
%    9 |  1     1     0        2
%   10 |  2     0     0        2
%      |
%   11 |  0     0     3        3
%   12 |  0     1     2        3
%   13 |  0     2     1        3
%   14 |  0     3     0        3
%   15 |  1     0     2        3
%   16 |  1     1     1        3
%   17 |  1     2     0        3
%   18 |  2     0     1        3
%   19 |  2     1     0        3
%   20 |  3     0     0        3
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the maximum degree.
%    0 <= N.
%
%    Input, integer X(M), the current monomial.
%    To start the sequence, set X = [ 0, 0, ..., 0, 0 ].
%
%    Output, integer X(M), the next monomial.
%    The last value in the sequence is X = [ N, 0, ..., 0, 0 ].
%
  if ( n < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MONO_UPTO_NEXT_GRLEX - Fatal error!\n' );
    fprintf ( 1, '  N < 0.\n' );
    error ( 'MONO_UPTO_NEXT_GRLEX - Fatal error!' );
  end

  if ( sum ( x(1:m) ) < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MONO_UPTO_NEXT_GRLEX - Fatal error!\n' );
    fprintf ( 1, '  Input X sums to less than 0.\n' );
    error ( 'MONO_UPTO_NEXT_GRLEX - Fatal error!' );
  end

  if ( n < sum ( x(1:m) ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MONO_UPTO_NEXT_GRLEX - Fatal error!\n' );
    fprintf ( 1, '  Input X sums to more than N.\n' );
    error ( 'MONO_UPTO_NEXT_GRLEX - Fatal error!' );
  end

  if ( n == 0 )
    return
  end

  if ( x(1) == n )
    x(1) = 0;
  else
    x = mono_next_grlex ( m, x );
  end

  return
end


text:
function [ var_node, var_num ] = ns_t6_var_count ( element_num, ...
  element_node, node_num )

%*****************************************************************************80
%
%% NS_T6_VAR_COUNT counts the Navier Stokes variables on a T6 grid.
%
%  Discussion:
%
%    We are given a mesh of T6 elements, and asked to count, in advance,
%    the number of Navier-Stokes variables associated with the grid.
%    In particular, every node has two velocity variables associated with
%    it, but only a node that is a vertex of the element will also have
%    an associated pressure variable.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM);
%    ELEMENT_NODE(I,J) is the global index of local node I in element J.
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Output, integer VAR_NODE(NODE_NUM+1), used to find the variables
%    associated with a given node, which are in VAR in locations
%    VAR_NODE(NODE) to VAR_NODE(NODE+1)-1.  Note that the last entry of
%    this array points to the location just after the last location in VAR.
%
%    Output, integer VAR_NUM, the number of variables.
%
  element_order = 6;
%
%  Our job is easy once we determine which nodes are vertices.
%  So to begin with, let VAR_NODE count the number of variables
%  associated with each node.
%
  var_node(1:node_num) = 2;

  for element = 1 : element_num
    for order = 1 : 3
      node = element_node(order,element);
      var_node(node) = 3;
    end
  end
%
%  Count them.
%
  var_num = sum ( var_node(1:node_num) );
%
%  Make pointers.
%
  total = 1;

  for node = 1 : node_num
    num = var_node(node);
    var_node(node) = total;
    total = total + num;
  end
  var_node(node_num+1) = total;

  return
end


text:
function out= proc_dBPercentiles(epo, varargin)
%PROC_DBAVERAGE - Classwise calculated averages for dB-scaled features
%
%This functions is exactly used as proc_average. It should be used
%for dB-scaled features (e.g. output of proc_power2dB; or
%proc_spectrum in the default setting 'scaling', 'dB').
if nargin==0,
  out= proc_average; return;
end

epo = misc_history(epo);
out= epo;
%% scale back
out.x= 10.^(epo.x/10);
out= rmfield(out, 'yUnit');  % otherwise we will enter an infinite recursion

%% average
out= proc_percentiles(out, varargin{:});

%% re-convert to dB
out.x= 10*log10(out.x);
out.yUnit= 'dB';

author: bbci


text:
function N = mtimes(K, M)
% N = mtimes(K, M)
%
%  kronMatrix multiplication;
%     multiply a kronMatrix by a matrix, a vector, or by
%     another kronMatrix (if possible),
%
%

% 9/2002 L. Perrone
% written for new kronMatrix class

% 6/2003 J. Nagy
% some minor modifications to incorporate new subsref.m capabilities

if (isa(K, 'kronMatrix'))
   if isa(M, 'double')
      N = left_mtimes(K, M);
   elseif isa(M,'kronMatrix')
      if length(M) == 1
         sizeMa = size(M.a{1});
         sizeMb = size(M.b{1});
         l = length(K);
         Anew = cell(l,1);
         Bnew = cell(l,1);
         for i = 1:l
           sizeKa = size(K.a{i});
           sizeKb = size(K.b{i});
           if sizeKa(2)==sizeMa(1) & sizeKb(2)==sizeMb(1)
             Anew{i} = K.a{i} * M.a{1};
             Bnew{i} = K.b{i} * M.b{1};
           else
             error('Kron factors must be of compatible sizes for multiplication')
           end
         end % end the for i=1:l loop
         N = kronMatrix(Anew,Bnew);
      elseif length(K) == 1
         sizeKa = size(K.a{1});
         sizeKb = size(K.b{1});
         l = length(M);
         Anew = cell(l,1);
         Bnew = cell(l,1);
         for i = 1:l
           sizeMa = size(M.a{i});
           sizeMb = size(M.b{i});
           if sizeKa(2)==sizeMa(1) & sizeKb(2)==sizeMb(1)
             Anew{i} = K.a{1} * M.a{i};
             Bnew{i} = K.b{1} * M.b{i};
           else
             error('Kron factors must be of compatible sizes for multiplication')
           end
         end % end the for i=1:l loop
         N = kronMatrix(Anew,Bnew);
      else
         error('This currently only works if one of the kronMatrix objects has one term in it''s sum.')
      end
   else % M isn't a 'double' or a 'kronMatrix'
      error('Wrong input arguments')
   end % end the case when K is a 'kronMatrix'

elseif isa(K,'double') & isa(M,'kronMatrix')
   N = right_mtimes(K, M);
else
   error('Wrong input arguments')
end


text:
function I = gcmi_model_cd(x, y, Ym)
% GCMI_MODEL_CD Gaussian-Copula Mutual Information between a continuous and a
%         discrete variable in bits based on ANOVA style model comparison.
%   I = gcmi_model_gd(x,y,Ym) returns the MI between the (possibly multidimensional)
%   continuous variable x and the discrete variable y.
%   For 1D x this is a lower bound to the mutual information.
%   Rows of x correspond to samples, columns to dimensions/variables.
%   (Samples first axis)
%   y should contain integer values in the range [0 Ym-1] (inclusive).
%   See also: GCMI_MIXTURE_CD

% ensure samples first axis for vectors
if isvector(x)
    x = x(:);
end
if ndims(x)~=2
    error('gcmi_model_cd: input array should be 2d')
end

if isvector(y)
    y = y(:);
else
    error('gcmi_model_cd: only univariate discrete variable supported');
end

Ntrl = size(x,1);
Nvar = size(x,2);

if size(y,1) ~= Ntrl
    error('gcmi_model_cd: number of trials do not match');
end

% check for repeated values
for xi=1:Nvar
    if length(unique(x(:,xi)))./Ntrl < 0.9
        warning('Input x has more than 10% repeated values.')
        break
    end
end

% check values of discrete variable
if min(y)~=0 || max(y)~=(Ym-1) || any(round(y)~=y)
    error('Values of discrete variable y are not correct')
end

% copula normalisation
cx = copnorm(x);
% parametric Gaussian MI
I = mi_model_gd(cx,y,Ym,true,true);


text:
% MatrixUser, a multi-dimensional matrix analysis software package
% https://sourceforge.net/projects/matrixuser/
%
% The MatrixUser is a matrix analysis software package developed under Matlab
% Graphical User Interface Developing Environment (GUIDE). It features
% functions that are designed and optimized for working with multi-dimensional
% matrix under Matlab. These functions typically includes functions for
% multi-dimensional matrix display, matrix (image stack) analysis and matrix
% processing.
%
% Author:
%   Fang Liu <leoliuf@gmail.com>
%   University of Wisconsin-Madison
%   Aug-30-2014

% --- Updata waitbar for processing DICOM files

function MU_update_waitbar(axes_handle,value,total_value)

axes(axes_handle);
cla(axes_handle);
axis(axes_handle,[0,total_value,0,1]);
patch([0,value,value,0],[0,0,1,1],'r');
% axis off;

if value==1 & value~=total_value
    set(axes_handle,'Visible','on');
elseif value==total_value
    set(axes_handle,'Visible','off');
    axes(axes_handle);
    cla(axes_handle);
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function varargout = rosenbrock(x)
%ROSENBROCK Rosenbrock's 'banana' function in 2D

if nargin == 0
    varargout{1} = [-2,-1];                 % LB
    varargout{2} = [2,3];                   % UB
    varargout{3} = [0,3000];                % Zlim
    varargout{4} = [-26.5,60];              % view
    varargout{5} = [1,1];                   % xmin
    varargout{6} = 'Rosenbrock function';   % name
    varargout{7} = [-0.5,2.5];              % default x0
    return;
end

if size(x,2) == 1; x = x'; end

varargout{1} = 100*(x(:,2)-x(:,1).^2).^2+(1-x(:,1)).^2;
if nargout > 1
    varargout{2} = [-400*(x(:,2)-x(:,1).^2).*x(:,1)-2*(1-x(:,1)), 200*(x(:,2)-x(:,1).^2)];
end

end

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function F = ismember(x,F)
% Internal class for constraint list

F = replace(F,recover(depends(F)),x);


text:
function val=jsonopt(key,default,varargin)
%
% val=jsonopt(key,default,optstruct)
%
% setting options based on a struct. The struct can be produced
% by varargin2struct from a list of 'param','value' pairs
%
% authors:Qianqian Fang (q.fang <at> neu.edu)
%
% $Id: loadjson.m 371 2012-06-20 12:43:06Z fangq $
%
% input:
%      key: a string with which one look up a value from a struct
%      default: if the key does not exist, return default
%      optstruct: a struct where each sub-field is a key
%
% output:
%      val: if key exists, val=optstruct.key; otherwise val=default
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of jsonlab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

val=default;
if(nargin<=2) return; end
opt=varargin{1};
if(isstruct(opt))
    if(isfield(opt,key))
       val=getfield(opt,key);
    elseif(isfield(opt,lower(key)))
       val=getfield(opt,lower(key));
    end
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
% This is a test script which demonstrates the usage of the "MofN_TrackInitiatorX" class.
% =========================================================================>

% Load the ground truth data
load('multiple-robot-tracking.mat');

% Plot settings
ShowPlots = 0;              % Set to 0 to hide plots
numTrueTracks = 3;

% Model parameter shortcuts
lambdaV = 10; % Expected number of clutter measurements over entire surveillance region
V = 10^2;     % Volume of surveillance region (10x10 2D-grid)
V_bounds = [0 10 0 10]; % [x_min x_max y_min y_max]
P_D = 1;    % Probability of detection
timestep_duration = duration(0,0,1);

%% Models
transition_model = ConstantVelocityX('VelocityErrVariance', 0.0001,...
                                     'NumDims', 2,...
                                     'TimestepDuration', timestep_duration);
measurement_model = LinearGaussianX('NumMeasDims', 2,...
                                    'NumStateDims', 4,...
                                    'MeasurementErrVariance', 0.02,...
                                    'Mapping', [1 3]);
clutter_model = PoissonRateUniformPositionX('ClutterRate',lambdaV,...
                                            'Limits',[V_bounds(1:2);...
                                                      V_bounds(3:4)]);
detection_model = ConstantDetectionProbabilityX('DetectionProbability',P_D);

% Compile the State-Space model
model = StateSpaceModelX(transition_model,measurement_model,'Clutter',clutter_model, 'Detection', detection_model);

%% Generate DataList
meas_simulator = MultiTargetMeasurementSimulatorX('Model',model);
% meas_simulator.DetectionProbability = 1;
DataList = meas_simulator.simulate(GroundTruthStateSequence);
N = numel(DataList);

%% Base Filter
obs_covar= measurement_model.covar();
PriorState = GaussianStateX(zeros(4,1), transition_model.covar() + blkdiag(obs_covar(1,1), 0, obs_covar(2,2),0));
base_filter = KalmanFilterX('Model', model, 'StatePrior', PriorState);

%% Data Associator
config.ClutterModel = clutter_model;
config.DetectionModel = detection_model;
config.Clusterer = NaiveClustererX();
config.Gater = EllipsoidalGaterX(2,'GateLevel',10)';
config.DetectionProbability = P_D;
jpdaf = JointProbabilisticDataAssocX(config);

%% Track Initiator

% Initiate Data Associator
config_pdaf.ClutterModel = clutter_model;
% config.Clusterer = NaiveClustererX();
config_pdaf.Gater = EllipsoidalGaterX(2,'GateLevel',5);
config_pdaf.DetectionProbability = P_D;
pdaf = ProbabilisticDataAssocX(config_pdaf);

% Initiate Tag Generator
tag_gen = RandSampleTagGeneratorX(1:10000);

% Prepare initiator parameters
config_ti.TagGenerator = tag_gen;
config_ti.InitFilter = base_filter;
config_ti.DataAssociator = pdaf;
config_ti.ConfirmThreshold = [8,10];
config_ti.DeleteThreshold = [5,10];
CovarThreshold = 4*PriorState.Covar;
config_ti.CustomDeleteConditionFcn = ...
    @(x,t) t.Filter.StatePosterior.Covar(1,1)>CovarThreshold(1,1) ...
           || t.Filter.StatePosterior.Covar(3,3)>CovarThreshold(3,3);

% Create the track initiator
myti = MofN_TrackInitiatorX(config_ti);

%% START OF SIMULATION
%  ===================>

% Create figure windows
if(ShowPlots)
    img = imread('maze.png');

    % set the range of the axes
    % The image will be stretched to this.
    min_x = 0;
    max_x = 10;
    min_y = 0;
    max_y = 10;

    % make data to plot - just a line.
    x = min_x:max_x;
    y = (6/8)*x;

    figure('units','normalized','outerposition',[0 0 .5 1])
    ax(1) = gca;
end

TrackList = [];
for k=2:N
    fprintf('Iteration = %d/%d\n================>\n',k,N);

    %% Extract DataList at time k
    MeasurementList = DataList(k);
    timestamp_km1 = DataList(k-1).Timestamp;
    timestamp_k = MeasurementList.Timestamp;
    dt = timestamp_k - timestamp_km1;
    transition_model.TimestepDuration = dt;
    fprintf('Timestamp = %s\n================>\n',timestamp_k);

    %% Process JPDAF
    jpdaf.MeasurementList = MeasurementList;
    jpdaf.TrackList = TrackList;
    jpdaf.predictTracks();
    jpdaf.associate();
    jpdaf.updateTracks();

    %% Perform Track initiation
    [TrackList, TentativeTrackList] = myti.initiateTracks(jpdaf.TrackList, MeasurementList, jpdaf.AssocWeightsMatrix);

    %% Plot update step results
    if(ShowPlots)

        cla(ax(1));
        %imagesc(ax(1),[min_x max_x], [min_y max_y], flipud(img));
        hold on;
        if(exist('data_plot','var'))
            delete(data_plot);
        end
        data_inv = measurement_model.finv(MeasurementList.Vectors);
        data_plot = plot(ax(1), data_inv(1,:), data_inv(3,:),'k*','MarkerSize', 10);

        % Plot confirmed tracks
        for j=1:numel(TrackList)
            means = [TrackList{j}.Trajectory.Mean];
            h2 = plot(ax(1), means(1,:),means(3,:),'-','LineWidth',1);
            h2 = plotgaussellipse(TrackList{j}.Filter.StatePosterior.Mean([1 3]),...
                                  TrackList{j}.Filter.StatePosterior.Covar([1 3],[1 3]),...
                                  'Color','r',...
                                  'Axis',ax(1));
        end

        % Plot tentative tracks
        for j=1:numel(TentativeTrackList)
            h2 = plot(ax(1), TentativeTrackList{j}.Filter.StatePosterior.Mean(1),TentativeTrackList{j}.Filter.StatePosterior.Mean(3),'.','LineWidth',1);
            if j==2
                set(get(get(h2,'Annotation'),'LegendInformation'),'IconDisplayStyle','off');
            end
            h2 = plotgaussellipse(TentativeTrackList{j}.Filter.StatePosterior.Mean([1 3]),...
                                  TentativeTrackList{j}.Filter.StatePosterior.Covar([1 3],[1 3]),...
                                  'Color','g',...
                                  'Axis',ax(1));
        end

        % set the y-axis back to normal.
        set(ax(1),'ydir','normal');
        str = sprintf('Robot positions (Update)');
        title(ax(1),str)
        xlabel('X position (m)')
        ylabel('Y position (m)')
        axis(ax(1),V_bounds)
        pause(0.01)
    end
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function area = polygon_area_2d ( n, v )

%*****************************************************************************80
%
%% POLYGON_AREA_2D computes the area of a polygon in 2D.
%
%  Discussion:
%
%    AREA = 1/2 * abs ( sum ( 1 <= I <= N ) X(I) * ( Y(I+1) - Y(I-1) ) )
%    where Y(0) should be replaced by Y(N), and Y(N+1) by Y(1).
%
%    If the vertices are given in counterclockwise order, the area
%    will be positive.  If the vertices are given in clockwise order,
%    the area will be negative.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of vertices of the polygon.
%
%    Input, real V(2,N), the vertices.
%
%    Output, real AREA, the area of the polygon.
%
  area = 0.0;

  for i = 1 : n

    im1 = i4_wrap ( i-1, 1, n );
    ip1 = i4_wrap ( i+1, 1, n );

    area = area + v(1,i) * ( v(2,ip1) - v(2,im1) );

  end

  area = 0.5 * area;

  return
end


text:
% Altitute chart containing relative temperature, pressure and density.
h=(0:100:80000)';
chart=h;
for i=2:4,
chart(:,i)=0;
end

author: Sable

repo: mcbench-benchmarks


text:
function symb = holdMarker(reset,symbList2)
%HOLDMARKER  Returns a new symbol for use in plotting or rotates through
% symb = holdMarker()
%   returns a new symbol for use in plotting
%
% holdMarker('reset')
%   will reset the rotation
%   (if there is no output, then the rotation is reset)
% holdMarker( {j} )
%   will reset the rotation to start at symbol position "j"
%
%   Use this with plot( ..., 'marker', holdMarker() )
%
% holdMarker('reset',symbList) sets the symbol list (should be cell array)
%
% holdMarker( h )
%   will apply this function to the handles in h
%   (if h=0, uses h = get(gca,'children') )
%   ex.
%       h = plot( myMatrix );
%       holdMarker();  % reset
%       holdMarker( h ); % apply to all the lines in the plot -- each
%                        % will have a unique symbol.
%
% see also:
% set(gcf,'DefaultAxesColorOrder',[1 0 0;0 1 0;0 0 1],...
%   'DefaultAxesLineStyleOrder','-|--|:')
%   (using a 0 instead of gcf will make it permanent for your matlab session)
% Stephen Becker, May 2011
%
%   See also add_few_markers.m

persistent currentLocation symbListP

if isempty( currentLocation)
    currentLocation = 1;
end

if isempty( symbListP )
    %symbList = { '+','o','*','.','x','s',...
    %'d','^','v','>','<','p','h' };

    % Exclude the "."
    symbList = { '+','o','*','x','s',...
        'd','^','v','>','<','p','h' };
else
%     symbList = symbListP;
    % May 2013
    if ~iscell( symbListP )
        symbList = cell(size(symbListP));
        for jj = 1:length( symbListP )
            symbList{jj} = symbListP(jj);
        end
    else
        symbList = symbListP;
    end
end

% April 2013
if nargin > 0 && isequal(reset,0)
    reset = get( gca, 'children' );
end

do_reset = false;
if nargin > 0
    if ischar(reset)
        do_reset = true;
    elseif iscell( reset )
        currentLocation = reset{1};
    elseif all( ishandle(reset) ) && ~( isscalar(reset) && ( reset == 1) )
        % we have handles
        handles     = reset;
        for h = handles(:)'
            symb    = get_symbol;
            set( h, 'Marker', symb );
        end
    elseif isscalar(reset)
        currentLocation = reset;
    else
        error('invalid input');
    end
end
% if nargout == 0 || do_reset
if false;
    do_reset  = true;
else
    symb = get_symbol;
end

if do_reset
    currentLocation = 1;
    if nargin >= 2
        symbListP = symbList2;
    else
        symbListP = [];
    end
    return;
end

% This subfunction has side-effects
function symb = get_symbol
    if currentLocation > length(symbList), currentLocation = 1; end
    symb    = symbList{ currentLocation };
    currentLocation = currentLocation + 1;
end

end % end of main function

author: stephenbeckr

repo: fastRPCA


text:
%########################################################################
%
%   - PPGI Toolbox -
%   A MATLAB toolbox for Photoplethysmography Imaging (PPGI)
%
% Author   : Christian S. Pilz
% Company  : The Nature of Space of Time
% Date     : 07.05.2019
%
% Contact  : cpi@partofthestars.com
% Web Page : www.partofthestars.com
%
% Version  : beta0.1
%
%########################################################################
%
%	bandpass_filter.m:
%
% Description:
%
%   A  MATLAB class to bandpass filter signals
%

classdef bandpass_filter
    properties
        fs_;
        low_;
        high_;
        low_a_;
        low_b_;
        high_a_;
        high_b_;
    end

   methods

       function obj = bandpass_filter(fs,low,high)

          if nargin > 0
             if isnumeric(fs)
                obj.fs_ = fs;
                obj.low_ = high;
                obj.high_ = low;

                fNorm = obj.low_ / (obj.fs_/2);
                [obj.low_b_,obj.low_a_] = butter(3, fNorm, 'low');

                fNorm = obj.high_ / (obj.fs_/2);
                [obj.high_b_,obj.high_a_] = butter(3, fNorm, 'high');
             else
                error('fs must be numeric')
             end
          end

       end

       function signal_f = get(obj,signal)
           signal_f = filtfilt(obj.low_b_, obj.low_a_, signal);
           signal_f = filtfilt(obj.high_b_, obj.high_a_, signal_f);
       end
   end
end

author: partofthestars

repo: PPGI-Toolbox


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% function B = getGradientNodal(omega,m);
%
% Builds the gradient operator B for a domain defined by omega
% and a nodal grid discretization defined by m:
%
%
% EXAMPLE FOR 2D:
%
% o--x--o--x--o--x--o
% |     |     |     |
% +     +     +     +
% |     |     |     |
% o--x--o--x--o--x--o
% |     |     |     |       WHERE:
% +     +     +     +       'o'  - uc       (nodal)
% |     |     |     |       'x'  - d_1 uc      (staggered-2)
% o--x--o--x--o--x--o       '+'  - d_2 uc      (staggered-1)
%
% Note: the matrices are simple, it is size that matters.
% =============================================================================

function B = getGradientNodal(omega,m)
if nargin == 0,
  help(mfilename);
  runMinimalExample;
  return;
end;

h      =  (omega(2:2:end)-omega(1:2:end))./m;
dim    =  size(omega,2)/2;
id     =  @(i) speye(m(i)+1);   % (m(i)+1) identity matrix
                             % short difference operator
dx    =  @(i) spdiags(ones(m(i),1)*[-1,1],[0,1],m(i),m(i)+1)/h(i);
switch dim
    case 2
      B = [kron(id(2),dx(1));kron(dx(2),id(1))];
      z = sparse(size(B,1),size(B,2));
      B = sparse([B z; z B]);
    case 3
        B = [
        kron(id(3),kron(id(2),dx(1)))
        kron(id(3),kron(dx(2),id(1)))
        kron(dx(3),kron(id(2),id(1)))
        ];
        z = sparse(size(B,1),size(B,2));
        B = sparse([B z z; z B z; z z B]);
    otherwise
        error('Dimension must be either 2 or 3.')
end
%------------------------------------------------------------------------------
function runMinimalExample

% 2D
omega = [0,3,0,5]; m = [10,13];
B2 = getGradientNodal(omega,m);
figure(1); clf;
subplot(1,2,1)
spy(B2);
title(sprintf('%s matrix (2D)',mfilename));

% 3D
omega = [0,1,0,2,0,3]; m = [4,5,6];
B2 = getGradientNodal(omega,m);
subplot(1,2,2);
spy(B2);
title(sprintf('%s matrix (3D)',mfilename));
%==============================================================================


text:
function example ( n )

%*****************************************************************************80
%
%% EXAMPLE is a small example of the return problem.
%
%  Discussion:
%
%    POTENTIAL1 and POTENTIAL2 are identical, except that POTENTIAL1
%    includes a RETURN statement.
%
%    Both functions are called many times.  Somehow, the overhead of
%    the RETURN statement dominates the computation time.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'Call potential1:\n' );
  fprintf ( 1, '\n' );

  tic
  for ii = 1 : n
    vx = potential1 ( 1, 1, 0.5, 0.5 );
  end
  toc

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Call potential2:\n' );
  fprintf ( 1, '\n' );

  tic
  for ii = 1 : n
    vx = potential2 ( 1, 1, 0.5, 0.5 );
  end
  toc

  return
end

function v = potential1 ( a, b, x, y )

%*****************************************************************************80
%
%% POTENTIAL1 evaluates the potential function V(X,Y).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters that define the ellipse.
%
%    Input, real X, Y, the coordinates of the point.
%
%    Output, real V, the value of the potential function.
%
  v = 2.0 * ( ( x / a^2 )^2 + ( y / b^2 )^2 ) + 1.0 / a^2 + 1.0 / b^2;

  return
end

function v = potential2 ( a, b, x, y )

%*****************************************************************************80
%
%% POTENTIAL2 evaluates the potential function V(X,Y).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters that define the ellipse.
%
%    Input, real X, Y, the coordinates of the point.
%
%    Output, real V, the value of the potential function.
%
  v = 2.0 * ( ( x / a^2 )^2 + ( y / b^2 )^2 ) + 1.0 / a^2 + 1.0 / b^2;

% return
end


text:
function distmsr=dtiFrenetDistance(curves1, curves2, samegroupflag)

%Computes distance between curves using frenet framework: curve matching
%method from Bakircioglu et al., HBM 6:329-333 (1998)

%Usage:
%dtiFrenetDistances(fibergroup1, fibergroup2, npoints)
%dtiFrenetDistances(fibergroup, npoints)

%ER 2008 04/2008 SCSNL

%1. Checks
if ~((samegroupflag==0)||(samegroupflag==1))
display('Same group flag should be 1 or zero');
    return;
end

if ~isequal(size(curves1,2), size(curves2,2))
    display('Curves must be resampled to the same number of nodes');
return
end

if (samegroupflag==1)&&~isequal(size(curves1,3), size(curves2,3))
    display('Are you sure your curves1 and curves2 are equivalent? curves2 will be ignored');
end

nfibers1=size(curves1, 3);

%actually if same group flag == 1 then curves2 is ignored

%2. Compute frenet representations for the curves
for i=1:nfibers1
[T(:, :, i),N(:, :, i),B(:, :,i ),v1(:, i),k1(:, i),t1(:, i)] = frenetAll(curves1(1, :, i),curves1(2, :, i),curves1(3, :, i));
end

if samegroupflag==0
nfibers2=size(curves2, 3);
    for i=1:nfibers2
[T(:, :, i),N(:, :, i),B(:, :,i ),v2(:, i),k2(:, i),t2(:, i)] = frenetAll(curves2(1, :, i),curves2(2, :, i),curves2(3, :, i));
end
else
    nfibers2=size(curves1, 3);
end

distmsr=zeros(nfibers1, nfibers2);

%CASE WHERE THE GROUPS OF FIBERS ARE NOT EQUIVALENT
if (samegroupflag==0)
display('2 different fiber groups');
    for i=1:nfibers1
        for j=1:nfibers2

            distmsr(i, j)=sum(3.*(k1(:, i).*v1(:, i)-k2(:, j).*v2(:, j)).^2+(t1(:, i).*v1(:, i)-t2(:, j).*v2(:, j)).^2);
        end
      end
elseif (samegroupflag==1)
display('2 equivalent fiber groups');
%A SHORTCUT-CASE, where the two fibergroups supplied are the same.

          for i=1:nfibers1
        for j=i:nfibers2
            distmsr(i, j)=sum(3.*(k1(:, i).*v1(:, i)-k1(:, j).*v1(:, j)).^2+(t1(:, i).*v1(:, i)-t1(:, j).*v1(:, j)).^2);
        end
          end

distmsr=distmsr+distmsr'- diag(diag(distmsr)); %To make a full matrix
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
filename='CantileverBeam_Triangle_Linear_Fine';
ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance';'perimeter'};
weights = [1 0.1];
constraint = {'volumeConstraint'};
optimizer = 'SLERP';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';

nsteps = 5;
Vfrac_final = 0.5;
Perimeter_target=1;
optimality_final =1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%  sandwich_solve solves the linear equation X=A*X*B+C
%
%  ::
%
%    [X,retcode]=sandwich_solve(A,B,C)
%
%  Args:
%     - A :
%     - B :
%     - C :
%
%  Returns:
%     :
%     - X :
%     - retcode :
%
%  Note:
%
%  Example:
%
%     See also:
%

author: jmaih


text:
% The vector corresponding to the second
% smallest eigenvalue of the Laplacian matrix
%
% INPUTs: adjacency matrix, nxn
% OUTPUTs: fiedler vector, nx1
%
% Other routines used: laplacianMatrix.m
% GB: last updated, Oct 10 2012

function fv=fiedlerVector(adj)

[V,D]=eig(laplacianMatrix(adj));
[~,Y]=sort(diag(D));
fv=V(:,Y(2));

author: aeolianine

repo: octave-networks-toolbox


text:
function [sol,m,Q,residuals,everything] = solvesos_find_blocks(F,obj,options,params,candidateMonomials)

tol = options.sos.numblkdg;
if tol > 1e-2
    disp(' ');
    disp('-> Are you sure you meant to have a tolerance in numblk that big!')
    disp('-> The options numblkdiag controls the tolerance, it is not a 0/1 switch.')
    disp(' ');
end
options.sos.numblkdg = 0;
[sol,m,Q,residuals,everything] = solvesos(F,obj,options,params,candidateMonomials);

% Save old structure to find out when we have stalled
for i = 1:length(Q)
    oldlengths{i} = length(Q{i});
end

go_on = (sol.problem == 0 | sol.problem == 4);
while go_on

    for sosfun = 1:length(Q)
        Qtemp = Q{sosfun};
        keep = diag(Qtemp)>tol;
        Qtemp(:,find(~keep)) = [];
        Qtemp(find(~keep),:) = [];

        m{sosfun} = m{sosfun}(find(keep));

        Qtemp(abs(Qtemp) < tol) = 0;
        [v1,dummy1,r1,dummy3]=dmperm(Qtemp+eye(length(Qtemp)));
        lengths{sosfun} = [];
        n{sosfun} = {};
        for blocks = 1:length(r1)-1
            i1 = r1(blocks);
            i2 = r1(blocks+1)-1;
            if i2>i1
                n{sosfun}{blocks} = m{sosfun}(v1(i1:i2));
            else
                n{sosfun}{blocks} = m{sosfun}(v1(i1));
            end
            lengths{sosfun} =  [lengths{sosfun}  length(n{sosfun}{blocks})];
        end
        lengths{sosfun} = sort(lengths{sosfun});
    end
    go_on = ~isequal(lengths,oldlengths);
    oldlengths = lengths;
    if go_on
        [sol,m,Q,residuals,everything] = solvesos(F,obj,options,params,n);
        go_on = go_on & (sol.problem == 0 | sol.problem == 4);
        if sol.problem == 1
            disp('-> Feasibility was lost during the numerical block-diagonalization.')
            disp('-> The setting sos.numblkdiag is probably too big')
        end
    end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function x = Inv_GUSFT_CG(b,lambda,MaxIts,ErrTol,Pre,x0,mu)
% Inv_UtU_CG: Invert GUSFT_Toeplitz.
%  Usage:
%    x = Inv_GUSFT_CG(b,lambda,MaxIts,ErrTol,Pre,x0,mu);
%  Inputs:
%    b       vector of length n
%    lambda  vector of Fourier multipliers; length 2*n-1
%    MaxIts  Maximum number of CG iterations. Default 6.
%    ErrTol  Error tolerance. Default 1.e-9.
%    Pre     Boolean variable. Pre = 1 if we apply a preconditioner
%    x0      Initial guess
%    mu      Variable for use with the preconditioner
%  Outputs:
%    x       Vector of length n
%  Description
%    Performs the inverse UtU where U is the nouniform FT. This is
%    an approximate inverse which uses a conjugate gradient
%    solver. Each iteration uses GUSFT_Toeplitz (which can be
%    applied by merely using only 2 FFTs).
% See Also
%   GUSFT_Toeplitz
%
% By Emmanuel candes, 2003-2004

if (nargin < 2)
   error('Not enough input arguments.');
end

if (nargin < 3) | isempty(MaxIts)
   MaxIts = 6;
end

if (nargin < 4) | isempty(ErrTol)
   ErrTol = 1.e-9;
end

if (nargin < 5) | isempty(Pre)
	Pre = 0;
end

if (nargin < 6) | isempty(x0)
	x0 = b;
end

if (nargin < 7) | isempty(mu)
	mu = ones(size(b));
end

% Initialization
xk = x0;

if MaxIts > 0,
	 temp = GUSFT_Toeplitz(xk,lambda);

	 rk=b-temp;
	 if Pre == 0,
	   dk = rk;
	   else
	     dk=ApplyCirculantPrecond(rk,mu);
	 end

	% Conjugate gradient iteration
	for j=1:MaxIts,
		 perr=max(abs(rk));
		 fprintf(['Iteration ',int2str(j-1),': Residual error=']);
		 fprintf('%f\n',perr);

		 if perr>ErrTol,
			  temp = GUSFT_Toeplitz(dk,lambda);
			  a0 = sum(conj(rk).*dk);
			  a1 = sum(conj(dk).*temp);

			  % Update x and residual
			  a=a0/a1;
			  xk=xk+a*dk;
			  rk=rk-a*temp;

			  % Update search direction
			  if Pre == 0,
			    b = sum(abs(rk.^2))/a0;
			    dk = rk+b*dk;
			  else
			    rkp = ApplyCirculantPrecond(rk,mu);
			    b = sum(conj(rkp).*rk)/a0;
			    dk = rkp+b*dk;
			  end
		 end
	end
end

x = xk;

%  References: Jonathan R. Shewchuk "An introduction to the Conjugate
%  Gradient Method Without the Agonizing Pain"
%
%  Amir Averbuch and David L. Donoho "BeamLab 200"

author: thomaskuestner


function y = vl_nnrelu(x,dzdy)
% VL_NNRELU  CNN rectified linear unit
%   Y = VL_NNRELU(X) applies the rectified linear unit to the data
%   X. X can have arbitrary size.
%
%   DZDX = VL_NNRELU(X, DZDY) computes the network derivative DZDX
%   with respect to the input X given the derivative DZDY with respect
%   to the output Y. DZDX has the same dimension as X.

% Copyright (C) 2014 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

if nargin <= 1 || isempty(dzdy)
  y = max(x, single(0)) ;
else
  y = dzdy .* (x > single(0)) ;
end


text:
% Plot simulation results with log x axis.
% Plot only + freq terms.
% Per common practice, do not show carrier term.
% Plot measured, plus specification.
% Dick Benson
% Copyright 2004-2013 The MathWorks, Inc.
index_start = 3;              % move out from carrier at DC and fft window effects
Frame_Length = length(Normalized_Spectrum);
Fvec_meas = Spectrum_dF*((index_start-1):Frame_Length/2);

h2  =  findobj('Tag','plot_VCO');
if isempty(h2)
   h2 =  figure('Tag','plot_VCO')
else
   figure(h2)
end

semilogx(Fvec_meas,10*log10(Normalized_Spectrum(index_start:(Frame_Length/2+1))),....
         Fvec,Lvec,'o','markersize',10,'linewidth',1,'markerfacecolor','g');
xlabel('Hz from Carrier'); ylabel('dBc/Hz'); grid on;
legend('Measured', 'specification');
axis([Fvec(1),Fvec_meas(end),-160,-60])


text:
function display(SO3F,varargin)
% standard output

if check_option(varargin,'skipHeader')
  disp('  bingham component');
else
  displayClass(SO3F,inputname(1),[],'moreInfo',symChar(SO3F),varargin{:});
  disp(' ');
end

if SO3F.antipodal, disp('  antipodal: true'); end

disp(['  kappa: ',xnum2str(SO3F.kappa)]);
disp(['  weight: ' xnum2str(mean(SO3F))])
disp(' ');

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function varargout = tpfp(varargin)
% VL_TPFP  Compute true positives and false positives
%  This is an helper function used by VL_PR(), VL_ROC(), VL_DET().
%
%  See also: VL_PR(), VL_ROC(), VL_DET(), VL_HELP().
[varargout{1:nargout}] = vl_tpfp(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function [ rgb ] = jmol_color( Z )
%JMOL_COLOR Get color assigned to an element.
%   rgb = JMOL_COLOR(Z) returns a vector containing the RGB color assigned
%   to the element with atomic number Z. See
%   http://jmol.sourceforge.net/jscolors/
%
%   rgb = JMOL_COLOR(symbol) returns a vector containing the RGB color
%   assigned to the element with the given chemical symbol

    if ischar(Z)
        Z = chemsym2number(Z);
    end

    colors = [...
     [255,255,255];... H
     [217,255,255];... He
     [204,128,255];... Li
     [194,255,0];... Be
     [255,181,181];... B
     [144,144,144];... C
     [48,80,248];... N
     [255,13,13];... O
     [144,224,80];... F
     [179,227,245];... Ne
     [171,92,242];... Na
     [138,255,0];... Mg
     [191,166,166];... Al
     [240,200,160];... Si
     [255,128,0];... P
     [255,255,48];... S
     [31,240,31];... Cl
     [128,209,227];... Ar
     [143,64,212];... K
     [61,255,0];... Ca
     [230,230,230];... Sc
     [191,194,199];... Ti
     [166,166,171];... V
     [138,153,199];... Cr
     [156,122,199];... Mn
     [224,102,51];... Fe
     [240,144,160];... Co
     [80,208,80];... Ni
     [200,128,51];... Cu
     [125,128,176];... Zn
     [194,143,143];... Ga
     [102,143,143];... Ge
     [189,128,227];... As
     [255,161,0];... Se
     [166,41,41];... Br
     [92,184,209];... Kr
     [112,46,176];... Rb
     [0,255,0];... Sr
     [148,255,255];... Y
     [148,224,224];... Zr
     [115,194,201];... Nb
     [84,181,181];... Mo
     [59,158,158];... Tc
     [36,143,143];... Ru
     [10,125,140];... Rh
     [0,105,133];... Pd
     [192,192,192];... Ag
     [255,217,143];... Cd
     [166,117,115];... In
     [102,128,128];... Sn
     [158,99,181];... Sb
     [212,122,0];... Te
     [148,0,148];... I
     [66,158,176];... Xe
     [87,23,143];... Cs
     [0,201,0];... Ba
     [112,212,255];... La
     [255,255,199];... Ce
     [217,255,199];... Pr
     [199,255,199];... Nd
     [163,255,199];... Pm
     [143,255,199];... Sm
     [97,255,199];... Eu
     [69,255,199];... Gd
     [48,255,199];... Tb
     [31,255,199];... Dy
     [0,255,156];... Ho
     [0,230,117];... Er
     [0,212,82];... Tm
     [0,191,56];... Yb
     [0,171,36];... Lu
     [77,194,255];... Hf
     [77,166,255];... Ta
     [33,148,214];... W
     [38,125,171];... Re
     [38,102,150];... Os
     [23,84,135];... Ir
     [208,208,224];... Pt
     [255,209,35];... Au
     [184,184,208];... Hg
     [166,84,77];... Tl
     [87,89,97];... Pb
     [158,79,181];... Bi
     [171,92,0];... Po
     [117,79,69];... At
     [66,130,150];... Rn
     [66,0,102];... Fr
     [0,125,0];... Ra
     [112,171,250];... Ac
     [0,186,255];... Th
     [0,161,255];... Pa
     [0,143,255];... U
     [0,128,255];... Np
     [0,107,255];... Pu
     [84,92,242];... Am
     [120,92,227];... Cm
     [138,79,227];... Bk
     [161,54,212];... Cf
     [179,31,212];... Es
     [179,31,186];... Fm
     [179,13,166];... Md
     [189,13,135];... No
     [199,0,102];... Lr
     [204,0,89];... Rf
     [209,0,79];... Db
     [217,0,69];... Sg
     [224,0,56];... Bh
     [230,0,46];... Hs
     [235,0,38];... Mt
     ];

    rgb = colors(Z,:);

end


text:
function [G]=erdosRenyi(nv,p,Kreg)
%Funciton [G]=edosRenyi(nv,p,Kreg) generates a random graph based on
%the Erdos and Renyi algoritm where all possible pairs of 'nv' nodes are
%connected with probability 'p'.
%
% Inputs:
%   nv - number of nodes
%   p  - rewiring probability
%   Kreg - initial node degree of for regular graph (use 1 or even numbers)
%
% Output:
%   G is a structure inplemented as data structure in this as well as other
%   graph theory algorithms.
%   G.Adj   - is the adjacency matrix (1 for connected nodes, 0 otherwise).
%   G.x and G.y -   are row vectors of size nv wiht the (x,y) coordinates of
%                   each node of G.
%   G.nv    - number of vertices in G
%   G.ne    - number of edges in G
%
%Created by Pablo Blinder. blinderp@bgu.ac.il
%
%Last update 25/01/2005

%build regular lattice
A=sparse(nv,nv);
Kreg=fix(abs(Kreg)/2);Kreg=(Kreg<1)+Kreg;

for k=1:Kreg
    A=sparse(A+diag(ones(1,length(diag(A,k))),k)+diag(ones(1,length(diag(A,nv-k))),nv-k));
end
ne0=nnz(A);
%find connected pairs
[v1,v2]=find(A);
% P=permPairs(nv);%my version is faster
Dis=(rand(length(v1),1)<=p);%pairs to disconnect
A(v1(Dis),v2(Dis))=0;
vDis=unique([v1(Dis),v2(Dis)]);%disconnected vertices
nDis=ne0-nnz(A);sum(Dis);

%cycle trough disconnected pairs
disconPairs=[v1(Dis),v2(Dis)];
for n=1:nDis
    %choose one of the vertices from the disconnected pair
    i=ceil(rand*size(disconPairs,1));
    j=logical(1+rand>0.5);
    vDisToRec=disconPairs(i,j);
    %find non adjacent vertices and reconnect
    adj=[find(A(:,vDisToRec)) ; find(A(vDisToRec,:))'];
    nonAdj=setdiff(1:nv,adj);
    vToRec=nonAdj(ceil(rand*length(nonAdj)));
    S=sort([vDisToRec vToRec]);
    A(S(1),S(2))=1;
end
[x,y]=getNodeCoordinates(nv);
%make adjacency matrix symetric
A=A+fliplr((flipud(triu(A))));
G=struct('Adj',A,'x',x','y',y','nv',nv,'ne',nnz(A));


text:
function map = haxby(m);
%HAXBY  Haxby color map
%   HAXBY(M) returns an M-by-3 matrix containing a colormap with Haxby's
%   colors, commonly used for displaying bathymetry data.
%   HAXBY, by itself, is the same length as the current colormap.
%
%   For example, to reset the colormap of the current figure:
%
%             colormap(haxby)
%
%   Use
%             colormap(flipud(haxby))
%
%   for bathymetry data (positive downward).
%
%   Colormap is based on the colors used by W. F. Haxby's Gravity
%   field of World's oceans, 1985, developed for geoid and gravity maps.
%   The version used here is formed from a linear interpolation of
%   the GMT color table used by MB-System by David W. Caress and Dale N. Chayes.
%   <http://www.ldeo.columbia.edu/res/pi/MB-System>
%
%   See also HSV, GRAY, PINK, COOL, BONE, COPPER, FLAG, HOT
%   COLORMAP, RGBPLOT.

% Kelsey Jordahl
% Marymount Manhattan College
% Time-stamp: <Fri Oct 30 12:45:12 EDT 2009>

if nargin < 1, m = size(get(gcf,'colormap'),1); end
% mbm_grdplot Haxby color pallette
ncolors=11;
c=[ 37    57   175;    40   127   251;    50   190   255;   106   235   255;
    138   236   174;   205   255   162;   240   236   121;   255   189    87;
    255   161    68;   255   186   133;   255   255   255];
pp=1:(m-1)/(ncolors-1):m;
r=interp1(pp,c(:,1),1:m);
g=interp1(pp,c(:,2),1:m);
b=interp1(pp,c(:,3),1:m);
map=[r' g' b']/255;


text:
function res_feat = extract_features(image, pos, scales, features, gparams, extract_info)

% Sample image patches at given position and scales. Then extract features
% from these patches.
% Requires that cell size and image sample size is set for each feature.

if ~iscell(features)
    error('Wrong input');
end

if ~isfield(gparams, 'use_gpu')
    gparams.use_gpu = false;
end
if ~isfield(gparams, 'data_type')
    gparams.data_type = zeros(1, 'single');
end
if nargin < 6
    % Find used image sample size
    extract_info = get_feature_extract_info(features);
end

num_features = length(features);
num_scales = length(scales);
num_sizes = length(extract_info.img_sample_sizes);

% Extract image patches
img_samples = cell(num_sizes,1);
for sz_ind = 1:num_sizes
    img_sample_sz = extract_info.img_sample_sizes{sz_ind};
    img_input_sz = extract_info.img_input_sizes{sz_ind};
    img_samples{sz_ind} = zeros(img_input_sz(1), img_input_sz(2), size(image,3), num_scales, 'uint8');
    for scale_ind = 1:num_scales
        img_samples{sz_ind}(:,:,:,scale_ind) = get_pixels(image, pos, round(img_sample_sz*scales(scale_ind)),img_input_sz);
    end
end

% Find the number of feature blocks and total dimensionality
num_feature_blocks = 0;
total_dim = 0;
for feat_ind = 1:num_features
    num_feature_blocks = num_feature_blocks + length(features{feat_ind}.fparams.nDim);
    total_dim = total_dim + sum(features{feat_ind}.fparams.nDim);
end

feature_map = cell(1, 1, num_feature_blocks);

% Extract feature maps for each feature in the list
ind = 1;
CNN_ind = -1;
for feat_ind = 1:num_features
    feat = features{feat_ind};

    % get the image patch index
    img_sample_ind = cellfun(@(sz) isequal(feat.img_sample_sz, sz), extract_info.img_sample_sizes);

    % do feature computation
    if feat.is_cell
        num_blocks = length(feat.fparams.nDim);
        feature_map(ind:ind+num_blocks-1) = feat.getFeature(img_samples{img_sample_ind}, feat.fparams, gparams);
        CNN_ind = ind;
    else
        num_blocks = 1;
        feature_map{ind} = feat.getFeature(img_samples{img_sample_ind}, feat.fparams, gparams);
    end

    ind = ind + num_blocks;
end

% Do CNN feature normalization per channel
if(CNN_ind ~= -1)
    CNN_feat = {feature_map{CNN_ind}};
    delta = 1e-7;

    for i = 1:numel(CNN_feat)
        min_ele = repmat(min(min(CNN_feat{i},[],2),[],1), [size(CNN_feat{i},1),size(CNN_feat{i},2),1,1]);
        max_ele = repmat(max(max(CNN_feat{i},[],2),[],1), [size(CNN_feat{i},1),size(CNN_feat{i},2),1,1]) + delta;
        CNN_feat{i} = ( CNN_feat{i} - min_ele ) ./ (max_ele - min_ele);
    end
end

% Combine the features. In particular, the hand-crafted features are concatnated together
% along the channels.
if(CNN_ind ~= -1)
    res_feat = cell(2, 1);
    res_feat{2} = CNN_feat{1};
end

% concatnate the grayscale and HOG features only when the input is the grayscale
% image.
if(CNN_ind ~=-1)
    if(num_features == 3)
        res_feat{1} = cat(3,feature_map{1}, feature_map{2});
    else
        res_feat{1} = cat(3,feature_map{1}, feature_map{2}, feature_map{3});
    end
else
    if(num_features == 3)
        res_feat{1} = cat(3,feature_map{1}, feature_map{2}, feature_map{3});
    else
        res_feat{1} = cat(3,feature_map{1}, feature_map{2});
    end
end

end


text:
function dbmute
%dbmute disables all breakpoints currently set in MATLAB.
% Copyright 2009 The MathWorks, Inc.

    % iterate over each entry in the result of dbstatus,
    % and disable each of the breakpoints.
    breakpoints = dbstatus('-completenames');
    for i=1 : length(breakpoints)
        muteDbStatusEntry(breakpoints(i));
    end

end

function muteDbStatusEntry(dbstatusEntry)
%muteDbStatusEntry disables each breapoint in the given entry.
    for i=1 : length(dbstatusEntry.line)
        file = dbstatusEntry.file;
        line = dbstatusEntry.line(i);
        anonymousIndex = dbstatusEntry.anonymous(i);
        expression = dbstatusEntry.expression{i};

        lineNumberString = [num2str(line) '@' num2str(anonymousIndex)];
        newExpression = createDisabledExpression(expression);

        dbstop(file, lineNumberString, 'if', newExpression);
    end
end

function newExpression = createDisabledExpression(expression)
%createDisabledExpression wraps the given expression in a disabling
%  expression if necessary.
    if (isDisabled(expression))
        newExpression = expression;
    elseif strcmp(expression, '')
        newExpression = 'false';
    else
        newExpression = ['false&&(' expression ')'];
    end
end


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

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym ellipticK (@var{m})
%% Complete elliptic integral of the first kind.
%%
%% The complete elliptic integral of the first kind
%% with parameter @var{m} is defined by:
%% @example
%% @group
%% syms m
%% ellipticK (m)
%%   @result{} ans = (sym) K(m)
%% @end group
%%
%% @group
%% rewrite (ans, 'Integral')         % doctest: +SKIP
%%   @result{} ans = (sym)
%%       π
%%       ─
%%       2
%%       ⌠
%%       ⎮         1
%%       ⎮ ────────────────── dα
%%       ⎮    _______________
%%       ⎮   ╱          2
%%       ⎮ ╲╱  1 - m⋅sin (α)
%%       ⌡
%%       0
%% @end group
%% @end example
%%
%% Examples:
%% @example
%% @group
%% diff (ellipticK (m), m)
%%   @result{} (sym)
%%       -(1 - m)⋅K(m) + E(m)
%%       ────────────────────
%%           2⋅m⋅(1 - m)
%% @end group
%%
%% @group
%% vpa (ellipticK (sym (pi)/4))
%%   @result{} (sym) 2.2252536839853959577044373301346
%% @end group
%% @end example
%%
%% There are other conventions for the inputs of elliptic integrals,
%% @pxref{@@sym/ellipticF}.
%%
%% @seealso{@@sym/ellipke, @@sym/ellipticF, @@sym/ellipticE, @@sym/ellipticPi}
%% @end defmethod

function y = ellipticK (m)
  if (nargin > 1)
    print_usage ();
  end

  % y = ellipticF (sym (pi)/2, m);
  y = elementwise_op ('elliptic_k', m);

end

%!error ellipticK (sym(1), 2)

%!assert (isequal (ellipticK (sym (0)), sym (pi)/2))
%!assert (isequal (ellipticK (sym (-inf)), sym (0)))

%!assert (double (ellipticK (sym (1)/2)), 1.854074677, 10e-10)
%!assert (double (ellipticK (sym (pi)/4)), 2.225253684, 10e-10)
%!assert (double (ellipticK (sym (-55)/10)), 0.9324665884, 10e-11)

%!test
%! % compare to double ellipke
%! m = 1/5;
%! ms = sym(1)/5;
%! [K, E] = ellipke (m);
%! assert (double (ellipticK (ms)), K, -1e-15)
%! assert (double (ellipticE (ms)), E, -1e-15)

%!test
%! % compare to double ellipke
%! if (exist ('OCTAVE_VERSION', 'builtin'))
%! m = -10.3;
%! ms = -sym(103)/10;
%! [K, E] = ellipke (m);
%! assert (double (ellipticK (ms)), K, -1e-15)
%! assert (double (ellipticE (ms)), E, -1e-15)
%! end

%!test
%! % compare to Maple
%! us = vpa (ellipticK (sym (7)), 40);
%! % > evalf(EllipticK(sqrt(7)), 40);
%! maple = vpa ('0.6168027921799632674669917683443602673441', 40) - ...
%!         vpa ('0.9114898734184488922164103102629560336918j', 40);
%! assert (abs (double (maple - us)), 0, 1e-39)

repo: octsympy


text:
function [X_Out, StructDat] = convertVecStruct(X_In,StructDat)
%
%[X_Out, StructDat] = convertVecStruct(X_In,StructDat)
%
%This function is used to flatten a state struct to a vector or vice versa.
%Each field of the struct must contain a matrix (or scalar) of doubles.
%
%INPUTS:
% X_In = either a vector or struct of data
% StructDat = a struct with information about the fields in X_In. On the
%     first call with X_In = Struct, this field should be ommitted or [].
%     After the first call, the Struct -> vector conversion is faster if
%     this argument is included.
%
%OUTPUTS:
% X_Out = either a vector or struct - opposite of the input
% StructDat = a struct with information about the fields in X
%
%
%METHOD:
% If the input is a struct, then this algorithm takes each field and uses
% the matlab reshape command to make the field data into a column vector.
% The algorithm then takes all of the column vectors and combines them into
% one bit long column vector that is returned as x
%
% If the input is a vector, then the algorithm breaks the vector into a
% column vector corresponding to each field, and then uses reshape to get
% each of these column vectors into the correct shape, which is then
% assigned to a field in the return struct with the correct name.
%

if isstruct(X_In)
  %Then the input is a struct - goal is to produce a vector
  if nargin==1 || isempty(StructDat)
      %Then we don't know the struct data yet - compute it
      StructDat.FieldNames = fieldnames(X_In);
      N_Fields = length(StructDat.FieldNames);
      StructDat.FieldSize = zeros(N_Fields,2);
      for i=1:N_Fields
        StructDat.FieldSize(i,:) = size(X_In.(StructDat.FieldNames{i}));
      end
  end

  %Get the number of elements in the data structure
  N_elements = StructDat.FieldSize(:,1).*StructDat.FieldSize(:,2);
  N_Fields= length(N_elements);

  X_Out = zeros(sum(N_elements),1);  %Initialize the vector

  %Store each field of the struct into the vector
  IdxLow = 1;
  for i=1:N_Fields
    IdxUpp = IdxLow-1+N_elements(i);
    fieldName = StructDat.FieldNames{i};
    X_Out(IdxLow:IdxUpp) = reshape(X_In.(fieldName),N_elements(i),1);
    IdxLow = IdxUpp+1;
  end

else
  %Then the input is a vector - goal is to produce a struct

  %Get the number of elements in the data structure
  N_elements = StructDat.FieldSize(:,1).*StructDat.FieldSize(:,2);
  N_Fields = length(N_elements);

  %Store each field of the struct into the vector
  IdxLow = 1;
  for i=1:N_Fields
    IdxUpp = IdxLow-1+N_elements(i);
    fieldName = StructDat.FieldNames{i};
    data = X_In(IdxLow:IdxUpp);
    N = StructDat.FieldSize(i,1);
    M = StructDat.FieldSize(i,2);
    X_Out.(fieldName) = reshape(data,N,M);
    IdxLow = IdxUpp+1;
  end

end

end

% An alternate forumlation that uses fewer lines of code:

% % % % function [X_Out, StructDat] = Convert_State(X_In,StructDat)
% % % % %
% % % % %[X_Out, StructDat] = Convert_State(X_In,StructDat)
% % % % %
% % % % %This function is used to flatten a state struct to a vector or vice versa.
% % % % %Each field of the struct must contain a matrix (or scalar) of doubles.
% % % % %
% % % % %INPUTS:
% % % % % X_In = either a vector or struct of data
% % % % % StructDat = a struct with information about the fields in X_In. On the
% % % % %     first call with X_In = Struct, this field should be ommitted or [].
% % % % %     After the first call, the Struct -> vector conversion is faster if
% % % % %     this argument is included.
% % % % %
% % % % %OUTPUTS:
% % % % % X_Out = either a vector or struct - opposite of the input
% % % % % StructDat = a struct with information about the fields in X
% % % % %
% % % %
% % % % if isstruct(X_In)
% % % %   %Then the input is a struct - goal is to produce a vector
% % % %   if nargin==1 || isempty(StructDat)
% % % %       %Then we don't know the struct data yet - compute it
% % % %       StructDat.FieldNames = fieldnames(X_In);
% % % %       N_Fields = length(StructDat.FieldNames);
% % % %       StructDat.FieldSize = zeros(N_Fields,2);
% % % %       StructDat.FieldIdx = zeros(N_Fields,2);
% % % %       IdxLow = 1;
% % % %       for i=1:N_Fields
% % % %         [N,M] = size(X_In.(StructDat.FieldNames{i}));
% % % %         StructDat.FieldSize(i,:) = [N,M];
% % % %         IdxUpp = IdxLow-1+N*M;
% % % %         StructDat.FieldIdx(i,:) = [IdxLow, IdxUpp];
% % % %         IdxLow = IdxUpp+1;
% % % %       end
% % % %   end
% % % %
% % % %   %Get the number of elements in the data structure
% % % %   N_elements = StructDat.FieldSize(:,1).*StructDat.FieldSize(:,2);
% % % %   N_Fields = length(N_elements);
% % % %
% % % %   X_Out = zeros(sum(N_elements),1);  %Initialize the vector
% % % %
% % % %   %Store each field of the struct into the vector
% % % %   for i=1:N_Fields
% % % %     Idx = StructDat.FieldIdx(i,1):StructDat.FieldIdx(i,2);
% % % %     X_Out(Idx) = reshape(X_In.(StructDat.FieldNames{i}),N_elements(i),1);
% % % %   end
% % % %
% % % % else
% % % %   %Then the input is a vector - goal is to produce a struct
% % % %
% % % %   %Get the number of elements in the data structure
% % % %   N_elements = StructDat.FieldSize(:,1).*StructDat.FieldSize(:,2);
% % % %   N_Fields = length(N_elements);
% % % %
% % % %   %Store each field of the struct into the vector
% % % %   for i=1:N_Fields
% % % %     Idx = StructDat.FieldIdx(i,1):StructDat.FieldIdx(i,2);
% % % %     N = StructDat.FieldSize(i,1);
% % % %     M = StructDat.FieldSize(i,2);
% % % %     X_Out.(StructDat.FieldNames{i}) = reshape(X_In(Idx),N,M);
% % % %   end
% % % %
% % % % end
% % % %
% % % %
% % % % end

author: MatthewPeterKelly

repo: dscTutorials


text:
% Copyright 2019 Jonas Koenemann, Moritz Diehl, University of Freiburg
% Redistribution is permitted under the 3-Clause BSD License terms. Please
% ensure the above copyright notice is visible in any derived work.
%
function deprecation(msg)
  ocl.utils.warning(msg);

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


totalNum = 1867290;

% [epsilon_soft_2_DP_test, time_DP_2_test] = DP_loopCompress(1, '/home/yinhuan/mapModel/yq/weightVector.txt', '/home/yinhuan/mapModel/yq/rows_cut/sparse/visMatrix/', '/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_rows_cuts_2_test.txt', totalNum, 1000, '/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_compress_2_test/');
%
% point_num_DP_2_test = salientNumCnt('/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_compress_2_test/', totalNum );
%
% min_cost_DP_2_test = get_min_cost('/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_compress_2_test/', '/home/yinhuan/mapModel/yq/weightVector.txt', 1, epsilon_soft_2_DP_test);

%% sparse test
% two nums: section_num search_domain
% [epsilon_soft_5_DP_80_max, time_DP_5_80_max] = DP_loopCompress(1, '/home/yinhuan/mapModel/yq/weightVector.txt', '/home/yinhuan/mapModel/yq/rows_cut/sparse/visMatrix/', '/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_rows_cuts_5_80_max.txt', totalNum, 1000, '/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_compress_5_80_max/');
%
% point_num_DP_5_80_max = salientNumCnt('/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_compress_5_80_max/', totalNum );
%
% min_cost_DP_5_80_max = get_min_cost('/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_compress_5_80_max/', '/home/yinhuan/mapModel/yq/weightVector.txt', 1, epsilon_soft_5_DP_80_max);

%
%
% [epsilon_soft_2_DP_50, time_DP_2_50] = DP_loopCompress(1, '/home/yinhuan/mapModel/yq/weightVector.txt', '/home/yinhuan/mapModel/yq/rows_cut/sparse/visMatrix/', '/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_rows_cuts_2_50.txt', totalNum, 1000, '/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_compress_2_50/');
%
% point_num_DP_2_50 = salientNumCnt('/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_compress_2_50/', totalNum );
%
% min_cost_DP_2_50 = get_min_cost('/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_compress_2_50/', '/home/yinhuan/mapModel/yq/weightVector.txt', 1, epsilon_soft_2_DP_50);
%
%
%
%
% [epsilon_soft_2_DP_20, time_DP_2_20] = DP_loopCompress(1, '/home/yinhuan/mapModel/yq/weightVector.txt', '/home/yinhuan/mapModel/yq/rows_cut/sparse/visMatrix/', '/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_rows_cuts_2_20.txt', totalNum, 1000, '/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_compress_2_20/');
%
% point_num_DP_2_20 = salientNumCnt('/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_compress_2_20/', totalNum );
%
% min_cost_DP_2_20 = get_min_cost('/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_compress_2_20/', '/home/yinhuan/mapModel/yq/weightVector.txt', 1, epsilon_soft_2_DP_20);
%
%
%
% [epsilon_soft_2_DP_10, time_DP_2_10] = DP_loopCompress(1, '/home/yinhuan/mapModel/yq/weightVector.txt', '/home/yinhuan/mapModel/yq/rows_cut/sparse/visMatrix/', '/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_rows_cuts_2_10.txt', totalNum, 1000, '/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_compress_2_10/');
%
% point_num_DP_2_10 = salientNumCnt('/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_compress_2_10/', totalNum );
%
% min_cost_DP_2_10 = get_min_cost('/home/yinhuan/mapModel/yq/rows_cut/sparse/DP_compress_2_10/', '/home/yinhuan/mapModel/yq/weightVector.txt', 1, epsilon_soft_2_DP_10);

%% 11.30

DP_save_rows('/home/yinhuan/mapModel/yq_500/greedy/DP/cut_costs.txt', 2, 4/5, '/home/yinhuan/mapModel/yq_500/greedy/DP/DP_rows_cuts_2_80.txt');

[epsilon_soft_2_80_max, time_DP_2_80_max] = DP_loopCompress(1, '/home/yinhuan/mapModel/yq_500/weightVector.txt', '/home/yinhuan/mapModel/yq_500/visMatrix/', '/home/yinhuan/mapModel/yq_500/greedy/DP/DP_rows_cuts_2_80.txt', totalNum, 1000, '/home/yinhuan/mapModel/yq_500/greedy/DP/DP_rows_cut_2_80/');

point_num_DP_2_80_max = salientNumCnt('/home/yinhuan/mapModel/yq_500/greedy/DP/DP_rows_cut_2_80/', totalNum );

min_cost_DP_2_80_max = get_min_cost('/home/yinhuan/mapModel/yq_500/greedy/DP/DP_rows_cut_2_80/', '/home/yinhuan/mapModel/yq_500/weightVector.txt', 1, epsilon_soft_2_80_max);


text:
classdef AcademicProblem < handle

    properties (Access = public)
        result
    end

    properties (Access = private)
        cost
        constraint
        optimizer
        filename
    end

    methods (Access = public)

        function obj = AcademicProblem(cParams)
            obj.init(cParams);
            obj.compute();
        end

        function compute(obj)
            obj.createOptimizer();
        end

    end

    methods (Access = private)

        function init(obj, cParams)
            obj.filename = cParams.filename;
        end

        function createOptimizer(obj)
            run(obj.filename);
            d                            = DesignVariableAcademic();
            d.init(x0);
            j.dV                         = d;
            c.dV                         = d;
            s.designVar                  = d;
            s.cost                       = AcademicCost(j);
            s.constraint                 = AcademicConstraint(c);
            s.constraint.nSF             = nConstr;
            s.nConstraints               = nConstr;
            s.dualVariable               = DualVariable(s);
            s.outputFunction.type        = "Academic";
            s.outputFunction.iterDisplay = "off";
            s.outputFunction.monitoring  = MonitoringManager(s);
            s.optimizerNames.primal     = 'PROJECTED GRADIENT';
            opt = Optimizer.create(s);
            opt.solveProblem();
            obj.result = d.value;
            close all;
        end

    end

end


text:
clear all

Nsamples = 41500;
EulerSaved = zeros(Nsamples, 3);

dt = 0.01;

for k = 1:Nsamples
  [p q r] = GetGyro();
  [phi theta psi] = EulerGyro(p, q, r, dt);

  EulerSaved(k, :) = [ phi theta psi ];
end

PhiSaved   = EulerSaved(:, 1) * 180/pi;
ThetaSaved = EulerSaved(:, 2) * 180/pi;
PsiSaved   = EulerSaved(:, 3) * 180/pi;

t = 0:dt:Nsamples*dt-dt;

figure
plot(t, PhiSaved)

figure
plot(t, ThetaSaved)

figure
plot(t, PsiSaved)

author: philbooks

repo: Kalman-Filter-for-Beginners