text:
function train_id_net_vgg16(varargin)
% -------------------------------------------------------------------------
% Part 4.1: prepare the data
% -------------------------------------------------------------------------

imdb = load('./dataset/MSCOCO-prepare/url_data.mat');
imdb = imdb.imdb;
load('./dataset/MSCOCO-prepare/coco_word2.mat');
%sort row
[imdb.images.label2,index] = sort(imdb.images.label2);
wordcnn = wordcnn(:,index);
imdb.charcnn = wordcnn;
%imdb.charmean = mean(imdb.charcnn(:,:,:,imdb.images.set==1),4);
% -------------------------------------------------------------------------
% Part 4.2: initialize a CNN architecture
% -------------------------------------------------------------------------
net = coco_word2_pool();
net.conserveMemory = true;
im_mean = imdb.rgbMean;
net.meta.normalization.averageImage = im_mean;
%net.meta.normalization.charmean = imdb.charmean;
% -------------------------------------------------------------------------
% Part 4.3: train and evaluate the CNN

% -------------------------------------------------------------------------
opts.train.averageImage = net.meta.normalization.averageImage;
opts.train.batchSize = 32;
opts.train.continue = true;
opts.train.gpus = 2;
opts.train.prefetch = false ;
opts.train.nesterovUpdate = true ;
opts.train.expDir = './data/res52_coco_batch32_pool_shift_both_drop0.5_no_w2v';
opts.train.derOutputs = {'objective_img',1,'objective_txt',1} ;
%opts.train.gamma = 0.9;
opts.train.momentum = 0.9;
%opts.train.constraint = 100;
opts.train.learningRate = [0.1*ones(1,90)] ;
opts.train.weightDecay = 0.0001;
opts.train.numEpochs = numel(opts.train.learningRate) ;
[opts, ~] = vl_argparse(opts.train, varargin) ;
% Call training function in MatConvNet
[net,info] = cnn_train_dag(net, imdb, @getBatch,opts) ;

% --------------------------------------------------------------------
function inputs = getBatch(imdb,batch,opts)
% --------------------------------------------------------------------
%-- img data
im_url = imdb.images.data(batch) ;
im = vl_imreadjpeg(im_url,'Pack','Resize',[224,224],'Flip',...
    'CropLocation','random','CropSize',[0.8,1],...
    'Interpolation', 'bicubic','NumThreads',16,... %'Brightness', double(0.1*imdb.rgbCovariance),...
    'SubtractAverage',imdb.rgbMean,...
    'CropAnisotropy',[3/4,4/3]);
oim = im{1}; %bsxfun(@minus,im{1},opts.averageImage);
label_img =  imdb.images.label(batch);

%-- txt data
batchsize = numel(batch);
txt_batch = zeros(1,batchsize);
for i=1:batchsize
  txt_batch(i) = rand_same_class_coco(imdb,label_img(i));  % train txt range 1~68126
end
%label_txt =  imdb.images.label2(txt_batch);
label_txt = label_img;
txt = single(imdb.charcnn(:,txt_batch));
txtinput = zeros(1,32,29972,batchsize,'single');
for i=1:batchsize
    len = sum(txt(:,i)>0);
    location = randi(33-len);
    for j=1:len
        v = txt(j,i);
       txtinput(1,location,v,i)=1;
       location = location+1;
    end
end
txtinput = gpuArray(txtinput);
%}
%--
inputs = {'data',gpuArray(oim),'data2',txtinput,'label_img',label_img,'label_txt',label_txt};

author: layumi

repo: Image-Text-Embedding


text:
function[fs]=morsespace(varargin)
%MORSESPACE  Logarithmically-spaced frequencies for generalized Morse wavelets.
%
%   F=MORSESPACE(GAMMA,BETA,N) generates a frequency array for the
%   generalized Morse wavelet transform of an N-point time series. The
%   wavelets are specified by GAMMA and BETA.
%
%   LOG(F) is uniformly spaced, following convention for wavelet analysis.
%
%   F has units of *radians* per sample point, and is a column vector with
%   the frequencies arranged in decending order.
%
%   In this usage, the frequencies F are determined using default settings,
%   described below, which should be appropriate for most applications.
%
%   Additional control over the frequency array F can be obtained using the
%   following alternate usages.
%
%   For details on this calculation, see Lilly (2016).
%   __________________________________________________________________
%
%   High- and low-frequency specification
%
%   F=MORSESPACE(GAMMA,BETA,HIGH,LOW) explicitly sets the high-frequency
%   and low-frequency cutoffs for the frequency array.
%
%   The first (largest) value of F is then just smaller than HIGH and the
%   smallest is just larger than LOW.
%
%   HIGH and LOW have units of *radian* per unit sample point.
%   __________________________________________________________________
%
%   High-frequency cutoff
%
%   The highest frequency can be set to be the minimum of a specified value
%   and a cutoff frequency based on a Nyquist overlap condition.
%
%   F=MORSESPACE(GAMMA,BETA,{ETA,HIGH},LOW) sets the highest frequency
%   to be the minimum of the specified value HIGH, and the largest
%   frequency for which the wavelet will satisfy the threshold level ETA.
%
%   Here ETA be a number between zero and one specifying the ratio of a
%   frequency-domain wavelet at the Nyquist frequency to its peak value.
%
%   Note that in this usage, {ETA,HIGH} is a cell array with two entries.
%
%   The simplified usage F=MORSESPACE(GAMMA,BETA,N) corresponds to the
%   choice ETA=0.1, so that by default, the highest-frequency wavelet
%   will decay to at least 10% of its peak value at the Nyquist frequency.
%   __________________________________________________________________
%
%   Low-frequency cutoff
%
%   The lowest frequency can be set to a cutoff frequency based on an
%   endpoint overlap condition.
%
%   F=MORSESPACE(GAMMA,BETA,HIGH,{P,N}) sets the lowest frequency such that
%   the lowest-frequency wavelet will reach some number P times its central
%   window width at the ends of the time series.
%
%   P is called the packing number.  A choice of P=1 corresponds to
%   roughly 95% of the time-domain wavelet energy being contained within
%   the time series endpoints for a wavelet at the center of the domain.
%
%   F=MORSESPACE(GAMMA,BETA,HIGH,{P,N,LOW}) alternately chooses the maximum
%   of the R-level cutoff frequency, and a specified low frequency LOW.
%
%   The simplified usage F=MORSESPACE(GAMMA,BETA,N) corresponds to the
%   default value P=5.  At the lowest frequency, five wavelets will then
%   fit into the time series, with 5% energy overlap between them.
%   __________________________________________________________________
%
%   Wavelet density
%
%   F=MORSESPACE(GAMMA,BETA,HIGH,LOW,D) controls the number of points in
%   the frequency array through the 'density' D.
%
%   Higher values of D mean more overlap in the frequency domain. The
%   default value of the density is D=4.
%
%   When D=1, the peak of one wavelet is located at the half-power points
%   of the adjacent wavelet. D=4 means that four other wavelets will occur
%   between the peak of one wavelet and its half-power point.
%   __________________________________________________________________
%
%   See also MORSEWAVE, WAVETRANS.
%
%   Usage: f=morsespace(ga,be,N);
%          f=morsespace(ga,be,high,low,D);
%          f=morsespace(ga,be,{eta,high},low,D);
%          f=morsespace(ga,be,{eta,high},{p,N},D);
%          f=morsespace(ga,be,{eta,high},{p,N,low},D);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2009--2016 J.M. Lilly --- type 'help jlab_license' for details

%   Note then when HIGH and LOW are explicitly input, the time series
%   length N is no longer input as the first argument.

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

if strcmpi(varargin{1},'--f')
  type makefigs_morsespace;
  makefigs_morsespace
  return
end

