text:
function test_tutorial_coherence

% MEM 5gb
% WALLTIME 00:20:00
% DEPENDENCY ft_freqanalysis ft_connectivityanalysis ft_multiplotER ft_singleplotER ft_topoplotER ft_sourceanalysis ft_sourceinterpolate ft_prepare_sourcemodel headsurface

addpath(dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/coherence')); % this contains the trial function

% find the interesting epochs of data
cfg = [];
cfg.trialfun                  = 'trialfun_left';
cfg.dataset                   = dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/SubjectCMC.ds');
cfg = ft_definetrial(cfg);

% detect EOG artifacts in the MEG data
cfg.continuous                = 'yes';
cfg.artfctdef.eog.padding     = 0;
cfg.artfctdef.eog.bpfilter    = 'no';
cfg.artfctdef.eog.detrend     = 'yes';
cfg.artfctdef.eog.hilbert     = 'no';
cfg.artfctdef.eog.rectify     = 'yes';
cfg.artfctdef.eog.cutoff      = 2.5;
cfg.artfctdef.eog.interactive = 'no';
cfg = ft_artifact_eog(cfg);

% detect jump artifacts in the MEG data
cfg.artfctdef.jump.interactive = 'no';
cfg.padding                    = 5;
cfg = ft_artifact_jump(cfg);

% detect muscle artifacts in the MEG data
cfg.artfctdef.muscle.cutoff   = 8;
cfg.artfctdef.muscle.interactive = 'no';
cfg = ft_artifact_muscle(cfg);

% reject the epochs that contain artifacts
cfg.artfctdef.reject          = 'complete';
cfg = ft_rejectartifact(cfg);

% preprocess the MEG data
cfg.demean                    = 'yes';
cfg.dftfilter                 = 'yes';
cfg.channel                   = {'MEG'};
cfg.continuous                = 'yes';
meg = ft_preprocessing(cfg);

% preprocess the EMG data
cfg              = [];
cfg.dataset      = meg.cfg.dataset;
cfg.trl          = meg.cfg.trl;
cfg.continuous   = 'yes';
cfg.demean       = 'yes';
cfg.dftfilter    = 'yes';
cfg.channel      = {'EMGlft' 'EMGrgt'};
cfg.hpfilter     = 'yes';
cfg.hpfreq       = 10;
cfg.rectify      = 'yes';
emg = ft_preprocessing(cfg);

% concatenate the two data-structures into one structure
data = ft_appenddata([], meg, emg);

% visualisation
figure
subplot(2,1,1);
plot(data.time{1},data.trial{1}(77,:));
axis tight;
legend(data.label(77));

subplot(2,1,2);
plot(data.time{1},data.trial{1}(152:153,:));
axis tight;
legend(data.label(152:153));

% spectral analysis: fourier
cfg            = [];
cfg.output     = 'fourier';
cfg.method     = 'mtmfft';
cfg.foilim     = [5 100];
cfg.tapsmofrq  = 5;
cfg.keeptrials = 'yes';
cfg.channel    = {'MEG' 'EMGlft' 'EMGrgt'};
freqfourier    = ft_freqanalysis(cfg, data);

% spectral analysis: powandcsd
cfg            = [];
cfg.output     = 'powandcsd';
cfg.method     = 'mtmfft';
cfg.foilim     = [5 100];
cfg.tapsmofrq  = 5;
cfg.keeptrials = 'yes';
cfg.channel    = {'MEG' 'EMGlft' 'EMGrgt'};
cfg.channelcmb = {'MEG' 'EMGlft'; 'MEG' 'EMGrgt'};
freq           = ft_freqanalysis(cfg, data);

% compute coherence
cfg            = [];
cfg.method     = 'coh';
cfg.channelcmb = {'MEG' 'EMG'};
fd             = ft_connectivityanalysis(cfg, freq);
fdfourier      = ft_connectivityanalysis(cfg, freqfourier);

% visualisation
cfg                  = [];
cfg.parameter        = 'cohspctrm';
cfg.xlim             = [5 80];
cfg.refchannel       = 'EMGlft';
cfg.layout           = 'CTF151.lay';
cfg.showlabels       = 'yes';
figure; ft_multiplotER(cfg, fd)

cfg.channel = 'MRC21';
figure; ft_singleplotER(cfg, fd);

cfg                  = [];
cfg.parameter        = 'cohspctrm';
cfg.xlim             = [15 20];
cfg.zlim             = [0 0.1];
cfg.refchannel       = 'EMGlft';
cfg.layout           = 'CTF151.lay';
figure; ft_topoplotER(cfg, fd)

% 2 Hz smoothing
cfg            = [];
cfg.output     = 'powandcsd';
cfg.method     = 'mtmfft';
cfg.foilim     = [5 100];
cfg.tapsmofrq  = 2;
cfg.keeptrials = 'yes';
cfg.channel    = {'MEG' 'EMGlft'};
cfg.channelcmb = {'MEG' 'EMGlft'};
freq2          = ft_freqanalysis(cfg,data);

cfg            = [];
cfg.method     = 'coh';
cfg.channelcmb = {'MEG' 'EMG'};
fd2            = ft_connectivityanalysis(cfg,freq2);

cfg               = [];
cfg.parameter     = 'cohspctrm';
cfg.refchannel    = 'EMGlft';
cfg.xlim          = [5 80];
cfg.channel       = 'MRC21';
figure; ft_singleplotER(cfg, fd, fd2);

% 10 Hz smoothing
cfg            = [];
cfg.output     = 'powandcsd';
cfg.method     = 'mtmfft';
cfg.foilim     = [5 100];
cfg.keeptrials = 'yes';
cfg.channel    = {'MEG' 'EMGlft'};
cfg.channelcmb = {'MEG' 'EMGlft'};
cfg.tapsmofrq = 10;
freq10        = ft_freqanalysis(cfg,data);

cfg            = [];
cfg.method     = 'coh';
cfg.channelcmb = {'MEG' 'EMG'};
fd10          = ft_connectivityanalysis(cfg,freq10);

cfg               = [];
cfg.parameter     = 'cohspctrm';
cfg.xlim          = [5 80];
cfg.ylim          = [0 0.2];
cfg.refchannel    = 'EMGlft';
cfg.channel       = 'MRC21';
figure;ft_singleplotER(cfg, fd, fd2, fd10);

% 50 trials
cfg            = [];
cfg.output     = 'powandcsd';
cfg.method     = 'mtmfft';
cfg.foilim     = [5 100];
cfg.tapsmofrq  = 5;
cfg.keeptrials = 'yes';
cfg.channel    = {'MEG' 'EMGlft'};
cfg.channelcmb = {'MEG' 'EMGlft'};
cfg.trials     = 1:50;
freq50         = ft_freqanalysis(cfg,data);

cfg            = [];
cfg.method     = 'coh';
cfg.channelcmb = {'MEG' 'EMG'};
fd50           = ft_connectivityanalysis(cfg,freq50);

cfg                  = [];
cfg.parameter        = 'cohspctrm';
cfg.xlim             = [5 100];
cfg.ylim             = [0 0.2];
cfg.refchannel       = 'EMGlft';
cfg.channel          = 'MRC21';
figure; ft_singleplotER(cfg, fd, fd50);

% source reconstruction
cfg            = [];
cfg.method     = 'mtmfft';
cfg.output     = 'powandcsd';
cfg.foilim     = [18 18];
cfg.tapsmofrq  = 5;
cfg.keeptrials = 'yes';
cfg.channelcmb = {'MEG' 'MEG';'MEG' 'EMGlft'};
freq           = ft_freqanalysis(cfg, data);

cfg           = [];
cfg.method    = 'dics';
cfg.refchan   = 'EMGlft';
cfg.frequency = 18;
cfg.headmodel = dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/SubjectCMC.hdm');
cfg.inwardshift = 1;
cfg.sourcemodel.resolution = 1;
cfg.sourcemodel.unit       = 'cm';
source        = ft_sourceanalysis(cfg, freq);

mri = ft_read_mri(dccnpath('/home/common/matlab/fieldtrip/data/ftp/test/ctf/SubjectCMC.mri'));
mri = ft_volumereslice([], mri);

cfg            = [];
cfg.parameter  = 'coh';
cfg.downsample = 2;
interp         = ft_sourceinterpolate(cfg, source, mri);

cfg              = [];
cfg.method       = 'ortho';
%cfg.interactive  = 'yes';
cfg.funparameter = 'coh';
figure; ft_sourceplot(cfg, interp);

%--------------------------------
% subfunction
function trl = trialfun_left(cfg)

% read in the triggers and create a trial-matrix
% consisting of 1-second data segments, in which
% left ECR-muscle is active.

event = ft_read_event(cfg.dataset);
trig  = [event(find(strcmp('backpanel trigger', {event.type}))).value];
indx  = [event(find(strcmp('backpanel trigger', {event.type}))).sample];

%left-condition
sel = [find(trig==1028):find(trig==1029)];

trig = trig(sel);
indx = indx(sel);

trl = [];
for j = 1:length(trig)-1
  trg1 = trig(j);
  trg2 = trig(j+1);
  if trg1<=100 && trg2==2080
    trlok      = [[indx(j)+1:1200:indx(j+1)-1200]' [indx(j)+1200:1200:indx(j+1)]'];
    trlok(:,3) = [0:-1200:-1200*(size(trlok,1)-1)]';
    trl        = [trl; trlok];
  end
end


text:
function x = p11_start ( )

%*****************************************************************************80
%
%% P11_START returns a starting point for optimization for problem 11.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real X, a starting point for the optimization.
%
  x = 0.75;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function rgb=hsl2rgb(hsl_in)
%Converts Hue-Saturation-Luminance Color value to Red-Green-Blue Color value
%
%Usage
%       RGB = hsl2rgb(HSL)
%
%   converts HSL, a M [x N] x 3 color matrix with values between 0 and 1
%   into RGB, a M [x N] X 3 color matrix with values between 0 and 1
%
%See also rgb2hsl, rgb2hsv, hsv2rgb

% (C) Vladimir Bychkovsky, June 2008
% written using:
% - an implementation by Suresh E Joel, April 26,2003
% - Wikipedia: http://en.wikipedia.org/wiki/HSL_and_HSV

hsl=reshape(hsl_in, [], 3);

H=hsl(:,1);
S=hsl(:,2);
L=hsl(:,3);

lowLidx=L < (1/2);
q=(L .* (1+S) ).*lowLidx + (L+S-(L.*S)).*(~lowLidx);
p=2*L - q;
hk=H; % this is already divided by 360

t=zeros([length(H), 3]); % 1=R, 2=B, 3=G
t(:,1)=hk+1/3;
t(:,2)=hk;
t(:,3)=hk-1/3;

underidx=t < 0;
overidx=t > 1;
t=t+underidx - overidx;

range1=t < (1/6);
range2=(t >= (1/6) & t < (1/2));
range3=(t >= (1/2) & t < (2/3));
range4= t >= (2/3);

% replicate matricies (one per color) to make the final expression simpler
P=repmat(p, [1,3]);
Q=repmat(q, [1,3]);
rgb_c= (P + ((Q-P).*6.*t)).*range1 + ...
        Q.*range2 + ...
        (P + ((Q-P).*6.*(2/3 - t))).*range3 + ...
        P.*range4;

rgb_c=round(rgb_c.*10000)./10000;
rgb=reshape(rgb_c, size(hsl_in));

author: Sable

repo: mcbench-benchmarks


text:
function [ttot, info, infoTransfdiff, imout] = transfdiffreg(transform, im, optReg, optDiff)
% TRANSFDIFFREG Transform diffusion registration of a sequence of images.
%
% TRANSFDIFFREG implements a registration algorithm for a sequence of
% images I = 1, 2, ..., N, such that each image I is aligned to its two
% neighbours I-1 and I+1. This algorithm is based on diffusion of the
% transforms between pairs of images.
%
%     |         |         |         |         |         |         |
%     |  <--->  |  <--->  |  <--->  |  <--->  |  <--->  |  <--->  |
%     |         |         |         |         |         |         |
%
% Note: Currently, only implemented for matched filter rigid registration
% (regmatchedfilt) and elastix B-spline transforms (BSplineTransform).
%
% [TTOT, INFO, INFO_TRANSFDIFF, IMOUT] = TRANSFDIFFREG(TRANSFORM, IM, OPTREG, OPTDIFF)
%
%   TRANSFORM is a string with the type of transform that will be computed
%   by the registration algorithm and diffused (more details below):
%
%     'regmatchedfilt': matched filter (i.e. cross-correlation) rigid
%        transforms.
%
%     'BSplineTransform': B-spline transforms.
%
%   IM is a sequence of 2D images, grayscale or colour. The images can be
%   in one of three formats:
%
%     * A plain array IM(R, C, I, CH), R=row, C=column, I=image index,
%       CH=colour channel.
%
%     * A vector of scimat structs.  See "help scimat" for details. IM(I)
%       has the I-th image.
%
%     * A cell vector. IM{I} has the I-th image.
%
%   OPTREG and OPTDIFF are structs with parameters for the registration and
%   diffusion parts of the algorithm, respectively (see below for details).
%
%   TTOT is a vector of N structs with the transforms (in elastix format)
%   that applied to each image forms an aligned volume.
%
%   INFO is a struct (see below for details).
%
%   INFO_TRANSFDIFF is a vector of structs. INFO_TRANSFDIFF(I) contains the
%   output information from the diffusion step in registration iteration I.
%
%   IMOUT is result of transforming each IM with the corresponding TTOT.
%   Its format (array or scimat struct) is the same as IM's.
%
% =========================================================================
%
% [TTOT, IMOUT, INFO] = TRANSFDIFFREG('regmatchedfilt', IM, OPTREG, OPTDIFF)
%
%   Diffusion of rigid transforms computed by matched filter or
%   cross-correlation registration. Global optimum.
%
%   OPTREG:
%
%     'verbose': (def false) Verbose output.
%
%     'Angle': (def (-45:45)/180*pi) Vector of angle values that the
%         matched filter will try to align adjacent slices.
%
%     'tp': (def {}) User provided intra-slice registrations. tp(I) is the
%           result of registering image I to image I+1. This is by default
%           computed internally, but with this option the user can skip the
%           slow registration process and go directly into diffusion. You
%           can obtain tp from a previous run of the algorithm in output
%           INFO.tp.
%
%     'outdir': (def tempname) If IM are provided as filenames, path to the
%        output directory. For other IM formats, it is ignored.
%
%   OPTDIFF:
%
%     'verbose': (def false) Verbose output.
%
%     'MaxIter': (def 100) Stopping criterion. The algorithm stops after
%                MaxIter diffusion iterations.
%
%     'Alpha':   (def 0.45) Diffusion coefficient. Values >0.45 may not
%                smooth high frequencies very well. Values >=0.5 cause
%                oscillations.
%
%   INFO is a struct with the output information from the top-level
%   registration algorithm.
%
%     'StopCondition': Reason why the algorithm stopped.
%
%     'tp', 'tm': Are vectors with the I -> I+1 and I -> I-1 registration
%                 results from the last registration iteration.
%
% =========================================================================
%
% [TTOT, IMOUT, INFO] = TRANSFDIFFREG('BSplineTransform', IM, OPTREG, OPTDIFF)
%
%   Diffusion of B-spline transforms. Because of B-spline non-linearity,
%   the algorithm needs to run several iterations of registration followed
%   by diffusion.
%
%   This method uses Choi and Lee (2000)'s injectivity conditions to
%   guarantee the diffused B-splines are injective and have no fold-overs.
%   Note that these conditions are conservative, so it may be possible to
%   diffuse a bit more without fold-overs.
%
%   OPTREG:
%
%     'verbose': (def false) Verbose output.
%
%     'SpatialSamplesRatio': (def 1.0) Scalar value in [0.0, 1.0], with the
%                            percentage of pixels sampled by the
%                            registration algorithm.
%
%     'MaxIter': (def 5) Stopping criterion. Maximum number of registration
%                iterations.
%
%     'MaxVal': (def 0) Stopping criterion. After a registration sweep and
%               diffusion, if the update of all B-spline coefficients is <=
%               MaxVal, then no more registration sweeps are computed and
%               the algorithm stops.
%
%     'RegParam': (def {}) A string with the path to a file, or a struct in
%                 elastix format with the registration parameters.
%
%     'mask': (def {}) An array of binary masks with the same dimensions as
%                      the images. Anything outside the binary masks will
%                      be ignored for registration. This can be used to
%                      avoid artifacts and to speed up the algorithm.
%                      Furthermore, note that masks don't need to be very
%                      precise, so performance can be improved a lot if
%                      they are provided as low resolution images (with the
%                      correct metainformation so that they overlap in real
%                      world coordinates).
%
%     't0': (def {}) A list of structs with initial transforms in elastix
%           format. Each t0(I) is applied to image I before registration.
%           This can be used to provide a pre-alignment of the images,
%           typically rigid.
%
%     'tp': See 'regmatchedfilt' above.
%
%     'tm': Like tp, but tm(I) is the result of registering image I to
%           image I-1. Note that for B-splines we cannot compute tm(I) from
%           tp(I).
%
%     'CacheFile': (def '') Name of a .mat file to save/read  internal
%                  variables to and from after each full registration
%                  sweep. This allows to kill the algorithm, and restart it
%                  without having to repeat previous registration sweeps.
%
%   OPTDIFF:
%
%     'verbose': (def false) Verbose output.
%
%     'MaxIter': (def 100) Stopping criterion. Maximum number of diffusion
%                iterations.
%
%     'Epsilon': (def 0.0) Stopping criterion. It stops when all
%                coefficient displacement components are small:
%                for all t, abs(t)<=Epsilon.
%
%     'Alpha':   (def 0.45) Diffusion coefficient. Values >0.45 may not
%                smooth high frequencies very well. Values >=0.5 cause
%                oscillations.
%
%   INFO is a struct with the output information from the top-level
%   registration algorithm.
%
%     'StopCondition': Reason why the algorithm stopped.
%
%     'tp', 'tm': Are vectors with the I -> I+1 and I -> I-1 registration
%                 results from the last registration iteration.
%
%     'MaxVal': MaxVal(I, J) is the maximum absolute value of the
%               coefficients of the I-th transform at registration
%               iteration J.
%
%     'MedVal': MedVal(I, J) is the median absolute value of the
%               coefficients of the I-th transform at registration
%               iteration J.
%
% =========================================================================

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2015-2016 University of Oxford
% Version: 0.7.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

%% Process input arguments

% check number of arguments
switch (transform)
    case {'regmatchedfilt', 'BSplineTransform'}
        narginchk(2, 4);
        nargoutchk(0, 4);
    otherwise
        error('Not implemented for this transform')
end

% defaults
if (nargin < 3 || isempty(optReg))
    optReg = struct;
end
if (nargin < 4 || isempty(optDiff))
    optDiff = struct;
end

% common diffusion parameters
if (~isfield(optReg, 'verbose'))
    optReg.verbose = false;
end

% list of transforms, the options they admit and their default values
optRegList = {
    'regmatchedfilt',   'verbose',             false
    'regmatchedfilt',   'Angle',               (-45:45)/180*pi
    'regmatchedfilt',   'tp',                  {}
    'regmatchedfilt',   'outdir',              tempname
    'BSplineTransform', 'verbose',             false
    'BSplineTransform', 'SpatialSamplesRatio', 1.0
    'BSplineTransform', 'MaxIter',             5
    'BSplineTransform', 'MaxVal',              0
    'BSplineTransform', 'RegParam',            {}
    'BSplineTransform', 'mask',                {}
    'BSplineTransform', 't0',                  {}
    'BSplineTransform', 'tp',                  {}
    'BSplineTransform', 'tm',                  {}
    'BSplineTransform', 'outdir',              tempname
    'BSplineTransform', 'CacheFile',           ''
    };
optDiffList = {
    'regmatchedfilt',   'Alpha',               0.45
    'regmatchedfilt',   'MaxIter',             100
    'BSplineTransform', 'Alpha',               0.45
    'BSplineTransform', 'MaxIter',             100
    'BSplineTransform', 'Epsilon',             0.0
    };
optReg = check_user_options(optRegList, transform, optReg);
optDiff = check_user_options(optDiffList, transform, optDiff);

% number of input images
if (iscell(im) || isstruct(im)) % filenames or scimats
    N = length(im);
elseif (ischar(im))
    error('IM cannot be a single filename')
else % array
    N = size(im, 3);
end

% transform-specific checks
switch (transform)

    case 'regmatchedfilt'

        % with affine transforms, registrations between images need to be
        % computed only once
        optReg.MaxIter = 1;

        % transform any input image into scimat format
        [im, imOutType] = imToScimat(im);

    case 'BSplineTransform'

        if (optReg.SpatialSamplesRatio < 0.0 ...
                || optReg.SpatialSamplesRatio > 1.0)
            error('optReg.SpatialSamplesRatio must be in [0, 1]')
        end

        % if registration parameters provided as file name, read the file
        % (we may need to change some of the parameters, e.g.
        % SpatialSamplesRatio, and it also makes it easier to check the
        % values of the parameters)
        if (isempty(optReg.RegParam))
            error('No registration parameters provided in optReg.RegParam')
        elseif (ischar(optReg.RegParam))
            optReg.RegParam = elastix_read_file2param(optReg.RegParam);
        end

        % hard-coded options for registration between two images
        optElastix.verbose = false;

        % if images are in memory instead of as files, we have to
        % create temp files for them
        [im, deleteImFiles, imOutType] = imToFiles(im);
        N = length(im);

        % if the user has provided masks, and they are in memory, they need to be
        % saved as temp files
        [optReg.mask, deleteMaskFiles, maskType, nVoxMask] = imToFiles(optReg.mask);
        Nmask = length(optReg.mask);
        if (Nmask > 0 && N ~= Nmask)
            error('There must be either no masks or the same number of masks (OPTREG.mask) as images (IM)')
        end

end

%% Main algorithm

% init number of registration iterations
RegIter = 0; % number of registration iterations

switch (transform)

    case 'BSplineTransform'

    % initialize total transforms applied to each image as empty
    if (isempty(optReg.t0))
        ttot = cell(1, N);
    else
        ttot = optReg.t0;
    end

    % if the user is running in cached mode, we set the internal state of
    % the algorithm to the latest full registration sweep
    if (exist(optReg.CacheFile, 'file') == 2)

        if (optReg.verbose)
            disp('***************************************************')
            disp('** READING CACHED INFO ****************************')
            disp('***************************************************')
        end

        % load intermediate results, except for optReg, im because the temp
        % files will have different random names
        load(optReg.CacheFile, ...
            'ttot', 'info', 'infoTransfdiff', 'transform', 'optDiff', ...
            'optRegList', 'optDiffList', ...
            'N', 'imOutType', 'optElastix', 'deleteImFiles', ...
            'deleteMaskFiles', 'nVoxMask', ...
            'maskType', 'Nmask', 'RegIter', 'J', 'imoutaux', ...
            'iterinfo', 'tp', 'tpMat',  't', 'tref')

    end

end

% loop registration steps
while (true)

    if (optReg.verbose)
        disp('***************************************************')
        disp('** REGISTRATION BLOCK *****************************')
        disp('***************************************************')
    end

    % increase iteration counter
    RegIter = RegIter + 1;
    if (optReg.verbose)
        disp(['** IterReg = ' num2str(RegIter) ...
            '/' num2str(optReg.MaxIter) ...
            '                                 **'])
    end

    %% registration block
    if (optReg.verbose)
        disp('** Register image i to image i+1                 **')
    end

    switch (transform)

        case 'regmatchedfilt'

            % register image to next neighbour
            if (isempty(optReg.tp))

                % register each image i to image i+1
                % images loop
                parfor J = 1:N-1

                    if (optReg.verbose)
                        fprintf('** ... %d/%d', J, N-1)
                        tic
                    end

                    % register each image J to its next neighbour J+1
                    aux(J) = regmatchedfilt(im(J+1), im(J), ...
                        optReg.Angle);

                    if (optReg.verbose)
                        fprintf(' (%0.2f sec)\n', toc)
                    end

                end % images loop
                info.tp = aux;

            % if the user has provided tp already, we are going to skip the
            % registration step. This allows for fast testing of diffusion
            % without having to recompute the registration again and again
            else % if (isempty(optReg.tp))

                if (optReg.verbose)
                    disp('** TP provided by user, skipping registration block')
                    info.tp = optReg.tp;
                end

            end % if (isempty(optReg.tp))

        case 'BSplineTransform'

            % register each image i to image i+1
            % images loop
            for J = 1:N-1

                if (optReg.verbose)
                    fprintf('** ... %d/%d', J, N-1)
                end

                % in case we only want to use a subset of the voxels within
                % the mask. We use the fixed image to set the number of
                % samped voxels for both images
                if (strcmp(optReg.RegParam.ImageSampler, 'RandomCoordinate'))

                    % number of samples we are going to take for
                    % registration
                    numSamples = round(nVoxMask(J) ...
                        * optReg.SpatialSamplesRatio);

                    % make sure that the number of pixels is not too low,
                    % or that we are taking more than what's in the mask
                    numSamples = max([numSamples 2000]);
                    numSamples = min([numSamples nVoxMask(J)]);

                    optReg.RegParam.NumberOfSpatialSamples = numSamples;
                end

                % if the user has provided tp already, we are going to skip
                % the registration step. This allows for fast testing of
                % diffusion without having to recompute the registration
                % again and again
                if (isempty(optReg.tp))

                    % register each image J to its next neighbour J+1
                    tic
                    [info.tp(J), imoutaux, iterinfo] = register_2_images(...
                        optReg.RegParam, ...         % registration param
                        im{J+1}, ...                 % fixed
                        im{J}, ...                   % moving
                        ttot(J+1), ...               % fixedT0
                        ttot(J), ...                 % movingT0
                        optReg.mask{J+1}, ...        % fixedMask
                        optReg.mask{J}, ...          % movingMask
                        optElastix);                 % elastix options

                    % DEBUG
                    if (optReg.verbose)
                        fprintf(' (tp %0.2f sec)', toc)
                    end

                else

                    info.tp(J) = optReg.tp(J);
                    fprintf(' (tp provided)')

                end

                % same as above but for tm
                if (isempty(optReg.tm))

                    % register each image J to its previous neighbour J-1
                    tic
                    [info.tm(J), imoutaux, iterinfo] = register_2_images(...
                        optReg.RegParam, ...         % registration param
                        im{J}, ...                   % fixed
                        im{J+1}, ...                 % moving
                        ttot(J), ...                 % fixedT0
                        ttot(J+1), ...               % movingT0
                        optReg.mask{J}, ...          % fixedMask
                        optReg.mask{J+1}, ...        % movingMask
                        optElastix);                 % elastix options

                    % DEBUG
                    if (optReg.verbose)
                        fprintf(' (tm %0.2f sec)\n', toc)
                    end

                else

                    info.tm(J) = optReg.tm(J);
                    fprintf(' (tm provided)\n')

                end

            end % images loop

            % average transforms so that the transform J->J+1 is the
            % inverse of J+1->J. This is important because otherwise
            % diffusion doesn't converge
            %
            % we only need to correct tp, because tm will be easily
            % computed as -tp inside of transfdiff()
            for J = 1:N-1

                tp(J) = elastix_colon(info.tp(J), 1);
                tm(J) = elastix_colon(info.tm(J), 1);

                tp(J).TransformParameters ...
                    = (tp(J).TransformParameters ...
                    - tm(J).TransformParameters) * .5;

            end
            clear tm

            % if the user provided tp and tm at the input, they were valid
            % only for the first registration iteration, and need to be
            % removed for subsequent iterations
            if (~isempty(optReg.tp))
                optReg.tp = [];
            end
            if (~isempty(optReg.tm))
                optReg.tm = [];
            end

    end % switch transform

    %% diffusion block

    if (optReg.verbose)
        disp('***************************************************')
        disp('** DIFFUSION BLOCK ********************************')
        disp('***************************************************')
    end

    switch (transform)

        case 'regmatchedfilt'

            % convert elastix parameter vectors to affine matrices
            tpMat = elastix_affine_struct2matrix(info.tp);

            % transform diffusion of affine parameters
            optDiff.Transform = 'EulerTransform';
            [t, infoTransfdiff] = transfdiff(optDiff, tpMat);

            % transfer diffused parameters to elastix structs
            ttot = elastix_affine_matrix2struct(t, info.tp(1));

        case 'BSplineTransform'

            % extract transform vector from each TP and concatenate into a
            % matrix
            tpMat = cat(1, tp(:).TransformParameters);

            % transform diffusion of B-spline parameters
            optDiff.Transform = transform; % transform = 'BSplineTransform'
            optDiff.tbsp = tp(1); % used to validate B-spline injectivity
            [t, infoTransfdiff(RegIter)] = transfdiff(optDiff, tpMat);

            % template for the current B-spline transforms to apply to each
            % image
            tref(1:N) = tp(1);

            % transfer diffused parameters to elastix structs
            for J = 1:N
                tref(J).TransformParameters = t(J, :);
            end

            % add the new level of B-splines to the series of transforms to
            % apply to each image
            ttot = elastix_cat(tref, ttot);

            % stopping criterion: all control point translations are
            % very small
            info.MaxVal(:, RegIter) = max(abs(t), [], 2);
            info.MedVal(:, RegIter) = median(abs(t), 2);

            % save intermediate results to cache file
            if (~isempty(optReg.CacheFile))
                save(optReg.CacheFile)
            end

            if (optReg.verbose)
                disp(['** max(info.MaxVal) = ' num2str(max(info.MaxVal(:, RegIter)))])
            end
            if (max(info.MaxVal(:, RegIter)) <= optReg.MaxVal)
                info.StopCondition = 'MaxVal';
                break
            end

    end

    % stop criterion
    if (RegIter == optReg.MaxIter)
        info.StopCondition = 'MaxIter';
        break
    end

end % end registration step loop

% apply transform results to the input images
if (nargout > 3)

    if (optReg.verbose)
        disp('***************************************************')
        disp('** OUTPUT IMAGE TRANSFORM BLOCK *******************')
        disp('***************************************************')
    end

    % create output directory if it doesn't exist already and the output is
    % going to be filenames
    if (iscell(im) && ~exist(optReg.outdir))
        mkdir(optReg.outdir);
    end

    % if the output is not filenames, then we will delete the files, as the
    % output will be provided in memory
    deleteImOutFiles = ~strcmp(imOutType, 'char') ...
        && ~strcmp(imOutType, 'filename');

    for I = 1:N

        if (optReg.verbose)
            fprintf('** ... %d/%d\n', I, N)
            tic
        end

        % if the output is going to be given as filenames, then we want
        % control over where those files are saved to
        if (strcmp(imOutType, 'filename') || strcmp(imOutType, 'char'))

            % filename we are going to provide to the user
            [~, filename, ext] = fileparts(im{I});
            optTransformix.outfile = [optReg.outdir filesep filename ext];

        else

            optTransformix.outfile = '';

        end

        % at this point, "im" can be:
        % * vector of scimat (regmatchedfilt)
        % * cell vector of filenames (BSplineTransform)
        if (isstruct(im))

            % transform input image
            aux = transformix(ttot(I), im(I));

        elseif (iscell(im))

            % transform input image
            aux = transformix(ttot(I), im{I}, optTransformix);

        else

            error('Internal check fail: IM should not have this type')

        end

        % make the output image match the format of the input image
        if (isstruct(im) && strcmp(imOutType, 'char'))

            error('TODO: Not implemented')

        elseif (iscell(im) && strcmp(imOutType, 'char'))

            imout = aux;

        elseif (isstruct(im) && strcmp(imOutType, 'filename'))

            imout{I} = optTransformix.outfile;
            scimat_save(imout{I}, aux);

        elseif (iscell(im) && strcmp(imOutType, 'filename'))

            imout{I} = aux;

        elseif (isstruct(im) && strcmp(imOutType, 'scimat'))

            imout(I) = aux;

        elseif (iscell(im) && strcmp(imOutType, 'scimat'))

            imout(I) = scimat_load(aux);
            delete(aux);

        elseif (isstruct(im) && strcmp(imOutType, 'array'))

            imout(:, :, I, :) = aux.data;

        elseif (iscell(im) && strcmp(imOutType, 'array'))

            error('TODO: Not implemented')

        end

    end

    % delete output files if the output is provided in memory, not files
    if (exist(optReg.outdir) && deleteImOutFiles)
        rmdir(optReg.outdir);
    end

end

% delete temp files
switch (transform)

    case 'BSplineTransform'

        % if temp files had to be created to pass images to elastix, we can
        % delete them now
        if (deleteImFiles)
            for I = 1:N
                delete(im{I})
            end
        end
        if (deleteMaskFiles && ~isempty(optReg.mask))
            for I = 1:N
                delete(optReg.mask{I})
            end
        end

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% nested functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% check_user_options
%
% check invalid options and set defaults
function opt = check_user_options(optList, transform, opt)

% list of options allowed for this transform
idx = cellfun(@(s) strcmp(s, transform), optList(:, 1));
optList = optList(idx, :);

% list of options in the user input
names = fieldnames(opt);

% check whether any option provided by the user is not used for this
% transform
for I = 1:length(names)

    if (~ismember(names{I}, optList(:, 2)))
        warning([transform ' does not use option ' inputname(3) ...
            '.' names{I}])
    end

end

% set options not provided by the user to their defaults
isProvided = ismember(optList(:, 2), names);
for I = find(~isProvided)'
    opt = setfield(opt, optList{I, 2}, optList{I, 3});
end

end

%% imToScimat()
%
% convert any input image format (filenames, arrays or scimat) to scimat
function [im, imType] = imToScimat(im)

% check the type and number of the input images
if (isempty(im)) % empty input

    % nothing to do
    return

elseif (ischar(im)) % single filename

    imType = 'char';
    im = scimat_load(im);

elseif (iscell(im)) % list of filenames

    imType = 'filename';
    for I = 1:length(im)
        aux(I) = scimat_load(im{I});
    end
    im = aux;

elseif (isstruct(im))

    imType = 'scimat';

else % plain array

    imType = 'array';

    for I = 1:size(im, 3)
        aux(I) = scimat_im2scimat(im(:, :, I, :));
    end
    im = aux;

end

end

%% imToFiles()
%
% convert any input format to filenames, creating temp files if necessary
%
% * Inputs
%
%    im: input images (filenames, scimats or arrays)
%
% * Outputs:
%
%    im: cell vector with paths and filenames to images
%    deleteFiles: whether we need to delete the image files when we finish
%                 registering them (i.e. whether they are temp files)
%    nVox: number of voxels ~= 0 (only useful for binary masks)
function [im, deleteFiles, imType, nVox] = imToFiles(im)

% check the type and number of the input images
if (isempty(im)) % empty input (typically, no masks provided)

    imType = '';
    deleteFiles = false;
    nVox = 0;
    return

elseif (ischar(im)) % single filename

    imType = 'char';
    deleteFiles = false;
    im = {im};

    % get metainformation
    if (nargout > 3) % user wants to know number of voxels ~= 0
        scimat0 = scimat_load(im{1});
        nVox = nnz(scimat0.data);
    end

elseif (iscell(im)) % list of filenames

    imType = 'filename';
    deleteFiles = false;
    N = length(im);

    % get metainformation (we assume that all slices have the same)
    if (nargout > 3) % user wants to know number of voxels ~= 0
        for I = 1:N
            scimat0 = scimat_load(im{I});
            nVox(I) = nnz(scimat0.data);
        end
    end

elseif (isstruct(im)) % list of scimat structs

    imType = 'scimat';

    % number of images
    N = length(im);

    % if user wants to know number of voxels ~= 0
    if (nargout > 3)
        for I = 1:N
            nVox(I) = nnz(im(I).data);
        end
    end

    % create a temp file for each input image
    deleteFiles = true;
    filename = cell(1, N);
    for I = 1:N
        filename{I} = [tempname '.mha'];
        scimat_save(filename{I}, im(I));
    end
    im = filename;

else % plain array

    imType = 'array';

    % number of images
    N = size(im, 3);

    % if user wants to know number of voxels ~= 0
    if (nargout > 3)
        if (size(im, 4) ~= 1)
            warning('User wants to count voxels~=0, so we assume binary mask, but it has more than one channel')
        end
        for I = 1:N
            nVox(I) = nnz(im(:, :, I, :));
        end
    end

    % create a temp file for each input image
    deleteFiles = true;
    filename = cell(1, N);
    for I = 1:N
        filename{I} = [tempname '.mha'];
        scimat_save(filename{I}, scimat_im2scimat(im(:, :, I, :)));
    end
    im = filename;

end

end

%% register_2_images: local function to register two images
function [t, imout, iterinfo] = register_2_images(regParam, fixed, moving, fixedT0, movingT0, fixedMask, movingMask, optElastix)
tic
if (iscell(fixedT0)) % no initial transform

    % the only way to pass empty initial transorms is as a cell with an
    % empty, but if we pass them to elastix, we need to remove the cell
    fixedT0 = fixedT0{1};
    optElastix.fMask = fixedMask;

else % initial transform

    % apply initial transform
    fixed = transformix(fixedT0, fixed);
    optElastix.fMask = transformix(fixedT0, fixedMask);

end
if (iscell(movingT0)) % no initial transform

    % the only way to pass empty initial transorms is as a cell with an
    % empty, but if we pass them to elastix, we need to remove the cell
    optElastix.mMask = movingMask;

else % initial transform

    % we apply the initial transform explictly to the moving image
    % and mask, because it's substantially faster than providing a list of
    % transforms in optElastix.t0
    optElastix.t0 = '';
    optElastix.mMask = transformix(movingT0, movingMask);
    moving = transformix(movingT0, moving);

end

% register both images
[t, imout, iterinfo] ...
    = elastix(regParam, fixed, moving, optElastix);

% concatenate new level of B-splines and previous initial transforms.
% Note: because of elastix internals, if we apply t0 explicitly before
% registration, then t has to be concatenated BEFORE t0, not after
t = elastix_cat(t, movingT0);

end

%% dbg_plot_reg_results
% debug function to plot the registration result
function dbg_plot_reg_results(imf, imm, fT0, mT0, mT)

% plot input images
subplot(2, 1, 1)
imagesc(imfuse(imf, imm))

% transform fixed image with its initial transform
imf = transformix(fT0, imf);

% apply the registration output transform to the moving image
immreg = transformix(mT, imm);

% transform moving image with its initial transform
imm = transformix(mT0, imm);

% plot registration results
subplot(2, 1, 2)
imagesc(imfuse(imf, immreg))

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [dxdB dxdD dxdS dxdalpha dxdgParam] = gpsimXGradient(model, i, j)

% GPSIMXGRADIENT ...
%
% COPYRIGHT : Pei Gao, 2008

% SHEFFIELDML

% i: ti
% j: jth gene

dxdB = 1./model.D(j);

endPoint = model.times_index(i);

dxdD = 0;
dxdS = 0;
dxdalpha = [];
dxdgParam = [];

if model.ngParam > 0
  ngParamk = model.ngParam/model.numGenes;
  dxdgParam = zeros(1, ngParamk);
  gInd = j;
else
  gInd = 1;
end

for m = 1:endPoint
    arg = model.t(i)-model.mapt(m);
    if arg >= 0
        dxdD = dxdD+model.g(m,gInd)*arg*exp(-model.D(j)*arg+log(model.step)+log(model.S(j)));
        dxdS = dxdS+exp(-model.D(j)*arg+log(model.step))*model.g(m,gInd); % g=Sf

        if model.ngParam > 0
          for gParamInd = 1:ngParamk
            dxdgParam(gParamInd) = dxdgParam(gParamInd) + exp(- ...
            model.D(j)*arg+log(model.step)+log(model.S(j)))*model.dg(m,gInd);
          end
        end
    end
end

dxdD = -model.B(j)/(model.D(j)*model.D(j))-dxdD;

% check if g(f) is different for each gene.
if isfield(model, 'isGroupNonlinearity') && strcmp(model.nonLinearity{j}, ...
                                                   'repression')
      dxdalpha = exp(-model.D(j)*model.t(i));
      dxdD = -model.t(i)*model.alpha(j)*exp(-model.D(j)*model.t(i))+dxdD;
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ M ] = simGaussian( M, sigma )
%SIMGAUSSIAN Calculates Gaussian similarity on matrix
%   simGaussian(M, sigma) returns a matrix of the same size as
%   the distance matrix M, which contains similarity values
%   that are computed by using a Gaussian similarity function
%   with parameter sigma.
%
%   Author: Ingo Buerk
%   Year  : 2011/2012
%   Bachelor Thesis

