text:
function [valScaled, valFactor, valUnits] = bst_getunits( val, DataType, FileName, DisplayUnits)
% BST_GETUNITS: Get in which units is expressed a value.
%
% USAGE:  [valScaled, valFactor, valUnits] = bst_getunits(val, DataType, FileName=[], DisplayUnits=[]);
%
% INPUT:
%    - val       : Value to analyze
%    - DataType  : Type of data in the value "val". Possible strings:
%                  'EEG', 'MEG', 'MEG MAG', 'MEG GRAD', 'ECOG', 'SEEG', '$MEG', '$EEG', '$ECOG', '$SEEG', 'results', 'sources', 'source', 'stat', ...
%    - fUnits : Units of the file (eg "" if unknown, "cm", "pA", "mml/l")
%
% OUTPUT:
%    - valScaled : value in the detected units (val * valFactor)
%    - valFactor : factor to convert val -> valScaled
%    - valUnits  : string that represents the units ('\muV', 'fT', 'pA.m', etc.)

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

% Parse inputs
if (nargin < 4) || isempty(DisplayUnits)
    DisplayUnits = [];
end

% Check if there is something special in the filename
if (nargin >= 3) && ~isempty(FileName)
    % Source files
    if ismember(lower(DataType), {'results', 'sources', 'source'})
        % Detect sLORETA source files
        if ~isempty(strfind(lower(FileName), 'sloreta'));
            DataType = 'sloreta';
        % Detect NIRS source files
        elseif ~isempty(strfind(lower(FileName), 'nirs'))
            DataType = 'nirs-src';
        end
    end
end

% Consider input data in absolute value
val = abs(val);
% If no modality (ex: surface mask, mri values...)
if isempty(DataType)
    DataType = 'none';
end

% If the display unit is already defined
if ~isempty(DisplayUnits)
    if ismember(lower(DataType), {'nirs', '$nirs','nirs-src'})
        if ~isempty(strfind(DisplayUnits, 'mol'))
            [valFactor, valUnits] = GetSIFactor(val, DisplayUnits);
        elseif ~isempty(strfind(DisplayUnits, 'cm'))
            valFactor = 1;
            valUnits  = 'cm';
        elseif ~isempty(strfind(DisplayUnits, 'delta'))
            [valFactor, valUnits] = GetExponent(val);
            valUnits = sprintf('%s(%s)',strrep(DisplayUnits,'delta ','\Delta'),valUnits);
        else
            [valFactor, valUnits] = GetExponent(val);
            iParent =  find(DisplayUnits == '(');
            if ~isempty(iParent)
                DisplayUnits = DisplayUnits(1:iParent-1);
            end
            valUnits = sprintf('%s(%s)',DisplayUnits,valUnits);
        end
    else
        valUnits = DisplayUnits;
        switch (DisplayUnits)
            case 'fT'
                valFactor = 1e15;
            case 'pA.m'
                valFactor = 1e12;
            case '\muV'
                valFactor = 1e6;
            case {'mV', 'mm'}
                valFactor = 1e3;
            otherwise
                valFactor = 1;
        end
    end

% Otherwise: Units depends= on the data modality
else
    switch lower(DataType)
        case {'meg', '$meg', 'meg grad', 'meg mag', '$meg grad', '$meg mag', 'meg grad2', 'meg grad3', 'meg gradnorm'}
            % MEG data in fT
            if (val < 1e-8)
                valFactor = 1e15;
                valUnits  = 'fT';
            % MEG data without units (zscore, stat...)
            else
                valFactor = 1;
                valUnits  = 'No units';
            end

        case {'eeg', '$eeg', 'ecog', '$ecog', 'seeg', '$seeg', 'eog', '$eog', 'ecg', '$ecg', 'emg', '$emg'}
            % EEG data in Volts, displayed in microVolts
            if (val < 0.01)
                valFactor = 1e6;
                valUnits = '\muV';
            % EEG data in Volts, displayed in milliVolts
            elseif (val < 1)
                valFactor = 1e3;
                valUnits = 'mV';
            % EEG data without units (zscore, stat...)
            else
                valFactor = 1;
                valUnits = 'No units';
            end
        case {'nirs', '$nirs','nirs-src'}
             [valFactor, valUnits] = GetExponent(val);
        case {'results', 'sources', 'source'}
            % Results in Amper.meter (display in picoAmper.meter)
            if (val < 1e-4)
                valFactor = 1e12;
                valUnits  = 'pA.m';
            % Results without units (zscore, stat...)
            else
                valFactor = 1;
                valUnits  = 'No units';
            end

        case 'sloreta'
            if (val < 1e-4)
                [valFactor, valUnits] = GetExponent(val);
            else
                valFactor = 1;
                valUnits  = 'No units';
            end

        case 'stat'
            if (val < 1e-13)
                [valFactor, valUnits] = GetExponent(val);
            elseif (val < 1e-8)
                valFactor = 1e12;
                valUnits  = 'pA.m';
            elseif (val < 1e-3)
                [valFactor, valUnits] = GetExponent(val);
            else
                valFactor = 1;
                valUnits  = 'No units';
            end

        case 'connect'
            valFactor = 1;
            valUnits  = 'score';

        case 'timefreq'
            [valFactor, valUnits] = GetExponent(val);

        case 'hlu'
            valFactor = 1e3;
            valUnits  = 'mm';

        case {'ica', 'ssp'}
            [valFactor, valUnits] = GetExponent(val);
            if (valFactor == 0)
                valUnits = 'a.u.';
            else
                valUnits = ['x' valUnits ' a.u.'];
            end

        otherwise
            if isempty(val) || ((val < 1000) && (val > 0.1))
                valFactor = 1;
                valUnits  = 'No units';
            else
                [valFactor, valUnits] = GetExponent(val);
            end
    end
end
% Scale input value
valScaled = val .* valFactor;

end

%% ===== GET EXPONENT =====
function [valFactor, valUnits] = GetExponent(val)
    if (val == 0)
        valFactor = 1;
        valUnits  = 'No units';
    else
        %exponent = round(log(val)/log(10) / 3) * 3;
        exponent = round(log(val)/log(10)) - 1;
        % Do not allow 10^-1 and 10^-2
        if (exponent == -1) || (exponent == -2)
            valFactor = 1;
            valUnits  = 'No units';
        else
            valFactor = 10 ^ -exponent;
            valUnits  = sprintf('10^{%d}', exponent);
        end
    end
end

function [valFactor, valUnits] = GetSIFactor(val, originalUnit)
    vpw = [    -24,    -21,   -18,    -15,   -12,    -9,     -6,   -3, 0   +3,    +6,    +9,   +12,   +15,  +18,    +21,    +24];
    pfn = {'yocto','zepto','atto','femto','pico','nano','micro','milli','','kilo','mega','giga','tera','peta','exa','zetta','yotta'};
    pfs = {'y'    ,'z'    ,'a'   ,'f'    ,'p'   ,'n'   ,'\mu'    ,'m' ,'','k'   ,'M'   ,'G'   ,'T'   ,'P'   ,'E'  ,'Z'    ,'Y'    };
    sgf = 5;
    dpw = mode(diff(vpw));

    [unit, modifier] = getUnit(originalUnit);
    if abs(val) > 10^-2
        valFactor = 1;
        valUnits = originalUnit;
        return
    end

    adj = n2pAdjust(log10(abs(val)),dpw);

    vec = val./10.^adj;
    % Determine the number of decimal places:
    p10 = 10.^(sgf-1-floor(log10(abs(vec))));
    % Round coefficients to decimal places:
    vec = round(vec.*p10)./p10;
    % Identify which prefix is required:
    idx = 1+any(abs(vec)==[10.^dpw,1]);
    pwr = 1+floor(log10(abs(vec(idx))));

    % Obtain the required prefix index:
    idp = find(adj(idx)==vpw);

    valFactor = 10^(- vpw(idp));
    valUnits  = sprintf('%s%s',pfs{idp-1}, unit);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%num2sip
function adj = n2pAdjust(pwr,dPw)
adj = dPw*((0:1)+floor(floor(pwr)/dPw));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%n2pAdjust
function [unit, modifier] = getUnit(data)
% return the unit and the modifier from a string
% getUnit('mol/l') should return mol/l and 0
% getUnit('mmol/l') should return mol/l and -1
% getUnit('pA') should return A and -4

pfs = {'y'    ,'z'    ,'a'   ,'f'    ,'p'   ,'n'   ,'\mu'    ,'m' ,'','k'   ,'M'   ,'G'   ,'T'   ,'P'   ,'E'  ,'Z'    ,'Y'    };
vpw = [    -24,    -21,   -18,    -15,   -12,    -9,     -6,   -3, 0   +3,    +6,    +9,   +12,   +15,  +18,    +21,    +24];

Units = {'A','mol.l-1'};

unit = Units{cellfun(@(x)contains(data,x), Units)};
modifier = find(strcmp(strcat(pfs,unit),data)) - 9 ;

end

author: brainstorm-tools