%/***************************************************************
%Sort out input arguments
% if ischar(varargin{end})
%     str=varargin{end};
%     varargin=varargin(1:end-1);
% else
%     str='peak';
% end

ga=varargin{1};
be=varargin{2};
D=4;

if length(varargin)==3
    N=varargin{3};

    %Default choices
    low={5,N};
    high={0.1,pi};
else
    high=varargin{3};
    low=varargin{4};
    if length(varargin)==5
        D=varargin{5};
    end
end

%\***************************************************************

if iscell(high)
    %Recall pi is Nyquist
    high=min(high{2},morsehigh(ga,be,high{1}));
    %high
end

if iscell(low)
    if length(low)==2
        low{3}=0;
    end
    low=max(low{3},morsespace_low(ga,be,low{1},low{2}));
end

r=1+frac(1,D*morseprops(ga,be));
N=floor(frac(log(frac(high,low)),log(r)));
fs=high*ones(N+1,1)./r.^[0:N]';

% if length(fs)>1000
%     warning('F is longer than 1000 points... ')
% end

% if strfind(str,'ene')
%      [fm,fe,fi] = morsefreq(ga,be);
%      fs=frac(fe,fm).*fs;
% end

function[fmin]=morsespace_low(ga,be,r,N)

p=morseprops(ga,be);
fmin=frac(2*sqrt(2)*p*r,N);

function[]=morsespace_test
morsespace_hightest;
morsespace_lowtest;

function[]=morsespace_hightest
%Test for high-frequncy cutoff

ga=3;
be=4;
eta=0.1;
ompeak=morsefreq(ga,be);
N=100;

dom=ompeak/N+zeros(4*N,1);
om=cumsum(dom,1);

%morse=frac(1,2)*morseafun(ga,be).*(om.^be).*exp(-om.^ga);

fhigh=morsehigh(ga,be,eta);
[psi,morse] = morsewave(10000,ga,be,fhigh);
reporttest('MORSESPACE high-frequency cutoff',aresame(morse(end/2)./2,eta,1e-3))

function[]=morsespace_lowtest
%Test for low-frequency cutoff

N=1001;
ga1=(1/3:1:11);
be1=(1:1:10);
%ga1=(1/3:.5:11);
%be1=(1:.5:10);

[ga,be]=meshgrid(ga1,be1);
vcolon(ga,be);
[a,sigt,sigo]=morsebox(ga,be);

psi=zeros(N,length(ga));
for i=1:length(ga)
    fs{i}=morsespace(ga(i),be(i),{0.95,pi},{1,N/10},4);
    psi(:,i)=morsewave(N,ga(i),be(i),fs{i}(end),'energy');
end
%figure,plot(sum(squared(psi((end+1)/2-50:(end+1)/2+50,:)),1))
bool=aresame(median(sum(squared(psi((end+1)/2-50:(end+1)/2+50,:)))),0.95,0.01);

reporttest('MORSESPACE low-frequency cutoff, P approximates 95% energy',bool)

%fs=morsespace(2,2,10000);
%psi=morsewave(10000,2,2,fs(end),'energy');
%Can make a figure with the wavelet shifted by 2000 points each time

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function [node,elem] = cvtuniformmesh(n)
% n is the number of intervales in one coordinate direction

h = 2/n;
[x y z] = meshgrid(-1:h:1,-1:h:1,-1:h:1);
[cx cy cz] = meshgrid(-1+h/2:h:1-h/2,-1+h/2:h:1-h/2,-1+h/2:h:1-h/2);
node(:,1) = [x(:); cx(:)];
node(:,2) = [y(:); cy(:)];
node(:,3) = [z(:); cz(:)];
elem = delaunayn(node);
elem = fixorder3(node,elem);

repo: ifem


text:
f=[-2; -3; 5];
a=[-2,5,-1;1,3,1]; b=[-10;12];
aeq=[1,1,1];
beq=7;
[x,y]=linprog(f,a,b,aeq,beq,zeros(3,1));
x, y=-y

author: Eurus-Holmes


text:
function [x, cost, info, options] = conjugategradient(problem, x, options)
% Conjugate gradient minimization algorithm for Manopt.
%
% function [x, cost, info, options] = conjugategradient(problem)
% function [x, cost, info, options] = conjugategradient(problem, x0)
% function [x, cost, info, options] = conjugategradient(problem, x0, options)
% function [x, cost, info, options] = conjugategradient(problem, [], options)
%
% Apply the conjugate gradient minimization algorithm to the problem
% defined in the problem structure, starting at x0 if it is provided
% (otherwise, at a random point on the manifold). To specify options whilst
% not specifying an initial guess, give x0 as [] (the empty matrix).
%
% The outputs x and cost are the best reached point on the manifold and its
% cost. The struct-array info contains information about the iterations:
%   iter : the iteration number (0 for the initial guess)
%   cost : cost value
%   time : elapsed time in seconds
%   gradnorm : Riemannian norm of the gradient
%   stepsize : norm of the last tangent vector retracted
%   beta : value of the beta parameter (see options.beta_type)
%   linesearch : information logged by options.linesearch
%   And possibly additional information logged by options.statsfun.
% For example, type [info.gradnorm] to obtain a vector of the successive
% gradient norms reached.
%
% The options structure is used to overwrite the default values. All
% options have a default value and are hence optional. To force an option
% value, pass an options structure with a field options.optionname, where
% optionname is one of the following and the default value is indicated
% between parentheses:
%
%   tolgradnorm (1e-6)
%       The algorithm terminates if the norm of the gradient drops below this.
%   maxiter (1000)
%       The algorithm terminates if maxiter iterations have been executed.
%   maxtime (Inf)
%       The algorithm terminates if maxtime seconds elapsed.
%   minstepsize (1e-10)
%       The algorithm terminates if the linesearch returns a displacement
%       vector (to be retracted) smaller in norm than this value.
%   beta_type ('H-S')
%       Conjugate gradient beta rule used to construct the new search
%       direction, based on a linear combination of the previous search
%       direction and the new (preconditioned) gradient. Possible values
%       for this parameter are:
%           'S-D', 'steep' for beta = 0 (preconditioned steepest descent)
%           'F-R' for Fletcher-Reeves's rule
%           'P-R' for Polak-Ribiere's modified rule
%           'H-S' for Hestenes-Stiefel's modified rule
%           'H-Z' for Hager-Zhang's modified rule
%           'L-S' for Sato's Liu-Storey rule
%       See Hager and Zhang 2006, "A survey of nonlinear conjugate gradient
%       methods" for a description of these rules in the Euclidean case and
%       for an explanation of how to adapt them to the preconditioned case.
%       The adaption to the Riemannian case is straightforward: see in code
%       for details. Modified rules take the max between 0 and the computed
%       beta value, which provides automatic restart, except for H-Z and L-S
%       which use a different modification. Sato's Liu-Storey rule is
%       described in Sato 2021, "Riemannian conjugate gradient methods:
%       General framework and specific algorithms with convergence analyses"
%   orth_value (Inf)
%       Following Powell's restart strategy (Math. prog. 1977), restart CG
%       (that is, make a -preconditioned- gradient step) if two successive
%       -preconditioned- gradients are "too" parallel. See for example
%       Hager and Zhang 2006, "A survey of nonlinear conjugate gradient
%       methods", page 12. An infinite value disables this strategy. See in
%       code formula for the specific criterion used.
%   linesearch (@linesearch_adaptive or @linesearch_hint)
%       Function handle to a line search function. The options structure is
%       passed to the line search too, so you can pass it parameters. See
%       each line search's documentation for info. Another available line
%       search in manopt is @linesearch, in /manopt/linesearch/linesearch.m
%       If the problem structure includes a line search hint, then the
%       default line search used is @linesearch_hint.
%   statsfun (none)
%       Function handle to a function that will be called after each
%       iteration to provide the opportunity to log additional statistics.
%       They will be returned in the info struct. See the generic Manopt
%       documentation about solvers for further information.
%   stopfun (none)
%       Function handle to a function that will be called at each iteration
%       to provide the opportunity to specify additional stopping criteria.
%       See the generic Manopt documentation about solvers for further
%       information.
%   verbosity (3)
%       Integer number used to tune the amount of output the algorithm
%       generates during execution (mostly as text in the command window).
%       The higher, the more output. 0 means silent.
%   storedepth (2)
%       Maximum number of different points x of the manifold for which a
%       store structure will be kept in memory in the storedb. If the
%       caching features of Manopt are not used, this is irrelevant. For
%       the CG algorithm, a store depth of 2 should always be sufficient.
%
%
% In most of the examples bundled with the toolbox (see link below), the
% solver can be replaced by the present one if need be.
%
% See also: steepestdescent trustregions manopt/solvers/linesearch manopt/examples

% An explicit, general listing of this algorithm, with preconditioning,
% can be found in the following paper:
%     @Article{boumal2015lowrank,
%       Title   = {Low-rank matrix completion via preconditioned optimization on the {G}rassmann manifold},
%       Author  = {Boumal, N. and Absil, P.-A.},
%       Journal = {Linear Algebra and its Applications},
%       Year    = {2015},
%       Pages   = {200--239},
%       Volume  = {475},
%       Doi     = {10.1016/j.laa.2015.02.027},
%     }

% This file is part of Manopt: www.manopt.org.
% Original author: Bamdev Mishra, Dec. 30, 2012.
% Contributors: Nicolas Boumal, Nick Vannieuwenhoven
% Change log:
%
%   March 14, 2013, NB:
%       Added preconditioner support : see Section 8 in
%       https://www.math.lsu.edu/~hozhang/papers/cgsurvey.pdf
%
%   Sept. 13, 2013, NB:
%       Now logging beta parameter too.
%
%   Nov. 7, 2013, NB:
%       The search direction is no longer normalized before it is passed
%       to the linesearch. This way, it is up to the designers of the
%       linesearch to decide whether they want to use the norm of the
%       search direction in their algorithm or not. There are reasons
%       against it, but practical evidence that it may help too, so we
%       allow it. The default linesearch_adaptive used does exploit the
%       norm information. The base linesearch does not. You may select it
%       by setting options.linesearch = @linesearch;
%
%   Nov. 29, 2013, NB:
%       Documentation improved: options are now explicitly described.
%       Removed the Daniel rule for beta: it was not appropriate for
%       preconditioned CG and I could not find a proper reference for it.
%
%   April 3, 2015 (NB):
%       Works with the new StoreDB class system.
%
%   Aug. 2, 2018 (NB):
%       Now using storedb.remove() to keep the cache lean.
%
%   Feb. 7, 2022 (NV):
%       Added support for Liu-Storey rule (L-S).

M = problem.M;

% Verify that the problem description is sufficient for the solver.
if ~canGetCost(problem)
    warning('manopt:getCost', ...
        'No cost provided. The algorithm will likely abort.');
end
if ~canGetGradient(problem) && ~canGetApproxGradient(problem)
    warning('manopt:getGradient:approx', ...
           ['No gradient provided. Using an FD approximation instead (slow).\n' ...
            'It may be necessary to increase options.tolgradnorm.\n' ...
            'To disable this warning: warning(''off'', ''manopt:getGradient:approx'')']);
    problem.approxgrad = approxgradientFD(problem);
end

% Set local defaults here
localdefaults.minstepsize = 1e-10;
localdefaults.maxiter = 1000;
localdefaults.tolgradnorm = 1e-6;
localdefaults.storedepth = 20;
% Changed by NB : H-S has the "auto restart" property.
% See Hager-Zhang 2005/2006 survey about CG methods.
% The auto restart comes from the 'max(0, ...)', not so much from the
% reason stated in Hager-Zhang I think. P-R also has auto restart.
localdefaults.beta_type = 'H-S';
localdefaults.orth_value = Inf; % by BM as suggested in Nocedal and Wright

% Depending on whether the problem structure specifies a hint for
% line-search algorithms, choose a default line-search that works on
% its own (typical) or that uses the hint.
if ~canGetLinesearch(problem)
    localdefaults.linesearch = @linesearch_adaptive;
else
    localdefaults.linesearch = @linesearch_hint;
end

% Merge global and local defaults, then merge w/ user options, if any.
localdefaults = mergeOptions(getGlobalDefaults(), localdefaults);
if ~exist('options', 'var') || isempty(options)
    options = struct();
end
options = mergeOptions(localdefaults, options);

timetic = tic();

% If no initial point x is given by the user, generate one at random.
if ~exist('x', 'var') || isempty(x)
    x = M.rand();
end

% Create a store database and generate a key for the current x
storedb = StoreDB(options.storedepth);
key = storedb.getNewKey();

% Compute cost-related quantities for x
[cost, grad] = getCostGrad(problem, x, storedb, key);
gradnorm = M.norm(x, grad);
Pgrad = getPrecon(problem, x, grad, storedb, key);
gradPgrad = M.inner(x, grad, Pgrad);

% Iteration counter (at any point, iter is the number of fully executed
% iterations so far)
iter = 0;

% Save stats in a struct array info and preallocate.
stats = savestats();
info(1) = stats;
info(min(10000, options.maxiter+1)).iter = [];

if options.verbosity >= 2
    fprintf(' iter\t               cost val\t    grad. norm\n');
end

% Compute a first descent direction (not normalized)
desc_dir = M.lincomb(x, -1, Pgrad);

% Start iterating until stopping criterion triggers
while true

    % Display iteration information
    if options.verbosity >= 2
        fprintf('%5d\t%+.16e\t%.8e\n', iter, cost, gradnorm);
    end

    % Start timing this iteration
    timetic = tic();

    % Run standard stopping criterion checks
    [stop, reason] = stoppingcriterion(problem, x, options, info, iter+1);

    % Run specific stopping criterion check
    if ~stop && abs(stats.stepsize) < options.minstepsize
        stop = true;
        reason = sprintf(['Last stepsize smaller than minimum '  ...
                          'allowed; options.minstepsize = %g.'], ...
                          options.minstepsize);
    end

    if stop
        if options.verbosity >= 1
            fprintf([reason '\n']);
        end
        break;
    end

    % The line search algorithms require the directional derivative of the
    % cost at the current point x along the search direction.
    df0 = M.inner(x, grad, desc_dir);

    % If we didn't get a descent direction: restart, i.e., switch to the
    % negative gradient. Equivalent to resetting the CG direction to a
    % steepest descent step, which discards the past information.
    if df0 >= 0

        % Or we switch to the negative gradient direction.
        if options.verbosity >= 3
            fprintf(['Conjugate gradient info: got an ascent direction '...
                     '(df0 = %2e), reset to the (preconditioned) '...
                     'steepest descent direction.\n'], df0);
        end
        % Reset to negative gradient: this discards the CG memory.
        desc_dir = M.lincomb(x, -1, Pgrad);
        df0 = -gradPgrad;

    end

    % Execute line search
    [stepsize, newx, newkey, lsstats] = options.linesearch( ...
                   problem, x, desc_dir, cost, df0, options, storedb, key);

    % Compute the new cost-related quantities for newx
    [newcost, newgrad] = getCostGrad(problem, newx, storedb, newkey);
    newgradnorm = M.norm(newx, newgrad);
    Pnewgrad = getPrecon(problem, newx, newgrad, storedb, newkey);
    newgradPnewgrad = M.inner(newx, newgrad, Pnewgrad);

    % Apply the CG scheme to compute the next search direction.
    %
    % This paper https://www.math.lsu.edu/~hozhang/papers/cgsurvey.pdf
    % by Hager and Zhang lists many known beta rules. The rules defined
    % here can be found in that paper (or are provided with additional
    % references), adapted to the Riemannian setting.
    %
    if strcmpi(options.beta_type, 'steep') || ...
       strcmpi(options.beta_type, 'S-D')              % Gradient Descent

        beta = 0;
        desc_dir = M.lincomb(newx, -1, Pnewgrad);

    else

        oldgrad = M.transp(x, newx, grad);
        orth_grads = M.inner(newx, oldgrad, Pnewgrad) / newgradPnewgrad;

        % Powell's restart strategy (see page 12 of Hager and Zhang's
        % survey on conjugate gradient methods, for example)
        if abs(orth_grads) >= options.orth_value
            beta = 0;
            desc_dir = M.lincomb(x, -1, Pnewgrad);

        else % Compute the CG modification

            desc_dir = M.transp(x, newx, desc_dir);

            switch upper(options.beta_type)

                case 'F-R'  % Fletcher-Reeves
                    beta = newgradPnewgrad / gradPgrad;

                case 'P-R'  % Polak-Ribiere+
                    % vector grad(new) - transported grad(current)
                    diff = M.lincomb(newx, 1, newgrad, -1, oldgrad);
                    ip_diff = M.inner(newx, Pnewgrad, diff);
                    beta = ip_diff / gradPgrad;
                    beta = max(0, beta);

                case 'H-S'  % Hestenes-Stiefel+
                    diff = M.lincomb(newx, 1, newgrad, -1, oldgrad);
                    ip_diff = M.inner(newx, Pnewgrad, diff);
                    beta = ip_diff / M.inner(newx, diff, desc_dir);
                    beta = max(0, beta);

                case 'H-Z' % Hager-Zhang+
                    diff = M.lincomb(newx, 1, newgrad, -1, oldgrad);
                    Poldgrad = M.transp(x, newx, Pgrad);
                    Pdiff = M.lincomb(newx, 1, Pnewgrad, -1, Poldgrad);
                    deno = M.inner(newx, diff, desc_dir);
                    numo = M.inner(newx, diff, Pnewgrad);
                    numo = numo - 2*M.inner(newx, diff, Pdiff)*...
                                     M.inner(newx, desc_dir, newgrad) / deno;
                    beta = numo / deno;

                    % Robustness (see Hager-Zhang paper mentioned above)
                    desc_dir_norm = M.norm(newx, desc_dir);
                    eta_HZ = -1 / ( desc_dir_norm * min(0.01, gradnorm) );
                    beta = max(beta, eta_HZ);

                case 'L-S' % Liu-Storey+ from Sato
                    diff = M.lincomb(newx, 1, newgrad, -1, oldgrad);
                    ip_diff = M.inner(newx, Pnewgrad, diff);
                    denom = -1*M.inner(x, grad, desc_dir);
                    betaLS = ip_diff / denom;
                    betaCD = newgradPnewgrad / denom;
                    beta = max(0, min(betaLS, betaCD));

                otherwise
                    error(['Unknown options.beta_type. ' ...
                           'Should be steep, S-D, F-R, P-R, H-S, H-Z, or L-S.']);
            end

            desc_dir = M.lincomb(newx, -1, Pnewgrad, beta, desc_dir);

        end

    end

    % Transfer iterate info.
    storedb.removefirstifdifferent(key, newkey);
    x = newx;
    key = newkey;
    cost = newcost;
    grad = newgrad;
    Pgrad = Pnewgrad;
    gradnorm = newgradnorm;
    gradPgrad = newgradPnewgrad;

    % iter is the number of iterations we have accomplished.
    iter = iter + 1;

    % Make sure we don't use too much memory for the store database.
    storedb.purge();

    % Log statistics for freshly executed iteration.
    stats = savestats();
    info(iter+1) = stats;

end

info = info(1:iter+1);

if options.verbosity >= 1
    fprintf('Total time is %f [s] (excludes statsfun)\n', info(end).time);
end

% Routine in charge of collecting the current iteration stats
function stats = savestats()
    stats.iter = iter;
    stats.cost = cost;
    stats.gradnorm = gradnorm;
    if iter == 0
        stats.stepsize = nan;
        stats.time = toc(timetic);
        stats.linesearch = [];
        stats.beta = 0;
    else
        stats.stepsize = stepsize;
        stats.time = info(iter).time + toc(timetic);
        stats.linesearch = lsstats;
        stats.beta = beta;
    end
    stats = applyStatsfun(problem, x, storedb, key, options, stats);
end

end

author: NicolasBoumal

repo: manopt


text:
function pC = prtUtilConfusion2PercentCorrect(confusionMat)
% prtUtilConfusion2PercentCorrect  Calculates the percent correct from a confunsion matrix.
%   This is done by summing the percentages allong the diagonal. If the
%   confusion matrix lists number of responses and not percentages, this
%   this is accounted for.
%
% Syntax:  pC = prtUtilConfusion2PercentCorrect(confusionMat)
%
% Inputs:
%   confusionMat - nClass x nClass matrix listing the number of responses
%       for a given truth. Truths are listed vertically moving downward.
%       Responses are listed horizontally moving right.
%
% Outputs:
%   pC - The percent correct for the confusion matrix
%
% Example 1:
% Example 2:
%   confusionMat = [4 1 0 0; 0 2 1 0; 0 1 2 1; 0 1 2 3];
%   prtUtilConfusion2PercentCorrect(confusionMat)
%
% Other m-files required: none
% Subfunctions: none
% MAT-files required: none
%

if any(~prtUtilIsNaturalNumber(confusionMat(:)))
    error('requires counting number - confusionCountMatrix, not confusionPercentMatrix');
end

assert(ndims(confusionMat)==2 && size(confusionMat,1)==size(confusionMat,2),'prt:prtUtilConfusion2PercentCorrect:BadInput','prtUtilConfusion2PercentCorrect requires a square confusion matrix. A non square confusion matrix possible implies a mismatch between the true targets and the assigned targets. This is ambiguous.')

% The confusion matrix lists the number of responses not percentages.
% We must change the matrix to a percentage matrix.
occurances = repmat(sum(confusionMat,2),1,size(confusionMat,2));

%For normalization, set 0 --> inf; this discounts rows where we had no
%examples in truth
normOccurances = occurances;
normOccurances(occurances == 0) = inf;
confusionMat = confusionMat./normOccurances;

pC = sum(diag(confusionMat).*occurances(:,1))./sum(occurances(:,1));

function is = prtUtilIsNaturalNumber(x)

is = (x == round(x) & x >= 0);

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function ExpFiles = file_expand_selection(FileFilter, Files)
% Expands a list of paths :
%    1) If it is a file => keep it unchanged
%    2) If they are .ds/ directory with .meg4/.res4 files,
%       or ERPCenter directories with .hdr/.erp files,
%       => keep them as "files to open"
%    3) Else : add all the data files they contains (subdirectories included)

% @=============================================================================
% 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, 2010-2017

% Returned variable
ExpFiles = {};
% Get file format
FileFormat = char(FileFilter.getFormatName());
% Loop on all the selected items
for i = 1:length(Files)
    % DIRECTORIES
    if isdir(Files{i})
        isDirFile = 0;
        % If CTF/ERP format
        if ismember(FileFormat, {'CTF', 'EEG-ERPCENTER', 'EEG-NEURALYNX', 'EEG-NEURONE', 'EEG-EGI-MFF'})
            switch (FileFormat)
                case 'CTF'
                    % Check if there are .meg4 or .res4 files in directory
                    dirFiles = [dir(bst_fullfile(Files{i}, '*.meg4'));
                                dir(bst_fullfile(Files{i}, '*.res4'))];
                case 'EEG-ERPCENTER'
                    % Check if there are .hdr file in directory
                    dirFiles = dir(bst_fullfile(Files{i}, '*.hdr'));
                case 'EEG-NEURALYNX'
                    % Check if there are .ncs or .nse files in directory
                    dirFiles = [dir(bst_fullfile(Files{i}, '*.ncs'));
                                dir(bst_fullfile(Files{i}, '*.nse'))];
                case 'EEG-NEURONE'
                    % Get binary files 1-9.bin
                    dirFiles = [dir(bst_fullfile(Files{i}, '1.bin')); dir(bst_fullfile(Files{i}, '2.bin')); dir(bst_fullfile(Files{i}, '3.bin')); dir(bst_fullfile(Files{i}, '4.bin')); dir(bst_fullfile(Files{i}, '5.bin'));
                                dir(bst_fullfile(Files{i}, '6.bin')); dir(bst_fullfile(Files{i}, '7.bin')); dir(bst_fullfile(Files{i}, '8.bin')); dir(bst_fullfile(Files{i}, '9.bin'))];
                case 'EEG-EGI-MFF'
                    % Check if there are .bin files in directory
                    dirFiles = dir(bst_fullfile(Files{i}, '*.bin'));
            end
            j = 1;
            while (~isDirFile && (j <= length(dirFiles)))
                if FileFilter.accept(java.io.File(bst_fullfile(Files{i}, dirFiles(j).name)))
                    isDirFile = 1;
                    ExpFiles = cat(2, ExpFiles, Files{i});
                else
                    j = j + 1;
                end
            end
        end

        % If directory is not a CTF/ERP dir to be opened
        if ~isDirFile
            % Get all files in this directory
            dirFiles = dir(bst_fullfile(Files{i}, '*'));
            % Build a new files list
            dirFullFiles = {};
            for j = 1:length(dirFiles)
                % Skip all the files that start with a '.'
                if (dirFiles(j).name(1) == '.')
                    continue;
                end
                % Skip all the 4D system files
                if strcmpi(FileFormat, '4D') && ismember(dirFiles(j).name, {'config', 'hs_file'})
                    continue;
                end
                dirFullFiles{end + 1} = bst_fullfile(Files{i}, dirFiles(j).name);
            end
            % And call function recursively
            ExpFiles = cat(2, ExpFiles, file_expand_selection(FileFilter, dirFullFiles));
        end

    % SINGLE FILE
    else
        % Test if is accepted by FileFilter
        if FileFilter.accept(java.io.File(Files{i}))
            % If single file is CTF : get the parent
            if strcmpi(FileFormat, 'CTF') || strcmpi(FileFormat, 'EEG-NEURALYNX')
                ExpFiles{end+1} = bst_fileparts(Files{i});
            else
                ExpFiles{end+1} = Files{i};
            end
        end
    end
end

author: brainstorm-tools


text:
function [dat] = read_nexstim_nxe(filename, begsample, endsample, chanindx)

% READ_NEXSTIM_NXE reads specified samples from a NXE continous datafile
%
% Use as
%   [hdr] = read_nexstim_nxe(filename)
% where
%    filename        name of the datafile, including the .bdf extension
% This returns a header structure with the following elements
%   hdr.Fs           sampling frequency
%   hdr.nChans       number of channels
%   hdr.nSamples     number of samples per trial
%   hdr.nSamplesPre  number of pre-trigger samples in each trial
%   hdr.nTrials      number of trials
%   hdr.label        cell-array with labels of each channel
%
% Or use as
%   [dat] = read_nexstim_nxe(filename, begsample, endsample, chanindx)
% where
%    filename        name of the datafile, including the .nxe extension
%    begsample       index of the first sample to read
%    endsample       index of the last sample to read
%    chanindx        index of channels to read (optional, default is all)
% This returns a Nchans X Nsamples data matrix

% Written by Vladimir Litvak based on functions provided by Nexstim
%
% Copyright (C) 2007, Vladimir Litvak
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% 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: read_nexstim_nxe.m 945 2010-04-21 17:41:20Z roboos $

if nargin==1
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % read the header
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    hdr.Fs          = 1450;
    hdr.nChans      = 64;
    hdr.label  = cell(64,1);
    hdr.label(1:4)= {'GATE', 'TRIG1', 'TRIG2','EOG'};
    hdr.label(5:64) = {
        'Fp1'
        'Fpz'
        'Fp2'
        'AF1'
        'AFz'
        'AF2'
        'F7'
        'F3'
        'F1'
        'Fz'
        'F2'
        'F4'
        'F8'
        'FT9'
        'FT7'
        'FC5'
        'FC3'
        'FC1'
        'FCz'
        'FC2'
        'FC4'
        'FC6'
        'FT8'
        'FT10'
        'T7'
        'C5'
        'C3'
        'C1'
        'Cz'
        'C2'
        'C4'
        'C6'
        'T8'
        'TP9'
        'TP7'
        'CP5'
        'CP3'
        'CP1'
        'CPz'
        'CP2'
        'CP4'
        'CP6'
        'TP8'
        'TP10'
        'P9'
        'P7'
        'P3'
        'P1'
        'Pz'
        'P2'
        'P4'
        'P8'
        'P10'
        'PO3'
        'POz'
        'PO4'
        'O1'
        'Oz'
        'O2'
        'Iz'};

    % it is continuous data, therefore append all records in one trial
    hdr.nTrials     = 1;

    fid=fopen(filename,'r','l');
    fseek(fid,0,'eof');
    numBytes = ftell(fid);
    hdr.nSamples = (numBytes/2)/hdr.nChans;

    hdr.nSamplesPre = 0;

    fclose(fid);

    % return the header
    dat = hdr;

else
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % read the data
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    sfEEG = (1/2000) * (10/65535) * 1000000;
    sfEOG = (1/400) * (10/65535) * 1000000;
    sfTRIG = 10 * (10/65535);

    numChannels = 64;

    fid = fopen(filename,'r','l');
    fseek(fid, 2*numChannels*(begsample-1),'bof');
    data = fread(fid,[numChannels endsample-begsample+1],'short');
    fclose(fid);

    data(1:3,:) = sfTRIG.*data(1:3,:);
    data(4,:) = sfEOG.*data(4,:);
    data(5:64,:) = sfEEG.*data(5:64,:);

    if nargin<4
        chanindx = 1:numChannels;
    end

    dat = data(chanindx,:);
end

author: PatternRecognition

repo: OpenBMI


text:
function filum_test06 ( )

%*****************************************************************************80
%
%% FILUM_TEST06 tests FILE_EXIST.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 June 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST06\n' );
  fprintf ( 1, '  FILE_EXIST is true if a file exists.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Exist?    FILENAME\n' );
  fprintf ( 1, '\n' );

  filename = 'filum_test.m';
  fprintf ( 1, '  %d  %s\n', file_exist ( filename ), filename );

  filename = 'filum.html';
  fprintf ( 1, '  %d  %s\n', file_exist ( filename ), filename );

  filename = 'rabbit.txt';
  fprintf ( 1, '  %d  %s\n', file_exist ( filename ), filename );

  filename = 'make.money.fast';
  fprintf ( 1, '  %d  %s\n', file_exist ( filename ), filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function disk_monte_carlo_test01 ( )

%*****************************************************************************80
%
%% DISK_MONTE_CARLO_TEST01 uses DISK01_SAMPLE with an increasing number of points.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2014
%
%  Author:
%
%    John Burkardt
%
  e_test = [ ...
    0, 0; ...
    2, 0; ...
    0, 2; ...
    4, 0; ...
    2, 2; ...
    0, 4; ...
    6, 0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'DISK_MONTE_CARLO_TEST01\n' );
  fprintf ( 1, '  Use DISK01_SAMPLE to estimate integrals in the unit disk.\n' );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '         N        1              X^2             Y^2        ' );
  fprintf ( 1, '     X^4             X^2Y^2           Y^4             X^6\n' );
  fprintf ( 1, '\n' );

  n = 1;

  while ( n <= 65536 )

    [ x, seed ] = disk01_sample ( n, seed );

    fprintf ( 1, '  %8d', n );
    for j = 1 : 7

      e(1:2) = e_test(1:2,j);

      value = monomial_value ( 2, n, e, x );

      result(j) = disk01_area ( ) * sum ( value(1:n) ) / n;
      fprintf ( 1, '  %14.6g', result(j) );
    end
    fprintf ( 1, '\n' );

    n = 2 * n;

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '     Exact' );
  for j = 1 : 7

    e(1:2) = e_test(1:2,j);

    result(j) = disk01_monomial_integral ( e );
    fprintf ( 1, '  %14.6g', result(j) );
  end

  fprintf ( 1, '\n' );

  return
end


text:
%   SCRIPT TO COMPUTE THE TORQUES AT EACH JOINT FOR DIFFERENT MOTION STATES OF
%   THE ARM.
%

function motor_selection

close all
global robot

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PARAMETERS SECTION
%   Feel free to change the values of q, maximum_speeds and
%   maximum_accels.
%
%   This script tries to allow the student to test any mechanism
%   at the worst case. In this sense, q should be adjusted
%   as the pose where each joint would (statically) be needing a
%   higher torque.
%   The maximum_speeds and maximum_acceleration define a trapezoidal
%   speed profile. This trapezoidal speed is used by most machines
%   to command changes in speed in any of their joints.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

echo on
%LOAD ANY ROBOT YOU WOULD LIKE TO ANALYZE
%EXPERIMENTA CAMBIANDO EL ROBOT Y VIENDO LOS RESULTADOS!
robot = load_robot('ABB', 'IRB2600');
echo off

% robot pose: experiment by changing the pose while observing the different
%             torques at each joint
q=[0 pi/2 -pi/2 0 0 0]; %rad

%Velocidad maxima en cada articulacion en rad/s.
% robot.velmax = [deg2rad(175); %Axis 1, rad/s
%                 deg2rad(175); %Axis 2, rad/s
%                 deg2rad(175); %Axis 3, rad/s
%                 deg2rad(360); %Axis 4, rad/s
%                 deg2rad(360); %Axis 5, rad/s
%                 deg2rad(500)];%Axis 6, rad/s

maximum_speeds=[3.5 3.5 3.5 6.5 6.5 8.5];%rad/second

%maximum acceleration/deceleration for each joint

maximum_accels=maximum_speeds*2; %rad/second^2
%maximum_accels=[5 5 6 7 8 9];

% time of the trapezoidal profile that the joint moves at maximum speed
time_at_constant_speed=0.4; %seconds

%load robot parameters. Just uncomment this line

drawrobot3d(robot, q)

%START BY COMPUTING TORQUES AT G
%    robot.motors.G = [1 1 1 1 1 1 ]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  FIRST, COMPUTE TRAPEZOIDAL PROFILES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%compute acceleration plus deceleration times for every joint
time_acc = 2*maximum_speeds./maximum_accels+time_at_constant_speed;

%compute the total time for the slowest joint
total_time=max(time_acc);

% Trapezoidal speed profiles for each joint
[input_speeds, input_accels, time]=build_trapezoidal_speed_profile(maximum_speeds, maximum_accels, total_time);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FINALLY, COMPUTE TORQUES FOR EACH MOTION STATE.
% Please note that we consider that the robot is placed at a fixed position and consider
% different motion situations when we change the acceleration and speed at
% each joint. For each motion state, the inverse dynamic model returns the
% torques at each joint that would bring the robot to that motion.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
compute_inverse_dynamics(q, input_speeds, input_accels, time);

  fprintf('\n\nOBSERVE THE PLOTS AND NOTE DOWN THE PEAK TORQUE, NOMINAL TORQUE AND MOTOR SPEEDS')
fprintf('\nNOW COMPUTE THE TORQUES FOR 5 DIFFERENT SELECTED MOTIONS STATES')
fprintf('\nPRESS ANY KEY TO CONTINUE...')

pause

% NOW COMPUTE THE WORST CASE CONSIDERING ONLY THE SPEED AND ACCELERATION AT
% 5 MOTIONS STATES
input_speeds = [zeros(6,1) maximum_speeds' maximum_speeds' maximum_speeds' zeros(6,1) ];
input_accels = [maximum_accels' maximum_accels' zeros(6,1) -maximum_accels' -maximum_accels' ];

compute_inverse_dynamics(q, input_speeds, input_accels, [1:5]);

%Computes a trapezoidal speed profile for every joint given maximum
%permitted accelerations and maximum joint speeds
function [input_speeds, input_accelerations, time]=build_trapezoidal_speed_profile(maximum_speeds, maximum_accels, total_time)

delta_time=0.01;

%build time vector: twice acceleration time plus time at constant speed
time = 0:delta_time:total_time;

input_speeds=[];
input_accelerations=[];

for j=1:length(maximum_speeds),
    vel_row=[];
    acc_row=[];
    for i=1:length(time),
        [vel acc] = compute_values(time(i), maximum_speeds(j), maximum_accels(j), total_time);
        vel_row = [vel_row vel];
        acc_row = [acc_row acc];
    end
    input_speeds = [input_speeds; vel_row];
    input_accelerations = [input_accelerations; acc_row];
end

%returns the values of velocity and speed corresponding to a given time
function [vel acc]=compute_values(time_i, vel_max, acc_max, total_time)

tacc = vel_max/acc_max;
tdec = total_time-tacc;

if time_i < tacc
    vel = time_i.*acc_max;
    acc = acc_max;
    return;
elseif (time_i >= tacc) & (time_i < tdec)
    vel = vel_max;
    acc = 0;
    return;
else % time_i> tdec
    vel = vel_max-(time_i-tdec)*acc_max;
    acc = -acc_max;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   COMPUTE THE INVERSE DYNAMICS FOR EACH MOTION STATE
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function compute_inverse_dynamics(q, input_speeds, input_accels, time)
global robot

%adjust_view(robot)
torques=[];
for j=1:length(time),
    fprintf('\nComputing time %d out of %d', j, length(time));
    % compute the torque to bring the robot instantaneously to this motion
    % state. change M=1  to add the effects of a 1kg mass load at the end effector
    M=20;
    %please note that the force due to the load acts on the z axis of
    tau=inversedynamic(robot, q, input_speeds(:,j), input_accels(:,j), [0  0 -9.81]', [ 0  0 -M 0 0 0 ]');
    torques=[torques tau];
end

%plot trapezoidal profiles
figure, hold, xlabel('time (s)'), ylabel('Input reference speeds (rad/s)')
plot(time, input_speeds(1,:), time, input_speeds(2,:), time, input_speeds(3,:),...
        time, input_speeds(4,:), time, input_speeds(5,:), time, input_speeds(6,:));
legend('Speed for joint 1 (qd1)','Speed for joint 2 (qd2)','Speed for joint 3 (qd3)',...
   'Speed for joint 4 (qd4)','Speed for joint 5 (qd5)','Speed for joint 6 (qd6)' )
%plot trapezoidal profiles, acceleration
figure, hold, xlabel('time (s)'), ylabel('Input reference acceleration (rad/s)')
plot(time, input_accels(1,:), time, input_accels(2,:), time, input_accels(3,:),...
        time, input_accels(4,:), time, input_accels(5,:), time, input_accels(6,:));
legend('Acceleration for joint 1 (qd1)','Acceleration for joint 2 (qd2)','Acceleration for joint 3 (qd3)',...
   'Acceleration for joint 4 (qd4)','Acceleration for joint 5 (qd5)','Acceleration for joint 6 (qd6)' )

% plot results. First, torques at each joint
figure, hold, xlabel('time (s)'), ylabel('Join Torques (N m)')
plot(time, torques(1,:), time, torques(2,:), time, torques(3,:),...
        time, torques(4,:), time, torques(5,:), time, torques(6,:));
legend('Torque for joint 1 ','Torque  for joint 2 ','Torque  for joint 3 ',...
   'Torque  for joint 4','Torque  for joint 5 ','Torque  for joint 6 ' )

%plot torques at each motor
figure, hold, xlabel('time (s)'), ylabel('Motor Torques (N m)')
plot(time, torques(1,:)/robot.motors.G(1), time, torques(2,:)/robot.motors.G(2), time, torques(3,:)/robot.motors.G(3),...
        time, torques(4,:)/robot.motors.G(4), time, torques(5,:)/robot.motors.G(5), time, torques(6,:)/robot.motors.G(6));
legend('Torque at motor 1 ','Torque at motor 2 ','Torque at motor 3 ',...
   'Torque at motor 4 ','Torque at motor 5 ','Torque at motor 6 ' )

%plot power needed by the motor at each time step, without considering the
%losses at the gears
figure, hold, xlabel('time (s)'), ylabel('Power needed by each motor (W)')
plot(time, torques(1,:).*input_speeds(1,:), time, torques(2,:).*input_speeds(2,:), time, torques(3,:).*input_speeds(3,:),...
        time, torques(4,:).*input_speeds(4,:), time, torques(5,:).*input_speeds(5,:), time, torques(6,:).*input_speeds(6,:));
legend('Power: motor 1','Power: motor 2','Power: motor 3',...
   'Power: motor 4','Power: motor 5','Power: motor 6' )

%plot motor speed in rpm for each motor
figure, hold, xlabel('time (s)'), ylabel('Speed in r.p.m of every motor (rev/min)')
plot(time, robot.motors.G(1)*input_speeds(1,:)*30/pi, time, robot.motors.G(2)*input_speeds(2,:)*30/pi, time, robot.motors.G(3)*input_speeds(3,:)*30/pi,...
        time, robot.motors.G(4)*input_speeds(4,:)*30/pi, time, robot.motors.G(5)*input_speeds(5,:)*30/pi, time, robot.motors.G(6)*input_speeds(6,:)*30/pi);
legend('Speed at motor 1 (qd1*G)','Speed at motor 2 (qd2*G)','Speed at motor 3 (qd3*G)',...
   'Speed at motor 4 (qd4*G)','Speed at motor 5 (qd5*G)','Speed at motor 6 (qd6*G)' )

%Now present results:
fprintf('\nMAIN RESULTS (referred to each motor): ')
fprintf('\n------------------------------------------------------------------------------------ ')
fprintf('\n                         Joint 1 - Joint 2 - Joint 3 - Joint 4  - Joint 5 - Joint 6: ')
fprintf('\nPeak Torque (N�m):        %.3f     %.3f     %.3f     %.3f      %.3f    %.3f ', max(abs(torques(1,:)/robot.motors.G(1))), max(abs(torques(2,:)/robot.motors.G(2))) , max(abs(torques(3,:)/robot.motors.G(3))) , max(abs(torques(4,:)/robot.motors.G(4))) , max(abs(torques(5,:)/robot.motors.G(5))) , max(abs(torques(6,:)/robot.motors.G(6))))
fprintf('\nNominal Torque (N�m):     %.3f     %.3f     %.3f     %.3f      %.3f    %.3f  ', abs(torques(1,round(length(torques)/2))/robot.motors.G(1)), abs(torques(2,round(length(torques)/2))/robot.motors.G(2)), abs(torques(3,round(length(torques)/2))/robot.motors.G(3))...
    , abs(torques(4,round(length(torques)/2))/robot.motors.G(4)), abs(torques(5,round(length(torques)/2))/robot.motors.G(5)), abs(torques(6,round(length(torques)/2))/robot.motors.G(6)))
fprintf('\nMax motor speed (r.p.m.):   %.1f     %.1f     %.1f     %.1f      %.1f    %.1f  ', max(abs(robot.motors.G(1)*input_speeds(1,:))*30/pi), max(abs(robot.motors.G(2)*input_speeds(2,:))*30/pi), max(abs(robot.motors.G(3)*input_speeds(3,:))*30/pi)...
    ,max(abs(robot.motors.G(4)*input_speeds(4,:))*30/pi), max(abs(robot.motors.G(5)*input_speeds(5,:))*30/pi), max(abs(robot.motors.G(6)*input_speeds(6,:))*30/pi))
fprintf('\n------------------------------------------------------------------------------------ ')

repo: ARTE


text:
function J2=numDiff2(x,f,fDim,N,epsilon)
%%NUMDIFF2 Numerically compute the second derivatives using central
%          difference formulae based on Lagrange interpolating polynomials.
%          The algorithm is made to alert (and not crash) in the event that
%          f fails. Specifically, if f returns an empty matrix or it any of
%          the components of f are NaNs, then the numDiff function will
%          terminate early, returning an empty matrix to indicate failure.
%          The function only computes second derivatives, not considering
%          cross terms (it does not find a Hessian matrix). Use numDiffHess
%          to get a Hessian matrix with cross terms.
%
%INPUTS: x The xDimX1 vector or scalar point at which the second derivative
%          of the (possibly vector) function is desired.
%        f The scalar or vector function that is to be differentiated. The
%          function f must take x as its parameter and its output should be
%          a scalar or a column vector.
%     fDim The dimensionality of the output of f.
%        N A number >=1 specifying the order of the second derivative
%          approximation. Values for n=1 through 3 are explicitly coded in.
%          For values 3 and above, the coefficients of the second
%          derivative of the Lagrange interpolating polynomial are
%          explicitly solved. If omitted, a value of N=1 is used.
%  epsilon A scalar or xDimX1 vector quantity specifying the finite step
%          size used for numerical differentiation. If a scalar value is
%          given, that value is used for differentiating with respect to
%          elements of xDim. If an xDimX1 value is given, then the
%          corresponding element of epsilon is used to differentiate each
%          element of x. If epsilon is omitted, then
%          epsilon=max(1e-5*x,1e-7); is used.
%
%OUTPUTS: J2 A fDimXxDim matrix of second derivatives. Each column is the
%            derivative vector of f with respect to the corresponding
%            element of x. If at any point the function f returned a NaN or
%            an empty matrix, J2 will be an empty matrix.
%
%The function is similar to the numDiff function. Central-difference
%numerical differentiation is discussed in terms of Lagrange interpolating
%polynomials in Chapter 4.1 of [1]. The Lagrange interpolating polynomials
%themselves are discussed in Chapter 3 of [1]. In general, the coefficients
%of the interpolating polynomial for first derivative central difference
%numerical differentiation come from Equation 4.2, which  expresses them in
%terms of the first derivative of a Lagrange interpolating polynomial. This
%function just extends the concept by using the second-derivative
%coefficients.
%
%REFERENCES:
%[1] R. L. Burden and J. D. Faires, Numerical Analysis, 9th ed. Boston, MA:
%    Brooks/ Cole, 2011.
%
%April 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4)
   N=1;
end

switch(N)
    case 1
        a=[1; -2; 1];
    case 2
        a=[-1/12; 4/3; -5/2; 4/3; -1/12];
    case 3
        a=[1/90; -3/20; 3/2; -49/18; 3/2; -3/20; 1/90];
    otherwise
        [~,~,d2Li]=LagrangeInterpPoly(-N:1:N);
        a=d2Li(end,:);
end

xDim=size(x,1);

if(nargin<5)
    %If epsilon is not specified, then use some ad-hoc default value
    epsilon=max(1e-5*x,1e-7);
end

if(isscalar(epsilon))
   epsilon=repmat(epsilon,[xDim,1]);
end

J2=zeros(fDim,xDim);
for curEl=1:xDim
    epsCur=epsilon(curEl);

    curIdx=1;
    for curP=-N:1:N
        xP=x;
        xP(curEl)=xP(curEl)+curP*epsCur;
        fxP=f(xP);
        if(isempty(fxP)||any(isnan(fxP)))
            J2=[];
            return;
        end
        J2(:,curEl)=J2(:,curEl)+a(curIdx)*fxP;
        curIdx=curIdx+1;
    end
    J2(:,curEl)=J2(:,curEl)/(epsCur^2);
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
gammaa=[];
deltaa=[];
alphaa=[];
for lc=1:25
    gamma=4*pi*(rand-0.5)
    delta=4*pi*(rand-0.5)
    alpha=4*pi*(rand-0.5)
    an=ax_an_bounding(gamma,delta,alpha)
    gamma=an(2);
    delta=an(3);
    alpha=an(4);
    gammaa=[gammaa gamma];
    deltaa=[deltaa delta];
    alphaa=[alphaa alpha];
end

figure;
plot(gammaa,zeros(1,length(gammaa)),'.b');
title('gamma');

figure;
plot(deltaa,zeros(1,length(deltaa)),'.r');
title('delta');

figure;
plot(alphaa,zeros(1,length(alphaa)),'.r');
title('alpha');

author: Sable

repo: mcbench-benchmarks


text:
function move_topo_layers_to_bottom(container)
    % move the topo layers to the bottom of an axes
    h=findobj(container,'Type','surface','-and','-regexp','Tag','topographic_map_.*');
    p=unique([h.Parent]);
    for i=1:numel(p)
        idx=startsWith(get(p(i).Children,'Tag'),'topographic_map_');
        h=p(i).Children(idx);
        others=p(i).Children(~idx);
        p(i).Children=[others;h];
    end

end

author: CelsoReyes


text:
function [ cost, grad, numTotal, pred_cell ] = drdae_discrim_obj( ...
    theta, eI, data_cell, targets_cell, fprop_only, pred_out)
% Copyright (c) 2014-present University of Illinois at Urbana-Champaign
% All rights reserved.
%
% Developed by:     Po-Sen Huang, Paris Smaragdis
%                   Department of Electrical and Computer Engineering
%                   Department of Computer Science
%
% discrim. training (not joint)
%
%PRNN_OBJ MinFunc style objective for Deep Recurrent Denoising Autoencoder
%   theta is the full parameter vector
%   eI contains experiment / network architecture
%   data_cell is a cell array of matrices. Each a distinct length is a cell
%             entry. Each matrix has a time series example in each column
%   targets_cell is parallel to data, but contains the labels for each time
%   fprop_only is a flag that only computes the cost, no gradient
%   numTotal is total number of frames evaluated
%   pred_out is a binary flag for whether pred_cell is populated
%            pred_cell only filled properly when utterances one per cell

%% Debug: Turns this into an identity-function for debugging rest of system
if isfield(eI, 'objReturnsIdentity') && eI.objReturnsIdentity
    cost = 0; grad = 0; numTotal = 0;
    for l = 1:numel(data_cell)
        numUtterances = size(data_cell{l}, 2);
        original_vector = reshape(data_cell{l}, eI.winSize*eI.featDim, []);
        midPnt = ceil(eI.winSize/2);
        original_vector = original_vector((midPnt-1)*14+1 : midPnt*14, :);
        pred_cell{l} = reshape(original_vector, [], numUtterances);
    end
    return;
end

return_activation = 0;

%% Load data from globals if not passed in (happens when run on RPC slave)
global g_data_cell;
global g_targets_cell;
isSlave = false;
if isempty(data_cell)
    data_cell = g_data_cell;
    targets_cell = g_targets_cell;
    isSlave = true;
end;
pred_cell = cell(1,numel(data_cell));
act_cell = cell(1,numel(data_cell));
%% default short circuits to false
if ~isfield(eI, 'shortCircuit')
    eI.shortCircuit = 0;
end;

%% default dropout to false
if ~isfield(eI, 'dropout')
  eI.dropout = 0;
end;

%% setup weights and accumulators
[stack, W_t] = rnn_params2stack(theta, eI);
cost = 0; numTotal = 0;
outputDim = eI.layerSizes(end);
%% setup structures to aggregate gradients
stackGrad = cell(1,numel(eI.layerSizes));
% W_t_grad = zeros(size(W_t));
if isfield(eI, 'fullRNN') && eI.fullRNN==1
   W_t_grad = cell(1,numel(eI.layerSizes)-1);
    for l = 1:numel(eI.layerSizes)-1
        W_t_grad{l}.W = zeros(size(W_t{l}.W));
    end
else
   W_t_grad = zeros(size(W_t));
end

for l = 1:numel(eI.layerSizes)
    stackGrad{l}.W = zeros(size(stack{l}.W));
    stackGrad{l}.b = zeros(size(stack{l}.b));
end
if eI.shortCircuit
    stackGrad{end}.W_ss = zeros(size(stack{end}.W_ss));
end;
%% check options
if ~exist('fprop_only','var')
    fprop_only = false;
end;
if ~exist('pred_out','var')
    pred_out = false;
end;

% DROPOUT: vector of length of hidden layers with 0 or 1
% (to drop or keep activation unit) with prob=0.5
hActToDrop = cell(numel(eI.layerSizes-1),1);
for i=1:numel(eI.layerSizes)-1
 if eI.dropout
    hActToDrop{i} = 1/eI.dropout * binornd(1,eI.dropout, eI.layerSizes(i),1);
 else
   hActToDrop{i} = ones(eI.layerSizes(i),1);
 end
end

%% loop over each distinct length
for c = 1:numel(data_cell)
    data = data_cell{c};
    targets = {};
    if ~isempty(targets_cell), targets = targets_cell{c}; end;
    uttPred = [];
    T =size(data,1) / eI.inputDim;
    % store hidden unit activations at each time instant
    hAct = cell(numel(eI.layerSizes)-1, T);
    for t = 1:T
        %% forward prop all hidden layers
        for l = 1:numel(eI.layerSizes)-1
            if l == 1
                hAct{1,t} = stack{1}.W * data((t-1)*eI.inputDim+1:t*eI.inputDim, :);
            else
                hAct{l,t} = stack{l}.W * hAct{l-1,t};
            end;
            hAct{l,t} = bsxfun(@plus, hAct{l,t}, stack{l}.b);
            % temporal recurrence. limited to single layer for now
            if l == eI.temporalLayer && t > 1
                hAct{l,t} = hAct{l,t} + W_t * hAct{l,t-1};
            end;
            % nonlinearity
            if strcmpi(eI.activationFn,'tanh')
                hAct{l,t} = tanh(hAct{l,t});
            elseif strcmpi(eI.activationFn,'logistic')
                hAct{l,t} = 1./(1+exp(-hAct{l,t}));
            elseif strcmpi(eI.activationFn,'RELU')
                hAct{l,t} = max(0,hAct{l,t});
            else
                error('unrecognized activation function: %s',eI.activationFn);
            end;
            %dropout (hActToDrop will be all ones if no dropout specified)
            hAct{1,t} = bsxfun(@times, hAct{1,t}, hActToDrop{l});
        end;
        % forward prop top layer not done here to avoid caching it
    end;
    %% compute cost and backprop through time
    if  eI.temporalLayer
        delta_t = zeros(eI.layerSizes(eI.temporalLayer),size(data,2));
    end;
    for t = T:-1:1
        l = numel(eI.layerSizes);
        %% forward prop output layer for this timestep
        curPred = bsxfun(@plus, stack{l}.W * hAct{l-1,t}, stack{l}.b);
        % add short circuit to regression prediction if model has it
        if eI.shortCircuit
            curPred = curPred + stack{end}.W_ss ...
                * data((t-1)*eI.inputDim+1:t*eI.inputDim, :);
        end;
        if pred_out, uttPred = [curPred, uttPred]; end;
        % skip loss computation if no targets given
        if isempty(targets), continue; end;

        curTargets = targets((t-1)*outputDim+1:t*outputDim, :);

        if eI.cleanonly==0,
        curTargets_neg = [curTargets(outputDim/2+1:outputDim,:); curTargets(1:outputDim/2,:)];
        end
        if eI.cleanonly==1,
            delta = curPred  - curTargets;
            cost = cost + 0.5 * sum( sum((curPred - curTargets).^2));
        else
            delta = (1- eI.r) * curPred + eI.r * curTargets_neg - curTargets;
            cost = cost + 0.5 * ( sum( sum((curPred - curTargets).^2))-  eI.r* sum( sum((curPred - curTargets_neg).^2)));
        end

        if fprop_only, continue; end;
        %% regression layer gradient and delta
        stackGrad{l}.W = stackGrad{l}.W + delta * hAct{l-1,t}';
        stackGrad{l}.b = stackGrad{l}.b + sum(delta,2);
        % short circuit layer
        if eI.shortCircuit
            stackGrad{end}.W_ss = stackGrad{end}.W_ss + delta ...
                * data((t-1)*eI.inputDim+1:t*eI.inputDim, :)';
        end;
        delta = stack{l}.W' * delta;
        %% backprop through hidden layers
        for l = numel(eI.layerSizes)-1:-1:1
            % aggregate temporal delta term if this is the recurrent layer
            if l == eI.temporalLayer
                delta = delta + delta_t;
            end;
            % push delta through activation function for this layer
            % tanh unit choice assumed
            if strcmpi(eI.activationFn,'tanh')
                delta = delta .* (1 - hAct{l,t}.^2);
            elseif strcmpi(eI.activationFn,'logistic')
                delta = delta .* hAct{l,t} .* (1 - hAct{l,t});
             elseif strcmpi(eI.activationFn,'RELU')
                delta = delta .* double(hAct{l,t}>0);
            else
                error('unrecognized activation function: %s',eI.activationFn);
            end;

            % gradient of bottom-up connection for this layer
            if l > 1
                stackGrad{l}.W = stackGrad{l}.W + delta * hAct{l-1,t}';
            else
                stackGrad{l}.W = stackGrad{l}.W + delta * data((t-1)*eI.inputDim+1:t*eI.inputDim, :)';
            end;
            % gradient for bias
            stackGrad{l}.b = stackGrad{l}.b + sum(delta,2);

            % compute derivative and delta for temporal connections
            if l == eI.temporalLayer && t > 1
                 W_t_grad = W_t_grad + delta * hAct{l,t-1}';
                 % push delta through temporal weights
                 delta_t = W_t' * delta;
            end;
            % push delta through bottom-up weights
            if l > 1
                delta = stack{l}.W' * delta;
            end;
        end
    end
    pred_cell{c} = uttPred;
    % Return the activations for this utterance.
    if return_activation,
      act_cell{c} = cell2mat(hAct);
    end
    % keep track of how many examples seen in total
    numTotal = numTotal + T * size(targets,2);
end

%% stack gradients into single vector and compute weight cost
wCost = numTotal * eI.lambda * sum(theta.^2);
grad = rnn_stack2params(stackGrad, eI, W_t_grad, true);
grad = grad + 2 * numTotal * eI.lambda * theta;

%% clipping
if isfield(eI,'clip') && eI.clip~=0, % if eI.clip==0, no clip
  if eI.clip > 0 % method one -clip the whole
      norm_grad = norm(grad);
      fprintf('norm_grad:%f\n', norm_grad);
      % avoid numerial problem
      if norm_grad <0 || norm_grad > 1e15 || isnan(norm_grad) || isinf(norm_grad),
          grad = zeros(size(grad));
          fprintf('set gradient to zeros\n');
      end
      if norm_grad > eI.clip
         grad = eI.clip * grad/ norm_grad;
      end
  else % method two - clip each entry
      clip_value = -1*eI.clip;
      grad(grad > clip_value)=clip_value;
      grad(grad < -clip_value)=-clip_value;
  end
end

%%
avCost = cost/numTotal;
avWCost = wCost/numTotal;
cost = cost + wCost;

% print output
if ~isSlave && ~isempty(targets_cell)
    fprintf('loss:  %f  wCost:  %f \t',avCost, avWCost);

    if isfield(eI, 'fullRNN') && eI.fullRNN==1
        fprintf('wNorm: %f  rNorm: %f  oNorm: %f\n',sum(stack{1}.W(:).^2),...
            sum(W_t{1}.W(:).^2), sum(stack{end}.W(:).^2));
    else
        fprintf('wNorm: %f  rNorm: %f  oNorm: %f\n',sum(stack{1}.W(:).^2),...
            sum(W_t(:).^2), sum(stack{end}.W(:).^2));
    end
end;

author: posenhuang

repo: deeplearningsourceseparation


text:
  function [x, ei_fun] = embed_in(x, mask, np)
%|function [x, ei_fun] = embed_in(x, mask, np)
%|
%| The matrix-vector versions of the linear forward models used herein
%| expect an input x of size [np *L] and produce an output y of size [nd *L],
%| where np = sum(mask(:)) is the number of estimated pixel values,
%| and nd is the number of data points.
%| However, for convenience it is also desirable for the overloaded
%| mtimes operation to map a [(N) (L)] input into a [(M) (L)] output,
%| where (N) is the d_in -dimensional input size
%| and (M) is the d_out -dimensional output size,
%| and (L) reflects the possibility of multiple inputs, e.g.: A * [u v w].
%|
%| Furthermore, most of the code for forward models is designed to
%| map input [(N) *L] to output [(M) *L].
%|
%| in
%|	x	[np (L)] or [(N) (L)]	input image(s), possibly as columns
%|	mask	[(N)]			logical support array
%|	np				sum(mask(:))
%| out
%|	x	[(N) *L]		output images, as array
%|	ei_fun	strum object with methods:
%|		y = ei_fun.shape(y)	reshape y from [(M) *L] to be either
%|					[nd (L)] or [(M) (L)], depending on x
%|
%| Copyright 2006-12-9, Jeff Fessler, University of Michigan

if nargin == 1 && streq(x, 'test'), embed_in_test, clear x, return, end
if nargin < 2, ir_usage, end
if nargin < 3, np = sum(mask(:)); end

% convert input to [(N) *L]
state.column = false;
if size(x,1) == np % convert [np (L)] to [(N) *L]
	state.column = true;
	dimx = size(x);
	state.diml = dimx(2:end); % [(L)]
	x = embed(x, mask, '*dim', 0); % [(N) *L]

else % convert [(N) (L)] to [(N) *L]
	dimi = size(x);
	if length(dimi) < ndims(mask), fail('bad image size'), end

	size1_mask = size(mask);
	ndims1_mask = ndims(mask);
	if size1_mask(end) == 1 % trick: handle '1d' mask well
		ndims1_mask = ndims1_mask - 1;
		size1_mask = size1_mask(1:end-1);
	end
	jf_equal(dimi(1:ndims1_mask), size1_mask)
%	jf_equal(dimi(1:length(size(mask))), size(mask)) % pre 2012-06-03
	state.diml = dimi((ndims1_mask+1):end); % (L)
%	x = reshape(x, [size(mask) prod(state.diml) 1]); % [(N) *L] % pre
	x = reshape(x, [size1_mask prod(state.diml) 1]); % [(N) *L]

%	if ~dims_same(x, mask, 'up_to_dim', ndims(mask))
%		error(['dimension mismatch.  x: ' num2str(size(x), ' %0d') ...
%			', mask: ' num2str(size(mask), ' %0d')])
%	end
end

ei_fun = strum(state, {'shape', @embed_in_shape, '(y)'});

%
% embed_in_shape()
%
function y = embed_in_shape(state, y)

diml = state.diml;

if state.column % column in yields column out, i.e., [(M) *L] to [*M (L)]
	if any(diml > 1)
		diml = num2cell(diml);
		y = reshape(y, [], diml{:}); % [*M (L)]
	else
		y = y(:); % [*M,1]
	end

else % [(M) *L] to [(M) (L)]
	if any(diml > 1)
		dimy = size(y);
		if dimy(end) ~= prod(diml), error 'diml bug', end
		dimm = dimy(1:end-1); % (M)
		y = reshape(y, [dimm diml]);
	end
end

function embed_in_test
ig = image_geom('nx', 10, 'ny', 8, 'dx', 1);
ig.mask = ig.circ > 0;

x = ig.unitv;
dl = [2 3];
c = repmat(x(ig.mask), [1 dl]); % [np (L)]
[x1 ei] = embed_in(c, ig.mask);
t = ei.shape(ones(4, 5, prod(dl)));
jf_equal(size(t), [4*5 dl])

x2 = repmat(x, [1 1 dl]); % [(N) (L)]
%size(x2)
[x3 ei] = embed_in(x2, ig.mask);
jf_equal(size(x3), [ig.dim prod(dl)])
t = ei.shape(ones(4, 5, prod(dl)));
jf_equal(size(t), [4 5 dl])

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [kernel,rawkernel] = calibrate(obj, AtA, nCoil, coil)

sampling = ones([obj.kernelSize(1:2),nCoil]);

dummyK = zeros(obj.kernelSize(1),obj.kernelSize(2),nCoil); dummyK((end+1)/2,(end+1)/2,coil) = 1;
idxY = find(dummyK);
sampling(idxY) = 0;
idxA = find(sampling);

Aty = AtA(:,idxY); Aty = Aty(idxA); % correlation values to target point, take complete neighbourhood and not just aquired ones
AtA = AtA(idxA,:); AtA =  AtA(:,idxA); % kick out the searched point

kernel = sampling*0;

lambda = norm(AtA,'fro')/size(AtA,1)*obj.calibTyk;

rawkernel = inv(AtA + eye(size(AtA))*lambda)*Aty; % grappa weighting values
kernel(idxA) = rawkernel;

end

author: thomaskuestner


text:
function x = vecpostproc(x, a)
% VECPOSTPROC is post-processing of a D-dimensional vector.
%
%   V = vecpostproc(V) outputs L2 normalized vector:
%     V = V ./ L2NORM(V);
%
%   V = vecpostproc(V, A) outputs L2 and power-law normalized vector:
%     V = SIGN(X) .* ABS(X) .^ A;
%     V = V ./ L2NORM(V);
%
% Authors: F. Radenovic, G. Tolias, O. Chum. 2017.

    if ~exist('a'), a = 1; end
    x = replacenan (l2_normalize (powerlaw (x, a)));

function x = l2_normalize(x)
    l = sqrt(sum(x.^2));
    x = bsxfun(@rdivide,x,l);
    x = replacenan(x);

function x = powerlaw (x, a)
	if a == 1, return; end
	x = sign (x) .* abs(x)  .^ a;

function y = replacenan (x, v)
	if ~exist ('v')
	  v = 0;
	end
	y = x;
	y(isnan(x)) = v;

author: filipradenovic

repo: cnnimageretrieval


text:
% Fig. 6.26   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%

clear all;
%close all;
clf

num=1;
den=conv([1 0],[1 2 1]);
w=logspace(-2,3,100);
[re,im]=nyquist(num,den,w);
plot(re,im,re,-im);
axis([-2 2 -2 2]);
axis('square')
xlabel('Re(G(s))');
ylabel('Im(G(s))');
title('Fig. 6.26 Nyquist plot');
nicegrid;
hold on
ii=[11 21 47 61];
plot(re(ii),im(ii),'*');
plot(-.5,0,'*')
hold off


text:
% APM Specify Options
%
% response = apm_option(server,app,name,value)
%
% Function apm_option sends an option specification
%   to the APM server with the following arguments:
%
%   server = address of server
%      app = application name
%     name = option name
%    value = option value
%
% A response is returned indicating whether the option
%   was successfully recorded
%
% Either global options (NLC.option_name) or parameter/variable
%   options (NAME.option_name) are able to be specified with this
%   function. A full list of available options with their default
%   values are available here:
%
% Global Options:
%    http://apmonitor.com/wiki/index.php/Main/DbsGlobal
%
% Variable Options:
%    http://apmonitor.com/wiki/index.php/Main/DbsVariable
%
function response = apm_option(server,app,name,value)
   app = lower(deblank(app));
   aline = ['option ' deblank(char(name)) ' = ' num2str(value)];
   response = apm(server,app,aline);

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
% RES = reconWpyr(PYR, INDICES, FILT, EDGES, LEVS, BANDS)
%
% Reconstruct image from its separable orthonormal QMF/wavelet pyramid
% representation, as created by buildWpyr.
%
% PYR is a vector containing the N pyramid subbands, ordered from fine
% to coarse.  INDICES is an Nx2 matrix containing the sizes of
% each subband.  This is compatible with the MatLab Wavelet toolbox.
%
% FILT (optional) can be a string naming a standard filter (see
% namedFilter), or a vector which will be used for (separable)
% convolution.  Default = 'qmf9'.  EDGES specifies edge-handling,
% and defaults to 'reflect1' (see corrDn).
%
% LEVS (optional) should be a vector of levels to include, or the string
% 'all' (default).  1 corresponds to the finest scale.  The lowpass band
% corresponds to wpyrHt(INDICES)+1.
%
% BANDS (optional) should be a vector of bands to include, or the string
% 'all' (default).   1=horizontal, 2=vertical, 3=diagonal.  This is only used
% for pyramids of 2D images.

% Eero Simoncelli, 6/96.

function res = reconWpyr(pyr, ind, filt, edges, levs, bands)

if (nargin < 2)
  error('First two arguments (PYR INDICES) are required');
end

%%------------------------------------------------------------
%% OPTIONAL ARGS:

if (exist('filt') ~= 1)
  filt = 'qmf9';
end

if (exist('edges') ~= 1)
  edges= 'reflect1';
end

if (exist('levs') ~= 1)
  levs = 'all';
end

if (exist('bands') ~= 1)
  bands = 'all';
end

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

maxLev = 1+wpyrHt(ind);
if strcmp(levs,'all')
  levs = [1:maxLev]';
else
  if (any(levs > maxLev))
    error(sprintf('Level numbers must be in the range [1, %d].', maxLev));
  end
  levs = levs(:);
end

if strcmp(bands,'all')
  bands = [1:3]';
else
  if (any(bands < 1) | any(bands > 3))
    error('Band numbers must be in the range [1,3].');
  end
  bands = bands(:);
end

if isstr(filt)
  filt = namedFilter(filt);
end

filt = filt(:);
hfilt = modulateFlip(filt);

%% For odd-length filters, stagger the sampling lattices:
if (mod(size(filt,1),2) == 0)
	stag = 2;
else
	stag = 1;
end

%% Compute size of result image: assumes critical sampling (boundaries correct)
res_sz = ind(1,:);
if (res_sz(1) == 1)
  loind = 2;
  res_sz(2) = sum(ind(:,2));
elseif (res_sz(2) == 1)
  loind = 2;
  res_sz(1) = sum(ind(:,1));
else
  loind = 4;
  res_sz = ind(1,:) + ind(2,:);  %%horizontal + vertical bands.
  hres_sz = [ind(1,1), res_sz(2)];
  lres_sz = [ind(2,1), res_sz(2)];
end

%% First, recursively collapse coarser scales:
if any(levs > 1)

  if (size(ind,1) > loind)
    nres = reconWpyr( pyr(1+sum(prod(ind(1:loind-1,:)')):size(pyr,1)), ...
	ind(loind:size(ind,1),:), filt, edges, levs-1, bands);
  else
    nres = pyrBand(pyr, ind, loind); 	% lowpass subband
  end

  if (res_sz(1) == 1)
    res = upConv(nres, filt', edges, [1 2], [1 stag], res_sz);
  elseif (res_sz(2) == 1)
    res = upConv(nres, filt, edges, [2 1], [stag 1], res_sz);
  else
    ires = upConv(nres, filt', edges, [1 2], [1 stag], lres_sz);
    res = upConv(ires, filt, edges, [2 1], [stag 1], res_sz);
  end

else

  res = zeros(res_sz);

end

%% Add  in reconstructed bands from this level:
if any(levs == 1)
  if (res_sz(1) == 1)
    upConv(pyrBand(pyr,ind,1), hfilt', edges, [1 2], [1 2], res_sz, res);
  elseif (res_sz(2) == 1)
    upConv(pyrBand(pyr,ind,1), hfilt, edges, [2 1], [2 1], res_sz, res);
  else
    if any(bands == 1) % horizontal
      ires = upConv(pyrBand(pyr,ind,1),filt',edges,[1 2],[1 stag],hres_sz);
      upConv(ires,hfilt,edges,[2 1],[2 1],res_sz,res);  %destructively modify res
    end
    if any(bands == 2) % vertical
      ires = upConv(pyrBand(pyr,ind,2),hfilt',edges,[1 2],[1 2],lres_sz);
      upConv(ires,filt,edges,[2 1],[stag 1],res_sz,res);  %destructively modify res
    end
    if any(bands == 3) % diagonal
      ires =  upConv(pyrBand(pyr,ind,3),hfilt',edges,[1 2],[1 2],hres_sz);
      upConv(ires,hfilt,edges,[2 1],[2 1],res_sz,res);  %destructively modify res
    end
  end
end

repo: SelfExSR


text:
function PlotBoxEdges
factor=1.00;
plot3(factor*[0 1],factor*[0 0],factor*[0 0],'k','LineWidth',5),hold on
plot3(factor*[0 0],factor*[0 0],factor*[0 1],'k','LineWidth',5)
plot3(factor*[0 0],factor*[0 1],factor*[0 0],'k','LineWidth',5)
plot3(factor*[0 1],factor*[1 1],factor*[1 1],'k','LineWidth',5)
plot3(factor*[1 1],factor*[1 1],factor*[1 0],'k','LineWidth',5)
plot3(factor*[1 1],factor*[0 1],factor*[1 1],'k','LineWidth',5)
plot3(factor*[0 0],factor*[0 1],factor*[1 1],'k','LineWidth',5)
plot3(factor*[1 1],factor*[0 1],factor*[0 0],'k','LineWidth',5)
plot3(factor*[0 1],factor*[1 1],factor*[0 0],'k','LineWidth',5)
plot3(factor*[0 1],factor*[0 0],factor*[1 1],'k','LineWidth',5)
plot3(factor*[0 0],factor*[1 1],factor*[0 1],'k','LineWidth',5)
plot3(factor*[1 1],factor*[0 0],factor*[1 0],'k','LineWidth',5)

path: github-repos/MATLAB/Sable-mcbench-benchmarks/mcbench-benchmarks-ba13b2f0296ef49491b95e3f984c7c41fccdb6d8/34985-monte-carlo-simulation-of-three-dimensional-grain-growth-code-version-no-1-basic/Monte Carlo Simulation Q-state Potts model 3D square-lattice - microstructure/PlotBoxEdges.m


text:
function M = fixedrankMNquotientfactory(m, n, k)
% Manifold of m-by-n matrices of rank k with two factor quotient geometry.
%
% function M = fixedrankMNquotientfactory(m, n, k)
%
% This follows the quotient geometry described in the following paper:
% P.-A. Absil, L. Amodei and G. Meyer,
% "Two Newton methods on the manifold of fixed-rank matrices endowed
%  with Riemannian quotient geometries", arXiv, 2012.
%
% Paper link: http://arxiv.org/abs/1209.0068
%
% A point X on the manifold is represented as a structure with two
% fields: M and N. The matrix M (mxk) is orthonormal, while the matrix N
% (nxk) is full-rank such that X = M*N';
%
% Tangent vectors are represented as a structure with two fields (M, N).
%
% Please cite the Manopt paper as well as the research paper:
%     @Article{absil2014fixedrank,
%       Title   = {Two Newton methods on the manifold of fixed-rank matrices endowed with Riemannian quotient geometries},
%       Author  = {Absil, P.-A. and Amodei, L. and Meyer, G.},
%       Journal = {Computational Statistics},
%       Year    = {2014},
%       Number  = {3-4},
%       Pages   = {569--590},
%       Volume  = {29},
%       Doi     = {10.1007/s00180-013-0441-6}
%     }

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:

    M.name = @() sprintf('MN'' quotient manifold of %dx%d matrices of rank %d', m, n, k);

    M.dim = @() (m+n-k)*k;

    % Choice of the metric is motivated by the symmetry present in the
    % space.
    M.inner = @(X, eta, zeta) eta.M(:).'*zeta.M(:) + eta.N(:).'*zeta.N(:);

    M.norm = @(X, eta) sqrt(M.inner(X, eta, eta));

    M.dist = @(x, y) error('fixedrankMNquotientfactory.dist not implemented yet.');

    M.typicaldist = @() 10*k;

    symm = @(X) .5*(X+X');
    stiefel_proj = @(M, H) H - M*symm(M'*H);

    M.egrad2rgrad = @egrad2rgrad;
    function eta = egrad2rgrad(X, eta)
        eta.M = stiefel_proj(X.M, eta.M);
    end

    M.ehess2rhess = @ehess2rhess;
    function Hess = ehess2rhess(X, egrad, ehess, eta)

        % Directional derivative of the Riemannian gradient.
        Hess.M = ehess.M - eta.M*symm(X.M'*egrad.M);
        Hess.M = stiefel_proj(X.M, Hess.M);

        Hess.N = ehess.N;

        % Projection onto the horizontal space.
        Hess = M.proj(X, Hess);
    end

    M.proj = @projection;
    function etaproj = projection(X, eta)

        % Start by projecting the vector from Rmp x Rnp to the tangent
        % space to the total space, that is, eta.M should be in the
        % tangent space to Stiefel at X.M and eta.N is arbitrary.
        eta.M = stiefel_proj(X.M, eta.M);

        % Now project from the tangent space to the horizontal space, that
        % is, take care of the quotient.

        % First solve a Sylvester equation (A symm., B skew-symm.)
        A = X.N'*X.N + eye(k);
        B = eta.M'*X.M + eta.N'*X.N;
        B = B-B';
        omega = lyap(A, -B);

        % And project along the vertical space to the horizontal space.
        etaproj.M = eta.M + X.M*omega;
        etaproj.N = eta.N + X.N*omega;

    end

    M.exp = @exponential;
    function Y = exponential(X, eta, t)
        if nargin < 3
            t = 1.0;
        end

        A = t*X.M'*eta.M;
        S = t^2*eta.M'*eta.M;
        Y.M = [X.M t*eta.M]*expm([A -S ; eye(k) A])*eye(2*k, k)*expm(-A);

        % re-orthonormalize (seems necessary from time to time).
        [Q R] = qr(Y.M, 0);
        Y.M = Q * diag(sign(diag(R)));

        Y.N = X.N + t*eta.N;

    end

    % Factor M lives on the Stiefel manifold, hence we will reuse its
    % random generator.
    stiefelm = stiefelfactory(m, k);

    M.retr = @retraction;
    function Y = retraction(X, eta, t)
        if nargin < 3
            t = 1.0;
        end

        Y.M = uf(X.M + t*eta.M); % This is a valid retraction
        Y.N = X.N + t*eta.N;
    end

    M.hash = @(X) ['z' hashmd5([X.M(:) ; X.N(:)])];

    M.rand = @random;
    function X = random()
        X.M = stiefelm.rand();
        X.N = randn(n, k);
    end

    M.randvec = @randomvec;
    function eta = randomvec(X)
        eta.M = randn(m, k);
        eta.N = randn(n, k);
        eta = projection(X, eta);
        nrm = M.norm(X, eta);
        eta.M = eta.M / nrm;
        eta.N = eta.N / nrm;
    end

    M.lincomb = @lincomb;

    M.zerovec = @(X) struct('M', zeros(m, k), 'N', zeros(n, k));

    M.transp = @(x1, x2, d) projection(x2, d);

end

% Linear combination of tangent vectors
function d = lincomb(x, a1, d1, a2, d2) %#ok<INMSL>

    if nargin == 3
        d.M = a1*d1.M;
        d.N = a1*d1.N;
    elseif nargin == 5
        d.M = a1*d1.M + a2*d2.M;
        d.N = a1*d1.N + a2*d2.N;
    else
        error('Bad use of fixedrankMNquotientfactory.lincomb.');
    end

end

function A = uf(A)
    [L, unused, R] = svd(A, 0);
    A = L*R';
end

author: MIT-SPARK

repo: GlobalOptimizationTutorial


text:
function order = tetra_unit_size ( rule )

%*****************************************************************************80
%
%% TETRA_UNIT_SIZE sizes quadrature weights and abscissas in the unit tetrahedron.
%
%  Integration region:
%
%      0 <= X,
%    and
%      0 <= Y,
%    and
%      0 <= Z,
%    and
%      X + Y + Z <= 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 April 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Hermann Engels,
%    Numerical Quadrature and Cubature,
%    Academic Press, 1980,
%    ISBN: 012238850X,
%    LC: QA299.3E5.
%
%    Patrick Keast,
%    Moderate Degree Tetrahedral Quadrature Formulas,
%    Computer Methods in Applied Mechanics and Engineering,
%    Volume 55, Number 3, May 1986, pages 339-348.
%
%    Olgierd Zienkiewicz,
%    The Finite Element Method,
%    Sixth Edition,
%    Butterworth-Heinemann, 2005,
%    ISBN: 0750663200,
%    LC: TA640.2.Z54
%
%  Parameters:
%
%    Input, integer RULE, the index of the rule.
%     1, order 1, precision 0, Newton Cotes formula #0, Zienkiewicz #1.
%     2, order 4, precision 1, Newton Cotes formula #1.
%     3, order 4, precision 2, Zienkiewicz #2.
%     4, order 10, precision 2, Newton Cotes formula #2
%     5, order 5, precision 3, Zienkiewicz #3.
%     6, order 8, precision 3, Newton Cotes formula #3.
%     7, order 35, precision 4, Newton Cotes formula #4.
%     8, order 11, precision 4, a Keast rule.
%
%    Output, integer  ORDER, the order of the rule.
%

%
%  Newton Cotes #0.
%
  if ( rule == 1 )

    order = 1;
%
%  Newton Cotes #1.
%
  elseif ( rule == 2 )

    order = 4;
%
%  Zienkiewicz #2.
%
  elseif ( rule == 3 )

    order = 4;
%
%  Newton Cotes #2.
%
  elseif ( rule == 4 )

    order = 10;
%
%  Zienkiewicz #3.
%
  elseif ( rule == 5 )

    order = 5;
%
%  Newton Cotes #3.
%  (This is actually formally a 20 point rule, but with 12 zero coefficients%)
%
  elseif ( rule == 6 )

    order = 8;
%
%  Newton Cotes #4.
%
  elseif ( rule == 7 )

    order = 35;
%
%  Keast Rule of order 11
%
  elseif ( rule == 8 )

    order = 11;

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'TETRA_UNIT_SIZE - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of RULE = %d\n', rule );
    error ( 'TETRA_UNIT_SIZE - Fatal error!' );

  end

  return
end


text:
function [ a, b ] = p17_lim ( dim_num )

%*****************************************************************************80
%
%% P17_LIM returns the integration limits for problem 17.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Output, real A(DIM_NUM), B(DIM_NUM), the lower and upper
%    limits of integration.
%
  a(1:dim_num) = 0.0;
  b(1:dim_num) = 1.0;

  return
end


text:
function pass = test_compose(pref)
% Check that CHEBFUN3T composition operations are working.

if ( nargin == 0 )
    pref = chebfunpref;
end
tol = 1000*pref.cheb3Prefs.chebfun3eps;

f = chebfun3t(@(x,y,z) cos(x.*y.*z) + sin(x.*y.*z) + y -.1);

% Multiplication
exact = @(x,y,z) (cos(x.*y.*z) + sin(x.*y.*z) + y -.1) .* ...
    sin((x-.1).*(y+.4).*(z+.8));
g = chebfun3t(@(x,y,z) exact(x,y,z));
pass(1) = norm(g - chebfun3t(@(x,y,z) f(x,y,z).*sin((x-.1).*(y+.4).*(z+.8)))) < tol;

% Sine
exact = @(x,y,z) sin(cos(x.*y.*z) + sin(x.*y.*z) + y -.1);
g = chebfun3t(@(x,y,z) exact(x,y,z));
pass(2) = norm(g - sin(f)) < tol;

% Cosine
exact = @(x,y,z) cos(cos(x.*y.*z) + sin(x.*y.*z) + y -.1);
g = chebfun3t(@(x,y,z) exact(x,y,z));
pass(3) = norm(g - cos(f)) < tol;

% Sinh
exact = @(x,y,z) sinh(cos(x.*y.*z) + sin(x.*y.*z) + y -.1);
g = chebfun3t(@(x,y,z) exact(x,y,z));
pass(4) = norm(g - sinh(f)) < tol;

% Cosh
exact = @(x,y,z) cosh(cos(x.*y.*z) + sin(x.*y.*z) + y -.1);
g = chebfun3t(@(x,y,z) exact(x,y,z));
pass(5) = norm( g - cosh(f) ) < tol;

% Tanh
exact = @(x,y,z) tanh(-(cos(x.*y.*z) + sin(x.*y.*z) + y -.1));
g = chebfun3t(@(x,y,z) exact(x,y,z));
pass(6) = norm(g - tanh(-f)) < tol;

% Exp
exact = @(x,y,z) exp(cos(x.*y.*z) + sin(x.*y.*z) + y -.1);
g = chebfun3t(@(x,y,z) exact(x,y,z));
pass(7) = norm(g - exp(f)) < tol;

% Multiple operations:
f = chebfun3t(@(x,y,z) sin(10*x.*y.*z), [-1 2 -1 1 -3 -1]);
pass(8) = norm(f+f+f-3*f) < 100*tol;

pass(9) = norm(f.*f-f.^2) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function timelock = fttimelock(this, chanind, timeind, trialind, freqind)
% Method for converting meeg object to Fieldtrip timelock/freq struct
% FORMAT timelock = fttimelock(this, chanind, timeind, trialind, freqind)
%
% The method support both time and TF data and outputs different variants
% of timelock or freq FT struct depending on the dataset type and requested
% data dimensions.
% _______________________________________________________________________
% Copyright (C) 2008-2012 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: fttimelock.m 6158 2014-09-09 12:23:49Z vladimir $

if ~islinked(this)
    error('There is no data linked to the object');
end

if nargin < 2 || ~isnumeric(chanind)
    chanind = 1:nchannels(this);
end

if nargin < 3 || ~isnumeric(timeind)
    timeind = 1:nsamples(this);
end

if nargin < 4 || ~isnumeric(trialind)
    trialind = 1:ntrials(this);
end

if strncmpi(transformtype(this),'TF',2) && ...
        (nargin < 5 || isempty(freqind))
     freqind = 1:nfrequencies(this);
end

timelock             = [];
timelock.label       = chanlabels(this, chanind)';

if isequal(transformtype(this), 'time')
   if isequal(type(this), 'continuous')
            error('For continuous data use ftraw method');
   end

   if isequal(type(this), 'single') || length(trialind)>1
       timelock.dimord  = 'rpt_chan_time';
       timelock.trial   =  permute(subsref(this, substruct('()', {chanind, timeind, trialind})), [3 1 2]);
   else
       timelock.dimord  = 'chan_time';
       timelock.avg     =  spm_squeeze(subsref(this, substruct('()', {chanind, timeind, trialind})), 3);
   end

   timelock.time       = time(this, timeind);

elseif strncmpi(transformtype(this),'TF',2)
    if length(timeind)>1
        if isequal(type(this), 'single') || length(trialind)>1
            timelock.dimord    = 'rpt_chan_freq_time';
            timelock.powspctrm = permute(subsref(this, substruct('()', {chanind, freqind, timeind, trialind})), [4 1 2 3]);
        else
            timelock.dimord     = 'chan_freq_time';
            timelock.powspctrm  =  spm_squeeze(subsref(this, substruct('()', {chanind, freqind, timeind, trialind})), 3);
        end

        timelock.time       = time(this, timeind);
    else
        if isequal(type(this), 'single') || length(trialind)>1
            timelock.dimord    = 'rpt_chan_freq';
            timelock.powspctrm = spm_squeeze(permute(subsref(this, substruct('()', {chanind, freqind, timeind, trialind})), [4 1 2 3]), 4);
        else
            timelock.dimord     = 'chan_freq';
            timelock.powspctrm  =  spm_squeeze(subsref(this, substruct('()', {chanind, freqind, timeind, trialind})), [3 4]);
        end
    end

    timelock.freq      = frequencies(this, freqind);
else
    error('Unknown transform type.');
end

if length(trialind)>1

    clist      =  condlist(this);
    condlabels = conditions(this, trialind);
    timelock.trialinfo = 0*trialind;

    for k = 1:numel(clist)
        fprintf('mapping condition label "%s" to condition code %d\n', clist{k}, k);
        sel = strcmp(clist{k}, condlabels);
        timelock.trialinfo(sel) = k;
    end

end

if ~isempty(sensors(this, 'MEG'))
    timelock.grad = sensors(this, 'MEG');
end

if ~isempty(sensors(this, 'EEG'))
    timelock.elec = sensors(this, 'EEG');
end

author: spm


text:
function example3 ( )

%*****************************************************************************80
%
%% EXAMPLE3 uses BVP4C to solve the EXAMPLE3 problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2013
%
%  Author:
%
%    Original MATLAB version by Shampine, Kierzenka, Reichelt.
%    This version by John Burkardt.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EXAMPLE3:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Use BVP4C to solve the following eigenvalue boundary value problem:\n' );
  fprintf ( 1, '  y" + (lambda - 2 q cos(2x) y = 0\n' );
  fprintf ( 1, '  y''(0) = 0, y''(pi) = 0, y(0) 1\n' );
  fprintf ( 1, '  The initial guess is set using the functional form for y_init.\n' );
%
%  We need an estimate for the unknown parameter lambda.
%
  lambda = 15.0;
%
%  Set SOLINIT, the structure defining the initial guess.
%
  x_init = linspace ( 0.0, pi, 11 );
%
%  There are two choices for specifying y_init:
%  1) y_init can be a vector, containing a constant value for each component of y.
%  2) y_init can be a function which returns the guess value for the components of y
%     at a point x.
%  It is NOT possible to specify an M by N vector of values of the M components
%  of Y at each of the N values X.
%
  solinit = bvpinit ( x_init, @example3_init, lambda );
%
%  Have BVP4C solve the problem.
%
  sol = bvp4c ( @example3_ode, @example3_bc, solinit );
%
%  Use DEVAL to evaluate the solution.
%
  x = linspace ( 0.0, pi, 101 );
  y = deval ( sol, x );
  lambda_computed = sol.parameters;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Computed eigenvalue lambda = %g\n', lambda_computed );
%
%  Display a plot of Y and Y'.
%
  plot ( x, y(1,:), 'r-', 'Linewidth', 2 );
  xlabel ( '<--- X --->', 'Fontsize', 16 );
  ylabel ( '<--- Y(X) --->', 'Fontsize', 16 );
  title ( 'EXAMPLE3: Eigenfunction 4 of Mathieu''s Equation', 'Fontsize', 16 )
  grid on
  filename = 'example3.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saving plot file as "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EXAMPLE3:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function dydx = example3_ode ( x, y, lambda )

%*****************************************************************************80
%
%% EXAMPLE3_ODE evaluates the right hand side of the ODE.
%
%  Discussion:
%
%    We assume that the differential equation has been rewritten as a
%    system of first order equations of the form
%
%      dydx = f(x,y)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the point at which the ODE is to be evaluated.
%
%    Input, real Y(M), the value of the solution at X.
%
%    Input, real LAMBDA, the estimated eigenvalue.
%
%    Output, real DYDX(M), the value of the right hand side given X and Y.
%
  q = 5.0;

  dydx(1) = y(2);
  dydx(2) = - ( lambda - 2.0 * q * cos ( 2.0 * x ) ) * y(1);

  return
end
function bc = example3_bc ( ya, yb, lambda )

%*****************************************************************************80
%
%% EXAMPLE3_BC evaluates the boundary conditions.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real YA(M), YB(M), the solution value at the left and right endpoints.
%
%    Input, real LAMBDA, the estimated eigenvalue.
%
%    Output, real BC(2), the value of the boundary conditions.
%
  bc(1) = ya(2);
  bc(2) = yb(2);
  bc(3) = ya(1) - 1.0;

  return
end
function y_init = example3_init ( x )

%*****************************************************************************80
%
%% EXAMPLE3_INIT evaluates the initial guess at a point X.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 September 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the evaluation point.
%
%    Output, real Y_INITC(2), the value of the initial guess.
%
  y_init(1) = cos ( 4.0 * x );
  y_init(2) = - 4.0 * sin ( 4.0 * x );

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end