M = exp(-M.^2 ./ (2*sigma^2));

end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
% CXSparse MATLAB demos.
%
%   cs_demo  - run all CXSparse demos.
%   cs_demo1 - MATLAB version of the CSparse/Demo/cs_demo1.c program.
%   cs_demo2 - MATLAB version of the CSparse/Demo/cs_demo2.c program.
%   cs_demo3 - MATLAB version of the CSparse/Demo/cs_demo3.c program.

% Example:
%   help cs_demo

%   Copyright 2006-2007, Timothy A. Davis

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [cov] = mne_read_noise_cov(fname)
%
% [cov] = mne_read_noise_cov(fname)
%
% Reads a noise-covariance matrix from a fiff file
%
% fname      - The name of the file
%
%

%
%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.5  2006/05/03 19:09:03  msh
%   Fixed even more compatibility issues.
%
%   Revision 1.4  2006/04/23 15:29:41  msh
%   Added MGH to the copyright
%
%   Revision 1.3  2006/04/20 21:49:38  msh
%   Added mne_read_inverse_operator
%   Changed some of the routines accordingly for more flexibility.
%
%   Revision 1.2  2006/04/18 20:44:46  msh
%   Added reading of forward solution.
%   Use length instead of size when appropriate
%
%   Revision 1.1  2006/04/12 17:09:28  msh
%   Added routines for reading noise-covariance matrices
%
%