text:
classdef BinaryDescriptorMatcher < handle
    %BINARYDESCRIPTORMATCHER  BinaryDescriptor matcher class
    %
    % Furnishes all functionalities for querying a dataset provided by user or
    % internal to class (that user must, anyway, populate) on the model of
    % Descriptor Matchers.
    %
    % Once descriptors have been extracted from an image (both they represent
    % lines and points), it becomes interesting to be able to match a
    % descriptor with another one extracted from a different image and
    % representing the same line or point, seen from a differente perspective
    % or on a different scale. In reaching such goal, the main headache is
    % designing an efficient search algorithm to associate a query descriptor
    % to one extracted from a dataset. In the following, a matching modality
    % based on *Multi-Index Hashing (MiHashing)* will be described.
    %
    % ### Multi-Index Hashing
    %
    % The theory described in this section is based on [MIH]. Given a dataset
    % populated with binary codes, each code is indexed `m` times into `m`
    % different hash tables, according to `m` substrings it has been divided
    % into. Thus, given a query code, all the entries close to it at least in
    % one substring are returned by search as *neighbor candidates*. Returned
    % entries are then checked for validity by verifying that their full codes
    % are not distant (in Hamming space) more than `r` bits from query code.
    % In details, each binary code `h` composed of `b` bits is divided into
    % `m` disjoint substrings `h^(1), ..., h^(m)`, each with length
    % `floor(b/m)` or `ceil(b/m)` bits. Formally, when two codes `h` and `g`
    % differ by at the most `r` bits, in at the least one of their `m`
    % substrings they differ by at the most `floor(r/m)` bits. In particular,
    % when `||h-g||_H <= r` (where `||.||_H` is the Hamming norm), there must
    % exist a substring `k` (with `1 <= k <= m`) such that:
    %
    %     || h^(k) - g^(k) ||_H <= floor(r/m)
    %
    % That means that if Hamming distance between each of the `m` substring is
    % strictly greater than `floor(r/m)`, then `||h-g||_H` must be larger that
    % `r` and that is a contradiction. If the codes in dataset are divided
    % into `m` substrings, then `m` tables will be built. Given a query `q`
    % with substrings `{q^(i)}_[i=1..m]`, `i`-th hash table is searched for
    % entries distant at the most `floor(r/m)` from `q^(i)` and a set of
    % candidates `N_i(q)` is obtained. The union of sets `N(q) = U_i {N_i(q)}`
    % is a superset of the `r`-neighbors of `q`. Then, last step of algorithm
    % is computing the Hamming distance between `q` and each element in `N(q)`,
    % deleting the codes that are distant more that `r` from `q`.
    %
    % ## References
    % [MIH]:
    % > Mohammad Norouzi, Ali Punjani, and David J Fleet. "Fast search in
    % > hamming space with Multi-Index Hashing". In Computer Vision and
    % > Pattern Recognition (CVPR), 2012 IEEE Conference on, pages 3108-3115.
    % > IEEE, 2012.
    %
    % See also: cv.BinaryDescriptor, cv.DescriptorMatcher, cv.drawLineMatches
    %

    properties (SetAccess = private)
        % Object ID
        id
    end

    methods
        function this = BinaryDescriptorMatcher()
            %BINARYDESCRIPTORMATCHER  Constructor
            %
            %     matcher = cv.BinaryDescriptorMatcher()
            %
            % The BinaryDescriptorMatcher constructed is able to store and
            % manage 256-bits long entries.
            %
            % See also: cv.BinaryDescriptorMatcher
            %
            this.id = BinaryDescriptorMatcher_(0, 'new');
        end

        function delete(this)
            %DELETE  Destructor
            %
            %     matcher.delete()
            %
            % See also: cv.BinaryDescriptorMatcher
            %
            if isempty(this.id), return; end
            BinaryDescriptorMatcher_(this.id, 'delete');
        end
    end

    %% Algorithm
    methods
        function clear(this)
            %CLEAR  Clear dataset and internal data
            %
            %     matcher.clear()
            %
            % See also: cv.BinaryDescriptorMatcher.empty
            %
            BinaryDescriptorMatcher_(this.id, 'clear');
        end

        function status = empty(this)
            %EMPTY  Returns true if there are no train descriptors in the collection
            %
            %     status = matcher.empty()
            %
            % ## Output
            % * __status__ boolean status
            %
            % See also: cv.BinaryDescriptorMatcher.clear
            %
            status = BinaryDescriptorMatcher_(this.id, 'empty');
        end

        function save(this, filename)
            %SAVE  Saves the algorithm parameters to a file
            %
            %     matcher.save(filename)
            %
            % ## Input
            % * __filename__ Name of the file to save to.
            %
            % This method stores the algorithm parameters in the specified
            % XML or YAML file.
            %
            % See also: cv.BinaryDescriptorMatcher.load
            %
            BinaryDescriptorMatcher_(this.id, 'save', filename);
        end

        function load(this, fname_or_str, varargin)
            %LOAD  Loads algorithm from a file or a string
            %
            %     matcher.load(fname)
            %     matcher.load(str, 'FromString',true)
            %     matcher.load(..., 'OptionName',optionValue, ...)
            %
            % ## Input
            % * __fname__ Name of the file to read.
            % * __str__ String containing the serialized model you want to
            %   load.
            %
            % ## Options
            % * __ObjName__ The optional name of the node to read (if empty,
            %   the first top-level node will be used). default empty
            % * __FromString__ Logical flag to indicate whether the input is a
            %   filename or a string containing the serialized model.
            %   default false
            %
            % This method reads algorithm parameters from the specified XML or
            % YAML file (either from disk or serialized string). The previous
            % algorithm state is discarded.
            %
            % See also: cv.BinaryDescriptorMatcher.save
            %
            BinaryDescriptorMatcher_(this.id, 'load', fname_or_str, varargin{:});
        end

        function name = getDefaultName(this)
            %GETDEFAULTNAME  Returns the algorithm string identifier
            %
            %     name = matcher.getDefaultName()
            %
            % ## Output
            % * __name__ This string is used as top level XML/YML node tag
            %   when the object is saved to a file or string.
            %
            % See also: cv.BinaryDescriptorMatcher.save, cv.BinaryDescriptorMatcher.load
            %
            name = BinaryDescriptorMatcher_(this.id, 'getDefaultName');
        end
    end

    %% BinaryDescriptorMatcher
    methods
        function add(this, descriptors)
            %ADD  Store locally new descriptors to be inserted in dataset, without updating dataset
            %
            %     matcher.add(descriptors)
            %
            % ## Input
            % * __descriptors__ cell array of matrices containing descriptors
            %   to be inserted into dataset. Each matrix `descriptors{i}`
            %   should contain descriptors relative to lines extracted from
            %   i-th image.
            %
            % See also: cv.BinaryDescriptorMatcher.clear
            %
            BinaryDescriptorMatcher_(this.id, 'add', descriptors);
        end

        function train(this)
            %TRAIN  Update dataset by inserting into it all descriptors that were stored locally by add function
            %
            %     matcher.train()
            %
            % NOTE: Every time this function is invoked, current dataset is
            % deleted and locally stored descriptors are inserted into
            % dataset. The locally stored copy of just inserted descriptors is
            % then removed.
            %
            % See also: cv.BinaryDescriptorMatcher.add
            %
            BinaryDescriptorMatcher_(this.id, 'train');
        end

        function matches = match(this, queryDescriptors, varargin)
            %MATCH  For every input query descriptor, retrieve the best matching one from a dataset provided from user or from the one internal to class
            %
            %     matches = matcher.match(queryDescriptors, trainDescriptors)
            %     matches = matcher.match(queryDescriptors)
            %     [...] = matcher.match(..., 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __queryDescriptors__ query descriptors.
            % * __trainDescriptors__ dataset of descriptors furnished by user.
            %
            % ## Output
            % * __matches__ vector to host retrieved matches.
            %
            % ## Options
            % * __Mask__ mask to select which input descriptors must be
            %   matched to one in dataset. In the second variant, a vector of
            %   masks (the i-th mask in vector indicates whether each input
            %   query can be matched with descriptors in dataset relative to
            %   i-th image). Not set by default.
            %
            % For every input descriptor, find the best matching one:
            %
            % - in the first variant, for a pair of images
            % - in the second variant, from one image to a set
            %
            % See also: cv.BinaryDescriptorMatcher.knnMatch,
            %  cv.BinaryDescriptorMatcher.radiusMatch
            %
            matches = BinaryDescriptorMatcher_(this.id, 'match', queryDescriptors, varargin{:});
        end

        function matches = knnMatch(this, queryDescriptors, varargin)
            %KNNMATCH  For every input query descriptor, retrieve the best k matching ones from a dataset provided from user or from the one internal to class
            %
            %     matches = matcher.knnMatch(queryDescriptors, trainDescriptors, k)
            %     matches = matcher.knnMatch(queryDescriptors, k)
            %     [...] = matcher.knnMatch(..., 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __queryDescriptors__ query descriptors.
            % * __trainDescriptors__ dataset of descriptors furnished by user.
            % * __k__ number of the closest descriptors to be returned for
            %   every input query.
            %
            % ## Output
            % * __matches__ vector to host retrieved matches.
            %
            % ## Options
            % * __Mask__ mask to select which input descriptors must be
            %   matched to ones in dataset. A vector of masks in the second
            %   variant (the i-th mask in vector indicates whether each input
            %   query can be matched with descriptors in dataset relative to
            %   i-th image). Not set by default.
            % * __CompactResult__ flag to obtain a compact result (if true, a
            %   vector that doesn't contain any matches for a given query is
            %   not inserted in final result). default false
            %
            % For every input descriptor, find the best k matching descriptors:
            %
            % - in the first variant, for a pair of images
            % - in the second variant, from one image to a set
            %
            % See also: cv.BinaryDescriptorMatcher.match,
            %  cv.BinaryDescriptorMatcher.radiusMatch
            %
            matches = BinaryDescriptorMatcher_(this.id, 'knnMatch', ...
                queryDescriptors, varargin{:});
        end

        function matches = radiusMatch(this, queryDescriptors, varargin)
            %RADIUSMATCH  For every input query descriptor, retrieve, from a dataset provided from user or from the one internal to class, all the descriptors that are not further than maxDist from input query
            %
            %     matches = matcher.radiusMatch(queryDescriptors, trainDescriptors, maxDistance)
            %     matches = matcher.radiusMatch(queryDescriptors, maxDistance)
            %     [...] = matcher.radiusMatch(..., 'OptionName', optionValue, ...)
            %
            % ## Input
            % * __queryDescriptors__ query descriptors.
            % * __trainDescriptors__ dataset of descriptors furnished by user.
            % * __maxDistance__ search radius.
            %
            % ## Output
            % * __matches__ vector to host retrieved matches.
            %
            % ## Options
            % * __Mask__ mask to select which input descriptors must be
            %   matched to ones in dataset. A vector of masks in the second
            %   variant (the i-th mask in vector indicates whether each input
            %   query can be matched with descriptors in dataset relative to
            %   i-th image). Not set by default.
            % * __CompactResult__ flag to obtain a compact result (if true, a
            %   vector that doesn't contain any matches for a given query is
            %   not inserted in final result). default false
            %
            % For every input desciptor, find all the ones falling in a
            % certaing matching radius:
            %
            % - in the first variant, for a pair of images
            % - in the second variant, from one image to a set
            %
            % See also: cv.BinaryDescriptorMatcher.match,
            %  cv.BinaryDescriptorMatcher.knnMatch
            %
            matches = BinaryDescriptorMatcher_(this.id, 'radiusMatch',...
                queryDescriptors, varargin{:});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Add local code directories to Matlab path
% addpaths;

imdir='../Images_resized/'; % directory with original images

workspcdir='../tempworkspace/'; % directory to save intermediate results
if ~exist(workspcdir,'dir')
    mkdir(workspcdir);
end

resdir='../Results/'; % This is where we will save final results using this demo script.
if ~exist(resdir,'dir')
    mkdir(resdir);
end

% You can run it on a single image as follows
imagename='indoor_0268.jpg';
[ boxlayout,surface_labels ] = getspatiallayout(imdir,imagename,workspcdir);

author: zouchuhang

repo: LayoutNet


text:
% Fig. 6.23   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%

clear all;
close all;

num=1;
den=[1 2 1];
w=logspace(-2,3,100);
[re,im]=nyquist(num,den,w);

plot(re,im,re,-im);
grid;
axis equal
xlabel('Re(G(s))');
ylabel('Im(G(s))');
title('Fig. 6.23 Nyquist plot');

author: Sable

repo: mcbench-benchmarks


text:
 function [T, reuse] = Gdsft_gram(ob, W, reuse)
%function [T, reuse] = Gdsft_gram(A, W, reuse)
%|
%| Construct Toeplitz gram matrix object T = A'WA for a Gdsft object.
%| This object performs T*x rapidly using an over-sampled FFT.
%|
%| in
%|	A	[M np]		Gdsft object (fatrix2)
%|	W	[M M]		W = Gdiag() for fatrix2 (often simply "1" or [])
%|				assume W=diag(wi) with real wi so T is Hermitian
%|	reuse	struct		stuff from the last call that can be reused
%|
%| out
%|	T	[np np]		fatrix2 object
%|	reuse	struct		stuff that can be reused on future calls
%|
%| Copyright 2012-06-05, Jeff Fessler, University of Michigan

if nargin == 1 && streq(ob, 'test'), Gdsft_gram_test, return, end
if nargin == 1 && streq(ob, 'test1'), Gdsft_gram_test1, return, end
if nargin == 1 && streq(ob, 'test2'), Gdsft_gram_test2, return, end
if nargin == 1 && streq(ob, 'test3'), Gdsft_gram_test3, return, end
if nargin < 2 || nargin > 3, ir_usage, end

if ~isvar('reuse'), reuse = []; end

if ~isa(ob, 'fatrix2')
	fail('A wrong class: %s', class(ob))
end

[nd np] = size(ob);
forw = @(arg, x) Gdsft_gram_mult(x, arg.fftkern, arg.mask);

% extract wi
if isnumeric(W) % includes case where W is empty
	wi = W;
elseif isa(W, 'fatrix2') && streq(W.caller, 'Gdiag')
	wi = W.arg.diag;
else
	error 'W must be diag_sp or Gdiag or wi array'
end
clear W

if isempty(wi)
	wi = ones(nd, 1); % the usual unweighted case
end

if isscalar(wi)
	wi = wi * ones(nd, 1);
end

if isreal(wi)
	back = forw; % trick: because Hermitian!
else
	fail('not implemented for complex wi')
end

arg = ob.arg;
arg.mask = ob.imask_array;
[arg.fftkern arg.reuse] = Gdsft_gram_init(arg, wi, reuse, ob.idim);

if any(~arg.mask(:))
	omask = arg.mask;
else
	omask = []; % trick: omask not needed in full case
end
T = fatrix2('arg', arg', ...
	'idim', ob.idim, 'imask', arg.mask, ...
	'odim', ob.idim, 'omask', omask, ...
	'forw', forw, 'back', back);

% Gdsft_gram_init()
% construct kernel of circulant matrix into which T is embedded
% and take its DFT to prepare for multiplication.
% out
%	fftkern [[2Nd]]		FFT of kernel of circulant matrix
%
function [fftkern, reuse] = Gdsft_gram_init(arg, wi, reuse, Nd)

switch numel(Nd)
case 1
	[fftkern reuse] = Gdsft_gram_init1(arg, wi, reuse, Nd, arg.show);
case 2
	[fftkern reuse] = Gdsft_gram_init2(arg, wi, reuse, Nd, arg.show);
case 3
	[fftkern reuse] = Gdsft_gram_init3(arg, wi, reuse, Nd, arg.show);
otherwise
	fail('dim %d not done', numel(Nd))
end

% Gdsft_gram_init1()
% 1d filter for circulant matrix
% note: only toeplitz kernel values from -(N-1) to (N-1) are relevant
% so the value at +/- N does not matter so we set it to zero.
function [fftkern, reuse] = Gdsft_gram_init1(arg, wi, reuse, N1, show)

if isempty(reuse)
	reuse.G1 = Gdsft(arg.om_t', arg.Nd); % with n_shift = 0
end

block1 = reuse.G1' * real(wi); % kludge

% kernel of Toeplitz matrix from -N to N-1 but with fftshift
% This is inherently Hermitian symmetric except for the middle [0] value.
% Use 0 for the irrelevant value at -N.
err1 = abs(imag(block1(1))) / abs(block1(1));
tol = 0;
if err1 > tol
	printm('removing imaginary h[0] part of relative size %g', err1)
	block1(1) = real(block1(1));
end
kern = [block1; 0; flipud(conj(block1(2:N1)))]; % [2*N1]
if show
	kern = dsft_gram_hermitify(kern, show); % need not due to block1(1) fix
end
fftkern = fft(kern); % [2*N1]

% Gdsft_gram_init2()
function [fftkern, reuse] = Gdsft_gram_init2(arg, wi, reuse, Nd, show)

if isempty(reuse)
	[reuse.G1 reuse.G2] = Gdsft_gram_setup2(arg);
end

N1 = Nd(1);
N2 = Nd(2);

block1 = reshape(reuse.G1' * real(wi), Nd); % kludge
block2 = reshape(reuse.G2' * real(wi), Nd);

z1 = zeros(N1,1);
z2 = zeros(N1-1,1);
kern = [
	[block1 z1 conj(fliplr([block1(1,2:N2); block2(2:N1,2:N2)]))];
	zeros(1,2*N2);
	[flipud(block2(2:N1,:)) z2 fliplr(flipud(conj(block1(2:N1, 2:N2))))]
]; % [(2Nd)]
kern = dsft_gram_hermitify(kern, show); % force Hermitian symmetry
fftkern = fftn_fast(kern);

% Gdsft_gram_init3()
function [fftkern reuse] = Gdsft_gram_init3(arg, wi, reuse, Nd, show)
if isempty(reuse)
	[reuse.G1 reuse.G2 reuse.G3 reuse.G4 ] = build_G1_G2_G3_G4(arg);
end
N1 = Nd(1);
N2 = Nd(2);
N3 = Nd(3);

filtblk1 = reshape(reuse.G1' * real(wi), Nd);
filtblk2 = reshape(reuse.G2' * real(wi), Nd);
filtblk3 = reshape(reuse.G3' * real(wi), Nd);
filtblk4 = reshape(reuse.G4' * real(wi), Nd);

tblk1 = filtblk1;
tblk2 = filtblk2(2:N1,:,:);	% remove the duplicated part with filtblk1
tblk3 = filtblk3(2:N1,2:N2,:);	% remove the duplicated part with block 1, 2, 4
tblk4 = filtblk4(:,2:N2,:);	% remove the duplicated part with block 1

z1 = zeros(N1,1,N3);
z2 = zeros(N1-1,1,N3);

% top half of the 3D filter
kern_top = [
	tblk1 z1 flipdim(tblk4,2);	% upper block
	zeros(1,2*N2,N3);		% zero padding in the middle
	flipdim(tblk2,1) z2 flipdim(flipdim(tblk3,1),2) % lower block
];

% construct the bottom half now
bblk1 = flipdim(conj(filtblk3(:,:,2:N3)),3);
bblk2 = flipdim(conj(filtblk4(:,:,2:N3)),3);
bblk2 = bblk2(2:N1,:,:);
bblk3 = flipdim(conj(filtblk1(:,:,2:N3)),3);
bblk3 = bblk3(2:N1,2:N2,:);
bblk4 = flipdim(conj(filtblk2(:,:,2:N3)),3);
bblk4 = bblk4(:,2:N2,:);

z4 = zeros(N1,1,N3-1);
z5 = zeros(N1-1,1,N3-1);
kern_bottom = [
	bblk1 z4 flipdim(bblk4,2);
	zeros(1,2*N2,N3-1);
	flipdim(bblk2,1) z5 flipdim(flipdim(bblk3,1),2)
];

kern = cat(3, kern_top, zeros(2*N1,2*N2,1));
kern = cat(3, kern, kern_bottom);

kern = dsft_gram_hermitify(kern, show); % force Hermitian symmetry
fftkern = fftn_fast(kern);

% Gdsft_gram_setup2()
% modified versions of A for 2D
% this routine is quite tricky.  nothing new is really rebuilt here,
% but some structure values are changed.
function [G1, G2] = Gdsft_gram_setup2(arg)
om = arg.om_t';
Nd = arg.Nd;
G1 = Gdsft(om, Nd);
om(:,1) = -om(:,1);
G2 = Gdsft(om, Nd);

% build_G1_G2_G3_G4()
% Created by Daehyun Yoon for 3D case
% build more Gnufft objects, with negative om.
function [G1, G2, G3, G4] = build_G1_G2_G3_G4(arg)
om = arg.om_t';
Nd = arg.Nd;
G1 = Gdsft(om, Nd);
om(:,1) = -om(:,1);
G2 = Gdsft(om, Nd);
om(:,2) = -om(:,2);
G3 = Gdsft(om, Nd);
om(:,1) = -om(:,1);
G4 = Gdsft(om, Nd);

% Gdsft_gram_mult()
% multiply an image x by a toeplitz matrix
% by embedding it into a circulant matrix of twice the size and using FFT.
% in
%	x	[(Nd)]
%	fftkern	[[2Nd]]
%	mask	[(Nd)]
% out
%	y	[(Nd)]
function y = Gdsft_gram_mult(x, fftkern, mask)

N2 = size(fftkern);
if numel(N2) == 2 && N2(2) == 1
	N2 = N2(1);
end
Nd = N2 / 2;

LL = size(x,1+numel(N2));
y = zeros([Nd LL]);

switch numel(N2)
case 1
	tmp = fft(x, N2); % [N2 L]
	tmp = tmp .* repmat(fftkern, [1 ncol(tmp)]);
	y = ifft(tmp);
	y = y(1:Nd,:); % [Nd L]
	y = y .* repmat(mask, [1 ncol(y)]); % fatrix2 requires mask

case 2
	for ll=1:LL
		tmp = ifftn_fast(fftkern .* fftn_fast(x(:,:,ll), N2));
		tmp = tmp(1:Nd(1), 1:Nd(2));
		tmp = tmp .* mask; % fatrix2 requires mask
		y(:,:,ll) = tmp;
	end

case 3
	for ll=1:LL
		tmp = ifftn_fast(fftkern .* fftn_fast(x(:,:,:,ll), N2));
		tmp = tmp(1:Nd(1), 1:Nd(2), 1:Nd(3));
		tmp = tmp .* mask; % fatrix2 requires mask
		y(:,:,:,ll) = tmp;
	end

otherwise
	fail('dim %d not done', numel(N2))
end

% Gdsft_gram_test1()
% test 1d version
function Gdsft_gram_test1
N = 16;
M = 21;
omega = 2*pi*rand(M,1);
wi = [1:size(omega,1)]';
mask = true(N,1); mask(end-3) = false;
A = Gdsft(omega, N, 'mask', mask, 'n_shift', N/2, 'show', 1);
T = build_gram(A, wi);
test_adjoint(T, 'complex', 1, 'tolre', 1e-7);
fatrix2_tests(A, 'complex', 1, 'tol_gram', 2e-7);

% Gdsft_gram_test2()
function Gdsft_gram_test2
N = [16 14];
fov = N;
ktype = 'spiral0';
[kspace omega wi] = mri_trajectory(ktype, {}, N, fov, {'voronoi'});
ig = image_geom('nx', N(1), 'ny', N(2), 'dx', 1);
mask = ellipse_im(ig, [0 0 14 15 0 1], 'oversample', 3) > 0;
A = Gdsft(omega, N, 'n_shift', N/2, 'show', 1);
fatrix2_tests(A, 'complex', 1, 'tol_gram', 2e-7);

T = build_gram(A, wi);
fatrix2_tests(T, 'complex', 1, 'tol_gram', 1e-7);
test_adjoint(T, 'complex', 1, 'tolre', 1e-7); % matches!

if 1
	x = ellipse_im(ig, 'shepplogan-emis');
	x1 = A' * (wi .* (A * x));
	x1 = embed(x1, mask);
	x2 = T * x;
	equivs(x1, x2, 'thresh', 2e-7)
end

tic
for ii=1:50, b1 = embed(A' * (wi .* (A * x(mask))), mask); end
t1 = toc;
tic
for ii=1:50, b2 = embed(T * x(mask), mask); end
t2 = toc;
equivs(b1, b2)
printm('time: A''Ax = %.3f, Tx=%.3f', t1, t2)

% Gdsft_gram_test3()
function Gdsft_gram_test3
N = [8 4 2];
rng(0)
omega = rand(21,3) * 2 * pi;
A = Gdsft(omega, N, 'n_shift', N/2, 'show', 1);
fatrix2_tests(A, 'complex', 1, 'tol_gram', 2e-7);
test_adjoint(A, 'complex', 1, 'big', 1);

wi = [1:size(omega,1)]';
T = build_gram(A, wi);
fatrix2_tests(T, 'complex', 1, 'tol_gram', 1e-7);
test_adjoint(T, 'complex', 1, 'tolre', 1e-7);

if 1
	x = cumsum(ones(N), 2);
	x1 = A' * (wi .* (A * x));
	x1 = iembed(A, x1);
	x2 = T * x;
	equivs(x1, x2)
end

% Gdsft_gram_test()
function Gdsft_gram_test
Gdsft_gram_test1
Gdsft_gram_test2
Gdsft_gram_test3

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function c=projkern(c,p2,p3,p4,p5);
%PROJKERN  Projection onto generating kernel space
%   Usage:  cout=projkern(cin,a);
%           cout=projkern(cin,g,a);
%           cout=projkern(cin,ga,gs,a);
%
%   Input parameters:
%         cin   : Input coefficients
%         g     : analysis/synthesis window
%         ga    : analysis window
%         gs    : synthesis window
%         a     : Length of time shift.
%   Output parameters:
%         cout  : Output coefficients
%
%   `cout=projkern(cin,a)` projects a set of Gabor coefficients *c* onto the
%   space of possible Gabor coefficients. This means that *cin* and *cout*
%   synthesize to the same signal. A tight window generated from a Gaussian
%   will be used for both analysis and synthesis.
%
%   The rationale for this function is a follows: Because the coefficient
%   space of a Gabor frame is larger than the signal space (since the frame
%   is redundant) then there are many coefficients that correspond to the
%   same signal.
%
%   Therefore, you might desire to work with the coefficients *cin*, but you
%   are in reality working with *cout*.
%
%   `cout=projkern(cin,g,a)` does the same, using the window *g* for analysis
%   and synthesis.
%
%   `cout=projkern(cin,ga,gs,a)` does the same, but for different analysis
%   *ga* and synthesis *gs* windows.
%
%   See also: dgt, idgt

complainif_argnonotinrange(nargin,2,4,mfilename);

M=size(c,1);
N=size(c,2);

if nargin==2
  a=p2;
  L=a*N;
  ga=gabtight(a,M,L);
  gs=ga;
end;

if nargin==3;
  ga=p2;
  gs=p2;
  a=p3;
  L=a*N;
end;

if nargin==4;
  ga=p2;
  gs=p3;
  a=p4;
  L=a*N;
end;

assert_squarelat(a,M,1,'PROJKERN');

c=dgt(idgt(c,gs,a),ga,a,M);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [Q, V, policy, mean_discrepancy] = mdp_MK_learning(P, R, discount, N)

% mdp_Q_learning   Evaluation of the matrix Q, using the Q learning algorithm
%
% Arguments
% -------------------------------------------------------------------------
% Let S = number of states, A = number of actions
%   P(SxSxA)  = transition matrix
%              P could be an array with 3 dimensions or
%              a cell array (1xA), each cell containing a sparse matrix (SxS)
%   R(SxSxA) or (SxA) = reward matrix
%              R could be an array with 3 dimensions (SxSxA) or
%              a cell array (1xA), each cell containing a sparse matrix (SxS) or
%              a 2D array(SxA) possibly sparse
%   discount  = discount rate in ]0; 1[
%   N(optional) = number of iterations to execute, default value: 10000.
%                 It is an integer greater than the default value.
% Evaluation --------------------------------------------------------------
%   Q(SxA) = learned Q matrix
%   V(S)   = learned value function.
%   policy(S) = learned optimal policy.
%   mean_discrepancy(N/100) = vector of V discrepancy mean over 100 iterations
%             Then the length of this vector for the default value of N is 100.

% MDPtoolbox: Markov Decision Processes Toolbox
% Copyright (C) 2009  INRA
% 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 the <ORGANIZATION> 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 HOLDER 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.

% check of arguments
if (discount <= 0 || discount >= 1)
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: Discount rate must be in ]0,1[')
    disp('--------------------------------------------------------')
elseif (nargin >= 4) && (N < 10000)
    disp('--------------------------------------------------------')
    disp('MDP Toolbox ERROR: N must be upper than 10000')
    disp('--------------------------------------------------------')
else

    % initialization of optional arguments
    if (nargin < 4); N=10000; end;

    % Find number of states and actions
    if iscell(P)
        S = size(P{1},1);
        A = length(P);
    else
        S = size(P,1);
        A = size(R,2);
    end;

    % Initialisations
    Q = zeros(S,A);
    dQ = zeros(S,A);
    mean_discrepancy = [];
    discrepancy = [];

    % Initial state choice
    s = randi([1,S]);
            prob=[NaN NaN];
            while nansum(prob)==0
             s = randi([1,S]);
            [s1, act, prob]=find(P(:,s,:));  %which actions are possible in that state?
            end

    h = waitbar(0,'Initializing waitbar...');

    for n=1:N

        % Reinitialisation of trajectories every 100 transitions
        if (mod(n,100)==0);
            prob=[NaN NaN];
            while nansum(prob)==0
             s = randi([1,S]);
            [s1, act, prob]=find(P(:,s,:));  %which actions are possible in that state?
            end
            waitbar(n/N,h,n/N*100);
        end;

        % Action choice : greedy with increasing probability
        % probability 1-(1/log(n+2)) can be changed

        pn = rand(1);
%         if (pn < (1-(1/log(n+2))))
%           [~,a] = max(Q(s,:));
%         else
          a = act(randi([1,numel(act)]));
%         end;

        % Simulating next state s_new and reward associated to <s,s_new,a>
        p_s_new = rand(1);
        p = 0;
        s_new = 0;
        while ((p < p_s_new) && (s_new < S))
            s_new = s_new+1;
            if iscell(P)
                p = p + P{a}(s,s_new);
            else
                p = p + P(s,s_new,a);
            end;
        end;
        if iscell(R)
            r = R{a}(s,s_new);
        elseif ndims(R) == 3
            r = R(s,s_new,a);
        else
            r = R(s,a);
        end;

        % Updating the value of Q
        % Decaying update coefficient (1/sqrt(n+2)) can be changed
        delta = r + discount*max(Q(s_new,:)) - Q(s,a);
        dQ = (1/sqrt(n+2))*delta;
        Q(s,a) = Q(s,a) + dQ;

        % Current state is updated
        s = s_new;

        % Computing and saving maximal values of the Q variation
        discrepancy(mod(n,100)+1) = abs(dQ);

        % Computing means all over maximal Q variations values
        if (length(discrepancy) == 100)
           mean_discrepancy = [ mean_discrepancy mean(discrepancy)];
           discrepancy = [];
        end;

    end;

    %compute the value function and the policy
    [V, policy] = max(Q,[],2);
close(h);
end;

author: matthieukomorowski


text:
function [xk,niter,residuals,outputData,opts] =NESTA(A,At,b,muf,delta,opts)
% [xk,niter,residuals,outputData] =NESTA(A,At,b,muf,delta,opts)
%
% Solves a L1 minimization problem under a quadratic constraint using the
% Nesterov algorithm, with continuation:
%
%     min_x || U x ||_1 s.t. ||y - Ax||_2 <= delta
%
% Continuation is performed by sequentially applying Nesterov's algorithm
% with a decreasing sequence of values of  mu0 >= mu >= muf
%
% The primal prox-function is also adapted by accounting for a first guess
% xplug that also tends towards x_muf
%
% The observation matrix A is a projector
%
% Inputs:   A and At - measurement matrix and adjoint (either a matrix, in which
%               case At is unused, or function handles).  m x n dimensions.
%           b   - Observed data, a m x 1 array
%           muf - The desired value of mu at the last continuation step.
%               A smaller mu leads to higher accuracy.
%           delta - l2 error bound.  This enforces how close the variable
%               must fit the observations b, i.e. || y - Ax ||_2 <= delta
%               If delta = 0, enforces y = Ax
%               Common heuristic: delta = sqrt(m + 2*sqrt(2*m))*sigma;
%               where sigma=std(noise).
%           opts -
%               This is a structure that contains additional options,
%               some of which are optional.
%               The fieldnames are case insensitive.  Below
%               are the possible fieldnames:
%
%               opts.xplug - the first guess for the primal prox-function, and
%                 also the initial point for xk.  By default, xplug = At(b)
%               opts.U and opts.Ut - Analysis/Synthesis operators
%                 (either matrices of function handles).
%               opts.normU - if opts.U is provided, this should be norm(U)
%                   otherwise it will have to be calculated (potentially
%                   expensive)
%               opts.MaxIntIter - number of continuation steps.
%                 default is 5
%               opts.maxiter - max number of iterations in an inner loop.
%                 default is 10,000
%               opts.TolVar - tolerance for the stopping criteria
%               opts.stopTest - which stopping criteria to apply
%                   opts.stopTest == 1 : stop when the relative
%                       change in the objective function is less than
%                       TolVar
%                   opts.stopTest == 2 : stop with the l_infinity norm
%                       of difference in the xk variable is less
%                       than TolVar
%               opts.TypeMin - if this is 'L1' (default), then
%                   minimizes a smoothed version of the l_1 norm.
%                   If this is 'tv', then minimizes a smoothed
%                   version of the total-variation norm.
%                   The string is case insensitive.
%               opts.Verbose - if this is 0 or false, then very
%                   little output is displayed.  If this is 1 or true,
%                   then output every iteration is displayed.
%                   If this is a number p greater than 1, then
%                   output is displayed every pth iteration.
%               opts.fid - if this is 1 (default), the display is
%                   the usual Matlab screen.  If this is the file-id
%                   of a file opened with fopen, then the display
%                   will be redirected to this file.
%               opts.errFcn - if this is a function handle,
%                   then the program will evaluate opts.errFcn(xk)
%                   at every iteration and display the result.
%                   ex.  opts.errFcn = @(x) norm( x - x_true )
%               opts.outFcn - if this is a function handle,
%                   then then program will evaluate opts.outFcn(xk)
%                   at every iteration and save the results in outputData.
%                   If the result is a vector (as opposed to a scalar),
%                   it should be a row vector and not a column vector.
%                   ex. opts.outFcn = @(x) [norm( x - xtrue, 'inf' ),...
%                                           norm( x - xtrue) / norm(xtrue)]
%               opts.AAtinv - this is an experimental new option.  AAtinv
%                   is the inverse of AA^*.  This allows the use of a
%                   matrix A which is not a projection, but only
%                   for the noiseless (i.e. delta = 0) case.
%               opts.USV - another experimental option.  This supercedes
%                   the AAtinv option, so it is recommended that you
%                   do not define AAtinv.  This allows the use of a matrix
%                   A which is not a projection, and works for the
%                   noisy ( i.e. delta > 0 ) case.
%                   opts.USV should contain three fields:
%                   opts.USV.U  is the U from [U,S,V] = svd(A)
%                   likewise, opts.USV.S and opts.USV.V are S and V
%                   from svd(A).  S may be a matrix or a vector.
%
%  Outputs:
%           xk  - estimate of the solution x
%           niter - number of iterations
%           residuals - first column is the residual at every step,
%               second column is the value of f_mu at every step
%           outputData - a matrix, where each row r is the output
%               from opts.outFcn, if supplied.
%           opts - the structure containing the options that were used
%
% Written by: Jerome Bobin, Caltech
% Email: bobin@acm.caltech.edu
% Created: February 2009
% Modified (version 1.0): May 2009, Jerome Bobin and Stephen Becker, Caltech
% Modified (version 1.1): Nov 2009, Stephen Becker, Caltech
%
% NESTA Version 1.1
%   See also Core_Nesterov

if nargin < 6, opts = []; end
if isempty(opts) && isnumeric(opts), opts = struct; end

%---- Set defaults
fid = setOpts('fid',1);
Verbose = setOpts('Verbose',true);
function printf(varargin), fprintf(fid,varargin{:}); end
MaxIntIter = setOpts('MaxIntIter',5,1);
TypeMin = setOpts('TypeMin','L1');
TolVar = setOpts('tolvar',1e-5);
[U,U_userSet] = setOpts('U', @(x) x );
if ~isa(U,'function_handle')
    Ut = setOpts('Ut',[]);
else
    Ut = setOpts('Ut', @(x) x );
end
xplug = setOpts('xplug',[]);
normU = setOpts('normU',[]);  % so we can tell if it's been set

residuals = []; outputData = [];
AAtinv = setOpts('AAtinv',[]);
USV = setOpts('USV',[]);
if ~isempty(USV)
    if isstruct(USV)
        Q = USV.U;  % we can't use "U" as the variable name
                    % since "U" already refers to the analysis operator
        S = USV.S;
        if isvector(S), s = S; %S = diag(s);
        else s = diag(S); end
        %V = USV.V;
    else
        error('opts.USV must be a structure');
    end
end

% -- We can handle non-projections IF a (fast) routine for computing
%    the psuedo-inverse is available.
%    We can handle a nonzero delta, but we need the full SVD
if isempty(AAtinv) && isempty(USV)
    % Check if A is a partial isometry, i.e. if AA' = I
    z = randn(size(b));
    if isa(A,'function_handle'), AAtz = A(At(z));
    else AAtz = A*(A'*z); end
    if norm( AAtz - z )/norm(z) > 1e-8
        error('Measurement matrix A must be a partial isometry: AA''=I');
    end
end

% -- Find a initial guess if not already provided.
%   Use least-squares solution: x_ref = A'*inv(A*A')*b
% If A is a projection, the least squares solution is trivial
if isempty(xplug) || norm(xplug) < 1e-12
    if ~isempty(USV) && isempty(AAtinv)
        AAtinv = Q*diag( s.^(-2) )*Q';
    end
    if ~isempty(AAtinv)
        if delta > 0 && isempty(USV)
            error('delta must be zero for non-projections');
        end
        if isa(AAtinv,'function_handle')
            x_ref = AAtinv(b);
        else
            x_ref = AAtinv * b;
        end
    else
        x_ref = b;
    end

    if isa(A,'function_handle')
        x_ref=At(x_ref);
    else
        x_ref = A'*x_ref;
    end

    if isempty(xplug)
        xplug = x_ref;
    end
    % x_ref itself is used to calculate mu_0
    %   in the case that xplug has very small norm
else
    x_ref = xplug;
end

% use x_ref, not xplug, to find mu_0
if isa(U,'function_handle')
    Ux_ref = U(x_ref);
else
    Ux_ref = U*x_ref;
end
switch lower(TypeMin)
    case 'l1'
        mu0 = 0.9*max(abs(Ux_ref));
    case 'tv'
        mu0 = ValMUTv(Ux_ref);
end

% -- If U was set by the user and normU not supplied, then calcuate norm(U)
if U_userSet && isempty(normU)
    % simple case: U*U' = I or U'*U = I, in which case norm(U) = 1
    z = randn(size(xplug));
    if isa(U,'function_handle'), UtUz = Ut(U(z)); else UtUz = U'*(U*z); end
    if norm( UtUz - z )/norm(z) < 1e-8
        normU = 1;
    else
        z = randn(size(Ux_ref));
        if isa(U,'function_handle')
            UUtz = U(Ut(z));
        else
            UUtz = U*(U'*z);
        end
        if norm( UUtz - z )/norm(z) < 1e-8
            normU = 1;
        end
    end

    if isempty(normU)
        % have to actually calculate the norm
        if isa(U,'function_handle')
            [normU,cnt] = my_normest(U,Ut,length(xplug),1e-3,30);
            if cnt == 30, printf('Warning: norm(U) may be inaccurate\n'); end
        else
            [mU,nU] = size(U);
            if mU < nU, UU = U*U'; else UU = U'*U; end
            % last resort is to call MATLAB's "norm", which is slow
            if norm( UU - diag(diag(UU)),'fro') < 100*eps
                % this means the matrix is diagonal, so norm is easy:
                normU = sqrt( max(abs(diag(UU))) );
            elseif issparse(UU)
                normU = sqrt( normest(UU) );
            else
                if min(size(U)) > 2000
                    % norm(randn(2000)) takes about 5 seconds on my PC
                    printf('Warning: calculation of norm(U) may be slow\n');
                end
                normU = sqrt( norm(UU) );
            end
        end
    end
    opts.normU = normU;
end

niter = 0;
Gamma = (muf/mu0)^(1/MaxIntIter);
mu = mu0;
Gammat= (TolVar/0.1)^(1/MaxIntIter);
TolVar = 0.1;

for nl=1:MaxIntIter

    mu = mu*Gamma;
    TolVar=TolVar*Gammat;    opts.TolVar = TolVar;
    opts.xplug = xplug;
    if Verbose, printf('\tBeginning %s Minimization; mu = %g\n',opts.TypeMin,mu); end
    [xk,niter_int,res,out,optsOut] = Core_Nesterov(...
        A,At,b,mu,delta,opts);

    xplug = xk;
    niter = niter_int + niter;

    residuals = [residuals; res];
    outputData = [outputData; out];

end
opts = optsOut;

%---- internal routine for setting defaults
function [var,userSet] = setOpts(field,default,mn,mx)
    var = default;
    % has the option already been set?
    if ~isfield(opts,field)
        % see if there is a capitalization problem:
        names = fieldnames(opts);
        for i = 1:length(names)
            if strcmpi(names{i},field)
                opts.(field) = opts.(names{i});
                opts = rmfield(opts,names{i});
                break;
            end
        end
    end
    if isfield(opts,field) && ~isempty(opts.(field))
        var = opts.(field);  % override the default
        userSet = true;
    else
        userSet = false;
    end
    % perform error checking, if desired
    if nargin >= 3 && ~isempty(mn)
        if var < mn
            printf('Variable %s is %f, should be at least %f\n',...
                field,var,mn); error('variable out-of-bounds');
        end
    end
    if nargin >= 4 && ~isempty(mx)
        if var > mx
            printf('Variable %s is %f, should be at least %f\n',...
                field,var,mn); error('variable out-of-bounds');
        end
    end
    opts.(field) = var;
end

%---- internal routine for setting mu0 in the tv minimization case
function th=ValMUTv(x)

    N = length(x);n = floor(sqrt(N));
    Dv = spdiags([reshape([-ones(n-1,n); zeros(1,n)],N,1) ...
            reshape([zeros(1,n); ones(n-1,n)],N,1)], [0 1], N, N);
        Dh = spdiags([reshape([-ones(n,n-1) zeros(n,1)],N,1) ...
            reshape([zeros(n,1) ones(n,n-1)],N,1)], [0 n], N, N);
        D = sparse([Dh;Dv]);

    Dhx = Dh*x;
    Dvx = Dv*x;

    sk = sqrt(abs(Dhx).^2 + abs(Dvx).^2);
    th = max(sk);

end

end %-- end of NESTA function

%%%%%%%%%%%% POWER METHOD TO ESTIMATE NORM %%%%%%%%%%%%%%%
% Copied from MATLAB's "normest" function, but allows function handles, not just sparse matrices
function [e,cnt] = my_normest(S,St,n,tol, maxiter)
%MY_NORMEST Estimate the matrix 2-norm via power method.
    if nargin < 4, tol = 1.e-6; end
    if nargin < 5, maxiter = 20; end
    if isempty(St)
        St = S;  % we assume the matrix is symmetric;
    end
    x = ones(n,1);
    cnt = 0;
    e = norm(x);
    if e == 0, return, end
    x = x/e;
    e0 = 0;
    while abs(e-e0) > tol*e && cnt < maxiter
       e0 = e;
       Sx = S(x);
       if nnz(Sx) == 0
          Sx = rand(size(Sx));
       end
       e = norm(Sx);
       x = St(Sx);
       x = x/norm(x);
       cnt = cnt+1;
    end
end

repo: REST


text:
function hdr = read_deymed_ini(filename)

% READ_DEYMED_INI reads EEG data from the Deymed Truescan file format
%
% Use as
%   hdr = read_deymed_ini(filename)
%
% See also READ_DEYMED_DAT

% Copyright (C) 2013, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

[p, f, x] = fileparts(filename);
headerfile = fullfile(p, [f '.ini']);
if ~exist(headerfile, 'file')
  headerfile = fullfile(p, [f '.Ini']);
end
datafile = fullfile(p, [f '.154576dat']);
if ~exist(datafile, 'file')
  datafile = fullfile(p, [f '.Dat']);
end

fid = fopen_or_error(datafile, 'rb');
hdr.id        = fread(fid, [1, 12], 'uint8=>char');
hdr.lastname  = fread(fid, [1, 12], 'uint8=>char');
hdr.firstname = fread(fid, [1, 12], 'uint8=>char');
hdr.date      = fread(fid, [1, 12], 'uint8=>char');
hdr.time      = fread(fid, [1, 12], 'uint8=>char');
hdr.Fs        = fread(fid, [1, 1],  'int16');
hdr.nChans    = fread(fid, [1, 1],  'int8');

dummy         = fread(fid, [1, 1],  'int8'); % this is 4 in the example file
if dummy~=4
  ft_warning('deviant header format detected');
end

% A numerical comparison between a deymed dat file and a BESA-exported edf version of the
% same file showed a calibration difference of 1/4. So this might be a calibration value.
hdr.calib = 1/dummy;

hdr.label = {};
for i=1:hdr.nChans
  str = fread(fid, [1, 6], 'uint8=>char');
  str(str==0)  = ' ';
  hdr.label{i} = strtrim(str);
end
fclose(fid);

% determine the number of samples
d = dir(datafile);
hdr.nSamples = (d.bytes-512-1024*hdr.nChans)/(2*hdr.nChans);

% also store some details from the ini file
try
  hdr.ini.patient.id        = inifile(headerfile, 'read', {'Pacient', '', 'ID'});
  hdr.ini.patient.lastname  = inifile(headerfile, 'read', {'Pacient', '', 'Prijmeni'});
  hdr.ini.patient.firstname = inifile(headerfile, 'read', {'Pacient', '', 'Jmeno'});
  hdr.ini.record.user       = inifile(headerfile, 'read', {'Record',  '', 'AcqUser'});
  hdr.ini.record.type       = inifile(headerfile, 'read', {'Record',  '', 'Type'});
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


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

function [img mat] = MU_load_nifti(fname)
% Load a NIfTI file in as a matlab marix.
% Supports a 4-D NIfTI File
% Format FUNCTION [img mat] = load_nifti(fname)
%
% Inputs:
%    fname - the name of the .nii file
%
% Outputs:
%    img - a 3d (or 4d) matlab image matrix
%    mat - affine matrix describing voxel size & position
%
% External Requirements: SPM (Tested with SPM8)
% Acknowledge to Samuel A. Hurley

% Check if its a .nii.gz
ftype = 0;

if strcmp(fname(length(fname)-2:length(fname)), '.gz');
  ftype = 1;  % gzipped image
  disp('Unzipping image...');
  eval(['!gunzip ' fname]);

  % Change the fname to reflect the un-zipped image
  fname = fname(1:length(fname)-3);
end

% Grab SPM header
D = spm_vol(fname);

% Preallocate img matrix
img = zeros([D(1).dim length(D)]);

% Grab the image
for ii = 1:length(D)
  img(:,:,:,ii) = spm_read_vols(D(ii));
end

% Re-orient the image
img = imflipud(imtranspose(img));

% Grab the mat
mat = D.mat;

% Check if the file needs to be re-zipped
if ftype == 1
  disp('Re-zipping image...');
  eval(['!gzip ' fname]);
end
% Done.

function imout = imtranspose(imin)

% Preallocate, but only for a square image
dims = size(imin);
if dims(1) == dims(2)
  imout = zeros(size(imin));
end

if ndims(imin) == 2
  imout = imin';

% 3-D Image Matrix
elseif ndims(imin) == 3
  zdim = size(imin, 3);

  if zdim == 1
    imout = imin';
  else
    for j = 1:zdim
      imout(:,:,j) = imin(:,:,j)';
    end
  end

% 4-D Image
elseif ndims(imin) == 4
  % Do some recursive magic!
  for ii = 1:size(imin, 4)
    imout(:,:,:,ii) = imtranspose(imin(:,:,:,ii));
  end

% Will we ever do a 5-D image?
else
  error('Does not support 5-D image');
end

% Flips a 3-d array in the up/down direction
function flipped = imflipud(im)

flipped = zeros(size(im));

% 3-D
if ndims(im) == 2
  flipped = flipud(im);

elseif ndims(im) == 3
  for j = 1:size(im,3)
    flipped(:,:,j) = flipud(im(:,:,j));
  end

% 4-D

elseif ndims(im) == 4
  for ii = 1:size(im, 4)
    flipped(:,:,:,ii) = imflipud(im(:,:,:,ii));
  end

else
  error('Does not support 5-D image');
end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function pass = test_iszero( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end

% Test with function 1-1
f = ballfun(ones(21,20,22));
g = f-f;
pass(1) = iszero(g);

% Test with function 10^-20
f = ballfun(10^(-20));
pass(2) = (iszero(f) == 0);

% Test with function 0
f = ballfun(0);
pass(3) = iszero(f);
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%drawCrackPattern : dessine un motif de fissures
%
%
%   ---------
%
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 15/08/2005.

%% Initialisations

% size of window
window = [0 100 0 100];

% box as polygon
box = [0 0;100 0;100 100;0 100];

% points density
lambda = .005;

%% compute diagram

% number of points
Np = round((window(2)-window(1))*(window(4)-window(3)) * lambda);

% points coordinate
x = rand(Np, 1)*(window(2)-window(1)) + window(1);
y = rand(Np, 1)*(window(4)-window(3)) + window(3);
points = [x y];

% random angles
angles = rand(Np, 1)*pi;

% compute crack pattern
edges = crackPattern(box, points, angles);

%% Draw Result

% create figure
figure(1); clf;
axis(window);hold on;

% draw diagram
drawEdge(edges);
drawPoint(points, 'bo');

% clean up window
set(gca, 'box', 'on');
set(gca, 'xtick', []);
set(gca, 'ytick', []);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [otf] = psf2otf_Dy(outSize)

psf = single(zeros(outSize));
psf(1, 1) = -1;
psf(end, 1) = 1;
otf = fft2(psf);

author: wliusjtu


text:
function [tUs, odrIdx, TXmean, Wgt]  = MPCA(TX,gndTX,testQ,maxK)
% MPCA: Multilinear Principle Component Analysis
%
% %[Prototype]%
% function [tUs, odrIdx, TXmean, Wgt]  = MPCA(TX,gndTX,testQ,maxK)
%
% %[Author Notes]%
% Author: Haiping LU
% Email : hplu@ieee.org   or   eehplu@gmail.com
% Affiliation: Department of Electrical and Computer Engineering
%              University of Toronto
% Release date: June 24, 2008 (Version 1.1)
% Please email me if you have any problem, question or suggestion
%
% %[Algorithm]%:
% This function implements the Multilinear Principal Component Analysis
% (MPCA) algorithm presented in the follwing paper:
%    Haiping Lu, K.N. Plataniotis, and A.N. Venetsanopoulos,
%    "MPCA: Multilinear Principal Component Analysis of Tensor Objects",
%    IEEE Transactions on Neural Networks,
%    Vol. 19, No. 1, Page: 18-39, January 2008.
% Please reference this paper when reporting work done using this code.
%
% %[Toolbox needed]%:
% This function needs the tensor toolbox available at
% http://csmr.ca.sandia.gov/~tgkolda/TensorToolbox/
%
% %[Syntax]%: [tUs, odrIdx, TXmean, Wgt]  = MPCA(TX,gndTX,testQ,maxK)
%
% %[Inputs]%:
%    TX: the input training data in tensorial representation, the last mode
%        is the sample mode. For Nth-order tensor data, TX is of
%        (N+1)th-order with the (N+1)-mode to be the sample mode.
%        E.g., 30x20x10x100 for 100 samples of size 30x20x10
%        If your training data is too big, resulting in the "out of memory"
%        error, you could work around this problem by reading samples one
%        by one from the harddisk, or you could email me for help.
%
%    gndTX: the ground truth class labels (1,2,3,...) for the training data
%           E.g., a 100x1 vector if there are 100 samples
%           If the class label is not available (unsupervised learning),
%           please set gndTX=-1;
%
%    testQ: the percentage of variation kept in each mode, suggested value
%           is 97, and you can try other values, e.g., from 95 to 100, to
%           see whether better performance can be obtained.
%
%    maxK: the maximum number of iterations, suggested value is 1, and you
%          can try a larger value if computational time is not a concern.
%
% %[Outputs]%:
%    tUs: the multilinear projection, consiting of N
%         projection matrices, one for each mode
%
%    odrIdx: the ordering index of projected features in decreasing
%            variance (if unsupervised) or discriminality (if supervised)
%            for vectorizing the projected tensorial features
%
%    TXmean: the mean of the input training samples TX
%
%    Wgt: the weight tensor for use in modified distance measures. Please
%         refer to Section IV.B and IV.C of the paper.
%
% %[Supported tensor order]%
% This function supports N=2,3,4, for other order N, please modify the
% codes accordingly or email hplu@ieee.org or eehplu@gmail.com for help
%
% %[Examples]%
%%%%%%%%%%%%%%%%%%%%%%%%%%Example on 2D face data%%%%%%%%%%%%%%%%%%%%%%%%%%
%       load FERETC80A45%each sample is a second-order tensor of size 32x32
%       N=ndims(fea2D)-1;%Order of the tensor sample
%       Is=size(fea2D);%32x32x320
%       numSpl=Is(3);%There are 320 face samples
%       testQ=97;%Keep 97% variation in each mode
%       maxK=1;%One iteration only
%       [tUs, odrIdx, TXmean, Wgt] = MPCA(fea2D,gnd,testQ,maxK);
%       fea2Dctr=fea2D-repmat(TXmean,[ones(1,N), numSpl]);%Centering
%       newfea=ttm(tensor(fea2Dctr),tUs,1:N);%MPCA projection
%       %Vectorization of the tensorial feature
%       newfeaDim=size(newfea,1)*size(newfea,2);
%       newfea=reshape(newfea.data,newfeaDim,numSpl)';%Note: Transposed
%       selfea=newfea(:,odrIdx(1:P));%Select the first "P" sorted features
%       %"P" is the dimension of the final feature vector to be fed into a
%       %standard classifier (e.g., nearest neighbor classifier), you may
%       %need to test different values of P for best performance
%       Wgt=reshape(Wgt,newfeaDim,1);%Vectorizing weight tensor
%       Wgt=Wgt(odrIdx);%Select the weights accordingly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%Example on 3D gait data%%%%%%%%%%%%%%%%%%%%%%%%%%
%       load USF17Gal %each sample is a third-order tensor of size 32x22x10
%       N=ndims(fea3D)-1;%Order of the tensor sample
%       Is=size(fea3D);%32x22x10x731
%       numSpl=Is(4);%There are 731 gait samples
%       testQ=97;%Keep 97% variation in each mode
%       maxK=1;%One iteration only
%       [tUs, odrIdx, TXmean, Wgt] = MPCA(fea3D,gnd,testQ,maxK);
%       fea3Dctr=fea3D-repmat(TXmean,[ones(1,N), numSpl]);%Centering
%       newfea = ttm(tensor(fea3Dctr),tUs,1:N);%MPCA projection
%       %Vectorization of the tensorial feature
%       newfeaDim=size(newfea,1)*size(newfea,2)*size(newfea,3);
%       newfea=reshape(newfea.data,newfeaDim,numSpl)';%Note: Transposed
%       selfea=newfea(:,odrIdx(1:P));%Select the first "P" sorted features
%       %"P" is the dimension of the final feature vector to be fed into a
%       %standard classifier (e.g., nearest neighbor classifier), you may
%       %need to test different values of P for best performance
%       Wgt=reshape(Wgt,newfeaDim,1);%Vectorizing weight tensor
%       Wgt=Wgt(odrIdx);%Select the weights accordingly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %[Notes]%:
% A. Developed using Matlab R2006a
% B. Revision history:
%       Version 1.0 released on March 1, 2008
%       Version 1.1 released on June 24, 2008
%           ---Example usage on 2D data is included
%       Version 1.2 released on March 12, 2011
%           ---No change to this code, minor change in MPCALDA.m
%           ---Inclusion of a survey paper for in-depth analysis
%           ---Inclusion of a BibTeX file with MPCA extensions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%TX: (N+1)-dimensional tensor of Tensor Sample Dimension x NumSamples
N=ndims(TX)-1;%The order of samples.
IsTX=size(TX);
Is=IsTX(1:N);%The dimensions of the tensor
numSpl=IsTX(N+1);%Number of samples

%%%%%%%%%%%%%Zero-Mean%%%%%%%%%%
TXmean=mean(TX,N+1);%The mean
TX=TX-repmat(TXmean,[ones(1,N), numSpl]);%Centering

%The full projection for initialization
Qs=ones(N,1)*testQ;
Us=cell(N,1);
tUs=cell(N,1);
Lmds=cell(N,1);
for n=1:N
    In=Is(n);Phi=zeros(In,In);
    for m=1:numSpl
        switch N
            case 2
                Xm=TX(:,:,m);
            case 3
                Xm=TX(:,:,:,m);
            case 4
                Xm=TX(:,:,:,:,m);
            otherwise
                error('Order N not supported. Please modify the code here or email hplu@ieee.org for help.')
        end
        tX=tensor(Xm);
        tXn=tenmat(tX,n);
        Xn=tXn.data;
        Phi=Phi+Xn*Xn';
    end
    [Un,Lmdn]=eig(Phi);
    Lmd=diag(Lmdn);
    [stLmd,stIdx]=sort(Lmd,'descend');
    Us{n}=Un(:,stIdx);
    tUs{n}=Us{n}';
    Lmds{n}=Lmd(stIdx);
end

%Cumulative distribution of eigenvalues
cums=cell(N,1);
for n=1:N
    In=length(Lmds{n});
    cumLmds=zeros(In,1);
    Lmd=Lmds{n};
    cumLmds(1)=Lmd(1);
    for in=2:In
        cumLmds(in)=cumLmds(in-1)+Lmd(in);
    end
    cumLmds=cumLmds./sum(Lmd);
    cums{n}=cumLmds;
end

%MPCA Iterations
if maxK>0
    tPs=cell(N,1);
    pUs=cell(N,1);
    %%%%%%%%%%%%%Determine Rn, the dimension of projected space%%%%
    for n=1:N
        cum=cums{n};
        idxs=find(cum>=Qs(n)/100);
        Ps(n)=idxs(1);
        tUn=tUs{n};
        tPn=tUn(1:Ps(n),:);
        tPs{n}=tPn;
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for iK=1:maxK
        for n=1:N
            In=Is(n);
            Phi=double(zeros(In,In));
            for m=1:numSpl
                switch N
                    case 2
                        Xm=TX(:,:,m);
                    case 3
                        Xm=TX(:,:,:,m);
                    case 4
                        Xm=TX(:,:,:,:,m);
                    otherwise
                        error('Order N not supported. Please modify the code here or email hplu@ieee.org for help.')
                end
                tX=ttm(tensor(Xm),tPs,-n);
                tXn=tenmat(tX,n);
                Xn=tXn.data;
                Phi=Phi+Xn*Xn';
            end
            Pn=Ps(n);
            Phi=double(Phi);
            if Pn<In
                option=struct('disp',0);
                [pUs{n},pLmdn]=eigs(Phi,Pn,'lm',option);
                pLmds{n}=diag(pLmdn);
            else
                [pUn,pLmdn]=eig(Phi);
                pLmd=diag(pLmdn);
                [stLmd,stIdx]=sort(pLmd,'descend');
                pUs{n}=pUn(:,stIdx(1:Pn));
                pLmds{n}=pLmd(stIdx(1:Pn));
            end
            tPs{n}=pUs{n}';
        end
    end
    Us=pUs;
    tUs=tPs;
    Lmds=pLmds;
    Is=Ps;
else
    if testQ<100
        error('At least one iteration is needed');
    end
end

%Calculate the weight tensor Wgt
Wgt=zeros(Is);
switch N
    case 2
        for i1=1:Is(1)
            for i2=1:Is(2)
                Wgt(i1,i2)=sqrt(Lmds{1}(i1)*Lmds{2}(i2));
            end
        end
    case 3
        for i1=1:Is(1)
            for i2=1:Is(2)
                for i3=1:Is(3)
                    Wgt(i1,i2,i3)=sqrt(Lmds{1}(i1)*Lmds{2}(i2)*Lmds{3}(i3));
                end
            end
        end
    case 4
        for i1=1:Is(1)
            for i2=1:Is(2)
                for i3=1:Is(3)
                    for i4=1:Is(4)
                        Wgt(i1,i2,i3,i4)=sqrt(Lmds{1}(i1)*Lmds{2}(i2)*Lmds{3}(i3)*Lmds{4}(i4));
                    end
                end
            end
        end
    otherwise
        error('Order N not supported. Please modify the code here or email hplu@ieee.org for help.')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Yps=ttm(tensor(TX),tUs,1:N);%MPCA projections of samples TX
vecDim=1;
for n=1:N, vecDim=vecDim*Is(n); end
vecYps=reshape(Yps.data,vecDim,numSpl); %vectorization of Yps
%%%%%%%%%%%%%%Now vecYps contains the feature vectors for training data

if max(gndTX)<0%%%%%%%%%%%%%%%%%%%%%%%%Sort by Variance%%%%%%%%%%%%%%%%%%%%
    TVars=diag(vecYps*vecYps');
    [stTVars,odrIdx]=sort(TVars,'descend');
else%%%%%%%%%%%%%%%Sort according to Fisher's discriminality%%%%%%%%%%%%%%%
    classLabel = unique(gndTX);
    nClass = length(classLabel);%Number of classes
    ClsIdxs=cell(nClass);
    Ns=zeros(nClass,1);
    for i=1:nClass
        ClsIdxs{i}=find(gndTX==classLabel(i));
        Ns(i)=length(ClsIdxs{i});
    end
    Ymean=mean(vecYps,2);
    TSW=zeros(vecDim,1);
    TSB=zeros(vecDim,1);
    for i=1:nClass
        clsYp=vecYps(:,ClsIdxs{i});
        clsMean=mean(clsYp,2);
        FtrDiff=clsYp-repmat(clsMean,1,Ns(i));
        TSW=TSW+sum(FtrDiff.*FtrDiff,2);
        meanDiff=clsMean-Ymean;
        TSB=TSB+Ns(i)*meanDiff.*meanDiff;
    end
    FisherRatio=TSB./TSW;
    [stRatio,odrIdx]=sort(FisherRatio,'descend');
end


text:
function [ fea, out ] = ex_navierstokes12( varargin )
%EX_NAVIERSTOKES12 3D Example flow over a backwards facing step
%
%   [ FEA, OUT ] = EX_NAVIERSTOKES12( VARARGIN ) Sets up and solves stationary and
%   laminar 3D flow over a backwards facing step. Accepts the following property/value pairs.
%
%       Input       Value/{Default}        Description
%       -----------------------------------------------------------------------------------
%       rho         scalar {1}             Density
%       miu         scalar {2/3/389}       Molecular/dynamic viscosity
%       uin         scalar {1}             Magnitude of inlet velocity
%       sf_u        string {sflag1}        Shape function for velocity
%       sf_p        string {sflag1}        Shape function for pressure
%       solver      string 'openfoam'/{''} Use OpenFOAM or default solver
%       iplot       scalar 0/{1}           Plot solution and error (=1)
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       fea         struct                 Problem definition struct
%       out         struct                 Output struct
%
%   See also EX_NAVIERSTOKES4

% Copyright 2013-2022 Precise Simulation, Ltd.

cOptDef = {   ...
  'rho',      1;
  'miu',      2/3/389;
  'uin',      1;
  'igrid',    1;
  'sf_u',     'sflag1';
  'sf_p',     'sflag1';
  'solver',   '';
  'iplot',    1;
  'tol',      0.55;
  'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Geometry.
h_step = 0.0049/0.0101;
l_inlet = 0.02/0.0101;
l_channel = 0.08/0.0101;
fea.sdim = { 'x', 'y', 'z' };
gobj1 = gobj_block( -l_inlet, l_channel, -0.5, 0.5, -h_step, 1-h_step, 'B1' );
gobj2 = gobj_block( -l_inlet, 0, -0.5, 0.5, -h_step, 0, 'B2' );
fea.geom.objects = { gobj1, gobj2 };
fea = geom_apply_formula( fea, 'B1-B2' );

% Grid generation.
if( opt.igrid>=1 )
  n = 4;
  fea.grid = rectgrid(10*n,n,[-l_inlet, l_channel; -0.5, 0.5]);
  fea.grid = delcells( fea.grid, selcells(fea.grid,'(y<=0).*(x<=0)') );
  ix = find( fea.grid.p(2,:) <= -0.5 + sqrt(eps) );
  fea.grid.p(2,ix) = -h_step;
  ix = find( fea.grid.p(2,:) >= 0.5 - sqrt(eps) );
  fea.grid.p(2,ix) = 1-h_step;

  fea.grid = gridextrude( fea.grid, n, 1, -2 );
  fea.grid.p(2,:) = fea.grid.p(2,:) + 0.5;
  fea.grid = assign_bdr( fea.grid, fea.geom );
  for i=1:opt.igrid
    fea.grid = gridrefine( fea.grid, fid );
  end
else
  fea.grid = gridgen( fea, 'hmax', 0.1, 'fid', fid );
  % fea.grid = gridsmooth( tet2hex( fea.grid ), 5 );
end

% Problem definition.
fea = addphys( fea, @navierstokes );
fea.phys.ns.eqn.coef{1,end} = { opt.rho };
fea.phys.ns.eqn.coef{2,end} = { opt.miu };
fea.phys.ns.sfun            = { opt.sf_u opt.sf_u opt.sf_u opt.sf_p };
% fea.phys.ns.prop.artstab.iupw = 4;
if( any(strcmp(opt.solver,{'openfoam','su2'})) )
  [fea.phys.ns.sfun{:}] = deal('sflag1');
end

% Boundary conditions.
i_inflow  = findbdr( fea, ['x<',num2str(-l_inlet+1e-3)] );   % Inflow boundary number.
i_outflow = findbdr( fea, ['x>',num2str( l_channel-1e-3)] );   % Outflow boundary number.
% s_inflow  = ['4*',num2str(umax),'*(y*(',num2str((1-y)*h),'-y))/',num2str((1-y)*h),'^2'];   % Definition of inflow profile.
s_inflow  = ['4*',num2str(opt.uin),'*(z*(',num2str(1-h_step),'-z))/(1-',num2str(1-h_step),')^2'];
u_init    = [s_inflow,'*(z>0)'];
fea.phys.ns.bdr.sel(i_inflow) = 2;
fea.phys.ns.bdr.sel(i_outflow) = 4;
fea.phys.ns.bdr.coef{2,end}{1,i_inflow} = s_inflow;
if( ~strcmp(opt.solver,'openfoam') )
  fea.phys.ns.eqn.coef{6,end} = { u_init };
end

% Parse and solve problem.
fea  = parsephys( fea );
fea  = parseprob( fea );
if( strcmp(opt.solver,'openfoam') )
  logfid = fid; if( ~got.fid ), fid = []; end
  fea.sol.u = openfoam( fea, 'fid', fid, 'logfid', logfid );
  fid = logfid;
elseif( strcmp(opt.solver,'su2') )
  logfid = fid; if( ~got.fid ), fid = []; end
  fea.sol.u = su2( fea, 'fid', fid, 'logfid', logfid );
  fid = logfid;
else
  fea.sol.u = solvestat( fea, 'maxnit',50, 'nlrlx',1, 'tolchg',1e-3, 'fid', fid );
end

% Postprocessing.
if( opt.iplot>0 )
  postplot( fea, 'sliceexpr', 'sqrt(u^2+v^2+w^2)' )
end

% Error checking.
[~,slen] = minmaxsubd( ['(u<-eps)*x/',num2str(h_step),'*(z<0)*(y<0.01)*(y>-0.01)'], fea );
if( ~isempty(fid) )
  fprintf(fid,'\nRecirculation zone length: %3f (Ref: 7.93)\n\n',slen)
  fprintf(fid,'\n\n')
end

out.slen = [slen, 7.93];
out.err  = abs(diff(out.slen))/out.slen(end);
out.pass = out.err<opt.tol;

if ( nargout==0 )
  clear fea out
end

author: precise-simulation

repo: featool-multiphysics


text:
function [idx dst] = prisearch(ann, query, k, epsl, asm)
%
% kNN search
%   Usage:
%     [idx dst] = prisearch(ann, query, k, asm, eps)
%
% Inputs:
%   ann - ann class object
%   query - (d)x(N) query points
%   k - number of nearest nieghbors (if points in ann < k than less than k
%                                    points are returned)
%   epsl - epsilon search precision
%   asm - allow self match flag, if false points with dst = 0 are ignored
%         (defualt is true)
%
if nargin == 4
    asm = true;
end

if ~asm
    k = k+1;
end

if ~isa(query, ann.ccls)
    query = ann.cfun(query);
end

[idx dst] = annmex(ann.modes.PRISEARCH, ann.kd_ptr, query, k, epsl);

if ~asm
    gsm = dst(1,:)==0;
    dst(1:end-1,gsm) = dst(2:end,gsm);
    idx(1:end-1,gsm) = idx(2:end,gsm);
    dst(end,:) = [];
    idx(end,:) = [];
end
idx = idx + 1; % fix zero indexing of ann

author: CHLNDDEV


text:
function [ know, x ] = p40_sol ( n )

%*****************************************************************************80
%
%% P40_SOL returns the solution for problem 40.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the problem.  This value
%    is only needed for those problems with variable N.
%
%    Output, integer KNOW.
%    If KNOW is 0, then the solution is not known.
%    If KNOW is positive, then the solution is known, and is returned in X.
%
%    Output, real X(N), the solution, if known.
%
  know = 1;

  x = zeros ( n, 1 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [net, info] = bcnn_train_simplenn(net, imdb, getBatch_train, getBatch_val, varargin)

% BNN_TRAIN_SW   training a symmetric BCNN
%    BCNN_TRAIN() is an example learner implementing stochastic gradient
%    descent with momentum to train a symmetric BCNN for image classification.
%    The net should built in simplenn structure

% INPUT
% net: a network with a bcnn layer
% imdb: imdb structure of a dataset
% getBatch_train: function to read a batch of training images
% getBatch_val: function to read a batch of validation images

% OUTPUT
% net: bcnn network after fine-tuning
% info: log of training and validation

% A symmetric BCNN consists of multiple layers of convolutions, pooling, and
% nonlinear activation. A bilinearpool layer is built on top followed by square-root
% and L2 normalization and a sofmaxloss layer.

% Copyright (C) 2015 Tsung-Yu Lin, Aruni RoyChowdhury, Subhransu Maji.
% All rights reserved.
%
% This file is part of the BCNN and is made available under
% the terms of the BSD license (see the COPYING file).

% This function is modified from CNN_TRAIN of MatConvNet

% basic setting
opts.train = [] ;
opts.val = [] ;
opts.numEpochs = 300 ;
opts.batchSize = 256 ;
opts.numSubBatches = 1 ;
opts.gpus = [] ; % which GPU devices to use (none, one, or more)
opts.learningRate = 0.001 ;
opts.continue = false ;
opts.expDir = fullfile('data','exp') ;
opts.conserveMemory = false ;
opts.backPropDepth = +inf ;
opts.sync = false ;
opts.prefetch = false ;
opts.cudnn = true ;
opts.weightDecay = 0.0005 ;
opts.momentum = 0.9 ;
opts.errorFunction = 'multiclass' ;
opts.errorLabels = {} ;
opts.plotDiagnostics = false ;
opts.memoryMapFile = fullfile(tempdir, 'matconvnet.bin') ;

opts = vl_argparse(opts, varargin) ;

if ~exist(opts.expDir, 'dir'), mkdir(opts.expDir) ; end
if isempty(opts.train), opts.train = find(imdb.images.set==1) ; end
if isempty(opts.val), opts.val = find(imdb.images.set==2) ; end
if isnan(opts.train), opts.train = [] ; end

% -------------------------------------------------------------------------
%                                                    Network initialization
% -------------------------------------------------------------------------
% set hyperparameters
evaluateMode = isempty(opts.train) ;

if ~evaluateMode
  for i=1:numel(net.layers)
    if isfield(net.layers{i}, 'weights')
      J = numel(net.layers{i}.weights) ;
      for j=1:J
        net.layers{i}.momentum{j} = zeros(size(net.layers{i}.weights{j}), 'single') ;
      end
      if ~isfield(net.layers{i}, 'learningRate')
        net.layers{i}.learningRate = ones(1, J, 'single') ;
      end
      if ~isfield(net.layers{i}, 'weightDecay')
        net.layers{i}.weightDecay = ones(1, J, 'single') ;
      end
    end
    % Legacy code: will be removed
    if isfield(net.layers{i}, 'filters')
      net.layers{i}.momentum{1} = zeros(size(net.layers{i}.filters), 'single') ;
      net.layers{i}.momentum{2} = zeros(size(net.layers{i}.biases), 'single') ;
      if ~isfield(net.layers{i}, 'learningRate')
        net.layers{i}.learningRate = ones(1, 2, 'single') ;
      end
      if ~isfield(net.layers{i}, 'weightDecay')
        net.layers{i}.weightDecay = single([1 0]) ;
      end
    end
  end
end

% -------------------------------------------------------------------------
%                                                Move network to GPU or CPU
% -------------------------------------------------------------------------

% setup GPUs
numGpus = numel(opts.gpus) ;
if numGpus > 1
  if isempty(gcp('nocreate')),
    parpool('local',numGpus) ;
    spmd, gpuDevice(opts.gpus(labindex)), end
  end
elseif numGpus == 1
  gpuDevice(opts.gpus)
end
if exist(opts.memoryMapFile), delete(opts.memoryMapFile) ; end

% setup error calculation function
if isstr(opts.errorFunction)
  switch opts.errorFunction
    case 'none'
      opts.errorFunction = @error_none ;
    case 'multiclass'
      opts.errorFunction = @error_multiclass ;
      if isempty(opts.errorLabels), opts.errorLabels = {'top1e', 'top5e'} ; end
    case 'binary'
      opts.errorFunction = @error_binary ;
      if isempty(opts.errorLabels), opts.errorLabels = {'bine'} ; end
    otherwise
      error('Uknown error function ''%s''', opts.errorFunction) ;
  end
end

modelPath = @(ep) fullfile(opts.expDir, sprintf('net-epoch-%d.mat', ep));
modelFigPath = fullfile(opts.expDir, 'net-train.pdf') ;

start = opts.continue * findLastCheckpoint(opts.expDir) ;
if start >= 1
  fprintf('resuming by loading epoch %d\n', start) ;
  load(modelPath(start), 'net', 'info') ;
end

for epoch=start+1:opts.numEpochs

  % train one epoch and validate
  learningRate = opts.learningRate(min(epoch, numel(opts.learningRate))) ;
  train = opts.train(randperm(numel(opts.train))) ; % shuffle
  val = opts.val ;
  if numGpus <= 1
    [net,stats.train] = process_epoch(opts, getBatch_train, epoch, train, learningRate, imdb, net) ;
    [~,stats.val] = process_epoch(opts, getBatch_val, epoch, val, 0, imdb, net) ;
  else
    spmd(numGpus)
      [net_, stats_train_] = process_epoch(opts, getBatch_train, epoch, train, learningRate, imdb, net) ;
      [~, stats_val_] = process_epoch(opts, getBatch_val, epoch, val, 0, imdb, net_) ;
    end
    net = net_{1} ;
    stats.train = sum([stats_train_{:}],2) ;
    stats.val = sum([stats_val_{:}],2) ;
  end

  % save
  if evaluateMode, sets = {'val'} ; else sets = {'train'} ; end
  if ~evaluateMode && ~isempty(opts.val)
      sets{2} = 'val';
  end
  for f = sets
    f = char(f) ;
    n = numel(eval(f)) ;
    info.(f).speed(epoch) = n / stats.(f)(1) * max(1, numGpus) ;
    info.(f).objective(epoch) = stats.(f)(2) / n ;
    info.(f).error(:,epoch) = stats.(f)(3:end) / n ;
  end
  if ~evaluateMode, save(modelPath(epoch), 'net', 'info') ; end

  figure(1) ; clf ;
  hasError = isa(opts.errorFunction, 'function_handle') ;
  subplot(1,1+hasError,1) ;
  if ~evaluateMode
    semilogy(1:epoch, info.train.objective, '.-', 'linewidth', 2) ;
    hold on ;
  end
  if ~isempty(opts.val)
      semilogy(1:epoch, info.val.objective, '.--') ;
      xlabel('training epoch') ; ylabel('energy') ;
  end
  grid on ;
  h=legend(sets) ;
  set(h,'color','none');
  title('objective') ;
  if hasError
    subplot(1,2,2) ; leg = {} ;
    if ~evaluateMode
      plot(1:epoch, info.train.error', '.-', 'linewidth', 2) ;
      hold on ;
      leg = horzcat(leg, strcat('train ', opts.errorLabels)) ;
    end

    if ~isempty(opts.val)
        plot(1:epoch, info.val.error', '.--') ;
        leg = horzcat(leg, strcat('val ', opts.errorLabels)) ;
    end
    set(legend(leg{:}),'color','none') ;
    grid on ;
    xlabel('training epoch') ; ylabel('error') ;
    title('error') ;
  end
  drawnow ;
  print(1, modelFigPath, '-dpdf') ;
end

% -------------------------------------------------------------------------
function err = error_multiclass(opts, labels, res)
% -------------------------------------------------------------------------
predictions = gather(res(end-1).x) ;
[~,predictions] = sort(predictions, 3, 'descend') ;

% be resilient to badly formatted labels
if numel(labels) == size(predictions, 4)
  labels = reshape(labels,1,1,1,[]) ;
end

% skip null labels
mass = single(labels(:,:,1,:) > 0) ;
if size(labels,3) == 2
  % if there is a second channel in labels, used it as weights
  mass = mass .* labels(:,:,2,:) ;
  labels(:,:,2,:) = [] ;
end

error = ~bsxfun(@eq, predictions, labels) ;
err(1,1) = sum(sum(sum(mass .* error(:,:,1,:)))) ;
err(2,1) = sum(sum(sum(mass .* min(error(:,:,1:5,:),[],3)))) ;

% -------------------------------------------------------------------------
function err = error_binaryclass(opts, labels, res)
% -------------------------------------------------------------------------
predictions = gather(res(end-1).x) ;
error = bsxfun(@times, predictions, labels) < 0 ;
err = sum(error(:)) ;

% -------------------------------------------------------------------------
function err = error_none(opts, labels, res)
% -------------------------------------------------------------------------
err = zeros(0,1) ;

% -------------------------------------------------------------------------
function  [net_cpu,stats,prof] = process_epoch(opts, getBatch, epoch, subset, learningRate, imdb, net_cpu)
% -------------------------------------------------------------------------

% move CNN to GPU as needed
numGpus = numel(opts.gpus) ;
if numGpus >= 1
  net = vl_simplenn_move(net_cpu, 'gpu') ;
else
  net = net_cpu ;
  net_cpu = [] ;
end

% validation mode if learning rate is zero
training = learningRate > 0 ;if training
  mode = 'train' ;
  evalMode = 'normal' ;
else
  mode = 'val' ;
  evalMode = 'test' ;
end
if nargout > 2, mpiprofile on ; end

numGpus = numel(opts.gpus) ;
if numGpus >= 1
  one = gpuArray(single(1)) ;
else
  one = single(1) ;
end
res = [] ;
mmap = [] ;
stats = [] ;
start = tic ;

for t=1:opts.batchSize:numel(subset)
  fprintf('%s: epoch %02d: batch %3d/%3d: ', mode, epoch, ...
          fix(t/opts.batchSize)+1, ceil(numel(subset)/opts.batchSize)) ;
  batchSize = min(opts.batchSize, numel(subset) - t + 1) ;
  numDone = 0 ;
  error = [] ;
  for s=1:opts.numSubBatches
    % get this image batch and prefetch the next
    batchStart = t + (labindex-1) + (s-1) * numlabs ;
    batchEnd = min(t+opts.batchSize-1, numel(subset)) ;
    batch = subset(batchStart : opts.numSubBatches * numlabs : batchEnd) ;
    [im, labels] = getBatch(imdb, batch) ;

    if opts.prefetch
      if s==opts.numSubBatches
        batchStart = t + (labindex-1) + opts.batchSize ;
        batchEnd = min(t+2*opts.batchSize-1, numel(subset)) ;
      else
        batchStart = batchStart + numlabs ;
      end
      nextBatch = subset(batchStart : opts.numSubBatches * numlabs : batchEnd) ;
      getBatch(imdb, nextBatch) ;
    end

    if numGpus >= 1
      im = gpuArray(im{1}) ;
    end

    % evaluate CNN
    net.layers{end}.class = labels ;
    if training, dzdy = one; else, dzdy = [] ; end
    if s==1, res = []; end
    res = vl_bilinearnn(net, im, dzdy, res, ...
                      'accumulate', s ~= 1, ...
                      'mode', evalMode, ...
                      'conserveMemory', opts.conserveMemory, ...
                      'backPropDepth', opts.backPropDepth, ...
                      'sync', opts.sync, ...
                      'cudnn', opts.cudnn) ;

    % accumulate training errors
    error = sum([error, [...
      sum(double(gather(res(end).x))) ;
      reshape(opts.errorFunction(opts, labels, res),[],1) ; ]],2) ;
    numDone = numDone + numel(batch) ;
  end

  % gather and accumulate gradients across labs
  if training
    if numGpus <= 1
      [net,res] = accumulate_gradients(opts, learningRate, batchSize, net, res) ;
    else
      if isempty(mmap)
        mmap = map_gradients(opts.memoryMapFile, net, res, numGpus) ;
      end
      write_gradients(mmap, net, res) ;
      labBarrier() ;
      [net,res] = accumulate_gradients(opts, learningRate, batchSize, net, res, mmap) ;
    end
  end

  % print learning statistics

  time = toc(start) ;
  stats = sum([stats,[0 ; error]],2); % works even when stats=[]
  stats(1) = time ;
  n = (t + batchSize - 1) / max(1,numlabs) ;
  speed = n/time ;
  fprintf('%.1f Hz%s\n', speed) ;

  fprintf(' obj:%.3g', stats(2)/n) ;
  for i=1:numel(opts.errorLabels)
    fprintf(' %s:%.3g', opts.errorLabels{i}, stats(i+2)/n) ;
  end
  fprintf(' [%d/%d]', numDone, batchSize);
  fprintf('\n') ;

  % debug info
  if opts.plotDiagnostics && numGpus <= 1
    figure(2) ; vl_simplenn_diagnose(net,res) ; drawnow ;
  end
end

if nargout > 2
  prof = mpiprofile('info');
  mpiprofile off ;
end

if numGpus >= 1
  net_cpu = vl_simplenn_move(net, 'cpu') ;
else
  net_cpu = net ;
end

% -------------------------------------------------------------------------
function [net,res] = accumulate_gradients(opts, lr, batchSize, net, res, mmap)
% -------------------------------------------------------------------------
if nargin >= 6
  numGpus = numel(mmap.Data) ;
else
  numGpus = 1 ;
end

for l=numel(net.layers):-1:1
  for j=1:numel(res(l).dzdw)

    % accumualte from multiple labs (GPUs) if needed
    if numGpus > 1
      tag = sprintf('l%d_%d',l,j) ;
      tmp = zeros(size(mmap.Data(labindex).(tag)), 'single') ;
      for g = setdiff(1:numGpus, labindex)
        tmp = tmp + mmap.Data(g).(tag) ;
      end
      res(l).dzdw{j} = res(l).dzdw{j} + tmp ;
    end

    if j == 3 && strcmp(net.layers{l}.type, 'bnorm')
      % special case for learning bnorm moments
      thisLR = net.layers{l}.learningRate(j) ;
      net.layers{l}.weights{j} = ...
        (1-thisLR) * net.layers{l}.weights{j} + ...
        (thisLR/batchSize) * res(l).dzdw{j} ;
    else
      % standard gradient training
      thisDecay = opts.weightDecay * net.layers{l}.weightDecay(j) ;
      thisLR = lr * net.layers{l}.learningRate(j) ;
      net.layers{l}.momentum{j} = ...
        opts.momentum * net.layers{l}.momentum{j} ...
        - thisDecay * net.layers{l}.weights{j} ...
        - (1 / batchSize) * res(l).dzdw{j} ;
      net.layers{l}.weights{j} = net.layers{l}.weights{j} + ...
        thisLR * net.layers{l}.momentum{j} ;
    end
  end
end

% -------------------------------------------------------------------------
function mmap = map_gradients(fname, net, res, numGpus)
% -------------------------------------------------------------------------
format = {} ;
for i=1:numel(net.layers)
  for j=1:numel(res(i).dzdw)
    format(end+1,1:3) = {'single', size(res(i).dzdw{j}), sprintf('l%d_%d',i,j)} ;
  end
end
format(end+1,1:3) = {'double', [3 1], 'errors'} ;
if ~exist(fname) && (labindex == 1)
  f = fopen(fname,'wb') ;
  for g=1:numGpus
    for i=1:size(format,1)
      fwrite(f,zeros(format{i,2},format{i,1}),format{i,1}) ;
    end
  end
  fclose(f) ;
end
labBarrier() ;
mmap = memmapfile(fname, 'Format', format, 'Repeat', numGpus, 'Writable', true) ;

% -------------------------------------------------------------------------
function write_gradients(mmap, net, res)
% -------------------------------------------------------------------------
for i=1:numel(net.layers)
  for j=1:numel(res(i).dzdw)
    mmap.Data(labindex).(sprintf('l%d_%d',i,j)) = gather(res(i).dzdw{j}) ;
  end
end

% -------------------------------------------------------------------------
function epoch = findLastCheckpoint(modelDir)
% -------------------------------------------------------------------------
list = dir(fullfile(modelDir, 'net-epoch-*.mat')) ;
tokens = regexp({list.name}, 'net-epoch-([\d]+).mat', 'tokens') ;
epoch = cellfun(@(x) sscanf(x{1}{1}, '%d'), tokens) ;
epoch = max([epoch 0]) ;

repo: DBCNN


text:
function y = vl_nntukeyloss(x, t, varargin)
% VL_NNTUKEYLOSS computes the Tukey Loss
%    This is a slight modification of the robust loss function
%    contained in the deep regression codebase
%    https://github.com/bazilas/matconvnet-deepReg
%    and described in the paper:
%
%    Robust Optimization for Deep Regression
%    V. Belagiannis, C. Rupprecht, G. Carneiro, and N. Navab,
%    ICCV 2015, Santiago de Chile.
%    Copyright (C) 2016 Visual Geometry Group, University of Oxford.
%
% (modified by Samuel Albanie, 2017)

  opts.scaleRes = 1 ;
  opts.adaptInliers = 0 ;
  opts.instanceWeights = ones(size(x)) ;
  [opts, dzdy] = vl_argparsepos(opts, varargin, 'nonrecursive') ;

  % residuals
  res = x - t ;

  % Median absolute deviation (MAD)
  MAD = 1.4826 * mad(res, 1) ;
  C = 4.685 ;

  % inliers (percentage of inliers)
  nonZer = round(100 * sum(abs(res(:)) < C) / numel(res)) ;

  % Big V says that this sometimes helps the convergence
  if opts.adaptInliers % as in the paper
    if nonZer < 70, MAD = MAD * opts.scaleRes ; end
  end

  res = bsxfun(@rdivide, res, MAD) ;

  if isempty(dzdy)
    scale = (C^2) / 6 ;
    yt = scale * (1 - (1 - (res ./ C).^2).^3) ;
    yt(abs(res) > C) = scale ;
    y = opts.instanceWeights(:)' * yt(:) ;
  else
    keep = boolean(abs(res) < C) ;
    tukDer = res .* ((1 - (res ./ C).^2).^2) ;
    res = tukDer * dzdy{1} .* keep .* opts.instanceWeights ;
    y = bsxfun(@rdivide, res, MAD) ;
  end

repo: MFT


text:
% test for radial basis function interpolation in 1D.

% number of points
m = 20;

n = 2024;

scaling = 2;
x = rand(m,1)*(scaling-1)+1;
x = rescale(cumsum(x));

xi = linspace(0,1,n)';

% random values
d = rescale( rand(m,1), .1,.9 );

type = {'abs', 'poly3', 'sqrt', 'gwenn'};

clf;
for i=1:length(type)

    options.rbf = type{i};
    v = perform_rbf_interpolation(xi,x,d, options);
    subplot(length(type),1,i);
    hold on;
    plot(xi, v); plot(x, d, 'r.');
    axis([0 1 0 1]);
    hold off;
    title(['RBF ' type{i}]);

end

author: gpeyre

repo: matlab-toolboxes


text:
function epitrochoid(r1,r2,len)
%animation of EPITROCHOID
%epitrochoid(r1,r2,len);
%r1=radius of fixed circle on which smaller circle rolls
%r2=radius of rolling circle
%len=line connected to centre of rolling circle
%example: epitrochoid(1,0.25,0.5);
%         epitrochoid(1,0.25,1);
%this code also trace epicycloid
%when r2=len then it becomes epicycloid
%example for epicycloid: epitrochoid(1,0.25,0.25);
%always keep r1>r2
if nargin==0
    r1=1;          %radius of fixed circle on which smaller circle rolls
    r2=0.25;        %radius of rolling circle
    len=0.5;       %line connected to centre of rolling circle
end
n=0;
dtheta=pi/20;  %increase in angle per iteration
roughx=[];     %variable to trace trochoid
roughy=[];     %variable to trace trochoid
upto=2*pi*r1/r2+dtheta; %rotate upto
t=linspace(0,2*pi,75);

[x,y]=pol2cart(t,r2);
h1=patch(x,y,'r');%draw rolling circle

[phil,rl]=cart2pol([-len 0],[0 0]);
h2=line(0,0);%draw line connected to centre of rolling circle

[x_fix,y_fix]=pol2cart(t,r1);
line(x_fix,y_fix);%draw fixed circle,only once because it is fixed

axis equal
axis(r1*[-2 2 -2 2]) %fix axis

while n<upto && ishandle(h1)
     [xl,yl]=pol2cart(phil-n,rl); %rotate line
     [x2,y2]=transform2d(xl,yl,0,-(r1+r2),-n*r2/r1,0,r1+r2);
     set(h2,'xdata',x2,'ydata',y2);%draw line

     [x,y]=pol2cart(t,r2); %rotate circle
     [x3,y3]=transform2d(x,y,0,-(r1+r2),-n*r2/r1,0,r1+r2);
     set(h1,'xdata',x3,'ydata',y3);%draw rolling circle

     line([roughx x2(1)],[roughy y2(1)]);%trace trochoid
     roughx=x2(1);
     roughy=y2(1);

     pause(0.025)
     n=n+dtheta;
end
 %function end