text:
function conf = spm_cfg_deformations
% Configuration file for deformation jobs
%_______________________________________________________________________
% Copyright (C) 2008-2016 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_cfg_deformations.m 7700 2019-11-21 17:09:15Z john $

hsummary = {
'Utility for working with deformation fields.',...
['They can be loaded, inverted, combined etc, and the results ',...
'either saved to disk, or applied to some image or surface file. ',...
'This utility was intended for imaging experts and may therefore ',...
'be a bit difficult for naive users. ',...
'It provides a great deal of flexibility, which may be confusing to some.']};

hinv = {[...
'Creates the inverse of a deformation field. ',...
'Deformations are assumed to be one-to-one, in which case they ',...
'have a unique inverse.  If y'':A->B is the inverse of y:B->A, then ',...
'y'' o y = y o y'' = Id, where Id is the identity transform.'],...
'',...
'Deformations are inverted using the method described in the appendix of:',...
['    * Ashburner J, Andersson JLR & Friston KJ (2000) ',...
 '"Image Registration using a Symmetric Prior - in Three-Dimensions." ',...
 'Human Brain Mapping 9(4):212-225']};

hcomp = {[...
'Deformation fields can be thought of as mappings. ',...
'These can be combined by the operation of "composition", which is ',...
'usually denoted by a circle "o". ',...
'Suppose x:A->B and y:B->C are two mappings, where A, B and C refer ',...
'to domains in 3 dimensions. ',...
'Each element a in A points to element x(a) in B. ',...
'This in turn points to element y(x(a)) in C, so we have a mapping ',...
'from A to C. ',...
'The composition of these mappings is denoted by yox:A->C. ',...
'Compositions can be combined in an associative way, such that zo(yox) = (zoy)ox.'],...
'',[...
'In this utility, the right-to-left order of the compositions is ',...
'from top to bottom (note that the rightmost deformation would ',...
'actually be applied first).']};

hsn = {[...
'Spatial normalisation, and the unified segmentation model of ',...
'SPM5 save a parameterisation of deformation fields.  These consist ',...
'of a combination of an affine transform, and nonlinear warps that ',...
'are parameterised by a linear combination of cosine transform ',...
'basis functions.  These are saved in *_sn.mat files, which can be ',...
'converted to deformation fields.']};

hvox = {[...
'Specify the voxel sizes of the deformation field to be produced. ',...
'Non-finite values will default to the voxel sizes of the template image',...
'that was originally used to estimate the deformation.']};

hbb = {[...
'Specify the bounding box of the deformation field to be produced. ',...
'Non-finite values will default to the bounding box of the template image',...
'that was originally used to estimate the deformation.']};

himgr = {[...
'Deformations can be thought of as vector fields, and represented ',...
'by three-volume images.  In SPM, deformation fields are saved in ',...
'NIfTI format, with dimensions xdim x ydim x zdim x 1 x 3. ',...
'Each voxel contains the x, y and z mm coordinates of where the deformation points.']};

himgw = {[...
'Save the result as a three-volume image.  "y_" will be prepended to the ',...
'filename.']};

hdetw = {[...
'Save the Jacobian determinants as an image.  "j_" will be prepended to the ',...
'filename.']};

happly = {[...
'Apply the resulting deformation field to some images. ',...
'The filenames will be prepended by "w".']};

hmatname = {...
'Specify the _sn.mat to be used.'};

himg = {...
'Specify the image file on which to base the dimensions, orientation etc.'};

hid = {[...
'This option generates an identity transform, but this can be useful for ',...
'changing the dimensions of the resulting deformation (and any images that ',...
'are generated from it).  Dimensions, orientation etc are derived from ',...
'an image.']};

hidbbvox = {[...
'This option generates an identity transform, but this can be useful for ',...
'changing the dimensions of the resulting deformation (and any images that ',...
'are generated from it).  Dimensions, orientation etc are derived from ',...
'a specified bounding box and voxel dimensions.']};

def          = files('Deformation Field','def','.*y_.*\.nii$',[1 1]);
def.help     = himgr;

matname      = files('Parameter File','matname','.*_sn\.mat$',[1 1]);
matname.help = hmatname;

vox          = entry('Voxel sizes','vox','r',[1 3]);
vox.val      = {[NaN NaN NaN]};
vox.help     = hvox;

bb           = entry('Bounding box','bb','r',[2 3]);
bb.val       = {[NaN NaN NaN;NaN NaN NaN]};
bb.help      = hbb;

sn2def       = branch('Imported _sn.mat','sn2def',{matname,vox,bb});
sn2def.help  = hsn;

img          = files('Image to base Id on','space','nifti',[1 1]);
img.help     = himg;

id           = branch('Identity (Reference Image)','id',{img});
id.help      = hid;

idbbvox      = branch('Identity (Bounding Box and Voxel Size)','idbbvox',{vox, bb});
idbbvox.help = hidbbvox;

ffield = files('Flow field','flowfield','nifti',[1 1]);
ffield.ufilter = '^u_.*';
ffield.help = {...
    ['The flow field stores the deformation information. '...
     'The same field can be used for both forward or backward deformations '...
     '(or even, in principle, half way or exaggerated deformations).']};
%------------------------------------------------------------------------
forbak = mnu('Forward/Backwards','times',{'Backward','Forward'},{[1 0],[0 1]});
forbak.val  = {[1 0]};
forbak.help = {[...
    'The direction of the Dartel flow.  '...
    'Note that a backward transform will warp an individual subject''s '...
    'to match the template (ie maps from template to individual). '...
    'A forward transform will warp the template image to the individual.']};
%------------------------------------------------------------------------
K = mnu('Time Steps','K',...
        {'1','2','4','8','16','32','64','128','256','512'},...
        {0,1,2,3,4,5,6,7,8,9});
K.val  = {6};
K.help = {...
    ['The number of time points used for solving the '...
     'partial differential equations.  A single time point would be '...
     'equivalent to a small deformation model. '...
     'Smaller values allow faster computations, '...
     'but are less accurate in terms '...
     'of inverse consistency and may result in the one-to-one mapping '...
     'breaking down.']};
% ---------------------------------------------------------------------
template        = cfg_files;
template.tag    = 'template';
template.name   = 'Dartel Template';
template.filter = 'nifti';
template.num    = [0 1];
template.val    = {{''}};
template.help   = {...
['Select the final Template file generated by Dartel. This will be affine '...
 'registered with a TPM file, such that the resulting spatially normalised '...
 'images are closer aligned to MNI space. Leave empty if you do not wish to '...
 'incorporate a transform to MNI space '...
 '(ie just click ``done'' on the file selector, without selecting any images).']};
%------------------------------------------------------------------------
drtl = branch('Dartel flow','dartel',{ffield,forbak,K,template});
drtl.help = {'Imported Dartel flow field.'};
%------------------------------------------------------------------------
other = {drtl,def,id,idbbvox,sn2def};

img          = files('Image to base inverse on','space','nifti',[1 1]);
img.help     = himg;

comp0        = repeat('Composition','comp',other);
comp0.help   = hcomp;

iv0          = branch('Inverse','inv',{comp0,img});
iv0.help     = hinv;

comp1        = repeat('Composition','comp',[other,{iv0},{comp0}]);
comp1.num    = [1 Inf];
comp1.help   = hcomp;

iv1          = branch('Inverse','inv',{comp1,img});
iv1.help     = hinv;

comp2        = repeat('Composition','comp',[other,{iv1},{comp1}]);
comp2.num    = [1 Inf];
comp2.help   = hcomp;

iv2          = branch('Inverse','inv',{comp2,img});
iv2.help     = hinv;

comp         = repeat('Composition','comp',[other,{iv2},{comp2}]);
comp.num     = [1 Inf];
comp.help    = hcomp;

saveas       = entry('Save as','ofname','s',[0 Inf]);
saveas.help  = himgw;

savedas       = entry('Save as','ofname','s',[0 Inf]);
savedas.help  = hdetw;

applyto      = files('Apply to','fnames','nifti',[0 Inf]);
applyto.help = happly;

savepwd      = cfg_const;
savepwd.name = 'Current directory';
savepwd.tag  = 'savepwd';
savepwd.val  = {1};
savepwd.help = {['All created files (deformation fields and warped images) ' ...
                 'are written to the current directory.']};

savesrc      = cfg_const;
savesrc.name = 'Source directories';
savesrc.tag  = 'savesrc';
savesrc.val  = {1};
savesrc.help = {['The combined deformation field is written into the ' ...
                 'directory of the first deformation field, warped images ' ...
                 'are written to the same directories as the source ' ...
                 'images.']};

savedef      = cfg_const;
savedef.name = 'Source directory (deformation)';
savedef.tag  = 'savedef';
savedef.val  = {1};
savedef.help = {['The combined deformation field and the warped images ' ...
                 'are written into the directory of the first deformation ' ...
                 'field.']};

saveusr      = files('Output directory','saveusr','dir',[1 1]);
saveusr.help = {['The combined deformation field and the warped images ' ...
                 'are written into the specified directory.']};

savedir      = cfg_choice;
savedir.name = 'Output destination';
savedir.tag  = 'savedir';
savedir.values = {savepwd savesrc saveusr};
savedir.val  = {savepwd};

savedir1      = cfg_choice;
savedir1.name = 'Output destination';
savedir1.tag  = 'savedir';
savedir1.values = {savepwd saveusr};
savedir1.val  = {savepwd};

interp      = cfg_menu;
interp.name = 'Interpolation';
interp.tag  = 'interp';
interp.labels = {'Nearest neighbour','Trilinear','2nd Degree B-spline',...
'3rd Degree B-Spline ','4th Degree B-Spline ','5th Degree B-Spline',...
'6th Degree B-Spline','7th Degree B-Spline','Categorical'};
interp.values = {0,1,2,3,4,5,6,7,-1};
interp.def  = @(val)spm_get_defaults('normalise.write.interp',val{:});
interp.help    = {
                  ['The method by which the images are sampled when ' ...
                  'being written in a different space. ' ...
                  '(Note that Inf or NaN values are treated as zero, ' ...
                  'rather than as missing data)']
                  '    Nearest Neighbour:'
                  '      - Fastest, but not normally recommended.'
                  '    Trilinear Interpolation:'
                  '      - OK for PET, realigned fMRI, or segmentations'
                  '    B-spline Interpolation:'
                  ['      - Better quality (but slower) interpolation' ...
                  '/* \cite{thevenaz00a}*/, especially with higher ' ...
                  'degree splines. Can produce values outside the ' ...
                  'original range (e.g. small negative values from an ' ...
                  'originally all positive image).']
                  '    Categorical:'
                  ['       - Slow (particularly when there are lots of '...
                  'categories). This is intended to warp categorical images ' ...
                  'such as label maps.']
}';

% ---------------------------------------------------------------------
fwhm         = cfg_entry;
fwhm.tag     = 'fwhm';
fwhm.name    = 'Gaussian FWHM';
fwhm.val     = {[0 0 0]};
fwhm.strtype = 'r';
fwhm.num     = [1 3];
fwhm.help    = {'Specify the full-width at half maximum (FWHM) of the Gaussian blurring kernel in mm. Three values should be entered, denoting the FWHM in the x, y and z directions.'};
% ---------------------------------------------------------------------

% ---------------------------------------------------------------------
prefix         = cfg_entry;
prefix.tag     = 'prefix';
prefix.name    = 'Filename Prefix';
prefix.val     = {''};
prefix.strtype = 's';
prefix.num     = [0 Inf];
prefix.help    = {'The name of the output file(s) will be the name of the input file(s) prefixed with this prefix. Leave empty to use SPM default prefixes.'};
% ---------------------------------------------------------------------

% ---------------------------------------------------------------------
mask         = cfg_menu;
mask.tag     = 'mask';
mask.name    = 'Masking';
mask.help    = {'Because of subject motion, different images are likely to have different patterns of zeros from where it was not possible to sample data. With masking enabled, the program searches through the whole time series looking for voxels which need to be sampled from outside the original images. Where this occurs, that voxel is set to zero for the whole set of images (unless the image format can represent NaN, in which case NaNs are used where possible).'};
mask.labels = {
               'Mask images'
               'Dont mask images'
}';
mask.values = {1 0};
mask.val    = {1};
% ---------------------------------------------------------------------