me='MNE:mne_read_noise_cov';

global FIFF;
if isempty(FIFF)
   FIFF = fiff_define_constants();
end
%
%   Open the file, create directory
%
[ fid, tree ] = fiff_open(fname);
try
   cov = mne_read_cov(fid,tree,FIFF.FIFFV_MNE_NOISE_COV);
catch
   fclose(fid);
   error(me,'%s',mne_omit_first_line(lasterr));
end

return;

end


text:
% THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION IS RELEASED "AS IS."  THE U.S. GOVERNMENT MAKES NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, CONCERNING THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION, INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.  IN NO EVENT WILL THE U.S. GOVERNMENT BE LIABLE FOR ANY DAMAGES, INCLUDING ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE, OR INABILITY TO USE, THIS SOFTWARE OR ANY ACCOMPANYING DOCUMENTATION, EVEN IF INFORMED IN ADVANCE OF THE POSSIBILITY OF SUCH DAMAGES.
%
% integrate_poisspdf.m
% integrates a poisson pdf with respect to lambda from a to b
%
% 022300 tdr created
% 030400 tdr replaced global variable approach to getting x and n to binopdf.
% 030600 tdr added adjustments to a and b to help quad8 converge quickly.
% 			quad8 has trouble with long stretches of near-zero values
% 			so we try to bound the interval to avoid these regions
% 031600 tdr created integrate_poisspdf.m from integrate_binopdf.m
% 012901 tdr changed quad8 to quadl per MatLab 6 recommendation.  Turned off warnings
%   for quadl call because it thinks intervals starting at zero or ending at one
%   involve singularities and generates a warning, although accuracy is okay.
% 031501 tdr changed thres2 from 0.001 to 1.

function y = integrate_poisspdf(x,a,b)

if nargin ~= 3,
    error('Requires three input arguments (x,a,b)');
end

TRACE = 0;
lambda_hat = x;
thres = 1e-6*poisspdf_0ok(x,x); % pdf values less than this are assumed to be negligible.
thres2 = 1; % resolution of adjustments to a and b.
tiny_a = (poisspdf_0ok(x,a) < thres);
tiny_b = (poisspdf_0ok(x,b) < thres);

if tiny_a & tiny_b & (a > lambda_hat | b < lambda_hat),
	%both a and b give results near zero, with no non-zero values between them
   y = 0;
   return
end;

% defaults
new_a = a;
new_b = b;

if tiny_a,
	% a needs to be increased, but not past b or lambda_hat
	ll = a; ul = min(b, lambda_hat);
   while (abs(ll - ul) > thres2),
   	test = (ll + ul)/2;
      if poisspdf_0ok(x,test) < thres,
      	ll = test;
      else
      	ul = test;
      end;
   end;
	new_a = ll;
end;

if tiny_b,
	% b needs to be decreased, but not past a or lambda_hat
	ll = max(a, lambda_hat); ul = b;
   while (abs(ll - ul) > thres2),
   	test = (ll + ul)/2;
      if poisspdf_0ok(x,test) < thres,
      	ul = test;
      else
      	ll = test;
      end;
   end;
	new_b = ul;
end;

% call quad with increased a
%   quadl thinks interval starting at zero or ending at one might involve
%   a singularity and produces a warning, despite being accurate.
warning off;
y=quadl('poisspdf_of_lambda',new_a,new_b,[],TRACE,x);
warning on;


text:
function [ n_data, n, area ] = sphere_unit_area_values ( n_data )

%*****************************************************************************80
%
%% SPHERE_UNIT_AREA_VALUES returns some areas of the unit sphere in ND.
%
%  Discussion:
%
%    The formula for the surface area of the unit sphere in N dimensions is:
%
%      Sphere_Unit_Area ( N ) = 2 * PI**(N/2) / Gamma ( N / 2 )
%
%    Some values of the function include:
%
%       N   Area
%
%       2    2        * PI
%       3  ( 4 /    ) * PI
%       4  ( 2 /   1) * PI**2
%       5  ( 8 /   3) * PI**2
%       6  ( 1 /   1) * PI**3
%       7  (16 /  15) * PI**3
%       8  ( 1 /   3) * PI**4
%       9  (32 / 105) * PI**4
%      10  ( 1 /  12) * PI**5
%
%    For the unit sphere, Area(N) = N * Volume(N)
%
%    In Mathematica, the function can be evaluated by:
%
%      2 * Pi^(n/2) / Gamma[n/2]
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.
%    On input, if N_DATA is 0, the first test data is returned, and
%    N_DATA is set to the index of the test data.  On each subsequent
%    call, N_DATA is incremented and that test data is returned.  When
%    there is no more test data, N_DATA is set to 0.
%
%    Output, integer N, the spatial dimension.
%
%    Output, real AREA, the area of the unit sphere
%    in that dimension.
%
  n_max = 20;

  area_vec = [ ...
     0.2000000000000000E+01, ...
     0.6283185307179586E+01, ...
     0.1256637061435917E+02, ...
     0.1973920880217872E+02, ...
     0.2631894506957162E+02, ...
     0.3100627668029982E+02, ...
     0.3307336179231981E+02, ...
     0.3246969701133415E+02, ...
     0.2968658012464836E+02, ...
     0.2550164039877345E+02, ...
     0.2072514267328890E+02, ...
     0.1602315322625507E+02, ...
     0.1183817381218268E+02, ...
     0.8389703410491089E+01, ...
     0.5721649212349567E+01, ...
     0.3765290085742291E+01, ...
     0.2396678817591364E+01, ...
     0.1478625959000308E+01, ...
     0.8858104195716824E+00, ...
     0.5161378278002812E+00 ];

  n_vec = [ ...
     1, ...
     2, ...
     3, ...
     4, ...
     5, ...
     6, ...
     7, ...
     8, ...
     9, ...
    10, ...
    11, ...
    12, ...
    13, ...
    14, ...
    15, ...
    16, ...
    17, ...
    18, ...
    19, ...
    20 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    area = 0.0;
  else
    n = n_vec(n_data);
    area = area_vec(n_data);
  end

  return
end


text:
function suborder_num = dunavant_suborder_num ( rule )

%*****************************************************************************80
%
%% DUNAVANT_SUBORDER_NUM returns the number of suborders for a Dunavant rule.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 December 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    David Dunavant,
%    High Degree Efficient Symmetrical Gaussian Quadrature Rules
%    for the Triangle,
%    International Journal for Numerical Methods in Engineering,
%    Volume 21, 1985, pages 1129-1148.
%
%    James Lyness, Dennis Jespersen,
%    Moderate Degree Symmetric Quadrature Rules for the Triangle,
%    Journal of the Institute of Mathematics and its Applications,
%    Volume 15, Number 1, February 1975, pages 19-32.
%
%  Parameters:
%
%    Input, integer RULE, the index of the rule.
%
%    Output, integer SUBORDER_NUM, the number of suborders of the rule.
%
  suborder = [ ...
     1,  1 , 2,  2,  3,  3,  4,  5,  6,  6, ...
     7,  8, 10, 10, 11, 13, 15, 17, 17, 19 ];

  if ( 1 <= rule & rule <= 20 )
    suborder_num = suborder(rule);
  else

    suborder_num = -1;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'DUNAVANT_SUBORDER_NUM - Fatal error!\n' );
    fprintf ( 1, '  Illegal RULE = %d\n', rule );
    error ( 'DUNAVANT_SUBORDER_NUM - Fatal error!\n' )

  end

  return
end


text:
clear all;close all;clc