% ---------------------------------------------------------------------
preserve         = cfg_menu;
preserve.tag     = 'preserve';
preserve.name    = 'Preserve';
preserve.help    = {
'Preserve Concentrations: Smoothed spatially normalised images (sw*) represent weighted averages of the signal under the smoothing kernel, approximately preserving the intensities of the original images. This option is currently suggested for eg fMRI.'
''
'Preserve Amount: Smoothed and spatially normalised images preserve the total amount of signal from each region in the images (smw*). Areas that are expanded during warping are correspondingly reduced in intensity. This option is suggested for VBM.'
''
'Preserve Labels: This is intended for warping label images. While it is quite slow to run, it is intended to give more accurately warped categorical data.'
}';
preserve.labels = {
                   'Preserve Concentrations (no "modulation")'
                   'Preserve Amount ("modulation")'
                   'Preserve Labels (categorical data)'
}';
preserve.values = {0 1 2};
preserve.val    = {0};
% ---------------------------------------------------------------------

% ---------------------------------------------------------------------
fromimage       = cfg_files;
fromimage.name   = 'Image Defined';
fromimage.tag    = 'file';
fromimage.filter = 'nifti';
fromimage.num    = [1 1];
fromimage.help   = {'Use the dimensions, orientation etc of some pre-existing image.'};
% ---------------------------------------------------------------------

% ---------------------------------------------------------------------
surfa        = cfg_files;
surfa.name   = 'Surface';
surfa.tag    = 'surface';
surfa.filter = 'mesh';
surfa.num    = [1 Inf];
surfa.help   = {'Select a GIFTI file to warp.'};
% ---------------------------------------------------------------------

bbvox         = cfg_branch;
bbvox.name    = 'User Defined';
bbvox.tag     = 'bbvox';
bbvox.val     = {bb,vox};
bbvox.help   = {[...
'The part of the deformation to use is specified by defining the bounding box and ',...
'voxel sizes that you would like to have. This is probably stating the obvious to many ',...
'but smaller voxels and a broader bounding box will take up more disk space, but may ',...
'give a little more accuracy.']};

deffov        = cfg_choice;
deffov.name   = 'Field of View';
deffov.tag    = 'fov';
deffov.values = {fromimage,bbvox};
deffov.help   = {[...
'The dimensions and voxel size of the resulting deformation may be defined from some image, ',...
'or by specifying voxel sizes and a bounding box.']};

savedef       = cfg_branch;
savedef.name  = 'Save Deformation';
savedef.tag   = 'savedef';
savedef.val   ={saveas,savedir1};
savedef.help  = {'The deformation may be saved to disk as a ``y_*.nii'''' file.'};

savedet       = cfg_branch;
savedet.name  = 'Save Jacobian Determinants';
savedet.tag   = 'savejac';
savedet.val   ={savedas,savedir1};
savedet.help  = {'The Jacobian determinants may be saved to disk as a ``j_*.nii'''' file.'};

pullback      = cfg_branch;
pullback.name = 'Pullback';
pullback.tag  = 'pull';
pullback.val  = {applyto,savedir,interp,mask,fwhm,prefix};
pullback.help = {[...
'This is the old way of warping images, which involves resampling images based on a mapping from ',...
'the new (warped) image space back to the original image.  ',...
'The deformation should be the inverse of the deformation that would be used for the pushforward procedure.']};

weight        = cfg_files;
weight.name   = 'Weight Image';
weight.tag    = 'weight';
weight.filter = 'nifti';
weight.num    = [0 1];
weight.help   = {'Select an image file to weight the warped data with.  This is optional, but the idea is the same as was used by JE Lee et al (2009) in their ``A study of diffusion tensor imaging by tissue-specific, smoothing-compensated voxel-based analysis'''' paper.  In principle, a mask of (eg) white matter could be supplied, such that the warped images contain average signal intensities in WM.'};
weight.val    = {{''}};

% add note on aliasing to fwhm.help for Pushforward
fwhm.help     = {[fwhm.help{1} ' Note that you can specify [0 0 0], ',...
    'but any "modulated" data will show aliasing, which occurs because of ',...
    'the way the warped images are generated.']};

pushfo        = cfg_branch;
pushfo.name   = 'Pushforward';
pushfo.tag    = 'push';
pushfo.val    = {applyto,weight,savedir,deffov,preserve,fwhm,prefix};
pushfo.help   = {[...
'This is a newer way of warping images (for SPM at least), and involves the ',...
'forward pushing of voxel values from the original image into the appropriate place in the warped image. ',...
'The deformation field should be the inverse of the one used for the pullback procedure.'],...
'',...
[...
'``Smoothed'''' (blurred) spatially normalised images are generated in such a ',...
'way that the original signal is preserved. Normalised images are ',...
'generated by a ``pushing'''' rather than a ``pulling'''' (the usual) procedure. ',...
'Note that a procedure related to trilinear interpolation is used, and no masking is done.  It ',...
'is therefore recommended that the images are realigned and resliced ',...
'before they are spatially normalised, in order to benefit from motion correction using higher order interpolation.  Alternatively, contrast images ',...
'generated from unsmoothed native-space fMRI/PET data can be spatially ',...
'normalised for a 2nd level analysis.'],[...
'Two ``preserve'''' options are provided.  One of them should do the ',...
'equavalent of generating smoothed ``modulated'''' spatially normalised ',...
'images.  The other does the equivalent of smoothing the modulated ',...
'normalised fMRI/PET, and dividing by the smoothed Jacobian determinants.']};

pushsurf      = cfg_branch;
pushsurf.name = 'Surface';
pushsurf.tag  = 'surf';
pushsurf.val  = {surfa,savedir};
pushsurf.help = {[...
'Surfaces may be warped using the resulting deformation. ',...
'Note that a procedure similar to the pushforward is used, so the deformation should ',...
'be the inverse of the one that would be used for spatially normalising images via the pullback procedure.']};