M_dir = {'F:\Janelia2014\Fish1_16states_30frames';
    'F:\Janelia2014\Fish2_20140714_2_4_16states_10frames';
    'F:\Janelia2014\Fish3_20140715_1_1_16_states_10frames';
    'F:\Janelia2014\Fish4_20140721_1_8_16states_20frames';
    'F:\Janelia2014\Fish5_20140722_1_2_16states_30frames';
    'F:\Janelia2014\Fish6_20140722_1_1_3states_30,40frames';
    'F:\Janelia2014\Fish7_20140722_2_3_3states_30,50frames';
    'F:\Janelia2014\Fish8_20141222_2_2_7d_PT_3OMR_shock_lowcut';
    'F:\Janelia2014\Fish9_20150120_2_1_photo_OMR_prey_blob_blue_cy74_6d_20150120_220356'};

poolobj=parpool(4);
%%
for i_fish = 9,%1:8,
    disp(['i_fish = ', num2str(i_fish)]);

    datadir = M_dir{i_fish};
    file = fullfile(datadir,['Fish' num2str(i_fish) '_direct_load.mat']);
    load(file);
    %     varList = {'CR_raw','nCells','CInfo','anat_yx','anat_yz','anat_zx','ave_stack','fpsec','frame_turn','periods'};

    %% detrend
    CR_dtr = zeros(size(CR_raw));
    tmax=size(CR_raw,2);
    ncells=size(CR_raw,1);

    tic
    parfor i=1:ncells,
        cr = CR_raw(i,:);
        crd = 0*cr;
        for j=1:100:tmax,
            if j<=150,
                tlim1 = 1;
                tlim2 = 300;
            elseif j>tmax-150,
                tlim1 = tmax-300;
                tlim2 = tmax;
            else
                tlim1 = j-150;
                tlim2 = j+150;
            end
            crr = cr(tlim1:tlim2);
            crd(max(1,j-50):min(tmax,j+50)) = prctile(crr,15);
        end
        if mod(i,100)==0,
            disp(num2str(i));
        end
        CR_dtr(i,:) = zscore(cr-crd);
    end
    toc

    CR_dtr = single(CR_dtr);

    %% save new version
    varList = {'CR_raw','CR_dtr','nCells','CInfo','anat_yx','anat_yz','anat_zx','ave_stack','fpsec','frame_turn','periods'};
    save(file,varList{:});
end

delete(poolobj);

author: xiuyechen

repo: FishExplorer


text:
function [trx,timestamps] = ReadMWTContours(contournames,ids,varargin)

[dotransposeimage] = myparse(varargin,'dotransposeimage',false);

ncontours = numel(contournames);
% read everthing in from the contour files
trx = [];
for i = 1:ncontours,

  trk = struct('id',ids(i),'xcontour',{{}},'ycontour',{{}});
  timestamps = [];
  fid = fopen(contournames{i},'r');
  j = 1;
  while true,
    s = fgetl(fid);
    if ~ischar(s),
      break;
    end
    d = sscanf(s,'%f');
    if numel(d) < 3,
      warning('Contour line %s length < 3, skipping this line',s);
      continue;
    end
    if mod(numel(d),2) ~= 1,
      warning('Contour line %s length is not odd, skipping this line',s);
      continue;
    end
    d = d';
    timestamp = d(1);
    x = d(2:2:end-1);
    y = d(3:2:end);
    if dotransposeimage,
      trk.xcontour{j} = x+1;
      trk.ycontour{j} = y+1;
    else
      trk.xcontour{j} = y+1;
      trk.ycontour{j} = x+1;
    end
    timestamps(j) = timestamp; %#ok<*AGROW>
    j = j + 1;
  end
  fclose(fid);
  trk.timestamp = timestamps;
  trx = structappend(trx,trk);
end

% make all the timestamps agree with each other
fnsmerge = {'xcontour','ycontour'};
fnscopy = {};
[trx,timestamps] = MergeMWTData(trx,fnsmerge,fnscopy);

author: kristinbranson

repo: JAABA


text:
function tokens = get_tokens(string, delimiter)
% GET_TOKENS	To extract all tokens from a string even if a token is null.
%
% Description:	GET_TOKENS extract all tokens from a string even if one or
%		more tokens are null.  It is a replacement to iteratively calling the
%		builtin function 'strtok.m' because it did not output null tokens as
%		would be normally expected.  Note that if there are two adjacent
%		delimiters then an empty token would be output representing the token
%		that could have been present between the delimiters.
%
% Usage:	tokens = get_tokens(string, delimiter)
%
% Input:
%  string: a string from which to extract tokens.
%  delimiter:	a string containing one or more delimiters. The default
%    delimiters are whitespace characters (i.e. the ASCII codes for TAB,
%    LF, VT, FF, CR and space). Note that each delimiter is assumed to be
%    a single character; so if the delimiter string has three characters,
%    each of the three characters is considered a delimiter.
%
% Output:
%		tokens: An cell array of strings; the tokens retrieved from the input
%			character array.
%
% Example:
%
%		>> s = 'what a great , day';
%
%  Note that the variable 's' is a string that uses two delimiters, a space
%  and a comma.
%
%		>> get_tokens(s, ' ,')
%
%		ans =
%
%       'what'
%       'a'
%       'great'
%       ''
%       ''
%       'day'
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Copyright 2007  Jeff Jackson (Ocean Sciences, DFO Canada)
%   Creation Date: Oct. 20, 2006
%   Last Updated:  Jan. 22, 2008
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Check the input arguments before continuing.

% Check to see if no input arguments were supplied.  If this is the case,
% stop execution and output an error message to the user.
if nargin == 0
  error('MATLAB:get_tokens:NrInputArguments', 'No input arguments were supplied. At least one is expected.');
% Check to see if the only input argument is a cell array.  If it isn't
% then stop execution and output an error message to the user. Also set the
% DIM value since it was not supplied.
elseif nargin == 1
	% Check to see if the input is a string.  Output a message if the input
	% value is a not a character array, then exit the function.
	if ~ischar(string)
		error('MATLAB:get_tokens:InvalidInputArgument', 'The input string must be a valid character array.');
	end
	% Set the default delimiter value to whitespace characters.
	delimiter = [9:13 32];
elseif nargin == 2
	% Check to see if the two input arguments are character strings.  If
	% either check fails then stop execution and output an error message to
	% the user.
	if ~ischar(string)
		error('MATLAB:get_tokens:InvalidInputArgument', 'The first input variable must be a valid character array.');
	end
	if ~ischar(delimiter)
		error('MATLAB:get_tokens:InvalidInputArgument', 'The second input variable must be a valid character array.');
	end
elseif nargin > 2
	% Check to see if too many arguments were input.  If there were then exit
	% the function issuing a error message to the user.
	error('MATLAB:get_tokens:TooManyInputArguments', 'Too many input arguments were supplied.  The maximum permitted is two.');
end

%% If the input arguments are valid continue processing.

% Get the sizes of the input parameters.
string_size = size(string);
delim_size = size(delimiter);

% Check to see if the character array is vertical instead of horizontal. If
% it is then change it. Do the opposite for the delimiter array.
if string_size(1) > 1
	s = string';
else
	s = string;
end
if delim_size(1) > 1
	d = delimiter;
else
	d = delimiter';
end

% Get the array sizes.
if isempty(d)
   error('MATLAB:get_tokens:InputArguments', 'Delimiter is invalid.');
end

% Surround the delimiter string with square brackets so the regexprep
% function works correctly.
delimiter2 = ['[' delimiter ']'];

% Change all delimiters to the character value of the number 1.
s2 = regexprep(s,delimiter2,char(1));

% Extract all tokens from this updated string.
tokens = strread(s2,'%s','delimiter',char(1));

% Check to see if the last character in the input string is a delimiter. If
% it was then add a cell with an empty string to the end of the tokens cell
% array.
if strcmp(s2(end), char(1)) == 1
	tokens(end+1,1) = {''};
end


text:
function J = lotkaObjectiveFCN_models(u,x,N,xref,u0,p,Q,R,Ru, select_model)
%% Cost function of nonlinear MPC for Lotka-Volterra system
%
% Inputs:
%   u:      optimization variable, from time k to time k+N-1
%   x:      current state at time k
%   Ts:     controller sample time
%   N:      prediction horizon
%   xref:   state references, constant from time k+1 to k+N
%   u0:     previous controller output at time k-1
%
% Output:
%   J:      objective function cost
%

%% Nonlinear MPC design parameters
% Q = diag([1,1,1]);
% R = 0.01;

%% Integrate system
if strcmp(select_model,'DelayDMDc')
    [xk,~] = lsim(p.sys,[p.udelay(1:N);u'],[0:N-1].*p.dt,[p.xdelay(:,1); x]-[p.xmean; p.xmean]);
    xk = xk(:,3:4);
    xk = xk + repmat(p.xmean',[N 1]); xk = xk';
elseif strcmp(select_model,'DMDc')
    [xk,~] = lsim(p.sys,[u' 0],[0:N].*p.dt,x-p.xmean);
    xk = xk(2:end,:) + repmat(p.xmean',[N 1]); xk = xk';
elseif strcmp(select_model,'SINDYc')
    Ns = size(x,1);
    xk = zeros(Ns,N+1); xk(:,1) = x;
    for ct=1:N
        % Obtain plant state at next prediction step.
        xk(:,ct+1) = rk4u(@sparseGalerkinControl_Discrete,xk(:,ct),u(ct),p.dt,1,[],p);
    end
    xk = xk(:,2:N+1);
elseif strcmp(select_model,'NARX')
    Hu = [u',0];
    Hx = zeros(2,length(Hu)); Hx(:,1) = x;
    [Us,Ui,Si] = preparets(p.net,con2seq(Hu),{},con2seq(Hx));
    xk = p.net(Us,Ui,Si);
    xk = cell2mat(xk);
end

%% Cost Calculation
% Set initial plant states, controller output and cost
uk = u(1);
J = 0;
% Loop through each prediction step
for ct=1:N
    % Obtain plant state at next prediction step
    xk1 = xk(:,ct);

    % Accumulate state tracking cost from x(k+1) to x(k+N)
    J = J + (xk1-xref)'*Q*(xk1-xref);
    % Accumulate MV rate of change cost from u(k) to u(k+N-1)
    if ct==1
        J = J + (uk-u0)'*R*(uk-u0) + uk'*Ru*uk;
    else
        J = J + (uk-u(ct-1))'*R*(uk-u(ct-1)) + uk'*Ru*uk;
    end
    % Update uk for the next prediction step
    if ct<N
        uk = u(ct+1);
    end
end

author: eurika-kaiser

repo: SINDY-MPC


text:
function val = straight_line_integral_2(xs,ys,zs,data3d,x1,y1,z1,x2,y2,z2)
%
%   Deshan Yang, PhD
%	Department of radiation oncology
%	Washington University in Saint Louis
%   01/16/2011, Saint Louis, MO, USA
%
%   val = straight_line_integral_2(xs,ys,zs,data3d,x1,y1,z1,x2,y2,z2)
%
%   Perform straight line integration on the 3D data from point 1 to point
%   2
%
%   Deshan Yang, PhD
%   01/16/2011, Saint Louis, MO, USA
%
%   This function assumes that the point positions are given in the same
%   coordinate system as the 3D data. Any coordination transformation or
%   scaling should be done priorly.
%
%   Both data and points should be in single, or double
%
dim = [size(data3d) 1];
% if numel(xs)~=length(xs) || length(xs) ~= dim(2) || ...
% 		numel(ys)~=length(ys) || length(ys) ~= dim(1) || ...
% 		numel(zs)~=length(zs) || length(zs) ~= dim(3)
% 	fprintf('Incorrect image data: length does not match to the image dimension.\n');
% 	return;
% end

% Make sure the image and the voxel coordinate be positively increasing. If not, then flip it
[xs,ys,zs,data3d] = make_sure_positive(xs,ys,zs,data3d);
psrc = [x1,y1,z1];
pcdet = [x2,y2,z2];

corner_points = [...
	xs(1) ys(1) zs(1);...
	xs(1) ys(1) zs(end);...
	xs(1) ys(end) zs(1);...
	xs(1) ys(end) zs(end);...
	xs(end) ys(1) zs(1);...
	xs(end) ys(1) zs(end);...
	xs(end) ys(end) zs(1);...
	xs(end) ys(end) zs(end)...
	];

corner_points_t = zeros(8,1);
for k =1:8
	p = corner_points(k,:);
	[q,corner_points_t(k)] = project_1_point_to_a_line(psrc,pcdet,p);
end

t_min = min(corner_points_t);
t_max = max(corner_points_t);

dist = norm(psrc-pcdet);

N = ceil(dist*(t_max-t_min));	% N planes
dt = (t_max-t_min)/N;
ts = t_min:dt:t_max;

xst = psrc(1)+(pcdet(1)-psrc(1))*ts;
yst = psrc(2)+(pcdet(2)-psrc(2))*ts;
zst = psrc(3)+(pcdet(3)-psrc(3))*ts;

val = sum(interp3(xs,ys,zs,data3d,xst,yst,zst,'linear',0))*dt;

% % Transform the point coordinate into the image voxel index
% x1 = (x1 - xs(1))/dx+0.5;
% x2 = (x2 - xs(1))/dx+0.5;
% y1 = (y1 - ys(1))/dy+0.5;
% y2 = (y2 - ys(1))/dy+0.5;
% z1 = (z1 - zs(1))/dz+0.5;
% z2 = (z2 - zs(1))/dz+0.5;
%
%
% nx = ceil(min(x1,x2)):floor(max(x1,x2));
% tx = (nx-x1)/(x2-x1);
%
% ny = ceil(min(y1,y2)):floor(max(y1,y2));
% ty = (ny-y1)/(y2-y1);
%
% nz = ceil(min(z1,z2)):floor(max(z1,z2));
% tz = (nz-z1)/(z2-z1);
%
%
%
% ts = sort([tx ty tz]);    % It may be slow here
% dist = sqrt((x2-x1).^2+(y2-y1).^2+(z2-z1).^2);
%
% nx = (x2-x1)*ts+x1;
% ny = (y2-y1)*ts+y1;
% nz = (z2-z1)*ts+z1;
%
% idxes = find(nx>=0 & nx<=dim(2) & ny>=0 & ny<=dim(1) & nz>=0 & nz<=dim(3));
% nx = nx(idxes);
% ny = ny(idxes);
% nz = nz(idxes);
% N_1 = length(nx)-1;
%
% xidx = max(ceil((nx(1:N_1)+nx(2:end))/2),1);
% yidx = max(ceil((ny(1:N_1)+ny(2:end))/2),1);
% zidx = max(ceil((nz(1:N_1)+nz(2:end))/2),1);
%
% if ~isempty(xidx)
%     if x2~=x1
%         lengths = diff(nx)/(x2-x1);
%     else
%         lengths = diff(ny)/(y2-y1);
%     end
%
%     val = sum(data3d(sub2ind(dim,yidx,xidx,zidx)).*lengths)*dist;
% 	val = val*dx*dy*dz;	% Convert the length back to the original image coordinate length
% else
%     val = 0;
% end
%
%
%
%


text:
clear all
close all
clc

load ../data/lc_data_20s_withpoints

lc_len = length(lc_data);

for i = 1: lc_len
    temp = lc_data{i};
    pp   = points(i,:);
    lat  = temp.veh_s.x;
    lon  = temp.veh_s.y;
    figure(1);hold off;plot(lon,lat,'b.');hold on;
    plot(lon(pp(1)),lat(pp(1)),'o','MarkerFaceColor','g');
    plot(lon(pp(2)),lat(pp(2)),'o','MarkerFaceColor','r');
end

author: donnydcy


text:
function r = sup(a)
%SUP          right bound of point matrix (only for completeness)
%
%   r = sup(a)
%

% written  10/16/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
%

  r = a;

repo: OpenMAS


text:
function P = Predict(Theta, X)

n = size(X, 1);

P = zeros(n, 1);

for i = 1 : n,

    if Sigmoid(X * Theta)(i) >= 0.5

        P(i) = 1;

    end

end

end

author: TheAlgorithms

repo: MATLAB-Octave


text:
function [pbest,perror,nchi2]=cnonlin(mfunc,x,y,sy,pt,v);
% function [pbest,perror,nchi2]=cnonlin(mfunc,x,y,sy,pt,v);
%
% Levenberg-Marquardt non-linear regression for complex functions.
% mfunc=name of the function file that calculates the function to be
%    fitted to the set of data. mfunc is a string variable, so
%    the name of the file must be put inside quotes, eg. 'model'.
%    The function file written to calculate your function should
%    be of the form:
%
%    function yfit=model(x,p)
%    yfit=p(1)+p(2)*x;
%
%    The example above implements y=a+bx.
%    The vector x contains the values of the independent variable.
%    The vector p contains whatever parameters you are using.
%
% x=vector of the independent variable (e.g. Time).
% y=vector of the dependent variable (e.g. Counts recorded).
% sy=vector of standard error in y values.
% pt=initial estimate of parameters to be fitted. All parameters
%   must be non-zero, and are assumed to be real.
% v=vector indicating which parameters are to be varied (fitted), and
%   which are to be held fixed. In most cases we will want all parameters
%   to be varied, but in some circumstances it is useful to hold certain
%   parameters fixed while the others are varied. v is a vector of the same
%   length as p, and should contain only ones and zeros; a one indicates that
%   the corresponding parameter should be varied (fitted), and a zero
%   indicates that the corresponding parameter should be held constant.
%   eg. [1,0,1,1] would be used to keep the second parameter fixed while
%   varying the rest (in the case of four parameters).
%
% OUTPUT VARIABLES
% pbest=fitted parameters
% perror=standard error in fitted parameters
% nchi2=normalised chi-squared parameter (=chi-Squared/degrees of freedom). This
%    is expected to be close to 1 for a good fit.

% Written 4/12/95 by Michael Fleming, University Of Auckland

chiCut=0.01;  lambda=0.001;  stepSize=0.001;  chiOld=Inf;
y=y(:);  x=x(:);  sy=sy(:);  pt=pt(:);  v=v(:);  maxIter=50;
m=length(pt); n=length(x);

evalstr=[mfunc '(x,['];      % Create string containing function call
count=0;
for i=1:m
 if v(i)==1         % Extract vector p of parameters which will be varied
  count=count+1;    % from vector pt which contains all input parameters
  p(count)=pt(i);
  evalstr=[evalstr 'p(' int2str(count) ') '];
 else
  evalstr=[evalstr num2str(pt(i)) ' '];
 end;
 if (pt(i)==0)   % Check all parameters are non-zero
  error('You have entered a zero parameter.All parameters must be nonzero');
 end;
end;
evalstr=[evalstr '])'];

dof=n-count;    p=p(:);   delta=p*stepSize;
% delta is the size of the small increment used for calculating numerical derivatives
% dof is the number of degrees of freedom

chiSqr=sum((abs(y-eval(evalstr))./sy).^2);     % Calc ChiSqr from initial parameters
if chiSqr/dof>5000;
 disp('You have made a bad choice of initial parameters');
end;

iter=0;
while (abs(chiOld-chiSqr)>chiCut)&(iter<maxIter);
 iter=iter+1;
 chiOld=chiSqr;
 [alpha,beta]=cnonlin2(evalstr,x,y,sy,p,delta,lambda);
 if det(alpha)==0;
  error('No convergence - try a different set of parameters');
 end;
 dp=alpha\beta;    % Evaluate parameter increments
 p=p+dp;
 chiSqr=sum((abs(y-eval(evalstr))./sy).^2);
 while (chiSqr>(chiOld+chiCut));
  p=p-dp;
  iter=iter+1;  lambda=lambda*10;
  [alpha,beta]=cnonlin2(evalstr,x,y,sy,p,delta,lambda);
  if det(alpha)==0;
   error('No convergence - try a different set of parameters');
  end;
  dp=alpha\beta;    % Evaluate parameter increments
  p=p+dp;
  chiSqr=sum((abs(y-eval(evalstr))./sy).^2);
 end;
 lambda=0.1*lambda;
end;

if iter==maxIter,
 disp('Maximum number of iterations exceeded - convergence not achieved');
end;

%xfine=linspace(x(1),x(n),201);    % Print graph of function
%yfitted=model(xfine,p);
%errorbar(x,y,sy);
% Calculate correlation coefficicent R^2
%f=eval(evalstr);
%r=corrcoef(y,f);
%r2=r(1,2).^2;

% Calculate standard errors in parameters
[alpha,beta]=cnonlin2(evalstr,x,y,sy,p,delta,0);	% Calculate final alpha matrix
sp=sqrt(diag(inv(alpha)));   % Evaluate standard errors in parameters.
disp('Parameters');
count=0;
for i=1:m
 if v(i)==1
  count=count+1;
  pbest(i)=p(count); perror(i)=sp(count);
  disp(['p(', int2str(i), ') = ' sprintf('%5.3e',p(count)), '+-' sprintf('%5.3e',sp(count))]);
 else
  disp(['p(', int2str(i), ') = ' sprintf('%5.3e',pt(i)), '+-0']);
  pbest(i)=pt(i); perror(i)=0;
 end;
end;

format compact
%disp('Correlation Coefficient R^2');
%disp(r2);
disp('Normalised chi-squared value');
nchi2=chiSqr/dof;
disp(nchi2);
format;

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
classdef SteeringModelEnum < matlab.mixin.SetGet
    %SteeringModelEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        RollPitchYawPoly('Roll/Pitch/Yaw Steering','RollPitchYawPolySteeringModel', true);
        AeroAnglesPoly('Body Fixed Aero Angles Steering','AeroAnglesPolySteeringModel', true);
        InertialAeroAnglesPoly('Inertial Aero Angles Steering','InertialAeroAnglesPolySteeringModel', true);
        GenericPoly('Generic Angles Steering','GenericPolySteeringModel', true);

        GenericQuatInterp('Generic Angles Quaternion Interp','GenericQuatInterpSteeringModel', false);
        GenericSelectable('Generic Selectable Model Steering','GenericSelectableSteeringModel', false);
    end

    enumeration(Hidden)
        GenericLinTan('Generic Linear Tangent Angles Steering','GenericLinearTangentSteeringModel', false);
        GenericSumOfSines('Generic Sum of Sines Angles Steering','GenericSumOfSinesSteeringModel', false);
    end

    properties
        nameStr char = '';
        classNameStr char = '';
        showInPolyUi(1,1) logical = false;
    end

    methods
        function obj = SteeringModelEnum(nameStr,classNameStr,showInPolyUi)
            obj.nameStr = nameStr;
            obj.classNameStr = classNameStr;
            obj.showInPolyUi = showInPolyUi;
        end
    end

    methods(Static)
        function steeringModelNameStrs = getSteeringModelTypeNameStrs(showInPylUiMustBe)
            [m,~] = enumeration('SteeringModelEnum');
            m = m([m.showInPolyUi] == showInPylUiMustBe);

            steeringModelNameStrs = {};
            for(i=1:length(m)) %#ok<*NO4LP>
                steeringModelNameStrs{end+1} = m(i).nameStr; %#ok<AGROW>
            end
        end

        function ind = getIndOfListboxStrsForSteeringModel(steeringModel)
            [m,~] = enumeration('SteeringModelEnum');
            inputClass = class(steeringModel);

            ind = -1;
            for(i=1:length(m))
                if(strcmpi(m(i).classNameStr,inputClass))
                    ind = i;
                    break;
                end
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [N,R]=repeatv(A)
%[number,repeat]=repeatv(input vector)
%
%these function calculate number of repeat all array of a input vector
%input must be a vector
%vector is var(1,:) or var(:,1)  var==>variable
%these function calculate repeats by deleting column or row
%
%for example:
%[n,r]=repeatv([1 2 3 4 9 6 5 2 3 1])
%n =
%
%     1     2     3     4     9     6     5
%
%
%r =
%
%     2     2     2     1     1     1     1

if isempty(A)
    N=[];
    R=[];
else
    l=length(A);
    i=1;
    while i<=l
        t=0;
        b(1,i)=A(i);
        j=i;
        while j<=l
            if A(i)==A(j)
                t=t+1;
                if i~=j
                    A(j)=[];
                    l=l-1;
                    j=j-1;
                end
            end
            j=j+1;
        end
        b(2,i)=t;
        i=i+1;
    end
    N=b(1,:);
    R=b(2,:);
end


text:
% Copyright 2011 Zdenek Kalal
%
% This file is part of TLD.
%
% TLD 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.
%
% TLD 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 TLD.  If not, see <http://www.gnu.org/licenses/>.

function [nX1,nX2,nEx1,nEx2] = tldSplitNegativeData(nX,nEx)
% Splits negative data to training and validation set

N    = size(nX,2);
idx  = randperm(N);
nX   = nX(:,idx);
nX1  = nX(:,1:N/2);
nX2  = nX(:,N/2+1:end);

N    = size(nEx,2);
idx  = randperm(N);
nEx  = nEx(:,idx);
nEx1 = nEx(:,1:N/2);
nEx2 = nEx(:,N/2+1:end);

repo: OpenTLD


text:
% created: Zoya Bylinskii, Aug 2014

% This finds the normalized scanpath saliency between two different
% saliency maps as the mean value of the normalized saliency map at
% fixation locations.

function score = NSS(saliencyMap, fixationMap)
% saliencyMap is the saliency map
% fixationMap is the human fixation map (binary matrix)
map = double(imresize(saliencyMap,size(fixationMap)));

% normalize saliency map
map = (map - mean(map(:)))/std(map(:));

% mean value at fixation locations
score = mean(map(logical(fixationMap)));

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency


text:
function P = agm_pi ( d )

%*****************************************************************************80
%
%% AGM_PI  Arithmetic-geometric mean for pi.
%
%  Discussion:
%
%    The Brent-Salamin algorithm is used to compute PI to D decimal digits.
%
%  Licensing:
%
%    Copyright (c) 2011, The MathWorks, Inc.
%    All rights reserved.
%
%    Redistribution and use in source and binary forms, with or without
%    modification, are permitted provided that the following conditions are
%    met:
%
%        * Redistributions of source code must retain the above copyright
%          notice, this list of conditions and the following disclaimer.
%        * Redistributions in binary form must reproduce the above copyright
%          notice, this list of conditions and the following disclaimer in
%          the documentation and/or other materials provided with the distribution
%        * Neither the name of the The MathWorks, Inc. nor the names
%          of its contributors may be used to endorse or promote products derived
%          from this software without specific prior written permission.
%
%    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%    AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%    IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%    ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%    LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%    POSSIBILITY OF SUCH DAMAGE.
%
%  Modified:
%
%    27 February 2012
%
%  Author:
%
%    Cleve Moler
%
%  Reference:
%
%    Cleve Moler,
%    Cleve's Corner, "Computing Pi",
%    http://www.mathworks.com/company/newsletters/news_notes/2011/
%
%  Parameters:
%
%    Input, integer D, the number of decimal digits desired.
%
%    Output, symbolic P, the value of pi to D digits.
%
  digits ( d )

  a = vpa(1,d);
  b = 1/sqrt(vpa(2,d));
  s = 1/vpa(4,d);
  p = 1;
  n = ceil(log2(d));

  for k = 1:n
    c = (a+b)/2;
    b = sqrt(a*b);
    s = s - p*(c-a)^2;
    p = 2*p;
    a = c;
  end

  P = a^2 / s;

  return
end


text:
function [coord,R]=calcgrip(Tr,Tc)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%EN FUNCION DE DONDE ESTE LA CAJA Y DONDE ESTE EL ROBOT,
%CALCULA QUE LADO DE LA CAJA ESTA A MENOS DISTANCIA DE LA
%POSICION ACTUAL, Y DEVUELVE POR DONDE TIENE QUE COGER LA CAJA
%Y CON QUE ORIENTACION

    Tc(3,4)=Tc(3,4)+0.05;
    Tc1=Tc(1:3,4);
    Tc2=Tc(1:3,4);
    Tc3=Tc(1:3,4);
    Tc4=Tc(1:3,4);
    Tc1(1,4)=Tc(1,4)+0.14;
    dist(1,1)=norm(Tc1-Tr(1:3,4));
    Tc2(1,4)=Tc(1,4)-0.14;
    dist(2,1)=norm(Tc2-Tr(1:3,4));
    Tc3(1,4)=Tc(2,4)+0.10;
    dist(3,1)=norm(Tc3-Tr(1:3,4));
    Tc4(1,4)=Tc(2,4)-0.10;
    dist(4,1)=norm(Tc4-Tr(1:3,4));
    [M,I]=min(dist);
    if I==1
        coord=[0.14;0;0.05];
        R=[0 0 1 ;0 1 0;1 0 0];
    end
    if I==2
        coord=[-0.14;0;0.05];

        R=[0 0 -1 ;0 1 0;-1 0 0];
    end
    if I==3
        coord=[0;0.10;0.05];
        R=[0 -1 0 ;0 0 1; -1 0 0];

    end
    if I==4
        coord=[0;-0.10;0.05];
        R=[0 1 0 ;0 0 -1; 1 0 0];
    end

end

repo: ARTE


text:
function outputImage = power(this, otherImage)
% Computes input image to the power of other image per pixel as new image
%
% NOTE: If voxel dimensions of 2nd image do not match - or a scalar is
% given as 2nd argument - data in the 2nd argument is automatically
% replicated to match this image geometry.
%
%
%   Y = MrImage()
%   outputImage = power(Y, otherImage, ...
%   functionHandle)
%
%   OR
%
%   outputImage = Y.^otherImage
%
% This is a method of class MrImage.
%
%
% IN
%   otherImage              image that will be subtracted from this one
%
% OUT
%   outputImage             new MrImage, difference of this and otherImage
%
% EXAMPLE
%
%   % Compute difference of 2 images
%		Y = MrImage();
%		Z = MrImage();
%		X = Y.power(Z);
%
%   % OR (cool overload!):
%       X = Y.^Z
%
%   e.g. take square of image
%       X = Y.^2
%
%   e.g. take square root of image
%       X = Y.^(1/2)
%
%   See also MrImage perform_binary_operation MrImage.plus

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2014-11-13
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.
%

outputImage = this.perform_binary_operation(otherImage, @power);

author: translationalneuromodeling

repo: tapas


text:
classdef CasadiSolver < handle

  properties
    timeMeasures
    stageList
    collocationList
    nlpData
    stats
  end

  properties (Access = private)
    controls_regularization
    controls_regularization_value
  end

  methods

    function self = CasadiSolver(stageList, transitionList, ...
                                 nlp_casadi_mx, ...
                                 controls_regularization, controls_regularization_value, ...
                                 casadi_options, verbose, problem_userdata, transition_type)

      ocl.utils.assert(length(stageList)==length(transitionList)+1, ...
                'You need to specify Ns-1 transitions for Ns stages.');

      self.controls_regularization = controls_regularization;
      self.controls_regularization_value = controls_regularization_value;

      constructTotalTic = tic;

      % create variables as casadi symbolics
      if nlp_casadi_mx
        casadi_sym = @casadi.MX.sym;
      else
        casadi_sym = @casadi.SX.sym;
      end

      vars = cell(length(stageList), 1);
      costs = cell(length(stageList), 1);
      constraints = cell(length(stageList), 1);
      constraints_LB = cell(length(stageList), 1);
      constraints_UB = cell(length(stageList), 1);

      v_stage = [];

      collocationList = cell(length(stageList), 1);

      for k=1:length(stageList)
        stage = stageList{k};

        x_struct = stage.x_struct;
        z_struct = stage.z_struct;
        u_struct = stage.u_struct;
        p_struct = stage.p_struct;
        x_order = stage.x_order;

        daefh = stage.daefh;
        pathcostsfh = stage.pathcostsfh;
        gridcostsfh = stage.gridcostsfh;
        gridconstraintsfh = stage.gridconstraintsfh;
        terminalcostfh = stage.terminalcostfh;

        userdata = stage.userdata;

        nx = stage.nx;
        nu = stage.nu;
        np = stage.np;
        d = stage.d;
        H_norm = stage.H_norm;
        T = stage.T;
        N = stage.N;

        if length(stageList) > 1
          name_suffix = ['_s',mat2str(k)];
        else
          name_suffix = '';
        end

        % create casadi primitives for model functions
        x = ocl.casadi.structToSym(x_struct, casadi_sym, name_suffix);
        z = ocl.casadi.structToSym(z_struct, casadi_sym, name_suffix);
        u = ocl.casadi.structToSym(u_struct, casadi_sym, name_suffix);
        p = ocl.casadi.structToSym(p_struct, casadi_sym, name_suffix);

        % casadi dae function
        [casadi_ode_sym, casadi_alg_sym] = ocl.model.dae( ...
          daefh, x_struct, z_struct, u_struct, p_struct, x_order, ...
          x, z, u, p, userdata);
        casadi_dae_fun = casadi.Function('odefun', {x,z,u,p}, {casadi_ode_sym, casadi_alg_sym});
        daefun = @(x,z,u,p) ocl.casadi.daefun(casadi_dae_fun,x,z,u,p);

        % casadi pathcost function
        pathcosts_sym = ocl.model.pathcosts(...
          pathcostsfh,x_struct, z_struct, u_struct, p_struct, ...
          x, z, u, p, userdata);
        casadi_pathcost_fun = casadi.Function('pathcosts', {x,z,u,p}, {pathcosts_sym});
        pathcostfun = @(x,z,u,p) ocl.casadi.pathcostfun(casadi_pathcost_fun,x,z,u,p);

        collocation = ocl.collocation.Collocation(x_struct, z_struct, u_struct, p_struct, x_order, daefun, pathcostfun, d);
        ni = collocation.num_i;

        x = casadi_sym(['x','_s',mat2str(k)], nx);
        vi = casadi_sym(['vi','_s',mat2str(k)], ni);
        u = casadi_sym(['u','_s',mat2str(k)], nu);
        h = casadi_sym(['h','_s',mat2str(k)]);
        p = casadi_sym(['p','_s',mat2str(k)], np);

        [xF, cost_integr, equations] = ocl.collocation.equations(collocation, x, vi, u, h, p);
        integrator_fun = casadi.Function('sys', {x,vi,u,h,p}, {xF, cost_integr, equations});

        integratormap = integrator_fun.map(N, 'serial');

        nv_stage = ocl.simultaneous.nvars(N, nx, ni, nu, np);
        v_last_stage = v_stage;
        v_stage = casadi_sym(['v','_s',mat2str(k)], nv_stage);

        gridcostfun = @(k,K,x,p) ocl.model.gridcosts(gridcostsfh, x_struct, p_struct, k, K, x, p, userdata);
        gridconstraintfun = @(k,K,x,p) ocl.model.gridconstraints(gridconstraintsfh, x_struct, p_struct, k, K, x, p, userdata);
        terminalcostfun = @(x,p) ocl.model.terminalcost(terminalcostfh, x_struct, p_struct, x, p, userdata);

        [costs_stage, constraints_stage, ...
         constraints_LB_stage, constraints_UB_stage] = ...
            ocl.simultaneous.equations(H_norm, T, nx, ni, nu, np, ...
                                       gridcostfun, gridconstraintfun, ...
                                       terminalcostfun, ...
                                       integratormap, v_stage, ...
                                       controls_regularization, controls_regularization_value);

        collocationList{k} = collocation;

        transition_eq = [];
        transition_lb = [];
        transition_ub = [];
        if k >= 2
          [x0_cur, p_cur] = ocl.simultaneous.getFirstState(stageList{k}, collocationList{k}, v_stage);
          [xF_prev, p_prev] = ocl.simultaneous.getLastState(stageList{k-1}, collocationList{k-1}, v_last_stage);
          transition_fun = transitionList{k-1};

          tansition_handler = ocl.Constraint(problem_userdata);

          x0_cur = ocl.Variable.create(stageList{k}.x_struct, x0_cur);
          xF_prev = ocl.Variable.create(stageList{k-1}.x_struct, xF_prev);

          p_cur = ocl.Variable.create(stageList{k}.p_struct, p_cur);
          p_prev = ocl.Variable.create(stageList{k-1}.p_struct, p_prev);

          if transition_type == 1
            transition_fun(tansition_handler,x0_cur,xF_prev);
          elseif transition_type == 2
            transition_fun(tansition_handler,x0_cur,xF_prev,p_cur,p_prev);
          else
            ocl.error('Transition type invalid.');
          end

          transition_eq = tansition_handler.values;
          transition_lb = tansition_handler.lowerBounds;
          transition_ub = tansition_handler.upperBounds;
        end

        vars{k} = v_stage;
        costs{k} = costs_stage;
        constraints{k} = vertcat(transition_eq, constraints_stage);
        constraints_LB{k} = vertcat(transition_lb, constraints_LB_stage);
        constraints_UB{k} = vertcat(transition_ub, constraints_UB_stage);

      end

      v = vertcat(vars{:});
      costs = sum([costs{:}]);
      constraints = vertcat(constraints{:});
      constraints_LB = vertcat(constraints_LB{:});
      constraints_UB = vertcat(constraints_UB{:});

      % get struct with nlp for casadi
      casadiNLP = struct;
      casadiNLP.x = v;
      casadiNLP.f = costs;
      casadiNLP.g = constraints;
      casadiNLP.p = [];

      if ~verbose
        casadi_options.ipopt.print_level = 0;
        casadi_options.print_time = 0;
      end

      constructSolverTic = tic;
      casadiSolver = casadi.nlpsol('my_solver', 'ipopt', casadiNLP, casadi_options);
      constructSolverTime = toc(constructSolverTic);

      nlpData = struct;
      nlpData.casadiNLP = casadiNLP;
      nlpData.constraints_LB = constraints_LB;
      nlpData.constraints_UB = constraints_UB;
      nlpData.solver = casadiSolver;

      timeMeasures.constructTotal = toc(constructTotalTic);
      timeMeasures.constructSolver = constructSolverTime;

      self.stageList = stageList;
      self.nlpData = nlpData;
      self.timeMeasures = timeMeasures;
      self.collocationList = collocationList;
    end

    function igMerged = getInitialGuessWithUserData(self)

      stage_list = self.stageList;
      colloc_list = self.collocationList;

      igMerged = cell(length(stage_list), 1);

      ig = self.getInitialGuess();
      for k=1:length(stage_list)
        stage = stage_list{k};
        colloc = colloc_list{k};

        igMerged{k} = ocl.simultaneous.getInitialGuessWithUserData(ig{k}, stage, colloc);
      end
    end

    function igList = getInitialGuess(self)
      stage_list = self.stageList;

      igList = cell(length(stage_list),1);
      for k=1:length(stage_list)
        stage = stage_list{k};

        colloc = self.collocationList{k};

        N = stage.N;
        x_struct = stage.x_struct;
        z_struct = stage.z_struct;
        u_struct = stage.u_struct;
        p_struct = stage.p_struct;

        nx = stage.nx;
        nu = stage.nu;
        np = stage.np;
        H_norm = stage.H_norm;
        T = stage.T;

        vi_struct = colloc.vars;
        ni = colloc.num_i;

        x_bounds = stage.x_bounds;
        x0_bounds = stage.x0_bounds;
        xF_bounds = stage.xF_bounds;
        z_bounds = stage.z_bounds;
        u_bounds = stage.u_bounds;
        p_bounds = stage.p_bounds;

        [x_lb, x_ub] = ocl.model.bounds(x_struct, x_bounds);
        [x0_lb, x0_ub] = ocl.model.bounds(x_struct, x0_bounds);
        [xF_lb, xF_ub] = ocl.model.bounds(x_struct, xF_bounds);

        [z_lb, z_ub] = ocl.model.bounds(z_struct, z_bounds);
        [u_lb_traj, u_ub_traj] = ocl.model.boundsTrajectory(u_struct, u_bounds, N);
        [p_lb, p_ub] = ocl.model.bounds(p_struct, p_bounds);

        varsStruct = ocl.simultaneous.variablesStruct(N, x_struct, vi_struct, u_struct, p_struct);

        ig = ocl.simultaneous.getInitialGuess(H_norm, T, nx, ni, nu, np, ...
                                    x0_lb, x0_ub, xF_lb, xF_ub, x_lb, x_ub, ...
                                    z_lb, z_ub, u_lb_traj, u_ub_traj, p_lb, p_ub, ...
                                    vi_struct);

        igList{k} = ocl.Variable.create(varsStruct, ig);
      end

    end

    function [sol,times,info] = solve(self,v0)
      % solve(initialGuess)

      solveTotalTic = tic;

      stage_list = self.stageList;
      collocation_list = self.collocationList;

      ig_list = v0;

      for k=1:length(v0)
        v0{k} = v0{k}.value;
      end

      lbv = cell(length(stage_list),1);
      ubv = cell(length(stage_list),1);
      for k=1:length(stage_list)

        stage = stage_list{k};
        colloc = collocation_list{k};

        nx = stage.nx;
        nu = stage.nu;
        np = stage.np;
        H_norm = stage.H_norm;
        T = stage.T;
        N = length(H_norm);

        x_struct = stage.x_struct;
        z_struct = stage.z_struct;
        u_struct = stage.u_struct;
        p_struct = stage.p_struct;

        x_bounds = stage.x_bounds;
        x0_bounds = stage.x0_bounds;
        xF_bounds = stage.xF_bounds;
        z_bounds = stage.z_bounds;
        u_bounds = stage.u_bounds;
        p_bounds = stage.p_bounds;

        vi_struct = colloc.vars;
        ni = colloc.num_i;

        [x_lb, x_ub] = ocl.model.bounds(x_struct, x_bounds);
        [x0_lb, x0_ub] = ocl.model.bounds(x_struct, x0_bounds);
        [xF_lb, xF_ub] = ocl.model.bounds(x_struct, xF_bounds);

        [z_lb, z_ub] = ocl.model.bounds(z_struct, z_bounds);
        [u_lb_traj, u_ub_traj] = ocl.model.boundsTrajectory(u_struct, u_bounds, N);
        [p_lb, p_ub] = ocl.model.bounds(p_struct, p_bounds);

        [vi_lb, vi_ub] = ocl.collocation.bounds(vi_struct, x_lb, x_ub, z_lb, z_ub);

        [lbv_stage, ubv_stage] = ocl.simultaneous.bounds(H_norm, T, nx, ni, nu, np, ...
                                      x_lb, x_ub, x0_lb, x0_ub, xF_lb, xF_ub, ...
                                      vi_lb, vi_ub, u_lb_traj, u_ub_traj, p_lb, p_ub);
        lbv{k} = lbv_stage;
        ubv{k} = ubv_stage;
      end

      v0 = vertcat(v0{:});
      lbv = vertcat(lbv{:});
      ubv = vertcat(ubv{:});

      args = struct;
      args.lbg = self.nlpData.constraints_LB;
      args.ubg = self.nlpData.constraints_UB;
      args.p = [];
      args.lbx = lbv;
      args.ubx = ubv;
      args.x0 = v0;

      % execute solver
      solveCasadiTic = tic;
      sol = self.nlpData.solver.call(args);
      solveCasadiTime = toc(solveCasadiTic);

      if strcmp(self.nlpData.solver.stats().return_status, 'NonIpopt_Exception_Thrown')
        ocl.utils.warning('Solver was interrupted by user.');
      end

      sol_values = sol.x.full();

      sol = cell(length(stage_list),1);
      times = cell(length(stage_list),1);
      i_stage = 1;
      nlpFunEvalTic = tic;
      for k=1:length(stage_list)

        stage = stage_list{k};
        colloc = collocation_list{k};

        nx = stage.nx;
        nu = stage.nu;
        np = stage.np;
        N = stage.N;
        H_norm = stage.H_norm;

        ni = colloc.num_i;
        nt = colloc.num_t;

        nv_stage = ocl.simultaneous.nvars(N, nx, ni, nu, np);

        % unpack solution of this stage to state/controls trajectories
        V = sol_values(i_stage:i_stage+nv_stage-1);
        sol_out = ocl.Variable.create(ig_list{k}.type, V);

        [~,~,~,~,H] = ocl.simultaneous.variablesUnpack(V, N, nx, ni, nu, np);
        colloc_times = ocl.simultaneous.times(H(1)*H_norm, colloc);
        times_struct = ocl.simultaneous.timesStruct(N, nt);
        times_out = ocl.Variable.create(times_struct, colloc_times);

        i_stage = i_stage + nv_stage;

        sol{k} = sol_out;
        times{k} = times_out;
      end
      nlpFunEvalTime = toc(nlpFunEvalTic);

      self.timeMeasures.solveTotal      = toc(solveTotalTic);
      self.timeMeasures.solveCasadi     = solveCasadiTime;
      self.timeMeasures.nlpFunEval      = nlpFunEvalTime;

      self.stats = self.nlpData.solver.stats();

      info = self.info();
    end

    function r = info(self)
      r = struct;
      r.timeMeasures = self.timeMeasures;
      r.success = self.stats.success;
      r.ipopt_stats = self.stats;
    end
  end

end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
%% Structured grid
% This is a sample code showing how to use ConstructProjector3D().
% This example is going to setup a transformation matrix based on forth
% order polynomials in 3D.
% Pay attention that no data is needed to construct the interpolator. Only
% the coordinates of the points on source and destination grid.
clear;clc;
%% Initializing Part I
disp('- Initializing')
xMin=0;
xMax=2*pi;
yMin=0;
yMax=2*pi;
zMin=0;
zMax=2*pi;

nx=20;
ny=20;
nz=20;

dx=(xMax-xMin)/(nx-1);
dy=(yMax-yMin)/(ny-1);
dz=(zMax-zMin)/(nz-1);

nPoly=1; % We gonna fit a fourth order polynomial,
         % i.e. sum_{n,m,t=0}^{n,m,t=nPoly} x^n*y^m*z^t
nInterp=8; % nPoly=4 requires 125 points. However,
            % we are going to use 150 points. This would be
            % the least-square fit of the surface.
            % Note that: min(nInterp)=(nPoly+1)^3. Otherwise it won't work.
%% Turning off the warnings
% Depending on the parameter that you choose the matrices might be nearly
% singular. In this example, despite being nearly singular, they are not
% affecting the final interpolation. Therefore, I set the warning to off.
% But generally do not turn it off and if you got this warning check if it
% is going to affect your results or not.
disp('- Turning the warnings off')
warning('OFF','MATLAB:nearlySingularMatrix');

%% Generating the Source Grid
disp('- Generating the source grid.')
[xn,yn,zn]=meshgrid(linspace(xMin,xMax,nx),linspace(yMin,yMax,ny),linspace(zMin,zMax,nz));

%% Finding the cell centers
disp('- Generating the destination grid')
xc=xn(1:ny-1,1:nx-1,1:nz-1)+dx/2;
yc=yn(1:ny-1,1:nx-1,1:nz-1)+dy/2;
zc=zn(1:ny-1,1:nx-1,1:nz-1)+dz/2;

figure
plot3(xn(:),yn(:),zn(:),'k.','MarkerSize',20);
hold on
plot3(xc(:),yc(:),zc(:),'b.','MarkerSize',20);
axis tight;
axis square;
title('Source Grid, black dots, & the destination grid, blue dots.', ...
      'FontName','Arial','FontSize',12,'FontWeight','Bold');

%% Constructing the Interpolant
% Note: that we do not need the data on the source grid to create the
% interpolant.
disp('- Constructing the interpolant')
P=ConstructProjector3D(xn(:),yn(:),zn(:),xc(:),yc(:),zc(:),nPoly,nInterp);

%% Generarting some Data
disp('- Generating some data and interpolating')
F1=@(x,y,z) (sin(sqrt(x.^2+y.^2+z.^2)));
Vn=F1(xn,yn,zn);
Vc_interp=reshape(P*Vn(:),size(xc));
Vc_Analytic=F1(xc,yc,zc);
RMSE=sqrt(mean((Vc_Analytic(:)-Vc_interp(:)).^2));
disp(['RMSE: ' num2str(RMSE)])

%% Generating some more data
disp('- Generating multiple data field on the source grid and interpolating them.')
F2=@(x,y,z) (sin(x).*cos(y).*sin(z));
F3=@(x,y,z) (exp(-sqrt(x.^2+y.^2+z.^2)));
F4=@(x,y,z,x0,y0,z0) (exp(-sqrt((x-x0).^2+(y-y0).^2+(z-z0).^2)));

Vn(:,:,:,1)=F1(xn,yn,zn);
Vn(:,:,:,2)=F2(xn,yn,zn);
Vn(:,:,:,3)=F3(xn,yn,zn);
Vn(:,:,:,4)=F4(xn,yn,zn,mean(mean(mean(xn))),mean(mean(mean(yn))),mean(mean(mean(zn))));

Vc_Analytic(:,:,:,1)=F1(xc,yc,zc);
Vc_Analytic(:,:,:,2)=F2(xc,yc,zc);
Vc_Analytic(:,:,:,3)=F3(xc,yc,zc);
Vc_Analytic(:,:,:,4)=F4(xc,yc,zc,mean(mean(mean(xn))),mean(mean(mean(yn))),mean(mean(mean(zn))));

% Now interpolating
% Note that the same interpolant is used for all data fields and they can
% be all interpolated with one sparse matrix multiplication.
Vc_interp=P*reshape(Vn,nx*ny*nz,4); % There are 4 data fields.
Vc_interp=reshape(Vc_interp,(ny-1),(nx-1),(nz-1),4); % these two commands can be combined in one.
                                              % They were separated for clarity.

% calculating the RMSE and plotting
RMSE=zeros(4,1);
for i=1:4
  RMSE(i)=sqrt(mean((reshape(Vc_Analytic(:,:,:,i),(nx-1)*(ny-1)*(nz-1),1)-reshape(Vc_interp(:,:,:,i),(nx-1)*(ny-1)*(nz-1),1)).^2));
  disp(['F' num2str(i) ', nPoly: ' num2str(nPoly) ', RMSE:' num2str(RMSE(i))])
end

repo: mcbench-benchmarks