output        = cfg_repeat;
output.name   = 'Output';
output.tag    = 'out';
output.values = {savedef,pullback, pushfo,pushsurf,savedet};
output.help = {[...
'Various output options are available.  ',...
'The deformation may be saved to disk as a ``y_*.nii'''' file.',...
'Images may be warped using the resulting deformation, either using a ``pullback'''' procedure, or a ``pushforward''''.',...
'The old style of spatial normalisation involved the pullback, whereas the pushforward requires ',...
'the inverse of the deformation used by the pullback.  ',...
'Finally, the deformation may be used to warp a GIFTI surface file.']};

conf         = exbranch('Deformations','defs',{comp,output});
conf.prog    = @spm_deformations;
conf.vout    = @vout;
conf.help    = hsummary;

%==========================================================================
function vo = vout(job)
vo = [];
savedef   = false;
saveimage = false;
savesurf  = false;
savejac   = false;
for i=1:numel(job.out)
    out = job.out{i};
    if isfield(out,'savedef') && ~savedef
        savedef = true;
        if isempty(vo), vo = cfg_dep; else vo(end+1) = cfg_dep; end
        vo(end).sname      = 'Deformation';
        vo(end).src_output = substruct('.','def');
        vo(end).tgt_spec   = cfg_findspec({{'filter','nifti'}});
    end
    if (isfield(out,'pull') || isfield(out,'push')) && ~saveimage
        saveimage = true;
        if isempty(vo), vo = cfg_dep; else vo(end+1) = cfg_dep; end
        vo(end).sname      = 'Warped Images';
        vo(end).src_output = substruct('.','warped');
        vo(end).tgt_spec   = cfg_findspec({{'filter','image'}});
    end
    if isfield(out,'surf') && ~savesurf
        savesurf = true;
        if isempty(vo), vo = cfg_dep; else vo(end+1) = cfg_dep; end
        vo(end).sname      = 'Warped Surfaces';
        vo(end).src_output = substruct('.','surf');
        vo(end).tgt_spec   = cfg_findspec({{'filter','mesh'}});
    end
    if isfield(out,'savejac') && ~savejac
        savejac = true;
        if isempty(vo), vo = cfg_dep; else vo(end+1) = cfg_dep; end
        vo(end).sname      = 'Jacobian';
        vo(end).src_output = substruct('.','jac');
        vo(end).tgt_spec   = cfg_findspec({{'filter','image'}});
    end
end

%==========================================================================
function entry_item = entry(name, tag, strtype, num)
entry_item         = cfg_entry;
entry_item.name    = name;
entry_item.tag     = tag;
entry_item.strtype = strtype;
entry_item.num     = num;

function files_item = files(name, tag, fltr, num)
files_item        = cfg_files;
files_item.name   = name;
files_item.tag    = tag;
files_item.filter = fltr;
files_item.num    = num;

function branch_item = branch(name, tag, val)
branch_item      = cfg_branch;
branch_item.name = name;
branch_item.tag  = tag;
branch_item.val  = val;

function exbranch_item = exbranch(name, tag, val)
exbranch_item      = cfg_exbranch;
exbranch_item.name = name;
exbranch_item.tag  = tag;
exbranch_item.val  = val;

function repeat_item = repeat(name, tag, values)
repeat_item        = cfg_repeat;
repeat_item.name   = name;
repeat_item.tag    = tag;
repeat_item.values = values;

function menu_item = mnu(name, tag, labels, values)
menu_item        = cfg_menu;
menu_item.name   = name;
menu_item.tag    = tag;
menu_item.labels = labels;
menu_item.values = values;

author: spm


text:
function [fixedstring]=fixstring(string)
%FIXSTRING utility function with several purposes:
% - clear possible initial spaces
% - turn possible multiple spaces into single spaces
% - replace irregular spaces (e.g. tab space) with regular spaces
% - suppress possible final spaces
% this guarantees good behaviour of the code

fixedstring = regexprep(strtrim(string), '[ ]{2,}',' ');

author: european-central-bank

repo: BEAR-toolbox


text:
function idx = findClosestCentroids(X, centroids)
%FINDCLOSESTCENTROIDS computes the centroid memberships for every example
%   idx = FINDCLOSESTCENTROIDS (X, centroids) returns the closest centroids
%   in idx for a dataset X where each row is a single example. idx = m x 1
%   vector of centroid assignments (i.e. each entry in range [1..K])
%

% Set K
K = size(centroids, 1);

% You need to return the following variables correctly.
idx = zeros(size(X,1), 1);

% ====================== YOUR CODE HERE ======================
% Instructions: Go over every example, find its closest centroid, and store
%               the index inside idx at the appropriate location.
%               Concretely, idx(i) should contain the index of the centroid
%               closest to example i. Hence, it should be a value in the
%               range 1..K
%
% Note: You can use a for-loop over the examples to compute this.
%
m = size(X,1);
for i = 1:m,
  minDistance = 10^6; % randomly chosen large enough value
  for j = 1:K,
    distance = sum((X(i,:) - centroids(j,:)).^2);
    if distance < minDistance,
      minDistance = distance;
      idx(i) = j;
    end
  end
end

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

end

author: anirudhjayaraman

repo: Machine-Learning


text:
function[blobSubY, blobSubX] = blobToImageSubs(blob)
% [blobSubY, blobSubX] = blobToImageSubs(blob)
%
% Convert a blob and an image to the indices of that blob in the image.
% These indices are relative to the size of the whole image. We don't need
% to know the image size for that!
% See "help blob" for more information.
%
% Copyright by Holger Caesar, 2014

% % For Matlab coder (not used)
% assert(isa(blob, 'struct'));
% assert(isa(blob.rect, 'double'));
% assert(all(size(blob.rect) == [1, 4]));
% assert(isa(blob.mask, 'logical'));
% assert(all(size(blob.mask) >= [1, 1]));
% assert(isa(blob.size, 'double'));
% assert(all(size(blob.size) == [1, 1]));

% Get all pix. coords for the mask
[blobSubY, blobSubX] = find(blob.mask);

% Make sure blobSub* are vectors (necessary when mask is a row vector)
blobSubY = blobSubY(:);
blobSubX = blobSubX(:);

% Translate from blob coords to image coords
blobSubY = blobSubY + blob.rect(1) - 1;
blobSubX = blobSubX + blob.rect(2) - 1;

author: nightrome

repo: matconvnet-calvin


text:
function pass = test_nonlinSysDampingBreaks_C2(pref)
% Test 2x2 system (sin/cos) where damping is required
%
% Asgeir Birkisson, May 2014.

if ( nargin == 0 )
    pref = cheboppref;
end

tol = 1e-10;

%% Piecewise:
d = [-pi 0 pi];
x = chebfun('x',d);
f = [ 0*x ; 0*x ];

%% CHEBCOLLOC2
pref.discretization = @chebcolloc2;

A = chebop(@(x,u,v) [u - diff(v,2); diff(u,2) + cos(v)], d);
A.lbc = @(u,v) [u-1/2; v+1/4];
A.rbc = @(u,v) [v-1/4; u+1/2];

[u1, u2, info] = solvebvp(A, f, pref);

% Want to check BCs as well.
bcFunLeft = chebfun(A.lbc(u1,u2));
bcFunRight = chebfun(A.rbc(u1,u2));

% And check that we're continuous over breakpoint
u1jump = jump(u1, 0);
u2jump = jump(u2, 0);

pass(1) = abs(info.error) < tol;
pass(2) = norm(bcFunLeft(d(1))) < tol && norm(bcFunRight(d(end))) < tol;
pass(3) = norm(u1jump) < tol && norm(u2jump) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function fMI = fusionMI(imgA,imgB,imgF)
% the overall mutual information between source images and fused image
%--------------------------------------------------------------------
[pAF, pA, pF] = estpab(imgA,imgF);
MIAF = estmutualinfo(pAF, pA, pF);

[pBF, pB, pF] = estpab(imgB,imgF);
MIBF = estmutualinfo(pBF, pB, pF);

fMI = MIAF + MIBF;
return

author: thfylsty


text:
function s = getAlgorithmsHeaderSize(this)
%
%
% I am sorry, but this is pathetic.
% BW.
%
s = 255 + 255 + 2 + 4;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function[sigma2, cv, mu, sig2inv, cvinv, S, lmu ] = estNoise_mask(I, w, mask, N, zoom);
%
%  Input parameters
%   I     : Input (noisy) image
%   w     : radius of analysis window
%   N     : Number of ellements to be considered (histogram) when
%           estimating sigma2 and cv2
%   zoom  : Conditional flag to improve the estimation of
%           sigma2 and cv2
%
%
%  Output parameters
%
%   sigma2  : estimation of the variance
%   cv2     : estimation of the coeficient of variation (cv2 = sigma^2 / mu^2)
%   mu      : estimation of the mean value
%   sig2inv : estimation of 1/sigma2
%   cvinv   : estimation of 1/cv2
%   S       : estimation of the local sample variance (same size as I)
%   lmu     : estimation of the local sample mean (same size as I)
%

% Based on
%   S. Aja-Fernandez, G. Vegas-Sanchez-Ferrero,
%   M. Martin-Fernandez, C. Alberola-Lopez
%   "Automatic noise estimation in images using local statistics.
%   Additive and multiplicative cases",
%   Image and Vision Computing 27 (2009) 756-770
%
% Legal:
%   estNoise.m is part of the INMD software
%   (http://sites.google.com/a/istec.net/prodrig/Home/sw). INMD is free
%   software, you can redistribute it and/or modify it under the terms of
%   the GNU General Public License (version 2).
%
% Authors
%   Paul Rodriguez    prodrig@pucp.edu.pe

if(nargin < 5)
  zoom = 0;
  if(nargin < 4)
    N = 100;
  end
end

  L=1;
  b = 0.5*( 1 - cos(pi*(0:2*L)/L) );

  [Nrows, Ncols, depth] = size(I);

for d = 1:depth,

  % Local noise estimation
  [S(:,:,d) lmu(:,:,d)] = localnoise_mask(I(:,:,d),w, mask);

  sigma2(d) = estimateMode(S(:,:,d), N, zoom, b, L, w, 0, 1);
  sig2inv(d) = estimateMode(1./S, N, zoom, b, L, w, 1, 1);

  mu(d) = estimateMode(lmu(:,:,d), N, zoom, b, L, w, 0, 1);

  cv(d) = estimateMode(S(:,:,d)./(lmu(:,:,d).*lmu(:,:,d)), N, zoom, b, L, w, 0, 0);
  cvinv(d) = estimateMode((lmu(:,:,d).*lmu(:,:,d))./S(:,:,d), N, zoom, b, L, w, 1, 0);

end

%    sigma2 = estimateMode(S, N, zoom, b, L, w, 0, 1);
%    sig2inv = estimateMode(1./S, N, zoom, b, L, w, 1, 1);
%
%    mu = estimateMode(lmu, N, zoom, b, L, w, 0, 1);
%
%    cv = estimateMode(S./(lmu.*lmu), N, zoom, b, L, w, 0, 0);
%    cvinv = estimateMode((lmu.*lmu)./S, N, zoom, b, L, w, 1, 0);

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

function[sigma2] = estimateMode(S, N, zoom, b, L, w, inv, factor)

v = S(:);

vmin = min(v);
vmax = max(v);
res = (vmax - vmin)/N;

[h bin] = hist(v, N);

if(inv == 1)
  tmp = conv(b, h) / (2*L+1);
  hs = tmp(L:N+L-1);
%    [dummy pos] = max(hs(end:-1:1));
  posvec = findmaxima( hs(end:-1:1) );
  if(length(posvec) >= 1)
    pos = posvec(1);
  else
    [dummy pos] = max(h(end:-1:1));
  end
else
  [dummy pos] = max(h(end:-1:1));
end

if(zoom)

  ind = find( abs( (v - bin(N-pos+1))  ) < 2*res );

  [h bin] = hist(v(ind), N);

  tmp = conv(b, h) / (2*L+1);
  hs = tmp(L:N+L-1);

  maxpos = findmaxima(hs);

  if( length(maxpos) >= 2 )
    sval = bin( maxpos(1) );
  else
    [dummy pos] = max(h(end:-1:1));
    sval = bin(N-pos+1);
  end

else
  sval = bin(N-pos+1);
end

if factor == 1,
  w_size = (2*w + 1)*(2*w + 1);
  sigma2 = (w_size - 3)/(w_size - 1)*sval;
else
  sigma2 = sval;
end

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

function minima = findminima(x)

minima = findmaxima(-x);

function maxima = findmaxima(x)

% Unwrap to vector
x = x(:);
% Identify whether signal is rising or falling
upordown = sign(diff(x));
% Find points where signal is rising before, falling after
maxflags = [upordown(1)<0; diff(upordown)<0; upordown(end)>0];
maxima   = find(maxflags);

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [mixparams, y, z, a] = mdnfwd(net, x)
%MDNFWD	Forward propagation through Mixture Density Network.
%
%	Description
%	 MIXPARAMS = MDNFWD(NET, X) takes a mixture density network data
%	structure NET and a matrix X of input vectors, and forward propagates
%	the inputs through the network to generate a structure MIXPARAMS
%	which contains the parameters of several mixture models.   Each row
%	of X represents one input vector and the corresponding row of the
%	matrices in MIXPARAMS  represents the parameters of a mixture model
%	for the conditional probability of target vectors given the input
%	vector.  This is not represented as an array of GMM structures to
%	improve the efficiency of MDN training.
%
%	The fields in MIXPARAMS are
%	  type = 'mdnmixes'
%	  ncentres = number of mixture components
%	  dimtarget = dimension of target space
%	  mixcoeffs = mixing coefficients
%	  centres = means of Gaussians: stored as one row per pattern
%	  covars = covariances of Gaussians
%	  nparams = number of parameters
%
%	[MIXPARAMS, Y, Z] = MDNFWD(NET, X) also generates a matrix Y of the
%	outputs of the MLP and a matrix Z of the hidden unit activations
%	where each row corresponds to one pattern.
%
%	[MIXPARAMS, Y, Z, A] = MLPFWD(NET, X) also returns a matrix A  giving
%	the summed inputs to each output unit, where each row  corresponds to
%	one pattern.
%
%	See also
%	MDN, MDN2GMM, MDNERR, MDNGRAD, MLPFWD
%

%	Copyright (c) Ian T Nabney (1996-2001)
%	David J Evans (1998)

% Check arguments for consistency
errstring = consist(net, 'mdn', x);
if ~isempty(errstring)
  error(errstring);
end

% Extract mlp and mixture model descriptors
mlpnet = net.mlp;
mixes = net.mdnmixes;

ncentres = mixes.ncentres;	% Number of components in mixture model
dim_target = mixes.dim_target;	% Dimension of targets
nparams = mixes.nparams;	% Number of parameters in mixture model

% Propagate forwards through MLP
[y, z, a] = mlpfwd(mlpnet, x);

% Compute the postion for each parameter in the whole
% matrix.  Used to define the mixparams structure
mixcoeff  = [1:1:ncentres];
centres   = [ncentres+1:1:(ncentres*(1+dim_target))];
variances = [(ncentres*(1+dim_target)+1):1:nparams];

% Convert output values into mixture model parameters

% Use softmax to calculate priors
% Prevent overflow and underflow: use same bounds as glmfwd
% Ensure that sum(exp(y), 2) does not overflow
maxcut = log(realmax) - log(ncentres);
% Ensure that exp(y) > 0
mincut = log(realmin);
temp = min(y(:,1:ncentres), maxcut);
temp = max(temp, mincut);
temp = exp(temp);
mixpriors = temp./(sum(temp, 2)*ones(1,ncentres));

% Centres are just copies of network outputs
mixcentres =  y(:,(ncentres+1):ncentres*(1+dim_target));

% Variances are exp of network outputs
mixwidths = exp(y(:,(ncentres*(1+dim_target)+1):nparams));

% Now build up all the mixture model weight vectors
ndata = size(x, 1);

% Return parameters
mixparams.type = mixes.type;
mixparams.ncentres = mixes.ncentres;
mixparams.dim_target = mixes.dim_target;
mixparams.nparams = mixes.nparams;

mixparams.mixcoeffs = mixpriors;
mixparams.centres   = mixcentres;
mixparams.covars      = mixwidths;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
% OtQueryOpen Open an Omnitrans query.
%
% Usage:
%    iret = OtQueryOpen(hquery)
%
% IMPORTANT: Always check the return value.
% The return value should be equal to 1, otherwise
% later calls to the OtQuery interface may cause application crashes.
%

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 2006-02-28 Created (Johan Meijdam, VORtech)

function iret = OtQueryOpen(hquery)
   if nargout < 1
      error('OtQueryOpen: always check the return value.');
   end
   iret = Ot('query.open', hquery);
end

author: Sable

repo: mcbench-benchmarks


text:
% CVX: Built-in operators and functions supported in CVX models.
%
%    The following operators and functions are included with MATLAB but
%    have either been verified to work properly within CVX models or have
%    been extended with additional code to do so. In many cases, typing
%       help cvx/<func>
%    where <func> is one of the names listed below will provide specific
%    help on the proper use of that item in CVX models---including any
%    restrictions imposed by the DCP and DGP rulesets.
%
%    For a list of new functions created specifically for CVX, type
%    "help cvx/functions".
%
% Computational operators:
%    plus (+), uplus (unary +), minus (-), uminus (unary -), times (.*),
%    mtimes (*), ldivide (.\), mldivide (\), rdivide (./), mrdivide (/),
%    power (.^), mpower (^), subsref/subsasgn/end (subscripting),
%    transpose (.'), ctranspose (')
% Relational operators:
%    eq (==), ge (>=), gt (>), le (<=), lt(<), ne (~=).
% Linear/affine functions:
%    blkdiag, cat, conj, conv, cumsum, diag, dot, find, flipdim, fliplr,
%    flipud, hankel, horzcat, imag, ipermute, kron, permute, polyval, real,
%    repmat, reshape, rot90, sparse, sum, toeplitz, tril, triu, vertcat
% Nonlinear functions:
%    abs, max, min, norm, prod, sqrt
% Query functions:
%    disp/display, end, isempty, isequal, length, isreal, ndims, nnz,
%    numel, size, spy

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
function Cb2n = ch_roty(theta)
% 3D初等旋转， theta为旋转角度，rad
Cb2n = [cos(theta), 0 sin(theta); 0 1 0; -sin(theta) 0 cos(theta)];

end

author: yandld


text:
classdef mme_branch_opf_acp < mp.mme_branch_opf_ac

%   MATPOWER
%   Copyright (c) 2021-2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%         name = 'branch';
%     end

    methods
        function obj = add_constraints(obj, mm, nm, dm, mpopt)
            %% call parent
            add_constraints@mp.mme_branch_opf_ac(obj, mm, nm, dm, mpopt);

            %% branch voltage angle difference limits
            [Aang, lang, uang, iang] = obj.ang_diff_params(...
                    dm, mpopt.opf.ignore_angle_lim);
            if length(iang)
                mm.add_lin_constraint('ang', Aang, lang, uang, {'Va'});
            end
            mm.userdata.ang_diff_constrained_branch_idx = iang;
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%% NUMERICAL LINEAR ALGEBRA WITH MATLAB
% This document looks at some important concepts in linear and shows how to
% solve some basic problems using MATLAB.

%% Matrix-Vector products
% Let's think about what a matrix-vector product actually does. It's more
% interesting that you think! Lets take a simple example. Consider the
% vector [1;0] multiplied by a matrix [1,1;1,1]
%%
figure('color',[1 1 1])
A=ones(2)
x=[1;0]
Ax=A*x
line([0 x(1)],[0 x(2)]), hold on, grid on
line([0,Ax(1)],[0,Ax(2)],'color',[1 0 0])
text(1.1,0,'x','fontweight','bold'), text(1.1,1.1,'A*x','fontweight','bold')
axis([-0.5 1.5 -0.5 1.5])
title('Action of a matrix A on a vector x')
%%
% We can see that the matrix multiplication resulted in a vector that was
% rotated by a certain angle and also stretched by a certain factor. This
% action (stretching and rotating) is extremely useful for switching from
% one co-ordinate system [x,y] to another [x', y']. Something that you will
% come across in special relativity and a number of other areas. To
% visualise what we mean by a co-ordinate transformation, consider the
% basis vectors x:=[1,0] and y:=[0,1], and a point p:=[0.2,0.1].
figure('color',[1 1 1])
xy=[1,0;0,1]
p=[0.5,0.3,0.1;0.4,0.2,0.0]
A=[-1 1;1 1]
xyp=A*xy; % Co-ordinate rotation step
Ap=A*p;
line([0 xy(1,1)],[0 xy(2,1)]), grid on, hold on
line([0,xy(1,2)],[0,xy(2,2)])
plot(p(1,:),p(2,:),'*','markersize',6,'color','b')
line([0, xyp(1,1)],[0 xyp(2,1)],'color',[1 0 0])
line([0, xyp(2,1)],[0,xyp(2,2)],'color',[1 0 0])
plot(Ap(1,:),Ap(2,:),'*','markersize',6,'color','r')
axis([-1.5 1.5 -1.5 1.5]), axis square
title('Changing the co-ordinate basis with matrix multiplication')

%% A rogue's gallery of useful matrices
% Matrices can be thought of as a 2-D array of data, or a collection of
% vectors. There are a number of important matrices that are extremely
% useful, and worth knowing about. We will use several of these later on in
% this lecture. Here are a few:
%%
% * Identity
% * Ones/Zeros
% * Upper/Lower Triangular
% * Symmetric
% * Vandermonde
% * Finite Difference
% * Rotation
%%
% The *Identity Matrix* or *I*
%%
I = eye(4)
%%
% A *matrix of ones or zeros*:
%%
O = ones(4)
Z = zeros(4)
%%
% *Upper and lower triangular matrices*. These are particularly important,
% as we shall see below, and consist of matrices that have zeros below, or
% above the _main diagonal_
%%
UT = triu(randn(4))
LT = tril(randn(4))
%%
% *Symmetric matrices* are matrices that have the same
% elements above he main diagonal as below the main diagonal:
%%
A = randn(4);
SYM = A*A'
%%
% The *Vandermonde matrix* is a matrix whose columns are polynomials in
% $$x$. This is useful in curve fitting and interpolation, for example, if
% I wanted to find the $$3^{rd}$ order interpolating polynomial for some
% data [x,y], then I would use the following matrix:
figure('color',[1 1 1])
x=linspace(0,1,5)';
IM = [x.^0, x.^1, x.^2, x.^3, x.^4]
plot(x,IM), axis([0 1 0 1.2]), grid on
title('First five powers of x')
%%
% Generally for interpolation, it's better to use a matrix whose columns
% are *orthogonal* polynomials, such as Legendre or Chebyshev, but this is
% beyond the remit of this course.
%%
% *Finite difference matrices* are extremely useful. They are matrices that
% approximate the derivative of a function on a grid. One way to calculate
% them is to fit an interpolating polynomial through the function, and then
% differentiate the polynomial to get what is called a finite difference
% stencil, while the other is to use a Taylor expansion to approximate
% the derivatives.
%%
% As an illustration, we will use the Taylor expansion.
%%
% <latex> {\Large
% $ f(x+\Delta x) = f(x) + \Delta x f'(x) + \frac{{\Delta x}^2f''(x)}{2!} +
% $}{\emph higher order terms} </latex>
%%
% Here the dash represents a derivative with respect to x. If we ignore
% terms above the first derivative, we see we can say that
%%
% <latex> {\Large
% $ f'(x) \approx \frac{f(x+ \Delta x) - f(x)}{\Delta x}$} </latex>
%%
% This is called the *forward difference* approximation. If we apply this
% formula to itself, we come up with an approximation for the second
% derivative of a function:
%%
% <latex>{\Large
% $ f''(x) \approx \frac{1}{\Delta x}\left(\frac{f(x+\Delta x+\Delta
% x)-f(x+\Delta x)-f(x+\Delta x)+f(x)}{\Delta x}\right)$}</latex>
%%
% or that, finally:
%%
% <latex> {\Large
% $ f''(x) \approx \frac{1}{\Delta x^2}\left(f(x+2\Delta x
% )-2f(x+\Delta x)+f(x)\right)$}</latex>
%%
% How to turn this into a matrix, however? If we imagine a function
% expressed on a grid of points $${x_1,x_2,\cdots,x_N}$ and replace the
% function evaluation $$f(x+\Delta x)$ with function evaluations at the
% grid points, then the second derivative can be expressed as the following
% matrix:
%%
N=5;
I=ones(N,1);
D2 = spdiags([I, -2*I, I],[-1:1],N,N);
full(D2)
%%
% Here we use the construct
%  spdiags
% to tell MATLAB that the matrix is sparse, i.e. is composed of mainly zero
% elements. In fact, this particular type of sparsity pattern is called
% *tridiagonal*. Consider the second row of this matrix multiplied by a
% vector which is the function evaluated on the grid of points $$x_i$,
% $$\left[f(x_0), f(x_1), f(x_2),\cdots,f(x_N)\right]^T$. We have the following
% result:
%%
% <latex>
% $ \left(\begin{array}{cccccc} 1&-2&1&0&\cdots&0 \end{array}\right)\times
% \left(\begin{array}{c} f(x_0) \\ f(x_1) \\ f(x_2) \\ \vdots \\
% f(x_N)\end{array}\right) = f(x_0)-2\times f(x_1)+f(x_2)$
% </latex>
%%
% That is, the result of multiplying the second row of the differentiation
% matrix by a vector representing our function on the grid, gives an
% approximation to the derivative of our function on the grid. As an
% example, let's calculate the second derivative of $$ f(x)=\exp^{\sin(x)}$
% on a grid of 21 points:
N=21; x=linspace(0,2*pi,N)'; % set up the grid
I=ones(N,1); dx=x(2)-x(1);
D2 = spdiags([I, -2*I, I],[-1:1],N,N)/dx^2; %2nd derivative matrix
f=@(x)(exp(sin(x))) % function
d2 = D2*f(x);       % numerical derivative
d2(1)=1; d2(N)=1;   % fix boundary points
figure('color',[1 1 1])
plot(x,f(x),x,d2), grid on, hold on
d2f=@(x)(cos(x).^2.*exp(sin(x))-sin(x).*exp(sin(x))) % actual second derivative
plot(x,d2f(x),'r-.')
title('Approximating a second derivative using a differentiation matrix');
legend('f(x)=exp^{sin(x)}','Approximate 2^{nd} deriv','Actual 2^{nd} deriv');
hold off
%% Systems of linear equations
% At its heart, the subject of linear algebra is concerned with how to
% solve the seemingly simple equation:
%
% $${\bf{A}}x=b$,
% where $${\bf{A}}$ is a matrix and $$x$ and $$b$ are vectors.
%%
%
% This might seem surprising, but this problem is right at the heart of
% almost all problems in scientific computing, from differential equations
% to eigenvalue problems to interpolation and curve fitting.

%%
% Of course, for a single scalar variable, x, we know what that means:
% $$ a\times x = y$
% from which we can deduce that
% $$ x = y/a$, for non-zero $$a$.
%%
% The question is how do we solve such a problem when we have more than one
% unknown variable, or, equivalently, when x is a vector? To investigate this
% question, let's consider the following electrical circuit:
%
[c,m]=imread('circuit2.gif','gif');

figure('color',[1 1 1]),
image(c), colormap(m), axis equal, axis off

%%
% We want to calculate the currents $$I_1,I_2,I_3$
%%
% Using Kirchoff's voltage laws we get the following three equations for
% the three unknown currents:
%%
% <latex>
% $\left\{\begin{array}{ccc}
% I_1 + 25(I_1-I_2)+50(I_1-I_3) & = & 10 \\
% 25(I_2-I_1) + 30I_2 + I_2 -I_3 & = & 0 \\
% 50(I_3-I_1)+I_3-I_2+55I_3 & = & 0 \end{array}\right . $
% </latex>
%%
% which, with a little re-arranging becomes:
%%
% <latex>
% $\left\{\begin{array}{ccc}
% 76I_1 - 25I_2-50I_3 & = & 10 \\
% -25I_1 +56I_2 -I_3 & = & 0 \\
% -50I_1-I_2+106I_3 & = & 0 \end{array}\right .$
% </latex>
%%
% Finally, we can gather all the coefficients into a matrix, and the
% unknowns into a vector, to write the single matrix linear equation:
%%
% <latex>
% $\left( {\begin{array}{ccc}
% 76 & -25 & -50 \\
% -25 & 56 & -1  \\
% -50& -1 & 106 \end{array}} \right)
% \times \left(\begin{array}{c}
% I_1\\
% I_2\\
% I_3 \end{array}\right) = \left( \begin{array}{c}
% 10 \\
% 0\\
% 0 \end{array}\right)$
% </latex>
%%
% In this way we have got the system of equations into the form that we
% described above, namely:
% $${\bf{A}}x=b$
%%
% Of course, for a small system (3 by 3 in this case) we can solve this by
% hand. For a larger system, however, we need a computer. The basic
% approach to solving such a system is an extension of that used by hand,
% namely we add and subtract multiples of one row to one another until we
% finally end up with an upper triangular matrix. This action
% of adding and subtrating multiples of one row to another is known as a
% basic row operation. Let us illustrate the process with MATLAB. Before we
% begin, we augment the right hand side vector to be an extra column in our
% matrix $$\bf A$, so it is now three rows by four columns.
A=[76 -25 -50 10;-25 56 -1 0; -50 -1 106 0]
rrefexample(A)
%%
% What is the advantage of reducing $${\bf A}$ to be upper triangular?
% Starting from the last row, we see that we now have an equation that
% depends only on $$I_3$, that is $$I_3 = 0.117$. The row above depends on
% $$I_3$ and $$I_2$ only, and we now know what $$I_3$ is! In this way, we
% work our way up through the rows substituting as we go, and in each row
% there is only a single unknown variable. This process is called _back
% substitution_. This entire process - that of reducing the matrix to a
% triangular form and solving via back substitution - is the way most
% computer packages solve these linear systems. In MATLAB, the shorthand
% way of solving such systems is with the backslash character ``\''
A=[76 -25 -50;-25 56 -1; -50 -1 106]
b =[10;0;0]
I= A\b
%%
% In the example above, we have performed a special case of what is a more
% general idea - that of transforming a matrix into a special form, so that
% the equations can be solved more easily. In general, the idea of
% transforming a matrix through a series of elementary row operations is a
% powerful one. May special types of _matrix factorisations_ exist, and two
% of the most important examples are:
%%
% *QR decomposition*
% In this case, a general matrix is written as the product of two matrices,
% one orthogonal and one upper triangular:
% $${\bf A = QR}$
% As the inverse of an orthogonal matrix is equal to its transpose, we find
% the solution to the general system $${\bf A}x=b$ can be calculated as:
%%
% <latex>
% $\begin{array}{l}
% {\bf A}x = b \\
% {\bf QR}x = b \\
% {\bf Q}({\bf R}x)=b \\
% {\bf R}x = {\bf Q^T} \times b \\
% x = {\bf R}\setminus{\bf Q^T} \times b \end{array}$
% </latex>
%%
% *LU Decomposition*
% In this case the matrix is decomposed into two, a lower and an upper
% triangular matrix. As both the matrices are triangular, the back (and
% forward) substitution process is very quick:
%%
% <latex>
% $\begin{array}{l}
% {\bf A}x = b \\
% {\bf LU}x = b \\
% {\bf L}({\bf U}x)=b \\
% {\bf L}x = {\bf U}\setminus b \\
% x = {\bf L}\setminus({\bf U}\setminus b) \end{array}$
% </latex>
%%
% As another example of using MATLAB to solve linear systems, let us
% perform a least squares fit on some experimental data. We have run an
% experiment to measure the distance fallen by an object in a given time.
% We know from the laws of motion that the distance is related to the
% initial velocity and position, and the acceleration through a quadratic
% relationship:
%%
% $$x(t) = x(0) + v(0)t + \frac{1}{2}gt^2$
%%
% We want to fit a quadratic line of best fit through the data, and hence
% work out a rough estimate for $$g$ the acceleration due to gravity.
load fall_data
t=data(:,1);
dist=data(:,2);
I=ones(size(t));
R=[I,t,t.^2];  % This is a VANDERMONDE matrix (see above)
coeffs = R\dist
g=2*coeffs(3)  % approximation for g (coeff is g/2)
plot(t,dist,'o',t,R*coeffs,'r'), grid on
title('Fitting a quadratic curve to experimental data')

%% Eigenvectors and Eigenvalues of a matrix
% As we have seen, the action of a mtrix on a vector is one of rotation and
% scaling. There are certain vectors, however, which in some sense
% _resonate_ with some fundamental property of the matrix, and do not get
% rotated. These special vectors (or directions in N-Dimensional space) are
% called *eigenvectors* and the scaling factor by which they are stretched
% are called the *eigenvalues*. Mathematically what we are saying is:
% $$\bf{A}x = \lambda x$
% Or visually we are saying:
evecPlot
%%
% Essentially, what we are saying is that (square) matrices have certain
% ``resonant modes'' or ``preferred'' co-ordinate directions, which we call
% _eigenvectors_ of the matrix. The eigenvectors are orthogonal to one
% another, and together form an orthogonal basis. If an $$N\times N$ matrix
% has N distinct eigenvalues, and correspondingly N orthogonal
% eigenvectors, then it is said to have *Rank* N, or to be *full rank*.

%%
% To make all this more concrete, let's consider an example from the
% mathematical theory of waves and vibration. In one dimension, the wave
% equation can be written as the following _partial differential equation_
% (don't worry about what this means, exactly, you'll learn more about all
% this later).
%%
% <latex>{\Large $\frac{\partial ^2 \Psi}{\partial t^2} = c\frac{\partial ^2
% \Psi}{\partial x^2}$}</latex>
%%
% We solve such a system by a process known as _separation of the
% variables_, and this involves writing the solution to the above equation
% as a product of two parts, one purely dependent on space and the other on
% time:
%%
% $${\Large\Psi\left(x,t\right) = \Theta\left(x\right)\Phi\left(t\right)}$
%%
% By substituting this form of the solution back into the differential
% equation, we get:
%%
% $$\ddot{\Phi}(t)\Theta(x)=c\Phi(t)\Theta''(x)$
%%
% or that
%%
% <latex> {\Large$\frac{\ddot{\Phi}(t)}{\Phi(t)} =
% c\frac{\Theta''(x)}{\Theta(x)} = -k^2 $}</latex>
%%
% The only way that the left hand side, which is completely independent of
% $$x$, can be equal to the right hand side, which is completely
% independent of $$t$, is if both sides are equal to a constant, $$-k^2$.
% We are interested in time independent solutions (standing waves) for this
% problem, so we can consider only the spatial parts of the equation, which
% can be written as:
%%
% <latex>
% {\Large$\frac{\partial ^2 \Theta(x)}{\partial x^2} = -k^2\Theta(x)
% $ }</latex>
%%
% This is an eigenvalue equation, exactly as we had above, and the
% solutions we are looking for are called *eigenfunctions*. This particular
% equation is important in physics and is known as the _Helmholtz
% equation_.
%%
% Using the second derivative differentiation matrix we described above
% (actually a more accurate version of it) and extending the notion to two
% dimensions, we transform a differential equation into a matrix eigenvalue
% equation. We can then use MATLAB to compute the eigenvalues and
% eigenvectors numerically to investigate the standing wave patterns, for
% example:
chladni

%   Copyright 2008-2009 The MathWorks, Inc.
%   $Revision: 35 $  $Date: 2009-05-29 15:27:34 +0100 (Fri, 29 May 2009) $


text:
function g = dnetOutputGrad(model, X)

% DNETOUTPUTGRAD Evaluate derivatives of dnet model outputs with respect to parameters.
% FORMAT
% DESC evaluates the derivates of a density network's
% outputs with respect to the parameters of the mapping function.
% ARG model : the model for which the derivatives are to be
% computed.
% ARG X : the input data locations where the gradients are to be
% computed.
% RETURN g : the gradient of the outputs of the density network
% perceptron with respect to each of the parameters. The size of
% the matrix is number of data x number of parameters x number of
% outputs of the model.
%
% SEEALSO : modelOutputGrad, dnetCreate
%
% COPYRIGHT : Neil D. Lawrence, 2008

% MLTOOLS

  g = modelOutputGrad(model.mapping, X);
  % Add output gradient wrt beta.
  g(:, end+1, :) = 0;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [y alpha eps eta] = sample_int(n, N, p, m, r, Znl, Tnl, Z, T, Hdyn, Zdyn, Tdyn, Rdyn, Qdyn, cdyn, Hmat, Rmat, Qmat, cmat, a1, P1)
% y is p*N*n
% alpha is m*N*n
% eps is p*N*n
% eta is r*N*n

%% Determine nonlinear functions %%
if Znl, Zdyn = false; else Zmat = getmat(Z); end
if Tnl, Tdyn = false; else Tmat = getmat(T); end

%% Draw from Gaussian distribution %%
alpha1  = randn(m, N);
eps     = randn(p, N, n);
eta     = randn(r, N, n);

%% Initialization for sampling %%
if Hdyn, eps(:,:, 1) = sigma_int(Hmat{1}, eps(:,:, 1));
else eps = reshape(sigma_int(Hmat, reshape(eps, p, N*n)), p, N, n); end
if Tdyn, T = Tmat{1}; elseif ~Tnl, T = Tmat; end
if Rdyn, R = Rmat{1}; else R = Rmat; end
if Qdyn, eta(:,:, 1) = sigma_int(Qmat{1}, eta(:,:, 1));
else eta = reshape(sigma_int(Qmat, reshape(eta, r, N*n)), r, N, n); end
if cdyn, c = repmat(cmat{1}, 1, N); else c = repmat(cmat, 1, N); end
P1(P1 == Inf) = 0;
y       = zeros(p, N, n);
alpha   = zeros(m, N, n);

%% Generate independent samples from the model %%
alpha(:,:, 1) = repmat(a1, 1, N) + sigma_int(P1, alpha1);
if Znl, y(:,:, 1) = getfunc(Z, alpha(:,:, 1), 1);
elseif Zdyn, y(:,:, 1) = Zmat{1}*alpha(:,:, 1);
end
for t = 2 : n
    if Hdyn, eps(:,:, t) = sigma_int(Hmat{t}, eps(:,:, t)); end
    if Qdyn, eta(:,:, t) = sigma_int(Qmat{t}, eta(:,:, t)); end
    if Tnl, alpha(:,:, t) = c + getfunc(T, alpha(:,:, t-1), t-1) + R*eta(:,:, t-1);
    else alpha(:,:, t) = c + T*alpha(:,:, t-1) + R*eta(:,:, t-1);
    end
    if Znl, y(:,:, t) = getfunc(Z, alpha(:,:, t), t); end
    if Zdyn, y(:,:, t) = Zmat{t}*alpha(:,:, t); end
    if Tdyn, T = Tmat{t}; end
    if Rdyn, R = Rmat{t}; end
    if cdyn, c = repmat(cmat{t}, 1, N); end
end
if ~Znl && ~Zdyn, y = reshape(Zmat*reshape(alpha, m, N*n), p, N, n) + eps;
else y = y + eps; end

%% Function for incorporating covariance into independent Gaussian samples %%
function x = sigma_int(Sigma, u)
dgSigma = diag(Sigma);
if isequal(Sigma, diag(dgSigma)), x = diag(sqrt(dgSigma))*u;
else % Sigma is not diagonal
    [U Lambda] = eig(full(Sigma));
    x = U*(diag(sqrt(diag(Lambda)))*u);
end

repo: REST


text:
% MOST
%   Version 1.2         13-Dec-2022
%
%   The MATPOWER Optimal Scheduling Tool (MOST) is framework for solving
%   generalized steady-state electric power scheduling problems.
%
%   MOST can be used to solve problems as simple as a deterministic,
%   single period economic dispatch problem with no transmission
%   constraints or as complex as a stochastic, security-constrained,
%   combined unit-commitment and multiperiod optimal power flow problem
%   with locational contingency and load-following reserves, ramping
%   costs and constraints, deferrable demands, lossy storage resources
%   and uncertain renewable generation.
%
%   While the problem formulation is general and incorporates a full
%   nonlinear AC network model, the current implementation is limited
%   to DC power flow modeling of the network. Some work has been done
%   on an AC implementation, but it is not yet ready for release.
%
%   The primary developers of MOST are Carlos E. Murillo-Sanchez and
%   Ray D. Zimmerman, with significant contributions from
%   Daniel Munoz-Alvarez and Alberto J. Lamadrid. It is built on top of
%   MATPOWER, a package of MATLAB/Octave M-files for solving power flow
%   and optimal power flow problems.
%
%   The latest version of MOST can be found on GitHub at:
%
%       https://github.com/MATPOWER/most
%
%   MOST is covered by the 3-clause BSD License (see LICENSE for details).

%   MOST
%   Copyright (c) 2010-2022, Power Systems Engineering Research Center (PSERC)
%   by Carlos E. Murillo-Sanchez, PSERC Cornell & Universidad Nacional de Colombia
%   and Ray Zimmerman, PSERC Cornell
%
%   This file is part of MOST.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/most for more info.


text:
%% rrt_star
% *RRT** is probabilistically optimal variant of RRT which converges
% to the optimal solution asymtotically. Due to the fact that RRT* takes
% cost of the path into consideration, the computational diffuculty
% increases with each node added to tree, a good example of this
% problem is finding neighbors.
%
%% Syntax
%  problem = rrt_star(map, max_iter, is_benchmark, rand_seed, variant)
%  function returns the object of the respective class with the result
%
%% Input
%  map           -- struct with appropriate fields (developer of
%                the class provides more information on this topic)
%  max_iter      -- number of iteration to solve the problem
%  is_benchmark  -- if true saves snapshots of the tree in a special directory
%                boolean variable
%  rand_seed     -- a random seed
%  variant       -- what class to choose, class used defines the problem space
%% Output
% *problem* is an object of an appropriate problem space representing class
%

author: olzhas


text:
function [U_final, V_final, dnorm, nIter_final, objhistory_final] = GNMF(X, k, W, options, U, V)
% Graph regularized Non-negative Matrix Factorization (GNMF)
%
% where
%   X
% Notation:
% X ... (mFea x nSmp) data matrix
%       mFea  ... number of words (vocabulary size)
%       nSmp  ... number of documents
% k ... number of hidden factors
% W ... weight matrix of the affinity graph
%
% options ... Structure holding all settings
%               options.alpha ... the regularization parameter.
%                                 [default: 100]
%                                 alpha = 0, GNMF boils down to the ordinary NMF.
%
%
% You only need to provide the above four inputs.
%
% X = U*V'
%
% References:
% [1] Deng Cai, Xiaofei He, Xiaoyun Wu, and Jiawei Han. "Non-negative
% Matrix Factorization on Manifold", Proc. 2008 Int. Conf. on Data Mining
% (ICDM'08), Pisa, Italy, Dec. 2008.
%
% [2] Deng Cai, Xiaofei He, Jiawei Han, Thomas Huang. "Graph Regularized
% Non-negative Matrix Factorization for Data Representation", IEEE
% Transactions on Pattern Analysis and Machine Intelligence, , Vol. 33, No.
% 8, pp. 1548-1560, 2011.
%
%
%   version 2.0 --April/2009
%   version 1.0 --April/2008
%
%   Written by Deng Cai (dengcai AT gmail.com)
%

if min(min(X)) < 0
    error('Input should be nonnegative!');
end

if ~isfield(options,'error')
    options.error = 1e-5;
end
if ~isfield(options, 'maxIter')
    options.maxIter = [];
end

if ~isfield(options,'nRepeat')
    options.nRepeat = 10;
end

if ~isfield(options,'minIter')
    options.minIter = 30;
end

if ~isfield(options,'meanFitRatio')
    options.meanFitRatio = 0.1;
end

if ~isfield(options,'alpha')
    options.alpha = 100;
end

nSmp = size(X,2);

if isfield(options,'alpha_nSmp') && options.alpha_nSmp
    options.alpha = options.alpha*nSmp;
end

if isfield(options,'weight') && strcmpi(options.weight,'NCW')
    feaSum = full(sum(X,2));
    D_half = X'*feaSum;
    X = X*spdiags(D_half.^-.5,0,nSmp,nSmp);
end

if ~isfield(options,'Optimization')
    options.Optimization = 'Multiplicative';
end

if ~exist('U','var')
    U = [];
    V = [];
end

[U_final, V_final, nIter_final, objhistory_final] = GNMF_Multi(X, k, W, options, U, V);
V_final = V_final';
dnorm = norm(X - U_final * V_final, 'fro');

author: jwyang


text:
function [feaNew] = SRDAtest(fea, model)
% SRDAtest: Spectral Regression Discriminant Analysis Testing
%               SRDAtest uses SRDA as a dimensionality reduction tool.
%
%       [feaNew,elapse] = SRDAtest(fea, model);
%
%             Input:
%
%               fea     - data matrix. Each row is a data point.
%             model     - model trained by SRKDAtrain.m
%
%             Output:
%
%             feaNew    - The data in the c-1 SRDA subspace, where c is the
%                         number of classes.
%
%    Examples:
%
%
% See also SRDAtrain, SRDApredict, SR, SR_caller
%
%Reference:
%
%   [1] Deng Cai, Xiaofei He and Jiawei Han, "SRDA: An Efficient Algorithm for
%   Large Scale Discriminant Analysis" IEEE Transactions on Knowledge and
%   Data Engineering, vol. 20, no. 1, pp. 1-12, January, 2008.
%
%   [2] Deng Cai, "Spectral Regression: A Regression Framework for
%   Efficient Regularized Subspace Learning", PhD Thesis, Department of
%   Computer Science, UIUC, 2009.
%
%
%   version 3.0 --Jan/2012
%   version 2.0 --December/2011
%   version 1.0 --May/2006
%
%   Written by Deng Cai (dengcai AT gmail.com)
%

if ~strcmp(model.TYPE,'SRDA')
    error('model does not match!');
end

if model.LARs
    feaNew = cell(length(model.LassoCardi),1);
    for i = 1:length(model.LassoCardi)
        feaNew{i} = fea*model.projection{i};
    end
else
    feaNew = fea*model.projection;
end

author: ZJULearning

repo: MatlabFunc


text:
function sim = gaussianKernel(x1, x2, sigma)
%RBFKERNEL returns a radial basis function kernel between x1 and x2
%   sim = gaussianKernel(x1, x2) returns a gaussian kernel between x1 and x2
%   and returns the value in sim

% Ensure that x1 and x2 are column vectors
x1 = x1(:); x2 = x2(:);

% You need to return the following variables correctly.
sim = 0;

% ====================== YOUR CODE HERE ======================
% Instructions: Fill in this function to return the similarity between x1
%               and x2 computed using a Gaussian kernel with bandwidth
%               sigma
%
%
temp = (x1-x2)'*(x1-x2);

sim = exp(-temp/(2*(sigma^2)));

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

end

repo: MachineLearningEx


text:
function [kPa] = dynpcm22kPa(dynpcm2)
% Convert pressure from dynes per square centimeter to kilopascals
% Chad Greene 2012
kPa = dynpcm2*0.000100000;


text:
function savenirfast(v,f,filestub, nodeseg, nodeprop, proptype)
%
% savenirfast(nirfaststruct,filestub)
%    or
% savenirfast(v,f,filestub, nodeseg, proptype, proptype)
%
% save a tetrahedral or surface mesh and associated properties to NIRFAST format
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%      nirfaststruct: a structure storing the NIRFAST mesh data, type
%         'help readnirfast' to read more; alternatively one can use:
%      v: input, node list, the first 3 columns are the x/y/z positions,
%         the remaining columns are combined with nodeprop as node-based
%         (optical) parameters
%      f: input, tetrahedral or surface element list, dimension (ne,3)
%      filestub: output file stub, output will include multiple files
%          filestub.node: node file
%          filestub.elem: element file to store the surface or tet mesh
%          filestub.param: parameter file
%          filestub.region: node label file
%      nodeseg: optional, an integer label field to group nodes into
%         segmentations, same length as v, number starting from 0; or empty
%      nodeprop: optional, additional nodal parameters, typically defined
%         as mua (1/mm), musp (1/mm) and refractive index (n)l; row number
%         equals to that of v, column number is user-defined
%      proptype: optional, the type of the node-property. by default it is
%         'stnd' - for standard properties; one can also define multi-row
%         header using a cell-array.
%
% example:
%    [node,face,elem]=meshabox([0 0 0],[10 10 10],0.3,1);
%    savenirfast(node,elem,'test', [], ones(size(node)), 'user');
%    mymesh=readnirfast('test')
%    plotmesh([mymesh.nodes mymesh.bndvtx], mymesh.elements,'x>5')
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(nargin<2)
    error('you must provide at least 2 inputs');
end

if(nargin==2)
    filestub=f;
    node=v.nodes;
    f=v.elements;
    proptype=v.type;
    if(isfield(v,'region'))
        nodeseg=v.region;
    end
    if(isfield(v,'mua'))
        nodeprop=[v.mua v.mus v.ri];
    end
    if(isfield(v,'bndvtx'))
        isboundary=v.bndvtx;
    end
    v=node;
end

if(size(v,2)>3)
    if(nargin>4)
        nodeprop=[v(:,4:end) nodeprop];
    else
        nodeprop=v(:,4:end);
    end
else
    if(nargin<5)
        nodeprop=[];
    end
end

if(nargin<6)
    proptype='stnd';
end

if(nargin<4 || isempty(nodeseg))
    nodeseg=zeros(size(v,1),1);
end

if(nargin<6)
    proptype='stnd';
end

if(size(f,2)>4)
    f(:,5:end)=[];
end

if(size(v,2)<3)
    error('v must contain at least 3 columns, and f must have at least 4 columns');
end

if(~exist('isboundary','var'))
    face=surfedge(f);
    isboundary=ismember(1:size(v,1), face(:));
end

fid=fopen([filestub,'.node'],'wt');
if(fid==-1)
    error('Saving node file failed, check permission or disk space.');
end
fprintf(fid,'%d\t%.16f\t%.16f\t%.16f\n',[isboundary(:) v(:,1:3)]');
fclose(fid);

if(size(f,2)<2 || size(f,2)>4)
    error('element list f must contain 3 or 4 columns');
end

fid=fopen([filestub,'.elem'],'wt');
if(fid==-1)
    error('Saving elem file failed, check permission or disk space.');
end
fprintf(fid,'%6d\t%6d\t%6d\t%6d\t\n',f');
fclose(fid);

if(~isempty(nodeseg))
    if(numel(nodeseg)~=size(v,1))
        error('nodeseg must have the same length as v');
    end
    fid=fopen([filestub,'.region'],'wt');
    if(fid==-1)
        error('Saving regin file failed, check permission or disk space.');
    end
    fprintf(fid,'%d\n',nodeseg(:));
    fclose(fid);
end

if(~isempty(nodeprop))
    if(size(nodeprop,1)~=size(v,1))
        error('nodeprop must have the same row number as v');
    end
    fid=fopen([filestub,'.param'],'wt');
    if(fid==-1)
        error('Saving param file failed, check permission or disk space.');
    end
    if(iscell(proptype))
        proptype=strjoin(proptype,'\n');
    end
    fprintf(fid,[proptype '\n']);
    fprintf(fid,[repmat('%.16f\t', 1, size(nodeprop,2)) '\n'],nodeprop');
    fclose(fid);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [stack] = fiff_read_mri(fname,read_data)
%
% [stack] = fiff_read_mri(fname,read_data)
%
% read_data argument is optional, if set to false the pixel data are
% not read. The default is to read the pixel data
%
% Read a fif format MRI description file
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%
%   Revision 1.5  2009/04/05 13:56:28  msh
%   Fixed a typo in the log messages
%
%   Revision 1.4  2009/04/05 13:50:08  msh
%   Added voxel -> MRI coordinate transform
%
%   Revision 1.3  2009/01/18 23:50:52  msh
%   Handle ushort embedded data correctly
%
%   Revision 1.2  2008/11/17 00:23:51  msh
%   Fixed error in the coordinate system transformation between voxel and RAS coordinates
%
%   Revision 1.1  2006/04/26 19:50:58  msh
%   Added fiff_read_mri
%

global FIFF;
if isempty(FIFF)
    FIFF = fiff_define_constants();
end

me='MNE:fiff_read_mri';

if nargin == 1
    read_data = true;
elseif nargin ~= 2
    error(me,'Incorrect number of arguments');
end
%
%   Try to open the file
%
[ fid, tree ] = fiff_open(fname);
%
%   Locate the data of interest
%   Pick the first MRI set within the first MRI data block
%
mri = fiff_dir_tree_find(tree,FIFF.FIFFB_MRI);
if length(mri) == 0
    fclose(fid);
    error(me,'Could not find MRI data');
end
mri = mri(1);
%
set = fiff_dir_tree_find(mri,FIFF.FIFFB_MRI_SET);
if length(set) == 0
    fclose(fid);
    error(me,'Could not find MRI stack');
end
set = set(1);
%
slices = fiff_dir_tree_find(set,FIFF.FIFFB_MRI_SLICE);
if length(slices) == 0
    fclose(fid);
    error(me,'Could not find MRI slices');
end
%
%   Data ID
%
tag = find_tag(mri,FIFF.FIFF_BLOCK_ID);
if ~isempty(tag)
    stack.id = tag.data;
end
%
%   Head -> MRI coordinate transformation
%   MRI (surface RAS) -> MRI (RAS) transformation
%
stack.trans = [];
stack.ras_trans = [];
for p = 1:set.nent
    kind = set.dir(p).kind;
    if kind == FIFF.FIFF_COORD_TRANS
        tag = fiff_read_tag(fid,set.dir(p).pos);
        if tag.data.to == FIFF.FIFFV_COORD_MRI && tag.data.from == FIFF.FIFFV_COORD_HEAD
            stack.trans = tag.data;
        elseif tag.data.to == FIFF.FIFFV_MNE_COORD_RAS && tag.data.from == FIFF.FIFFV_COORD_MRI
            stack.ras_trans = tag.data;
        end
    end
end
if isempty(stack.trans)
    stack.trans.from  = FIFF.FIFFV_COORD_HEAD;
    stack.trans.to    = FIFF.FIFFV_COORD_MRI;
    stack.trans.trans = eye(4,4);
end
if isempty(stack.ras_trans)
    stack.ras_trans.from   = FIFF.FIFFV_COORD_MRI;
    stack.ras_trans.to     = FIFF.FIFFV_MNE_COORD_RAS;
    stack.ras_trans.trans = eye(4,4);
end
stack.voxel_trans = [];
stack.nslice = length(slices);
if read_data
    fprintf(1,'\tReading slice information and pixel data.');
else
    fprintf(1,'\tReading slice information.');
end
for k = 1:stack.nslice
    try
        stack.slices(k) = read_slice(slices(k));
    catch
        fclose(fid);
        error(me,'%s',mne_omit_first_line(lasterr));
    end
    if mod(k,50) == 0
        fprintf(1,'.%d.',k);
    end
end
fprintf(1,'.%d..[done]\n',k);

add_voxel_transform

fclose(fid);

return;

    function add_voxel_transform

        if stack.nslice < 2
            fprintf(1,'\tOnly one slice in the stack. Voxel transformation will not be included\n');
            return;
        end
        %
        %   Check that the slices really form a stack
        %
        for k = 1:stack.nslice-1
            d = stack.slices(k+1).trans.trans(1:3,4)-stack.slices(k).trans.trans(1:3,4);
            d = sqrt(d'*d);
            if k == 1
                d0 = d;
                r0 = stack.slices(k).trans.trans(1:3,4);
                n0 = stack.slices(k).trans.trans(1:3,3);
            else
                if abs(stack.slices(k).pixel_width-stack.slices(1).pixel_width) > 1e-4 || abs(stack.slices(k).pixel_height-stack.slices(1).pixel_height) > 1e-4
                    fprintf(1,'\tPixel sizes are not equal. Voxel transformation will not be included\n');
                    return;
                end
                if stack.slices(k).width ~= stack.slices(1).width || stack.slices(k).height ~= stack.slices(1).height
                    fprintf(1,'\tImage sizes are not equal. Voxel transformation will not be included\n');
                    return;
                end
                if abs(d-d0) > 1e-4
                    fprintf(1,'\tThe slices are not equally spaced. Voxel transformation will not be included\n');
                    return;
                end
                %
                %   Rectangular volume?
                %
                r1 = r0 + (k-1)*d0*n0;
                d = stack.slices(k).trans.trans(1:3,4) - r1;
                d = sqrt(d'*d);
                if abs(d) > 1e-4
                    fprintf(1,'\tThe slices do not form a rectangular volume. Voxel transformation will not be included\n');
                    return;
                end
            end
        end
        %
        %   Ready to proceed
        %
        voxel_trans = stack.slices(1).trans;
        t = eye(4);
        t(1,1) = stack.slices(1).pixel_width;
        t(2,2) = stack.slices(1).pixel_height;
        d = stack.slices(2).trans.trans(1:3,4)-stack.slices(1).trans.trans(1:3,4);
        t(3,3) = sqrt(d'*d);
        voxel_trans.trans = voxel_trans.trans*t;
        voxel_trans.from  = FIFF.FIFFV_MNE_COORD_MRI_VOXEL;
        stack.voxel_trans = voxel_trans;

        fprintf(1,'\tVoxel transformation added\n');

    end

    function [slice] = read_slice(node)
        %
        %   Read all components of a single slice
        %
        tag = find_tag(node,FIFF.FIFF_COORD_TRANS);
        if isempty(tag)
            error(me,'Could not find slice coordinate transformation');
        end
        slice.trans = tag.data;
        if slice.trans.from ~= FIFF.FIFFV_COORD_MRI_SLICE || ...
                slice.trans.to ~= FIFF.FIFFV_COORD_MRI
            error(me,'Illegal slice coordinate transformation');
        end
        %
        %   Change the coordinate transformation so that
        %   ex is right
        %   ey is down (up in the fif file)
        %   ez steps to the next slice in the series
        %
        slice.trans.trans(1:3,2) = -slice.trans.trans(1:3,2);
        %
        %   Pixel data info
        %
        tag = find_tag(node,FIFF.FIFF_MRI_PIXEL_ENCODING);
        if isempty(tag)
            error(me,'Pixel encoding tag missing');
        end
        slice.encoding = tag.data;
        %
        %    Offset in the MRI data file if not embedded
        %
        tag = find_tag(node,FIFF.FIFF_MRI_PIXEL_DATA_OFFSET);
        if isempty(tag)
            slice.offset = -1;
        else
            slice.offset = tag.data;
        end
        %
        %    Format of the MRI source file
        %
        tag = find_tag(node,FIFF.FIFF_MRI_SOURCE_FORMAT);
        if isempty(tag)
            slice.source_format = 0;
        else
            slice.source_format = tag.data;
        end
        %
        %   Suggested scaling for the pixel values
        %   (not applied here)
        %
        tag = find_tag(node,FIFF.FIFF_MRI_PIXEL_SCALE);
        if isempty(tag)
            slice.scale = 1.0;
        else
            slice.scale = tag.data;
        end
        %
        %   Width and height in pixels
        %
        tag = find_tag(node,FIFF.FIFF_MRI_WIDTH);
        if isempty(tag)
            error(me,'Slice width missing');
        end
        slice.width = tag.data;
        %
        tag = find_tag(node,FIFF.FIFF_MRI_HEIGHT);
        if isempty(tag)
            error(me,'Slice height missing');
        end
        slice.height = tag.data;
        %
        %   Pixel sizes
        %
        tag = find_tag(node,FIFF.FIFF_MRI_WIDTH_M);
        if isempty(tag)
            error(me,'Pixel width missing');
        end
        slice.pixel_width = double(tag.data)/double(slice.width);
        %
        tag = find_tag(node,FIFF.FIFF_MRI_HEIGHT_M);
        if isempty(tag)
            error(me,'Pixel height missing');
        end
        slice.pixel_height = double(tag.data)/double(slice.height);
        %
        %   Are the data here or in another file?
        %
        tag = find_tag(node,FIFF.FIFF_MRI_SOURCE_PATH);
        if isempty(tag)
            slice.offset = -1;
            slice.source = [];
            %
            %   Pixel data are embedded in the fif file
            %
            if read_data
                tag = find_tag(node,FIFF.FIFF_MRI_PIXEL_DATA);
                if isempty(tag)
                    error(me,'Embedded pixel data missing');
                end
                if slice.encoding == FIFF.FIFFV_MRI_PIXEL_WORD
                    if tag.type ~= slice.encoding && tag.type ~= FIFF.FIFFT_USHORT
                        error(me,'Embedded data is in wrong format (expected %d, got %d)',...
                            slice.encoding,tag.type);
                    end
                else
                    if tag.type ~= slice.encoding
                        error(me,'Embedded data is in wrong format (expected %d, got %d)',...
                            slice.encoding,tag.type);
                    end
                end
                if length(tag.data) ~= slice.width*slice.height
                    error(me,'Wrong length of pixel data');
                end
                %
                %   Reshape into an image
                %
                slice.data = reshape(tag.data,slice.width,slice.height)';
            end
        else
            if slice.offset < 0
                error(me,'Offset to external file missing');
            end
            slice.source = tag.data;
            %
            %   External slice reading follows
            %
            if read_data
                pname = search_pixel_file(slice.source,fname);
                if isempty(pname)
                    error(me,'Could not locate pixel file %s',slice.source);
                else
                    try
                        slice.data   = read_external_pixels(pname,...
                            slice.offset,slice.encoding,slice.width,slice.height);
                    catch
                        error(me,'%s',mne_omit_first_line(lasterr));
                    end
                end
            end
        end
    end

    function [name] = search_pixel_file(pname,sname)
        %
        %   First try the file name as it is
        %
        if exist(pname,'file') == 2
            name = pname;
        else
            %
            %   Then <set file dir>/../slices/<slice file name>
            %
            a = findstr(sname,'/');
            if isempty(a)
                d = pwd;
            else
                d = sname(1:a(length(a))-1);
            end
            a = findstr(pname,'/');
            if ~isempty(a)
                pname = pname(a(length(a))+1:length(pname));
            end
            pname = sprintf('%s/../slices/%s',d,pname);
            if exist(pname,'file') == 2
                name = pname;
            else
                name = [];
            end
        end
        return;

    end

    function [pixels] = read_external_pixels(pname,offset,encoding,width,height)
        %
        %   Read pixel data from an external file
        %
        if (encoding == FIFF.FIFFV_MRI_PIXEL_SWAP_WORD)
            sfid = fopen(pname,'rb','ieee-le');
        else
            sfid = fopen(pname,'rb','ieee-be');
        end
        if sfid < 0
            error(me,'Could not open pixel data file : %s',pname);
        end
        try
            fseek(sfid,double(offset),'bof');
        catch
            error(me,'Could not position to pixel data @ %d',offset);
            fclose(sfid);
        end
        %
        %   Proceed carefully according to the encoding
        %
        switch encoding

            case FIFF.FIFFV_MRI_PIXEL_BYTE
                pixels = fread(sfid,double(width*height),'uint8=>uint8');
            case FIFF.FIFFV_MRI_PIXEL_WORD
                pixels = fread(sfid,double(width*height),'int16=>int16');
            case FIFF.FIFFV_MRI_PIXEL_SWAP_WORD
                pixels = fread(sfid,double(width*height),'int16=>int16');
            case FIFF.FIFFV_MRI_PIXEL_FLOAT
                pixels = fread(sfid,double(width*height),'single=>double');
            otherwise
                fclose(sfid);
                error(me,'Unknown pixel encoding : %d',encoding);
        end
        fclose(sfid);
        %
        %   Reshape into an image
        %
        pixels = reshape(pixels,width,height)';
    end

    function [tag] = find_tag(node,findkind)

        for p = 1:node.nent
            kind = node.dir(p).kind;
            pos  = node.dir(p).pos;
            if kind == findkind
                tag = fiff_read_tag(fid,pos);
                return;
            end
        end
        tag = [];
        return

    end

end


text:
function [ weights CompMat fuzzyTFN ] = FuzzyAHP( CompMat )
%FUZZYAHP Fuzzy AHP
%   Fuzzy AHP selection algorithm
%
% AUTHOR:
%           F. Ozgur CATAK
% CREATED:
%           October, 2011

% fuzzy tfn and inverse fuzzy tfn constants
fuzzyTFN = {[1     1     1  ] 	[1      1    1  ]
            [1/2   3/4   1  ] 	[1      4/3  2  ]
            [2/3   1     3/2] 	[2/3    1    3/2]
            [1     3/2   2  ] 	[1/2    2/3  1  ]
            [3/2   2     5/2] 	[2/5    1/2  2/3]
            [2     5/2   3  ] 	[1/3    2/5  1/2]
            [5/2   3     7/2] 	[2/7    1/3  2/5]
            [3     7/2   4  ] 	[1/4    2/7  1/3]
            [7/2   4     9/2] 	[2/9    1/4  2/7]};

fuzzyCompMatCell={};

%%
% convert ordinal numbers to
% triangular fuzzy number using fuzzyTFN matrix
[m n] = size(CompMat);

for i=1:m
    for j=i+1:m
       CompMat(j,i) = 1 / CompMat(i,j);
    end
end

for i=1:m
    for j=1:n
        criteria = CompMat(i,j);
        if criteria >= 1
           fuzzyCompMatCell{i,j} = fuzzyTFN{ criteria ,1 };
        else
           fuzzyCompMatCell{i,j} = fuzzyTFN{ round(criteria^-1) ,2 };
        end
    end
end

%%
% find sum of every l,m,u values for triangular fuzzy number
for i=1:m
    vec = [fuzzyCompMatCell{i,:}];
    mExtendAnalysis{1,i} = sum(reshape(vec,3,[])');
end

vec = [mExtendAnalysis{1,:}];
mExtendAnalysisSum = sum(reshape(vec,3,[])');

for i=1:m
    vec = [mExtendAnalysis{1,i}];
    for j=1:3
        val = mExtendAnalysisSum(1,j);
        %valSum(1,j) = val*vec(1,j);
        valSum(1,j) = (vec(1,j))*(1/val);
    end
    mExtendAnalysis{1,i} = valSum;
end

%%
% degree of possibility calculation
%              /---
%              | 1    if m2>=m1
%              |
%              | 0    if l1>=l2
% V(M2>=M1) = <
%              |     l1-u2
%              | --------------- otherwise
%              | (m1-u2)-(m1-l1)
%              \---
degreeOfPossibility = zeros(m*(m-1),3);
rowIndex = 1;
for i=1:m
    for j=1:m
        if i~=j
            degreeOfPossibility(rowIndex,[1 2]) = [i j];
            M1 = mExtendAnalysis{1,i};
            M2 = mExtendAnalysis{1,j};
            if M1(1,2) >= M2(1,2)
                degreeOfPossibility(rowIndex,3) = 1;
            elseif M2(1,1) >= M1(1,3)
                degreeOfPossibility(rowIndex,3) = 0;
            else
                degreeOfPossibility(rowIndex,3) = (M2(1,1)-M1(1,3))/((M1(1,2)-M1(1,3))-(M2(1,2)-M2(1,1)));
            end
            rowIndex = rowIndex + 1;
        end
    end
end
%%
% normalized weight calculation
weights = zeros(1,m);
for i=1:m,
    weights(1,i) = min(degreeOfPossibility([find(degreeOfPossibility(:,1) == i)], [3]));
end
weights = weights/sum(weights);

end


text:
function viewHotKeyAcc(figNo)
	% Allows hotkeys to be added to the current figure for manipulation
	% Biafra Ahanonu
	% started: 2017.04.12
	% inputs
	    %
	% outputs
	    %
	% changelog
	    %
	% TODO
	    % add options for how much to offset

	% allows hotkeys to be added to the current figure for manipulation
	% A figure
	% fig = figure('menubar','none');
	fig = figure(figNo);
	% Add menus with Accelerators
	mymenu = uimenu('Parent',fig,'Label','Hot Keys');
	uimenu('Parent',mymenu,'Label','Data Cursor','Accelerator','d','Callback',@(src,evt)datacursormode(fig,'on'));
	uimenu('Parent',mymenu,'Label','Zoom','Accelerator','f','Callback',@(src,evt)zoom(fig,'on'));
	uimenu('Parent',mymenu,'Label','Pan','Accelerator','g','Callback',@(src,evt)pan(fig,'on'));
	% Some plot
	%plot(rand(1,50))
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function field = field_composition(field, N)
% function field = field_composition(field, N)
%
% N compositions of the provided field.
%
% INPUT ARGUMENTS
% field 		- Field to compose
% N 			- Number of compositions to perform
%
% OPTIONAL INPUT ARGUMENTS
% N/A
%
% OUTPUT ARGUMENTS
% field 		- Composed field

% Copyright (c) 2012 Daniel Forsberg
% danne.forsberg@outlook.com
%
% 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/>.

for l = 1 : N
    for k = 1 : length(field)
        newField{k} = deformation(field{k},field,'linear');
        newField{k} = newField{k}+field{k};
    end
    field = newField;
end

author: fordanic

repo: image-registration


text:
function y=fst(x,a,b)
if x>b
    y=1;
elseif x<a
    y=0;
else
    y=(x-a)/(b-a);
end

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
function visualizeDiff(inputA, inputB, diff, hueA, hueB)
% images are inputA and inputB in rgb format or grayscale. Both has
% the same size with values in range [0-255].
% diff image contains any real values, zero means zero difference

numOfPixels = size(inputA,1)*size(inputA,2);

if (size(inputA,3) ~= 1) % image is in RGB format
    % convert it into grayscale
    inputA = rgb2gray(inputA);
end

if (size(inputB,3) ~= 1) % input B is in RGB format
    % convert it into grayscale
    inputB = rgb2gray(inputB)
end

% scale input values from range 0/255 into [0,1]
inputA = double(inputA)/255;
inputB = double(inputB)/255;

% if no output color is defined, use default
if (nargin < 4)
    hueA = pi/3;
    hueB = 4*pi/3;
end

% if difference is not in double values - convert it also
diff = double(diff);
% normalize values of difference into range [-pi/4, pi/4]
diff = diff/max(max(abs(diff)))*pi/4

% here magic begins
%  - Lsh (lightness, saturation, hue) conversion
%  - M represents length of color vector
M = double(reshape(inputA + inputB,numOfPixels,1))*50;
hue = reshape(diff > 0, numOfPixels, 1);
hue = (hue * hueA) - ((hue - 1) * hueB);
saturation = reshape(diff, numOfPixels, 1);

% from hue and saturation extract standard values of Lab color space
L = M .* cos(saturation);
a = M .* abs(sin(saturation)) .* cos(hue);
b = M .* abs(sin(saturation)) .* sin(hue);

% create output image in Lab color format
output(:,:,1) = reshape(L, size(inputA,1), size(inputA,2), 1);
output(:,:,2) = reshape(a, size(inputA,1), size(inputA,2), 1);
output(:,:,3) = reshape(b, size(inputA,1), size(inputA,2), 1);

% convert image into RGB
image = Lab2RGB(output);

% image preview
imshow(image);

end


text:
function [ pc_deg, pc ] = plymul2 ( add, mul, pa_deg, pa, pb_deg, pb )

%*****************************************************************************80
%
%% PLYMUL2 multiplies two polynomials in the field of order 2.
%
%  Discussion:
%
%    Polynomials stored as arrays have the coefficient of degree N in POLY(N+1).
%
%    A polynomial which is identically 0 is given degree -1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2003
%
%  Author:
%
%    Original FORTRAN77 version by Paul Bratley, Bennett Fox, Harald Niederreiter.
%    MATLAB version by John Burkardt
%
%  Parameters:
%
%    Input, integer ADD(2,2), MUL(2,2), the addition and multiplication
%    tables, mod 2.
%
%    Input, integer PA_DEG, the degree of polynomial A.
%
%    Input, integer PA(PA_DEG+1), the coefficients of polynomial A.
%
%    Input, integer PB_DEG, the degree of polynomial B.
%
%    Input, integer PB(PB_DEG+1), the coefficients of polynomial B.
%
%    Output, integer PC_DEG, the degree of polynomial C.
%
%    Output, integer PC(PC_DEG+1), the product polynomial, C = A * B mod 2.
%
  if ( pa_deg == -1 | pb_deg == -1 )
    pc_deg = -1;
  else
    pc_deg = pa_deg + pb_deg;
  end

  for i = 0 : pc_deg
    term = 0;
    for j = max ( 0, i-pa_deg ) : min ( pb_deg, i )
      term = add ( term+1, mul ( pa(i-j+1)+1, pb(j+1)+1 ) + 1 );
    end
    pt(i+1) = term;
  end

  for i = 0 : pc_deg
    pc(i+1) = pt(i+1);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function err = test_factorize (A)
%TEST_FACTORIZE test the accuracy of the factorization object
%
% Example
%   test_factorize (A) ;    % where A is a square matrix (sparse or dense)
%
% See also test_all, factorize1, factorize, inverse, mldivide

% Copyright 2009, Timothy A. Davis, University of Florida

if (nargin < 1)
    A = rand (100) ;
end

[m n] = size (A) ;
err = 0 ;
if (min (m,n) > 0)
    anorm = norm (A,1) ;
else
    anorm = 1 ;
end

for nrhs = 1:3

    for bsparse = 0:1

        % fprintf ('n %4d nrhs %d bsparse %d : ', n, nrhs, bsparse) ;

        b = rand (m,nrhs) ;
        if (bsparse)
            b = sparse (b) ;
        end

        %------------------------------------------------------------------
        % test backslash and related methods
        %------------------------------------------------------------------

        % method 0:
        x = A\b ;           err = check_resid (err, anorm, A, x, b) ;

        % method 1:
        S = inverse (A) ;
        x = S*b ;           err = check_resid (err, anorm, A, x, b) ;

        % method 2:
        if (m == n)
            F = factorize1 (A) ;
            x = F\b ;       err = check_resid (err, anorm, A, x, b) ;
        end

        % method 3:
        if (m == n)
            S = inverse (F) ;
            x = S*b ;       err = check_resid (err, anorm, A, x, b) ;
        end

        % method 4:
        F = factorize (A) ;
        x = F\b ;           err = check_resid (err, anorm, A, x, b) ;

        % method 5:
        S = inverse (F) ;
        x = S*b ;           err = check_resid (err, anorm, A, x, b) ;

        % method 6:
        if (m == n)
            [L,U,p] = lu (A, 'vector') ;
            x = U \ (L \ (b (p,:))) ;
            err = check_resid (err, anorm, A, x, b) ;
        end

        % method 7 (ack!)
        if (m == n)
            S = inv (A) ;
        else
            S = pinv (full (A)) ;
            if (issparse (A))
                S = sparse (S) ;
            end
        end
        x = S*b ;           err = check_resid (err, anorm, A, x, b) ;

        %------------------------------------------------------------------
        % test mtimes
        %------------------------------------------------------------------

        S = inverse (F) ;
        d = rand (n,1) ;
        x = F*d ;
        y = A*d ;
        z = S\d ;
        e = max (norm (x-y,1), norm (x-z,1)) ;
        if (e > 0)
            error ('mtimes error') ;
        end

        if (m == n)
            F = factorize1 (A) ;
            S = inverse (F) ;
            d = rand (n,1) ;
            x = F*d ;
            y = A*d ;
            z = S\d ;
            e = max (norm (x-y,1), norm (x-z,1)) ;
            if (e > 0)
                error ('mtimes error') ;
            end
        end

        %------------------------------------------------------------------
        % test slash and related methods
        %------------------------------------------------------------------

        b = rand (nrhs,n) ;
        if (bsparse)
            b = sparse (b) ;
        end

        % method 0:
        x = b/A ;           err = check_resid (err, anorm, A, x, b, 1) ;

        % method 1:
        S = inverse (A) ;
        x = b*S ;           err = check_resid (err, anorm, A, x, b, 1) ;

        % method 2:
        if (m == n)
            F = factorize1 (A) ;
            x = b/F ;       err = check_resid (err, anorm, A, x, b, 1) ;
        end

        % method 3:
        if (m == n)
            S = inverse (F) ;
            x = b*S ;       err = check_resid (err, anorm, A, x, b, 1) ;
        end

        % method 4:
        F = factorize (A) ;
        x = b/F ;           err = check_resid (err, anorm, A, x, b, 1) ;

        % method 5:
        S = inverse (F) ;
        x = b*S ;           err = check_resid (err, anorm, A, x, b, 1) ;

        % method 6:
        if (m == n)
            [L,U,p] = lu (A, 'vector') ;
            x = (b / U) / L ; x (:,p) = x ;
            err = check_resid (err, anorm, A, x, b, 1) ;
        end

        % method 7 (ack!)
        if (m == n)
            S = inv (A) ;
        else
            S = pinv (full (A)) ;
            if (issparse (A))
                S = sparse (S) ;
            end
        end
        x = b*S ;           err = check_resid (err, anorm, A, x, b, 1) ;

        %------------------------------------------------------------------
        % test double
        %------------------------------------------------------------------

        Y = double (inverse (A)) ;
        if (m == n)
            Z = inv (A) ;
        else
            Z = pinv (full (A)) ;
        end
        e = norm (Y-Z,1) ;
        if (n > 0)
            e = e / norm (Z,1) ;
        end
        err = max (e, err) ;

        %------------------------------------------------------------------
        % test subsref
        %------------------------------------------------------------------

        F = factorize (A) ;
        Y = inverse (A) ;
        if (numel (A) > 1)
            if (F (end) ~= A (end))
                error ('factorization subsref error') ;
            end
            if (F.A (end) ~= A (end))
                error ('factorization subsref error') ;
            end
        end
        if (n > 0)
            if (F (1,1) ~= A (1,1))
                error ('factorization subsref error') ;
            end
            if (F.A (1,1) ~= A (1,1))
                error ('factorization subsref error') ;
            end
            e = abs (Y (1,1) - Z (1,1)) ;
            err = max (e,err) ;
            if (m > 1 && n > 1)
                e = norm (Y (1:2,1:2) - Z (1:2,1:2), 1) ;
                err = max (e,err) ;
            end
        end
        if (m > 3 && n > 1)
            if (any (F (2:end, 1:2) - A (2:end, 1:2)))
                error ('factorization subsref error') ;
            end
            if (any (F (2:4, :) - A (2:4, :)))
                error ('factorization subsref error') ;
            end
            if (any (F (:, 1:2) - A (:, 1:2)))
                error ('factorization subsref error') ;
            end
        end

        %------------------------------------------------------------------
        % test update/downdate
        %------------------------------------------------------------------

        if (isa (F, 'factorization_dense_chol'))
            w = rand (n,1) ;
            b = rand (n,1) ;
            % update
            G = F + w ;
            x = G\b ;       err = check_resid (err, anorm, A+w*w', x, b) ;
            % downdate
            G = G - w ;
            x = G\b ;       err = check_resid (err, anorm, A, x, b) ;
            clear G
        end

        %------------------------------------------------------------------
        % test size
        %------------------------------------------------------------------

        [m1 n1] = size (F) ;
        [m n] = size (A) ;
        if (m1 ~= m || n1 ~= n)
            error ('size error') ;
        end
        [m1 n1] = size (Y) ;
        if (m1 ~= n || n1 ~= m)
            error ('pinv size error') ;
        end
        if (size (Y,1) ~= n || size (Y,2) ~= m)
            error ('pinv size error') ;
        end
        if (size (F,1) ~= m || size (F,2) ~= n)
            error ('size error') ;
        end

        %------------------------------------------------------------------
        % test mtimes
        %------------------------------------------------------------------

        d = rand (1,m) ;
        x = d*F ;
        y = d*A ;
        z = d/Y ;
        e = max (norm (x-y,1), norm (x-z,1)) ;
        if (e > 0)
            error ('mtimes error') ;
        end

        if (m == n)
            F = factorize1 (A) ;
            Y = inverse (F) ;
            d = rand (1,m) ;
            x = d*F ;
            y = d*A ;
            z = d/Y ;
            e = max (norm (x-y,1), norm (x-z,1)) ;
            if (e > 0)
                error ('mtimes error') ;
            end
        end

        %------------------------------------------------------------------
        % test inverse
        %------------------------------------------------------------------

        Y = double (inverse (inverse (A))) ;
        e = norm (A-Y,1) ;
        if (e > 0)
            error ('inverse error') ;
        end

    end
end

fprintf ('.') ;

if (err > 1e-8)
    fprintf ('error: %8.3e\n', err) ;
    error ('error is too high!') ;
end

%--------------------------------------------------------------------------

function err = check_resid (err, anorm, A, x, b, transposed)
if (nargin < 6)
    transposed = 0 ;
end
[m n] = size (A) ;

if (transposed)
    if (m >= n)
        e = norm (A'*x'-b',1) / (anorm + norm (x,1)) ;
    else
        e = norm (A*(A'*x')-A*b',1) / (anorm + norm (x,1)) ;
    end
else
    if (m <= n)
        e = norm (A*x-b,1) / (anorm + norm (x,1)) ;
    else
        e = norm (A'*(A*x)-A'*b,1) / (anorm + norm (x,1)) ;
    end
end

if (min (m,n) > 1)
    if (issparse (A) && issparse (b))
        if (~issparse (x))
            error ('x must be sparse') ;
        end
    else
        if (issparse (x))
            error ('x must be full') ;
        end
    end
end
err = max (err, e) ;


text:
function h5ButtonCallback(blockHandle)
% Callback function for H5button
% Allows user to specify what .h5 file to use for body

    [filename,filepath] = uigetfile('.h5');                % Return filename, filepath of user-chosen file
    % Don't set value if no file is chosen, or prompt canceled.
    if ~isequal(filename,0) && ~isequal(filepath,0)
        mask = Simulink.Mask.get(blockHandle);
        fileParam = mask.getParameter('h5File');
        fileParam.Value = [filepath,filename];             % Update h5File with new filename
    end
end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim