text:
classdef MimLinkedDatasetChooser < CoreBaseClass
    % MimLinkedDatasetChooser. Part of the internal framework of the TD MIM Toolkit.
    %
    %     You should not use this class within your own code. It is intended to
    %     be used internally within the framework of the TD MIM Toolkit.
    %
    %     MimLinkedDatasetChooser is used to select between linked datasets.
    %     By default, each dataset acts independently, but you can link datasets
    %     together (for example, if you wanted to register images between two
    %     datasets). When datasets are linked, one is the primary dataset, and
    %     linked results are stored in the primary cache. The primary dataset
    %     may access results from any of its linked datasets (but not vice
    %     versa). Linking can be nested.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Access = private)
        LinkedRecorderSingleton
        DatasetCache
        PrimaryDatasetResults % Handle to the MimDatasetResults object for this dataset
        LinkedDatasetChooserList % Handles to MimLinkedDatasetChooser objects for all linked datasets, including this one
        PrimaryDatasetUid     % The uid of this dataset
    end

    events
        % This event is fired when a plugin has been run for this dataset, and
        % has generated a new preview thumbnail.
        PreviewImageChanged
    end

    methods
        function obj = MimLinkedDatasetChooser(framework_app_def, context_def, image_info, dataset_disk_cache, linked_recorder_singleton, plugin_cache, reporting)
            obj.LinkedRecorderSingleton = linked_recorder_singleton;
            obj.DatasetCache = dataset_disk_cache;
            primary_dataset_results = MimDatasetResults(framework_app_def, context_def, image_info, obj, obj, dataset_disk_cache, plugin_cache, reporting);
            obj.PrimaryDatasetUid = primary_dataset_results.GetImageInfo.ImageUid;
            obj.PrimaryDatasetResults = primary_dataset_results;
            obj.LinkedDatasetChooserList = containers.Map;
            obj.LinkedDatasetChooserList(obj.PrimaryDatasetUid) = obj;
        end

        function AddLinkedDataset(obj, linked_name, linked_dataset_chooser, reporting)
            % Links a different dataset to this one, using the specified name.
            % The name exists only within the scope of this dataset, and is used
            % to identify the linked dataset from which results should be
            % obtained.

            linked_uid = linked_dataset_chooser.PrimaryDatasetUid;
            obj.LinkedDatasetChooserList(linked_uid) = linked_dataset_chooser;
            obj.LinkedDatasetChooserList(linked_name) = linked_dataset_chooser;

            obj.LinkedRecorderSingleton.AddLink(obj.PrimaryDatasetUid, linked_uid, linked_name, reporting);
        end

        function dataset_results = GetDataset(obj, reporting, varargin)
            % Returns a handle to the DatasetResults object for a particular linked dataset.
            % The dataset is identified by its uid in varargin, or an empty
            % input will return the primary dataset.

            if nargin < 3
                dataset_name = [];
            else
                dataset_name = varargin{1};
            end
            if isempty(dataset_name)
                dataset_name = obj.PrimaryDatasetUid;
            end
            if ~obj.LinkedDatasetChooserList.isKey(dataset_name)
                reporting.Error('MimLinkedDatasetChooser:DatasetNotFound', 'No linked dataset was found with this name. Did you add the dataset with LinkDataset()?');
            end
            linked_dataset_chooser = obj.LinkedDatasetChooserList(dataset_name);
            dataset_results = linked_dataset_chooser.PrimaryDatasetResults;
        end

        function is_linked_dataset = IsLinkedDataset(obj, linked_name_or_uid, reporting)
            % Returns true if another dataset has been linked to this one, using
            % the name or uid specified

            is_linked_dataset = obj.LinkedDatasetChooserList.isKey(linked_name_or_uid);
        end

        function ClearMemoryCacheInAllLinkedDatasets(obj)
            % Clears the temporary memory cache of this and all linked
            % datasets
            for linker = obj.LinkedDatasetChooserList.values
                if linker{1} == obj
                    obj.DatasetCache.ClearTemporaryMemoryCache;
                else
                    linker{1}.ClearMemoryCacheInAllLinkedDatasets;
                end
            end
        end

        function NotifyPreviewImageChanged(obj, plugin_name)
            notify(obj,'PreviewImageChanged', CoreEventData(plugin_name));
        end

        function primary_uid = GetPrimaryDatasetUid(obj)
            primary_uid = obj.PrimaryDatasetUid;
        end

    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function save_txt(key,outfilename)
% Saves the Key to a text file.
% Inputs:
%   key: An object of type Key.
%   outfilename: The name for the output text file.

assert(nargin==2)
assert(isa(key,'Key'))
assert(isa(outfilename,'char'))
assert(key.validate())

outfile = fopen(outfilename,'w');

[x,y] = ind2sub(size(key.tar),find(key.tar));
numentries = length(x);
for ii=1:numentries
    modelname = key.modelset{x(ii)};
    segname = key.segset{y(ii)};
    fprintf(outfile,'%s %s target\n',modelname,segname);
end

[x,y] = ind2sub(size(key.non),find(key.non));
numentries = length(x);
for ii=1:numentries
    modelname = key.modelset{x(ii)};
    segname = key.segset{y(ii)};
    fprintf(outfile,'%s %s nontarget\n',modelname,segname);
end

fclose(outfile);

author: nesl


text:
%SerialLink.fkine Forward kinematics
%
% T = R.fkine(Q, OPTIONS) is the pose of the robot end-effector as an SE3
% object for the joint configuration Q (1xN).
%
% If Q is a matrix (KxN) the rows are interpreted as the generalized joint
% coordinates for a sequence of points along a trajectory.  Q(i,j) is the
% j'th joint parameter for the i'th trajectory point.  In this case T is a
% an array of SE3 objects (K) where the subscript is the index along the path.
%
% [T,ALL] = R.fkine(Q) as above but ALL (N) is a vector of SE3 objects describing
% the pose of the link frames 1 to N.
%
% Options::
%  'deg'    Assume that revolute joint coordinates are in degrees not
%           radians
%
% Note::
% - The robot's base or tool transform, if present, are incorporated into the
%   result.
% - Joint offsets, if defined, are added to Q before the forward kinematics are
%   computed.
% - If the result is symbolic then each element is simplified.
%
% See also SerialLink.ikine, SerialLink.ikine6s.

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function [t,allt] = fkine(robot, q, varargin)

%
% evaluate fkine for each point on a trajectory of
% theta_i or q_i data
%

n = robot.n;

opt.deg = false;

opt = tb_optparse(opt, varargin);

if opt.deg
    % in degrees mode, scale the columns corresponding to revolute axes
    q = robot.todegrees(q);
end

if nargout > 1
    allt(n) = SE3;
end

L = robot.links;
if numel(q) == n
    % single configuration
    t = SE3(robot.base);
    for i=1:n
        t = t * L(i).A(q(i));

        if nargout > 1
            allt(i) = t; % intermediate transformations
        end
    end
    t = t * SE3(robot.tool);
else
    if numcols(q) ~= n
        error('q must have %d columns', n)
    end
    t(numrows(q)) = SE3;  % preallocate storage
    for k=1:numrows(q)		% for each trajectory point
        qk = q(k,:);
        tt = SE3(robot.base);
        for i=1:n
            tt = tt * L(i).A(qk(i));
        end
        t(k) = tt * SE3(robot.tool);
    end
end

if issym(t)
    t = simplify(t);
else
    t = trnorm(t);
end

%robot.T = t;
%robot.notify('Moved');

author: petercorke

repo: robotics-toolbox-matlab


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 ntuple = ntuples(varargin)
% Computes all possible ntupples.

x = varargin;

ntuple = x{1};
for i = 2:length(x)
    num_col  = size(ntuple,2);
    num_item = length(x{i});
    ntuple   = repcel(ntuple, 1, num_item);
    newline  = repmat(x{i},1,num_col);
    ntuple   = [ntuple;newline];
end

repo: OpenTLD


text:
function test_suite = test_createBasisTransform3d
%TEST_CREATEBASISTRANSFORM3D  Test case for the file createBasisTransform3d
%
%   Test case for the file createBasisTransform3d

%   Example
%   test_createBasisTransform3d
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2011-10-13,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function test_Translate(testCase) %#ok<*DEFNU>
% Basic test to check the function runs

p1 = [3 4 5];
p2 = [10 20 30];
basis1 = [p1 1 0 0 0 1 0];
basis2 = [p2 1 0 0 0 1 0];

dp = p1-p2;
exp = [eye(3) dp' ; 0 0 0 1];

trans = createBasisTransform3d(basis1, basis2);
testCase.assertEqual(exp, trans, 'AbsTol', .01);

function test_TransformPointTranslation(testCase) %#ok<*DEFNU>
% Basic test to check the function runs

% two bases with different origins and same directions
basis1 = [0 0 0     1 0 0  0 1 0];
basis2 = [10 20 30  1 0 0  0 1 0];

trans = createBasisTransform3d(basis1, basis2);

% pt1 = [0 0 0];
% pt1T = transformPoint3d(pt1, trans);
% exp1 = [0 0 0];
% testCase.assertEqual(exp1, pt1T, 'AbsTol', .01);

pt2 = [10 20 30];
pt2T = transformPoint3d(pt2, trans);
exp2 = [0 0 0];
testCase.assertEqual(exp2, pt2T, 'AbsTol', .01);

function test_Rotate(testCase) %#ok<*DEFNU>
% Basic test to check the function runs

p1 = [3 4 5];
basis1 = [p1  1 0 0   0 1 0];
basis2 = [p1  0 1 0  -1 0 0];

trans = createBasisTransform3d(basis1, basis2);

% restrict the test to the linear part
exp = [0 1 0; -1 0 0; 0 0 1];
testCase.assertEqual(exp, trans(1:3, 1:3), 'AbsTol', .01);

function test_TransformedPoints(testCase) %#ok<*DEFNU>
% test a combination of translation and axis permutation

% origin basis
basis1 = [10 10 10   1 0 0   0 1 0];
% translation and permutation
basis2 = [0 0 0   0 1 0   0 0 1];

trans = createBasisTransform3d(basis1, basis2);

pt1 = [0 0 0];
pt1T = transformPoint3d(pt1, trans);

exp1 = [10 10 10];
testCase.assertEqual(exp1, pt1T, 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [Nm] = cal2Nm(cal)
% Convert energy or work from calories to newtons-meters.
% Note: these calories are different from the capitalized Calories found on
% American cereal boxes.  American consumers are to healthy to eat a
% 250,000 calorie candy bar.
% Chad A. Greene 2012
Nm = cal*4.1868;

author: Sable

repo: mcbench-benchmarks


text:
function [Z]=rgt(X,Y)
%
% Replicating greater than
%
% Does element by element operations on X and Y where non-same sized
% dimensions are implicity wrapped round to match the size of the larger
% to give a result matrix Z with size max(size(X),size(Y));
%
% In this case returns logical array with true where X>Y and 0 otherwise.
%
% N.B. for complex inputs this compares the *norms* of the values, *not*
%      the real parts as GE does!
% See also repops, gt
%
% Copyright 2006- by Jason D.R. Farquhar (jdrf@zepler.org)
% Inspired by code by Douglas M. Schwarz & Aki Vehtari.
Z=repop(X,'>',Y);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef IMODE < ALGORITHM
% <single> <real/integer> <large/none> <constrained/none>
% Improved multi-operator differential evolution
% minN  ---   4 --- Minimum population size
% aRate --- 2.6 --- Ratio of archive size to population size

%------------------------------- Reference --------------------------------
% K. M. Sallam, S. M. Elsayed, R. K. Chakrabortty, and M. J. Ryan, Improved
% multi-operator differential evolution algorithm for solving unconstrained
% problems, Proceedings of the IEEE Congress on Evolutionary Computation,
% 2020.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            [minN,aRate] = Algorithm.ParameterSet(4,2.6);

            %% Generate random population
            Population = Problem.Initialization();
            Archive    = [];
            MCR = zeros(20*Problem.D,1) + 0.2;
            MF  = zeros(20*Problem.D,1) + 0.2;
            k   = 1;
            MOP = ones(1,3)/3;

            %% Optimization
            while Algorithm.NotTerminated(Population)
                % Reduce the population size
                N          = ceil((minN-Problem.N)*Problem.FE/Problem.maxFE) + Problem.N;
                [~,rank]   = sort(FitnessSingle(Population));
                Population = Population(rank(1:N));
                Archive    = Archive(randperm(end,min(end,ceil(aRate*N))));
                % Generate parents, CR, F, and operator for each offspring
                Xp1 = Population(ceil(rand(1,N).*max(1,0.25*N))).decs;
                Xp2 = Population(ceil(rand(1,N).*max(2,0.5*N))).decs;
                Xr1 = Population(randi(end,1,N)).decs;
                Xr3 = Population(randi(end,1,N)).decs;
                P   = [Population,Archive];
                Xr2 = P(randi(end,1,N)).decs;
                CR  = randn(N,1).*sqrt(0.1) + MCR(randi(end,N,1));
                CR  = sort(CR);
                CR  = repmat(max(0,min(1,CR)),1,Problem.D);
                F   = min(1,trnd(1,N,1).*sqrt(0.1) + MF(randi(end,N,1)));
                while any(F<=0)
                    F(F<=0) = min(1,trnd(1,sum(F<=0),1).*sqrt(0.1) + MF(randi(end,sum(F<=0),1)));
                end
                F  = repmat(F,1,Problem.D);
                OP = arrayfun(@(S)find(rand<=cumsum(MOP),1),1:N);
                OP = arrayfun(@(S)find(OP==S),1:length(MOP),'UniformOutput',false);
                % Generate offspring
                PopDec = Population.decs;
                OffDec = PopDec;
                OffDec(OP{1},:) = PopDec(OP{1},:) + F(OP{1},:).*(Xp1(OP{1},:)-PopDec(OP{1},:)+Xr1(OP{1},:)-Xr2(OP{1},:));
                OffDec(OP{2},:) = PopDec(OP{2},:) + F(OP{2},:).*(Xp1(OP{2},:)-PopDec(OP{2},:)+Xr1(OP{2},:)-Xr3(OP{2},:));
                OffDec(OP{3},:) = F(OP{3},:).*(Xr1(OP{3},:)+Xp2(OP{3},:)-Xr3(OP{3},:));
                if rand < 0.4
                    Site = rand(size(CR)) > CR;
                    OffDec(Site) = PopDec(Site);
                else
                    p1 = randi(Problem.D,N,1);
                    p2 = arrayfun(@(S)find([rand(1,Problem.D),2]>CR(S,1),1),1:N);
                    for i = 1 : N
                        Site = [1:p1(i)-1,p1(i)+p2(i):Problem.D];
                        OffDec(i,Site) = PopDec(i,Site);
                    end
                end
                Offspring = Problem.Evaluation(OffDec);
                % Update the population and archive
                delta   = FitnessSingle(Population) - FitnessSingle(Offspring);
                replace = delta > 0;
                Archive = [Archive,Population(replace)];
                Archive = Archive(randperm(end,min(end,ceil(aRate*N))));
                Population(replace) = Offspring(replace);
                % Update CR, F, and probabilities of operators
                if any(replace)
                    w      = delta(replace)./sum(delta(replace));
                    MCR(k) = (w'*CR(replace,1).^2)./(w'*CR(replace,1));
                    MF(k)  = (w'*F(replace,1).^2)./(w'*F(replace,1));
                    k      = mod(k,length(MCR)) + 1;
                else
                    MCR(k) = 0.5;
                    MF(k)  = 0.5;
                end
                delta = max(0,delta./abs(FitnessSingle(Population)));
                if any(cellfun(@isempty,OP))
                	MOP = ones(1,3)/3;
                else
                    MOP = cellfun(@(S)mean(delta(S)),OP);
                    MOP = max(0.1,min(0.9,MOP./sum(MOP)));
                end
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function f = p27_f ( n, x )

%*****************************************************************************80
%
%% P27_F evaluates the objective function for problem 27.
%
%  Discussion:
%
%    F can be regarded as a function of R = SQRT ( X(1)^2 + X(2)^2 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Zbigniew Michalewicz,
%    Genetic Algorithms + Data Structures = Evolution Programs,
%    Third Edition,
%    Springer Verlag, 1996,
%    ISBN: 3-540-60676-9,
%    LC: QA76.618.M53.
%
%  Parameters:
%
%    Input, integer N, the number of variables.
%
%    Input, real X(N), the argument of the objective function.
%
%    Output, real F, the value of the objective function.
%
  r = sqrt ( x(1)^2 + x(2)^2 );

  a = ( 1.0 + 0.001 * r^2 )^( -2 );

  b = ( sin ( r ) )^2 - 0.5;

  f = 0.5 + a * b;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [success,msg] = ConvertCtrax2JAABA(varargin)

success = false;
msg = {};

[inmoviefile,intrxfile,annfile,...
  expdir,moviefilestr,trxfilestr,perframedirstr,...
  arenatype,arenacenterx,arenacentery,...
  arenaradius,arenawidth,arenaheight,roi2,...
  pxpermm,fps,overridefps,overridearena,...
  dosoftlink,...
  dofliplr,doflipud,dotransposeimage,...
  inperframedir,frameinterval] = myparse(varargin,...
  'inmoviefile','','intrxfile','','annfile','',...
  'expdir','','moviefilestr','movie.ufmf','trxfilestr','trx.mat','perframedirstr','perframe',...
  'arenatype','None','arenacenterx',0,'arenacentery',0,...
  'arenaradius',123,'arenawidth',123,'arenaheight',123,'roi2',[],...
  'pxpermm',1,'fps',30,...
  'overridefps',false,'overridearena',false,...
  'dosoftlink',false,...
  'fliplr',false,...
  'flipud',false,...
  'dotransposeimage',false,...
  'inperframedir','',...
  'frameinterval',[1,inf]);

% note that if flipping, transposing trajectories and input arena is given,
% the input arena is NOT flipped

wingunits = struct(...
  'nwingsdetected',parseunits('unit'),...
  'wing_areal',parseunits('px^2'),...
  'wing_arear',parseunits('px^2'),...
  'wing_trough_angle',parseunits('rad'));

%% check that required inputs are given
ismoviefile = ~isempty(inmoviefile);
if isempty(intrxfile),
  msg = 'Input trx mat file is empty';
  return;
end
% if isempty(annfile),
%   msg = 'Input ann file is empty';
%   return;
% end

% output file locations
if ismoviefile,
  moviefile = fullfile(expdir,moviefilestr);
else
  moviefile = '';
end
trxfile = fullfile(expdir,trxfilestr);
perframedir = fullfile(expdir,perframedirstr);

%% load in the trx
try
  [trx,~,success1,timestamps] = load_tracks(intrxfile,moviefile,...
    'dosave',false,'annname',annfile,'verbose',false);
  if ~success1,
    msg = sprintf('Could not load trx from file %s',intrxfile);
    return;
  end
catch ME,
  msg = sprintf('Could not load trx from file %s: %s',intrxfile,getReport(ME));
  return;
end

if all(isnan(timestamps)),
  timestamps = [];
end

iswings = isfield(trx,'xwingl');

if iswings,
  if isempty(inperframedir),
    msg = 'Per-frame directory containing wing data is not set';
    return;
  end
  fns = fieldnames(wingunits);
  for i = 1:numel(fns),
    fn = fns{i};
    filename = fullfile(inperframedir,[fn,'.mat']);
    if ~exist(filename,'file'),
      msg = sprintf('Wing feature file %s does not exist.',filename);
      return;
    end
  end
end

if dofliplr || doflipud,

  if ~ismoviefile,
    msg = 'Cannot flip without input movie file';
    return;
  end

  try
    % get image size
    [readframe,~,fid] = get_readframe_fcn(inmoviefile);
    im = readframe(1);
    [imheight,imwidth,~] = size(im);
    if fid > 1,
      fclose(fid);
    end
  catch ME,
    msg = getReport(ME);
    return;
  end
end

% flip if necessary
if dofliplr,
  for i = 1:numel(trx),
    trx(i).x = imwidth-trx(i).x+1;
    trx(i).theta = modrange(pi - trx(i).theta,-pi,pi);
    if iswings,
      trx(i).xwingl = imwidth-trx(i).xwingl+1;
      trx(i).xwingr = imwidth-trx(i).xwingr+1;
    end
  end
  if isfield(trx,'arena') && isfield(trx(1).arena,'x'),
    for i = 1:numel(trx),
      trx(i).arena.x = imwidth-trx(i).arena.x+1;
    end
  end
  msg{end+1} = 'Flipped the trajectories left-right.';
end
if doflipud,
  for i = 1:numel(trx),
    trx(i).y = imheight-trx(i).y+1;
    trx(i).theta = modrange(-trx(i).theta,-pi,pi);
    if iswings,
      trx(i).ywingl = imheight-trx(i).xwingl+1;
      trx(i).ywingr = imheight-trx(i).xwingr+1;
    end
  end
  if isfield(trx,'arena') && isfield(trx(1).arena,'y'),
    for i = 1:numel(trx),
      trx(i).arena.y = imheight-trx(i).arena.y+1;
    end
  end
  msg{end+1} = 'Flipped the trajectories up-down.';
end
if dotransposeimage,
  for i = 1:numel(trx),
    tmp = trx(i).x;
    trx(i).x = trx(i).y;
    trx(i).y = tmp;
    c = cos(trx(i).theta);
    s = sin(trx(i).theta);
    trx(i).theta = atan2(c,s);
    if iswings,
      tmp = trx(i).xwingl;
      trx(i).xwingl = trx(i).ywingl;
      trx(i).ywingl = tmp;
      tmp = trx(i).xwingr;
      trx(i).xwingr = trx(i).ywingr;
      trx(i).ywingr = tmp;
    end
  end
  if isfield(trx,'arena') && isfield(trx(1).arena,'x') && isfield(trx(1).arena,'y'),
    for i = 1:numel(trx),
      tmp = trx(i).arena.x;
      trx(i).arena.x = trx(i).arena.y;
      trx(i).arena.y = tmp;
    end
  end
  msg{end+1} = 'Switched x and y in trajectories.';
end

% Clip the tracks if necessary
outofrange = [];
clipranges = zeros(length(trx),2);
if (frameinterval(1)>1) || ~isinf(frameinterval(end)),
  fields2clip = {'x','y','a','b','theta','timestamps',...
                 'wing_anglel','wing_angler',...
                 'xwingl','ywingl','xwingr','ywingr'};
  for fly = 1:length(trx),
    if trx(fly).endframe<frameinterval(1) || ...
        trx(fly).firstframe>frameinterval(2),
      outofrange(end+1) = fly;  %#ok<AGROW>
      continue;
    end
    tstart = max(1,frameinterval(1)-trx(fly).firstframe+1);
    tend = min(trx(fly).nframes,frameinterval(2)-trx(fly).endframe+trx(fly).nframes);
    clipranges(fly,1) = tstart;
    clipranges(fly,2) = tend;
    trx(fly).firstframe = max(frameinterval(1),trx(fly).firstframe);
    trx(fly).endframe = min(frameinterval(2),trx(fly).endframe);
    trx(fly).nframes = trx(fly).endframe-trx(fly).firstframe+1;
    for fieldndx = 1:numel(fields2clip)
      if isfield(trx(fly),fields2clip{fieldndx});
        trx(fly).(fields2clip{fieldndx}) = trx(fly).(fields2clip{fieldndx})(tstart:tend);
      end
    end
  end

end
trx(outofrange) = [];

% compute frame rate if not overriding and possible
if ~overridefps && ~isempty(timestamps),
  fps = 1/nanmedian(diff(timestamps));
  msg{end+1} = sprintf('Computed frame rate = %f fps from timestamps.',fps);
else
  msg{end+1} = sprintf('Using input frame rate %f fps to set timestamps',fps);
end

if ~overridearena && isfield(trx,'pxpermm'),
  tmp = nanmean([trx.pxpermm]);
  if ~isnan(tmp),
    pxpermm = tmp;
    msg{end+1} = sprintf('Read pxpermm = %f from trxfile.',pxpermm);
  end
else
  msg{end+1} = sprintf('Using input pixels-per-mm scaling = %f.',pxpermm);
end

if ~overridearena && isfield(trx,'arena') && all(isfield(trx(1).arena,{'x','y','r'})) && ...
  ~isempty(trx(1).arena.r) && ~isempty(trx(1).arena.x) && ~isempty(trx(1).arena.y),
  arenaradius = trx(1).arena.r;
  arenacenterx = trx(1).arena.x;
  arenacentery = trx(1).arena.y;
  arenatype = 'Circle';
  msg{end+1} = sprintf('Read circle arena parameters from trxfile. Center = (%.1f, %.1f) px, radius = %.1f px.',arenacenterx,arenacentery,arenaradius);
elseif strcmpi(arenatype,'None'),
  arenacenterx = 0;
  arenacentery = 0;
end

if ~isempty(roi2),
  roi2(:,1) = (roi2(:,1)-arenacenterx) / pxpermm;
  roi2(:,2) = (roi2(:,2)-arenacentery) / pxpermm;
  roi2(:,3) = roi2(:,3) / pxpermm;
  roi2(:,4) = roi2(:,4) / pxpermm;
end

% apply to trajectories
for fly = 1:length(trx),

  % scale and translate position
  trx(fly).x_mm = (trx(fly).x - arenacenterx) / pxpermm;
  trx(fly).y_mm = (trx(fly).y - arenacentery) / pxpermm;
  trx(fly).a_mm = trx(fly).a / pxpermm;
  trx(fly).b_mm = trx(fly).b / pxpermm;
  trx(fly).theta_mm = trx(fly).theta;

  % add in dt
  if overridefps || isempty(timestamps),
    trx(fly).dt = repmat(1/fps,[1,trx(fly).nframes-1]);
  else
    trx(fly).dt = diff(timestamps(trx(fly).firstframe:trx(fly).endframe));
    trx(fly).dt(isnan(trx(fly).dt)) = 1/fps;
  end

  trx(fly).fps = 1/median(trx(fly).dt);
  trx(fly).pxpermm = pxpermm;

end

msg{end+1} = sprintf('Read trx for %d flies, frame range [%d,%d]',numel(trx),min([trx.firstframe]),max([trx.endframe]));

if isempty(timestamps) || overridefps,
  timestamps = (0:max([trx.endframe]))*1/fps; %#ok<NASGU>
end

% set landmark parameters
arenaradius_mm = arenaradius / pxpermm;
arenawidth_mm = arenawidth / pxpermm;
arenaheight_mm = arenaheight / pxpermm;
arenacenterx_mm = 0;
arenacentery_mm = 0;
trx = SetLandmarkParameters(trx,arenatype,arenacenterx_mm,arenacentery_mm,...
  arenaradius_mm,arenawidth_mm,arenaheight_mm,roi2);
trx = SetLandmarkParameters_px(trx,arenatype,arenacenterx,arenacentery,arenaradius);

if strcmpi(arenatype,'Circle')
  msg{end+1} = sprintf('Set circular arena parameters. Center = (0,0) mm, radius = %.1f mm.',arenaradius_mm);
elseif strcmpi(arenatype,'Rectangle'),
  msg{end+1} = sprintf('Set rectangular arena parameters. Center = (0,0) mm, width= %.1f mm, height = %.1f mm',arenawidth_mm,arenaheight_mm);
end
msg{end+1} = sprintf('x_mm ranges within [%.1f,%.1f], y_mm ranges within [%.1f,%.1f]',...
  min([trx.x_mm]),max([trx.x_mm]),min([trx.y_mm]),max([trx.y_mm]));

% create the experiment directory
if ~exist(expdir,'dir'),
  [success1,msg1] = mkdir(expdir);
  if ~success1,
    msg = msg1;
    return;
  end
  msg{end+1} = sprintf('Created experiment directory %s.',expdir);
end

% save the trx file
try
  save(trxfile,'trx','timestamps');
catch ME,
  msg = sprintf('Could not save to file %s: %s',trxfile,getReport(ME));
  return;
end
if ~exist(trxfile,'file'),
  msg = sprintf('Failed to save trx to file %s',trxfile);
  return;
end
msg{end+1} = sprintf('Saved trx to file %s.',trxfile);

% copy/soft-link movie
if ismoviefile,
  if strcmp(fullfile(inmoviefile),fullfile(moviefile)),
    msg{end+1} = 'Input and out movie files are the same, not copying/linking.';
  else

    if dosoftlink,
      if exist(moviefile,'file'),
        delete(moviefile);
      end
      if isunix,
        cmd = sprintf('ln -s %s %s',inmoviefile,moviefile);
        unix(cmd);
        % test to make sure it worked
        [status,result] = unix(sprintf('readlink %s',moviefile));
        result = strtrim(result);
        if status ~= 0 || ~strcmp(result,inmoviefile),
          res = questdlg(sprintf('Failed to make soft link. Copy %s to %s instead?',inmoviefile,moviefile));
          if ~strcmpi(res,'Yes'),
            msg = sprintf('Failed to make soft link from %s to %s.',inmoviefile,moviefile);
            return;
          end
          dosoftlink = false;
        end
      elseif ispc,
        if exist([moviefile,'.lnk'],'file'),
          delete([moviefile,'.lnk']);
        end
        cmd = sprintf('mkshortcut.vbs /target:"%s" /shortcut:"%s"',inmoviefile,moviefile);
        fprintf('Making a Windows shortcut file at "%s" with target "%s"\n',inmoviefile,moviefile);
        system(cmd);
        % test to make sure that worked
        [equalmoviefile,didfind] = GetPCShortcutFileActualPath(moviefile);
        if ~didfind || ~strcmp(equalmoviefile,inmoviefile),
          res = questdlg(sprintf('Failed to make shortcut. Copy %s to %s instead?',inmoviefile,moviefile));
          if ~strcmpi(res,'Yes'),
            msg = sprintf('Failed to make shortcut from %s to %s.',inmoviefile,moviefile);
            return;
          end
          dosoftlink = false;
        end
      else
        res = questdlg(sprintf('Unknown OS, cannot soft-link movie file %s. Copy instead?',inmoviefile));
        if ~strcmpi(res,'Yes'),
          msg = sprintf('Failed to make softlink from %s to %s.',inmoviefile,moviefile);
          return;
        end
        dosoftlink = false;
      end
      if dosoftlink,
        msg{end+1} = sprintf('Made a link to movie file %s at %s',inmoviefile,moviefile);
      end
    end

    if ~dosoftlink,
      if ispc,
        if exist([moviefile,'.lnk'],'file'),
          delete([moviefile,'.lnk']);
        end
      end
      if exist(moviefile,'file'),
        delete(moviefile);
      end
      [success1,msg1] = copyfile(inmoviefile,moviefile);
      if ~success1,
        msg = msg1;
        success = false;
        return;
      end
      msg{end+1} = sprintf('Copied movie file %s to %s',inmoviefile,moviefile);
    end

  end

end

% make per-frame directory
if ~exist(perframedir,'dir'),
  [success1,msg1] = mkdir(perframedir);
  if ~success1,
    msg = msg1;
    return;
  end
end

% copy over wing features
if iswings,
  fns = fieldnames(wingunits);
  doclip = false;
  docopy = false;
  fn = fns{1};
  filenamein = fullfile(inperframedir,[fn,'.mat']);
  filenameout = fullfile(perframedir,[fn,'.mat']);
  if (frameinterval(1)>1) || ~isinf(frameinterval(end)),
    if strcmp(filenamein,filenameout),
      qstr = sprintf('Input and output file locations for %s/%s.mat (and possibly others) are the same. The original file will be overwritten with newer version with fewer frames. Overwrite?',...
        perframedir,fn);
      res = questdlg(qstr,'Overwrite perframe file','Yes','No','No');
      if strcmpi(res,'Yes'),doclip = true; end
    else
      doclip = true;
    end
  else
    if strcmp(filenamein,filenameout),
      msg{end+1} = fprintf('%s/%s.mat input and output file locations are the same, not copying.',perframedir,fn);
    else
      docopy = true;
    end
  end
  for i = 1:numel(fns),
    fn = fns{i};
    filenamein = fullfile(inperframedir,[fn,'.mat']);
    filenameout = fullfile(perframedir,[fn,'.mat']);

    if docopy
      if exist(filenameout,'file'),
        delete(filenameout);
      end
      [success1,msg1] = copyfile(filenamein,filenameout);
      if ~success1,
        msg = msg1;
        return;
      end
    end
    if doclip
      Q = load(filenamein);
      for fly = 1:numel(Q.data)
        if ismember(fly,outofrange),
          continue;
        end
        Q.data{fly} = Q.data{fly}(clipranges(fly,1):clipranges(fly,2));
      end
      Q.data(outofrange) = [];

      save(filenameout,'-struct','Q');

    end
  end
end

success = true;

author: kristinbranson

repo: JAABA


text:
function mov = rmPlotGUI_predictionMovie(M, aviFile, fps);
% Animate a movie of the stimulus sweeping across the pRF, and the
% predicted response, for the selected voxel in the rm Plot GUI.
%
%  mov = rmPlotGUI_predictionMovie([M], [aviFile], [fps=6]);
%
% INPUTS:
%	M: rmPlotGUI struct. [default: get from cur fig]
%
%	aviFile: optional .avi file to save the movie. [if not path provided,
%	won't save an .avi file.]
%
%	fps: frames per second for the avi movie (default: 6)
%
% OUTPUTS:
%	mov: movie structure. See GETFRAME.
%
% ras, 05/2009.
if notDefined('M'),		M = get(gcf, 'UserData');			end
if notDefined('aviFile'),	aviFile = '';					end
if notDefined('fps'),		fps = 6;						end

%% get the prediction for this voxel
v = get(M.ui.voxel.sliderHandle, 'Value');
coords = M.coords(:,v);
if isequal(M.roi.viewType, 'Gray')  % convert coords into an index
    coords = M.coords(v);
end
[pred RF rfParams variance_explained] = rmPlotGUI_makePrediction(M, coords);

%% create the animation, get the frames
nFrames = size(M.tSeries, 1);

for n = 1:nFrames
	%% plot the time series up to this frame
	axes(M.ui.tsAxes); cla; hold on;

    % We can either plot a static time series and growing prediction...
    % 	hTSeries = plot(M.x, M.tSeries(:,v), 'k--', 'LineWidth', 1.5);
    % 	hFit = plot(M.x(1:n), pred(1:n,1), 'b', 'LineWidth', 2);

    % Or a static prediction and a growing time series (my preference - jw)
    hTSeries = plot(M.x, M.tSeries(:,v), 'k--', 'LineWidth', 1.5);
	hFit = plot(M.x(1:n), pred(1:n,1), 'b', 'LineWidth', 2);

    % Don't plot the residual - too much clutter (jw)
    % 	hResidual = plot(M.x(1:n), M.tSeries(1:n,v)-pred(:,1), 'r:', 'LineWidth', 1);

	allPlotted = [M.tSeries(:,v); pred(:,1); M.tSeries(:,v)-pred(:,1)];
	axis([min(M.x) max(M.x) min(allPlotted) max(allPlotted)]);
	h = axis;
% 	for n=1:numel(M.sepx),
% 		plot([1 1].*M.sepx(n), [h(3) h(4)], 'k:', 'LineWidth', 2);
% 	end;
	xlabel('Time (sec)');
	ylabel('BOLD signal change (%)');

	%% show the stimulus position over the pRF for this frame
	TR = M.params.stim(1).framePeriod;  % TODO: deal w/ diff't TRs in diff't scans
	t = get(M.ui.time.sliderHandle, 'Value') * TR;

	% get the stimulus image for this frame
	[stimImage RFresampled] = getCurStimImage(M, n, RF);

	% overlay and display
	axes(M.ui.rfAxes); cla
% 	RF_img(:,:,1) = RF;
% 	RF_img(:,:,2) = 1-RF;
% 	RF_img(:,:,3) = stimImage;
	RF_img(:,:,1) = stimImage;
	RF_img(:,:,2) = RFresampled;
	RF_img(:,:,3) = RFresampled;

    [x,y] = prfSamplingGrid(M.params);
    x = x(:); y = y(:);
    imagesc(x, -y, RF_img); hold on;
    plot([min(x) max(x)], [0 0], 'k-');
    plot( [0 0], [min(y) max(y)], 'k-');

    axis image  xy;
    ylabel('y (deg)');
    xlabel('x (deg)');

	%% grab the current movie frame
	mov(n) = getframe(gcf);
end

%% save the .avi file, if a path is provided.
if ~isempty(aviFile)
	quality = 100;
	description = aviFile;
	keyframe = fps;
	if ispc
		compression = 'Indeo5';
	else
		compression = 'None';
	end

	ensureDirExists( fileparts(aviFile) );

	movie2avi(mov, aviFile, 'FPS', fps, 'Compression', compression, ...
		  'Quality', quality, 'Videoname', description, 'Keyframe', keyframe);
	fprintf('[%s]: Saved movie as %s.\n', mfilename, aviFile);

end

return
%--------------------------------------

%--------------------------------------
function [stimImage RF] = getCurStimImage(M, f, RFvals)
% Get a stimulus image matching the sampling positions as the RF.
% Also returns the RF resampled into a square grid.
x = prfSamplingGrid(M.params);

% account for the different stimuli that are shown next to each other
% f originally refers to the frame in the combined time series across scans:
% we want to break this down into scan n, frame f within that scan.
n = 1;
nStimScans = numel(M.params.stim);
while n <= nStimScans,
    tmp = f + M.params.stim(n).prescanDuration;
    if tmp > size(M.params.stim(n).images_org,2),
        f = tmp - size(M.params.stim(n).images_org,2);
        n = n + 1;
    else
        f = tmp;
        break;
    end
end

% stim image
stimImage     = NaN(size(x));
stimImage(M.params.stim(1).instimwindow) = M.params.stim(n).images_org(:,f);
stimImage     = reshape(stimImage, size(x));

% RF
RF     = NaN(size(x));
RF(M.params.stim(1).instimwindow) = normalize(RFvals, 0, 1);
RF     = reshape(RF, size(x));

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text: %MINUS Dataset overload

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
classdef iTreeNode < handle
    % iTreeNode A class to represent an internal (non-leaf) node of an
    % isolation tree
    %   An isolation tree is a type of binary tree.

    properties
        SplitAttribute
        SplitValue
    end

    properties(SetAccess = private)
      Left = iTreeNode.empty;
      Right = iTreeNode.empty;
   end

    methods
        function node = iTreeNode(splitatt, splitvalue)
            % Construct an iTreeNode object.
            if nargin > 0
                node.SplitAttribute = splitatt;
                node.SplitValue = splitvalue;
            end
        end

        function insertLeft(parentNode, newNode)
            parentNode.Left = newNode;
        end

        function insertRight(parentNode, newNode)
            parentNode.Right = newNode;
        end

        function leftorright = compare(node, x)
            leftorright = x(node.SplitAttribute) >= node.SplitValue;
        end
    end

end

author: dsmi-lab-ntust

repo: AnomalyDetectionToolbox


text:
%% Scene Text Detection using CNN
%
% References:
%
% * Minghui Liao, Baoguang Shi, Xiang Bai, Xinggang Wang, Wenyu Liu.
%   "TextBoxes: A Fast Text Detector with a Single Deep Neural Network",
%   AAAI 2017. <https://arxiv.org/abs/1611.06779>
%
% Sources:
%
% * <https://github.com/opencv/opencv_contrib/blob/3.4.1/modules/text/samples/textbox_demo.cpp>
% * <https://github.com/opencv/opencv_contrib/blob/3.4.1/modules/text/samples/deeptextdetection.py>
%

%%
% load model
[net, blobOpts] = TextBoxes();

%%
% detect text
im = which('handicapSign.jpg');
assert(~isempty(im), 'specify input image');
img = cv.imread(im, 'Color',true);
[bboxes, confs] = detectText(img, net, blobOpts);
fprintf('%d text detections\n', size(bboxes,1));

%%
% post-processing
if true
    % non-maximum suppression
    ind = cv.Net.NMSBoxes(bboxes, confs, 0.3, 0.4);
    ind = ind + 1;  % zero to one-based indices
    bboxes = bboxes(ind,:);
    confs = confs(ind);
else
    % confidence thresholding
    idx = confs > 0.6;
    bboxes = bboxes(idx,:);
    confs = confs(idx);
end
fprintf('%d text detections after filtering\n', size(bboxes,1));

%%
% show results
for i=1:size(bboxes,1)
    img = insertAnnotation(img, bboxes(i,:), sprintf('%.2f', confs(i)), ...
        'Color',[255 0 255], 'TextColor',[255 255 255], 'Thickness',2, 'FontScale',0.3);
    fprintf('text box: [%3d %3d %3d %3d], confidence: %3.0f%%\n', ...
        bboxes(i,:), confs(i)*100);
end
imshow(img), title('TextBox Demo')

%%
% Helper functions

function [rects, confs] = detectText(img, net, blobOpts)
    %DETECTTEXT  Detect text on input image

    % detect text
    net.setInput(cv.Net.blobFromImages(img, blobOpts{:}));
    dets = net.forward();

    % SSD output is 1-by-1-by-ndetections-by-7
    % d = [img_id, class_id, confidence, left, top, right, bottom]
    dets = permute(dets, [3 4 2 1]);

    % filter
    dets = dets(dets(:,2) == 1,:);  % 0: background, 1: text

    % adjust relative coordinates to image size
    sz = [size(img,2) size(img,1)];
    dets(:,4:7) = bsxfun(@times, dets(:,4:7), [sz sz]);

    % output detections (clamp coords)
    rects = cv.Rect.from2points(dets(:,4:5), dets(:,6:7));
    rects = cv.Rect.intersect(rects, [0 0 sz]);
    rects = round(rects);
    confs = dets(:,3);
end

function img = insertAnnotation(img, rect, str, varargin)
    % See also: insertObjectAnnotation, insertShape, insertText
    p = inputParser();
    p.addParameter('Alpha', 0.6);
    p.addParameter('Thickness', 1);
    p.addParameter('Color', [255 255 0]);
    p.addParameter('TextColor', [0 0 0]);
    p.addParameter('FontFace', 'HersheySimplex');
    p.addParameter('FontScale', 0.4);
    p.addParameter('AntiAlias', true);
    p.addParameter('Shape', 'rectangle');
    p.parse(varargin{:});
    opts = p.Results;
    opts.Shape = validatestring(opts.Shape, {'rectangle','circle'});
    thick = 1;

    [sz,b] = cv.getTextSize(str, 'Thickness',thick, ...
        'FontFace',opts.FontFace, 'FontScale',opts.FontScale);
    txt_rect = [rect(1), rect(2)-sz(2)-b, sz(1), sz(2)+b];
    txt_orig = [rect(1), rect(2)-b];

    if opts.AntiAlias
        alias = {'LineType','AA'};
    else
        alias = {'LineType',8};
    end

    overlay = img;
    if strcmp(opts.Shape, 'rectangle')
        overlay = cv.rectangle(overlay, rect, ...
            'Color',opts.Color, 'Thickness',opts.Thickness, alias{:});
    else
        c = rect(1:2) + rect(3:4)/2;
        r = max(rect(3:4)/2);
        overlay = cv.circle(overlay, c, r, ...
            'Color',opts.Color, 'Thickness',opts.Thickness, alias{:});
    end
    overlay = cv.rectangle(overlay, txt_rect, ...
        'Color',opts.Color, 'Thickness','Filled', alias{:});
    if opts.Thickness > 1
        overlay = cv.rectangle(overlay, txt_rect, ...
            'Color',opts.Color, 'Thickness',opts.Thickness, alias{:});
    end
    overlay = cv.putText(overlay, str, txt_orig, ...
        'FontFace',opts.FontFace, 'FontScale',opts.FontScale, ...
        'Color',opts.TextColor, 'Thickness',thick, alias{:});

    img = cv.addWeighted(img,1-opts.Alpha, overlay,opts.Alpha, 0);
end

%%
% Pretrained models

function dname = get_dnn_dir(dname)
    %GET_DNN_DIR  Path to model files, and show where to get them if missing

    dname = fullfile(mexopencv.root(), 'test', 'dnn', dname);
    b = isdir(dname);
    if ~b
        % display help of calling function
        % (assumed to be a local function in current file)
        st = dbstack(1);
        help([mfilename() filemarker() st(1).name])
    end
    assert(b, 'Missing model: %s', dname);
end

function [net, blobOpts] = TextBoxes()
    %TEXTDETECTOR  Text detector model [Caffe]
    %
    % homepage = https://github.com/MhLiao/TextBoxes
    %
    % ## Model
    %
    % file = test/dnn/TextBoxes/textbox.prototxt
    % url  = https://github.com/opencv/opencv_contrib/raw/3.4.1/modules/text/samples/textbox.prototxt
    % hash = c294416fe6d156b9383342d62b9158ab707170c0
    %
    % ## Weights
    %
    % file = test/dnn/TextBoxes/TextBoxes_icdar13.caffemodel
    % url  = https://www.dropbox.com/s/g8pjzv2de9gty8g/TextBoxes_icdar13.caffemodel?dl=0
    % hash = e4b32aef3db3d66fec0630c96006e10999f785c4
    % size = 90.68 MB
    %

    dname = get_dnn_dir('TextBoxes');
    net = cv.Net('Caffe', ...
        fullfile(dname, 'textbox.prototxt'), ...
        fullfile(dname, 'TextBoxes_icdar13.caffemodel'));
    blobOpts = {'SwapRB',true, 'Crop',false, 'Size',[300 300], 'Mean',[104 117 123]};
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
% Copyright (C) 2008 Eric Chassande-Mottin, CNRS (France) <ecm@apc.univ-paris7.fr>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING. If not, see
% <https://www.gnu.org/licenses/>.

% -*- texinfo -*-
% @deftypefn  {Function File} {[@var{y}, @var{h}] =} resample (@var{x}, @var{p}, @var{q})
% @deftypefnx {Function File} {@var{y} =} resample (@var{x}, @var{p}, @var{q}, @var{h})
% Change the sample rate of @var{x} by a factor of @var{p}/@var{q}.  This is
% performed using a polyphase algorithm.  The impulse response @var{h} of
% the antialiasing filter is either specified or either designed with a
% Kaiser-windowed sinecard.
%
% Ref [1] J. G. Proakis and D. G. Manolakis,
% Digital Signal Processing: Principles, Algorithms, and Applications,
% 4th ed., Prentice Hall, 2007. Chap. 6
%
% Ref [2] A. V. Oppenheim, R. W. Schafer and J. R. Buck,
% Discrete-time signal processing, Signal processing series,
% Prentice-Hall, 1999
% @end deftypefn

function  [y, h] = resample( x, p, q, h )

if nargin<3
  error('at least 3 inputs are needed');
end

if any([p q]<=0) || any([p q]~=floor([p q]))
  error('resample.m: p and q must be positive integers');
end

% simplify decimation and interpolation factors

great_common_divisor=gcd(p,q);
if (great_common_divisor>1)
  p = double (p) / double (great_common_divisor);
  q = double (q) / double (great_common_divisor);
else
  p = double (p);
  q = double (q);
end

% filter design if required

if (nargin < 4)

  % properties of the antialiasing filter

  log10_rejection = -3.0;
  stopband_cutoff_f = 1 / (2 * max (p, q));
  roll_off_width = stopband_cutoff_f / 10.0;

  % determine filter length
  % use empirical formula from [2] Chap 7, Eq. (7.63) p 476

  rejection_dB = -20.0*log10_rejection;
  L = ceil((rejection_dB-8.0) / (28.714 * roll_off_width));

  % ideal sinc filter

  t=(-L:L)';

  % JMS: The matlab implementation uses a firls filter
  %ideal_filter=2*p*stopband_cutoff_f*sinc(2*stopband_cutoff_f*t);
  Lx = 2*10*max(p,q) + 1; % this is according to the matlab default
  ideal_filter=firls( Lx-1, [0 2*stopband_cutoff_f 2*stopband_cutoff_f 1], [1 1 0 0]);

  % determine parameter of Kaiser window
  % use empirical formula from [2] Chap 7, Eq. (7.62) p 474

  if ((rejection_dB>=21) && (rejection_dB<=50))
    beta = 0.5842 * (rejection_dB-21.0)^0.4 + 0.07886 * (rejection_dB-21.0);
  elseif (rejection_dB>50)
    beta = 0.1102 * (rejection_dB-8.7);
  else
    beta = 0.0;
  end

  % JMS: beta = 5 by default in the matlab version, this overrules the
  % above
  beta = 5;

  % apodize ideal filter response

  %h=kaiser(2*L+1,beta).*ideal_filter;
  h=kaiser(Lx,beta).*ideal_filter(:);

end

% check if input is a row vector
isrowvector=false;
if size(x,1)==1 && size(x,2)>1
  x=x(:);
  isrowvector=true;
end

% check if filter is a vector
if ~isvector(h)
  error('resample.m: the filter h should be a vector');
end

Lx = size(x,1);
Lh = length(h);
L = ( Lh - 1 )/2.0;
Ly = ceil(Lx*p/q);

% pre and postpad filter response

nz_pre = floor(q-mod(L,q));
hpad = [zeros(nz_pre,1);h];

offset = floor((L+nz_pre)/q);
nz_post = 0;
while ceil( ( (Lx-1)*p + nz_pre + Lh + nz_post )/q ) - offset < Ly
  nz_post = nz_post+1;
end
hpad = [hpad; zeros(nz_post,1)];

% filtering
if p==1 && q==1
  % no resampling needed.
  y = x;
else
  xfilt = upfirdn(x,hpad,p,q);
  y = xfilt(offset+1:offset+Ly,:);
end

if isrowvector
  y=y.';
end

%!test
%! N=512;
%! p=3; q=5;
%! r=p/q;
%! NN=ceil(r*N);
%! t=0:N-1;
%! tt=0:NN-1;
%! err=zeros(N/2,1);
%! for n = 0:N/2-1,
%!   phi0=2*pi*rand;
%!   f0=n/N;
%!   x=sin(2*pi*f0*t' + phi0);
%!   [y,h]=resample(x,p,q);
%!   xx=sin(2*pi*f0/r*tt' + phi0);
%!   t0=ceil((length(h)-1)/2/q);
%!   idx=t0+1:NN-t0;
%!   err(n+1)=max(abs(y(idx)-xx(idx)));
%! endfor;
%! rolloff=.1;
%! rejection=10^-3;
%! idx_inband=1:ceil((1-rolloff/2)*r*N/2)-1;
%! assert(max(err(idx_inband))<rejection);

%!test
%! N=512;
%! p=3; q=5;
%! r=p/q;
%! NN=ceil(r*N);
%! t=0:N-1;
%! tt=0:NN-1;
%! reject=zeros(N/2,1);
%! for n = 0:N/2-1,
%!   phi0=2*pi*rand;
%!   f0=n/N;
%!   x=sin(2*pi*f0*t' + phi0);
%!   [y,h]=resample(x,p,q);
%!   xx=sin(2*pi*f0/r*tt' + phi0);
%!   t0=ceil((length(h)-1)/2/q);
%!   idx=t0+1:NN-t0;
%!   reject(n+1)=max(abs(y(idx)));
%! endfor;
%! rolloff=.1;
%! rejection=10^-3;
%! idx_stopband=ceil((1+rolloff/2)*r*N/2)+1:N/2;
%! assert(max(reject(idx_stopband))<=rejection);

%!test
%! N=1024;
%! p=2; q=7;
%! r=p/q;
%! NN=ceil(r*N);
%! t=0:N-1;
%! tt=0:NN-1;
%! err=zeros(N/2,1);
%! for n = 0:N/2-1,
%!   phi0=2*pi*rand;
%!   f0=n/N;
%!   x=sin(2*pi*f0*t' + phi0);
%!   [y,h]=resample(x,p,q);
%!   xx=sin(2*pi*f0/r*tt' + phi0);
%!   t0=ceil((length(h)-1)/2/q);
%!   idx=t0+1:NN-t0;
%!   err(n+1)=max(abs(y(idx)-xx(idx)));
%! endfor;
%! rolloff=.1;
%! rejection=10^-3;
%! idx_inband=1:ceil((1-rolloff/2)*r*N/2)-1;
%! assert(max(err(idx_inband))<rejection);

%% Test integer-type arguments
%!test
%! N = 512;
%! f = 0.1;
%! x = sin (2*pi*f*[0:N-1]);
%! y1 = resample (x, 3, 2);
%! y2 = resample (x, uint8 (3), 2);
%! assert (y1, y2);


text:
function ADA = getada2(ADA, DAt,Aord, K) %#ok
% ADA = getada2(ADA, DAt,Aord, K)
%
% GETADA2  Compute ADA += DAt.q'*DAt.q
%   IMPORTANT: Updated ADA only on triu(ADA(Aord.qperm,Aord.qperm)).
%     Remaining entries are not affected.
%
% ******************** INTERNAL FUNCTION OF SEDUMI ********************
%
% See also sedumi, getada1, getada3

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% 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 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

%Indicate to the user Matlab cannot find the SeDuMi binaries
sedumi_binary_error();

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [J, grad] = lrCostFunction(theta, X, y, lambda)
%LRCOSTFUNCTION Compute cost and gradient for logistic regression with
%regularization
%   J = LRCOSTFUNCTION(theta, X, y, lambda) computes the cost of using
%   theta as the parameter for regularized logistic regression and the
%   gradient of the cost w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Hint: The computation of the cost function and gradients can be
%       efficiently vectorized. For example, consider the computation
%
%           sigmoid(X * theta)
%
%       Each row of the resulting matrix will contain the value of the
%       prediction for that example. You can make use of this to vectorize
%       the cost function and gradient computations.
%
% Hint: When computing the gradient of the regularized cost function,
%       there're many possible vectorized solutions, but one solution
%       looks like:
%           grad = (unregularized gradient for logistic regression)
%           temp = theta;
%           temp(1) = 0;   % because we don't add anything for j = 0
%           grad = grad + YOUR_CODE_HERE (using the temp variable)
%

htheta = sigmoid(X * theta);
J = 1 / m * sum(-y .* log(htheta) - (1 - y) .* log(1 - htheta)) + lambda / (2 * m) * sum(theta(2:end) .^ 2);
temp = theta;
temp(1) = 0;
grad = 1 / m * (X' * (htheta - y) + lambda * temp);

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

end

author: merwan

repo: ml-class

save_path: github-repos/MATLAB/merwan-ml-class


text:
function [y,swp]=mvk2(a,x,eps,nswp,z,rmax,varargin)
%Two-sided DMRG fast matrix-by-vector product
%   [Y,SWP]=MVK2(A,X,EPS,[NSWP],[Y],[RMAX],[OPTIONS]) Two-sided DMRG (mvk
%   is one-sided). Matrix-by-vector product of a TT-matrix A
%   by a TT-tensor X with accuracy EPS. Also, one can specify the number of
%   sweeps NSWP, initial approximation Z and the maximal TT-rank RMAX (if
%   they become too large)
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------
n=a.n;
m=a.m;
att=a.tt;
corea=att.core;
psa=att.ps;
ra=att.r;
d=att.d;
%start_val='fast_mv';
start_val='rough_mv';
%start_val='random';
rmax_loc=8; %For the "rough" matrix-by-vector product
if ( nargin <= 5 || isempty(rmax) )
   rmax=1000;
end
if ( nargin <= 4 || isempty(z) )

    if ( strcmp(start_val,'random') )
    rz1=rank(a,1)*rank(x,1);
    rzd=rank(a,d+1)*rank(x,d+1);
    kf=5;
    rz=[rz1;kf*ones(d-1,1);rzd];
    z=tt_rand(n,ndims(x),rz);
    elseif (strcmp(start_val,'rough_mv'))
        rmax_loc=8; %
        xloc=round(x,0,rmax_loc);
        aloc=round(a,0,rmax_loc);
        z=round(aloc*xloc,eps);
        %z=aloc*xloc;
    elseif (strcmp(start_val,'fast_mv') )
         rz1=rank(a,1)*rank(x,1);
      rzd=rank(a,d+1)*rank(x,d+1);
      kf=5;
      rz=[rz1;kf*ones(d-1,1);rzd];
      z=tt_rand(n,ndims(x),rz);
       z=mvk2(a,x,max(eps,1e-2),10,z,rmax); %First, do it with bad accuracy
    end
end
if ( nargin <= 3 || isempty(nswp) )
  nswp = 40;
end
y=z;

%Parameters section
kick_rank=6;
verb=false;

%Warmup is to orthogonalize Y from right-to-left and compute psi-matrices
%for Ax
psi=cell(d+1,1); %Psi-matrices
%psi{d+1}=1; psi{1}=1;
psi{d+1}=eye(rank(y,d+1));
psi{1}=eye(rank(y,1));
%Here we will add convergence test
%Warmup: right-to-left QR + computation of psi matrices

corex=x.core;
psx=x.ps;
rx=x.r;

swp=1;
converged=false;
while (swp <= nswp && ~converged)
psy=y.ps;
  ry=y.r;
  corey=y.core;
     pos1=psy(d+1);
cr1=corey(psy(d):psy(d+1)-1);

for i=d:-1:2
   cr2=corey(psy(i-1):psy(i)-1);
   cr1=reshape(cr1,[ry(i),n(i)*ry(i+1)]);
   cr2=reshape(cr2,[ry(i-1)*n(i-1),ry(i)]);
   cr1=cr1.';
   [q,rm]=qr(cr1,0); rn=size(q,2); rm=rm.';
   q=q.';
   ry(i)=rn;
   corey(pos1-ry(i+1)*n(i)*ry(i):pos1-1)=q(:);
   %Convolution is now performed for psi(i) using psi(i+1) and corea, and
   %(new) core q
   cra=corea(psa(i):psa(i+1)-1); cra=reshape(cra,[ra(i),n(i),m(i),ra(i+1)]);
   cry=reshape(conj(q),[ry(i),n(i),ry(i+1)]);
   crx=corex(psx(i):psx(i+1)-1); crx=reshape(crx,[rx(i),m(i),rx(i+1)]);
   pscur=psi{i+1}; pscur=reshape(pscur,[ra(i+1),rx(i+1),ry(i+1)]); %ra,rx,ry
   %First, convolve over rx(i+1)
   crx=reshape(crx,[rx(i)*m(i),rx(i+1)]);
   pscur=permute(pscur,[2,1,3]); pscur=reshape(pscur,[rx(i+1),ra(i+1)*ry(i+1)]);
   pscur=crx*pscur; %pscur is now rx(i)*m(i)*ra(i+1)*ry(i+1)
   %Convolve over m(i),ra(i+1),n(i),ry(i+1)
    pscur=reshape(pscur,[rx(i),m(i)*ra(i+1),ry(i+1)]);
    pscur=permute(pscur,[1,3,2]);
    pscur=reshape(pscur,[rx(i)*ry(i+1),m(i)*ra(i+1)]);
    cra=reshape(cra,[ra(i)*n(i),m(i)*ra(i+1)]); cra=cra.';
    pscur=pscur*cra;
    %pscur is now rx(i)*ry(i+1)*ra(i)*n(i), it is left to convolve over
    %n(i)*ry(i+1)
    pscur=reshape(pscur,[rx(i),ry(i+1),ra(i),n(i)]);
    pscur=permute(pscur,[3,1,4,2]);
    pscur=reshape(pscur,[rx(i)*ra(i),n(i)*ry(i+1)]);
    cry=reshape(cry,[ry(i),n(i)*ry(i+1)]); cry=cry.';
    pscur=pscur*cry;
    psi{i}=pscur;
   %End of psi-block
   pos1=pos1-ry(i+1)*n(i)*ry(i);
   cr1=cr2*rm;

end
corey(pos1-ry(2)*n(1)*ry(1):pos1-1)=cr1(:);
pos1=pos1-ry(2)*n(1)*ry(1);
corey=corey(pos1:numel(corey)); %Truncate unused elements

  %left-to-right dmrg sweep
  pos1=1;
  cry_old=corey;
  psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
  converged=true;
  ermax=0;
  for i=1:d-1
     %We care for two cores, with number i & number i+1, and use
     %psi(i) and psi(i+2) as a basis; also we will need to recompute
     %psi(i+1)
     ps1=psi{i}; ps2=psi{i+2};
     cra1=corea(psa(i):psa(i+1)-1); cra2=corea(psa(i+1):psa(i+2)-1);
     crx1=corex(psx(i):psx(i+1)-1); crx2=corex(psx(i+1):psx(i+2)-1);
     %our convolution is
     %ps1(ra(i),rx(i),ry(i))*cra1(ra(i),n(i),m(i),ra(i+1))*
     %cra2(ra(i+1),n(i+1),m(i+1),ra(i+2))*
     %*ps2(ra(i+2),rx(i+2),ry(i+2))
     %*crx1(rx(i),m(i),rx(i+1))*cr2x(rx(i+1)*m(i+1)*rx(i+2))
     %Scheme ps1*crx1 over rx(i)

     ps1=reshape(ps1,[ra(i),rx(i),ry(i)]);
     ps1=permute(ps1,[1,3,2]); ps1=reshape(ps1,[ra(i)*ry(i),rx(i)]);
     crx1=reshape(crx1,[rx(i),m(i)*rx(i+1)]);
     ps1=ps1*crx1; %ps1 is now ra(i)*ry(i)*m(i)*rx(i+1)
     %Now convolve with matrix A over ra(i)*m(i)
     ps1=reshape(ps1,[ra(i),ry(i),m(i),rx(i+1)]);
     ps1=permute(ps1,[2,4,1,3]);
     ps1=reshape(ps1,[ry(i)*rx(i+1),ra(i)*m(i)]);
     cra1=reshape(cra1,[ra(i),n(i),m(i),ra(i+1)]);
     cra1=permute(cra1,[1,3,2,4]);
     cra1=reshape(cra1,[ra(i)*m(i),n(i)*ra(i+1)]);
     ps1=ps1*cra1; %ps1 is now ry(i)*rx(i+1)*n(i)*ra(i+1)
     %Then the ``same'' convolution is carried over for second pair of
     %cores
     ps2=reshape(ps2,[ra(i+2),rx(i+2),ry(i+2)]);
     ps2=permute(ps2,[2,1,3]);
     ps2=reshape(ps2,[rx(i+2),ra(i+2)*ry(i+2)]);
     crx2=reshape(crx2,[rx(i+1)*m(i+1),rx(i+2)]);
     ps2=crx2*ps2; %ps2 is now rx*(i+1)*m(i+1)*ra(i+2)*ry(i+2)
     %Convolve over m(i+1)*ra(i+2)
     ps2=reshape(ps2,[rx(i+1),m(i+1)*ra(i+2),ry(i+2)]);
     ps2=permute(ps2,[2,1,3]);
     ps2=reshape(ps2,[m(i+1)*ra(i+2),rx(i+1)*ry(i+2)]);
     cra2=reshape(cra2,[ra(i+1)*n(i+1),m(i+1)*ra(i+2)]);
     ps2=cra2*ps2;
     %ps2 is now ra(i+1)*n(i+1)*rx(i+1)*ry(i+2)
     %Now form superblock by contraction ps1 & ps2
     %over ra(i+1)*rx(i+1)
     ps0=reshape(ps1,[ry(i),rx(i+1),n(i),ra(i+1)]);
     ps0=permute(ps0,[1,3,2,4]);
     ps0=reshape(ps0,[ry(i)*n(i),rx(i+1)*ra(i+1)]);
     ps2=reshape(ps2,[ra(i+1),n(i+1),rx(i+1),ry(i+2)]);
     ps2=permute(ps2,[3,1,2,4]);
     ps2=reshape(ps2,[rx(i+1)*ra(i+1),n(i+1)*ry(i+2)]);
     super_core=ps0*ps2; %super_core is ry(i)*n(i)*n(i+1)*ry(i+2)
    if ( i > 1 )
     %Compute previous supercore
        cr1=corey(pos1:pos1+ry(i)*n(i)*ry(i+1)-1);
        cr2=cry_old(psy(i+1):psy(i+2)-1);
        cr1=reshape(cr1,[ry(i)*n(i),ry(i+1)]);
        cr2=reshape(cr2,[ry(i+1),n(i+1)*ry(i+2)]);
         super_core_old=cr1*cr2;
         er=norm(super_core_old(:)-super_core(:))/norm(super_core(:));

        if ( er > eps )
            converged=false;
        end
        ermax=max(er,ermax);
     %if ( verb )
     %  fprintf('i=%d er=%3.2e \n',i,er);
     %end
     end
     [u,s,v]=svd(super_core,'econ');
     s=diag(s);
     r=my_chop2(s,eps/sqrt(d-1)*norm(s)); r=min(r,rmax);
     u=u(:,1:r); s=s(1:r); v=v(:,1:r); v=v*diag(s);

     %Kick rank

     ur=randn(size(u,1),kick_rank);
     %Orthogonalize ur to u by Golub-Kahan reorth
     u=reort(u,ur);
     radd=size(u,2)-r;
     if ( radd > 0 )
        vr=zeros(size(v,1),radd);
        v=[v,vr];
     end
     r=size(u,2);

     ry(i+1)=r;
     %u is ry(i)*n(i)*ry(i+1)
     %core_new(pos1:pos1+ry(i)*n(i)*ry(i+1)-1)=u(:);
     corey(pos1:pos1+ry(i)*n(i)*ry(i+1)-1)=u(:);

     u=reshape(u,[ry(i),n(i),ry(i+1)]);

     %Compute new psi
     %ps1 is ry(i)*rx(i+1)*n(i)*ra(i+1) with u over ry(i)*n(i)
     u=conj(u);
     ps1=reshape(ps1,[ry(i),rx(i+1),n(i),ra(i+1)]);
     ps1=permute(ps1,[4,2,1,3]);
     ps1=reshape(ps1,[ra(i+1)*rx(i+1),ry(i)*n(i)]);
     u=reshape(u,[ry(i)*n(i),ry(i+1)]);
     ps1=ps1*u;
     psi{i+1}=ps1;
     %Compute (?) new v
     pos1=pos1+ry(i)*n(i)*ry(i+1);
     v=v';
     %v=reshape(v,[ry(i+1),n(i+1),ry(i+2)]);
     corey(pos1:pos1+ry(i+1)*n(i+1)*ry(i+2)-1)=v(:);
  end
  psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
y.core=corey;
y.r=ry;
y.ps=psy;
swp=swp+1;
if ( verb )
fprintf('swp=%d er=%3.2e trunk=%3.2e \n',swp,ermax,eps/sqrt(d-1));
end
end
if ( swp == nswp )
  fprintf('mvk2 warning: error is not fixed for maximal number of sweeps %d\n', swp);
end%end
%p1=tt_mvdot(core(a),core(x),y0);
%p2=dot(y,tt_tensor(y0));
%abs(p1-p2)/abs(p1)
%keyboard;

author: oseledets

repo: TT-Toolbox


text:
function lik = lik_qgp(varargin)
%LIK_QGP  Create a Quantile Gaussian Process likelihood (utility) structure
%
%  Description
%    LIK = LIK_QGP('PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    creates a quantile gp likelihood structure in which the named
%    parameters have the specified values. Any unspecified
%    parameters are set to default values.
%
%    LIK = LIK_QGP(LIK,'PARAM1',VALUE1,'PARAM2,VALUE2,...)
%    modify a likelihood function structure with the named
%    parameters altered with the specified values.
%
%    Parameters for QGP likelihood function [default]
%      sigma2       - variance [0.1]
%      sigma2_prior - prior for sigma2 [prior_logunif]
%      quantile     - Quantile of interest [0.5]
%
%    Note! If the prior is 'prior_fixed' then the parameter in
%    question is considered fixed and it is not handled in
%    optimization, grid integration, MCMC etc.
%
%    The likelihood is defined as follows:
%                            __ n
%      p(y|f, sigma2, tau) = || i=1 tau*(1-tau)/sigma*exp(-(y-f)/sigma*
%                                 (tau - I(t <= f)))
%
%    where tau is the quantile of interest, sigma is the standard deviation
%    of the distribution and I(t <= f) = 1 if t <= f, 0 otherwise.
%
%    Note that because the form of the likelihood, second order derivatives
%    with respect to latent values are 0. Because this, EP should be used
%    instead of Laplace approximation.
%
%  See also
%    GP_SET, PRIOR_*, LIK_*
%
%   Reference
%     Boukouvalas et al. (2012). Direct Gaussian Process Quantile
%     Regression Using Expectation Propagation. In Proceedings of the
%     29th International Conference on Machine Learning, Edinburgh,
%     Scotland, UK, 2012.
%
%
% Copyright (c) 2012 Ville Tolvanen

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

  ip=inputParser;
  ip.FunctionName = 'LIK_QGP';
  ip.addOptional('lik', [], @isstruct);
  ip.addParamValue('sigma2',0.1, @(x) isscalar(x) && x>0);
  ip.addParamValue('sigma2_prior',prior_logunif(), @(x) isstruct(x) || isempty(x));
  ip.addParamValue('quantile',0.5, @(x) isscalar(x) && x>0 && x<1);
  ip.parse(varargin{:});
  lik=ip.Results.lik;

  if isempty(lik)
    init=true;
    lik.type = 'QGP';
  else
    if ~isfield(lik,'type') || ~isequal(lik.type,'QGP')
      error('First argument does not seem to be a valid likelihood function structure')
    end
    init=false;
  end

  % Initialize parameters
  if init || ~ismember('sigma2',ip.UsingDefaults)
    lik.sigma2 = ip.Results.sigma2;
  end
  if init || ~ismember('quantile',ip.UsingDefaults)
    lik.quantile = ip.Results.quantile;
  end
  % Initialize prior structure
  if init
    lik.p=[];
  end
  if init || ~ismember('sigma2_prior',ip.UsingDefaults)
    lik.p.sigma2=ip.Results.sigma2_prior;
  end
  if init
    % Set the function handles to the subfunctions
    lik.fh.pak = @lik_qgp_pak;
    lik.fh.unpak = @lik_qgp_unpak;
    lik.fh.lp = @lik_qgp_lp;
    lik.fh.lpg = @lik_qgp_lpg;
    lik.fh.ll = @lik_qgp_ll;
    lik.fh.llg = @lik_qgp_llg;
    lik.fh.llg2 = @lik_qgp_llg2;
    lik.fh.llg3 = @lik_qgp_llg3;
    lik.fh.tiltedMoments = @lik_qgp_tiltedMoments;
    lik.fh.siteDeriv = @lik_qgp_siteDeriv;
    lik.fh.predy = @lik_qgp_predy;
    lik.fh.invlink = @lik_qgp_invlink;
    lik.fh.recappend = @lik_qgp_recappend;
  end

end

function [w s h] = lik_qgp_pak(lik)
%LIK_QGP_PAK  Combine likelihood parameters into one vector.
%
%  Description
%    W = LIK_QGP_PAK(LIK) takes a likelihood structure LIK
%    and combines the parameters into a single row vector W.
%    This is a mandatory subfunction used for example in
%    energy and gradient computations.
%
%       w = [ log(lik.sigma2)
%             (hyperparameters of lik.magnSigma2)]'
%
%  See also
%    LIK_QGP_UNPAK

  w = []; s = {}; h=[];
  if ~isempty(lik.p.sigma2)
    w = [w log(lik.sigma2)];
    s = [s; 'log(qgp.sigma2)'];
    h = [h 0];
    % Hyperparameters of sigma2
    [wh, sh, hh] = lik.p.sigma2.fh.pak(lik.p.sigma2);
    w = [w wh];
    s = [s; sh];
    h = [h hh];
  end

end

function [lik, w] = lik_qgp_unpak(lik, w)
%LIK_QGP_UNPAK  Extract likelihood parameters from the vector.
%
%  Description
%    W = LIK_QGP_UNPAK(W, LIK) takes a likelihood structure
%    LIK and extracts the parameters from the vector W to the LIK
%    structure. This is a mandatory subfunction used for example
%    in energy and gradient computations.
%
%    Assignment is inverse of
%       w = [ log(lik.sigma2)
%             (hyperparameters of lik.magnSigma2)]'
%
%  See also
%    LIK_QGP_PAK

  if ~isempty(lik.p.sigma2)
    lik.sigma2 = exp(w(1));
    w = w(2:end);

    % Hyperparameters of sigma2
    [p, w] = lik.p.sigma2.fh.unpak(lik.p.sigma2, w);
    lik.p.sigma2 = p;
  end
end

function lp = lik_qgp_lp(lik)
%LIK_QGP_LP  Evaluate the log prior of likelihood parameters
%
%  Description
%    LP = LIK_QGP_LP(LIK) takes a likelihood structure LIK and
%    returns log(p(th)), where th collects the parameters. This
%    subfunction is needed when there are likelihood parameters.
%
%  See also
%    LIK_QGP_PAK, LIK_QGP_UNPAK, LIK_QGP_G, GP_E

  lp = 0;

  if ~isempty(lik.p.sigma2)
    likp=lik.p;
    lp = likp.sigma2.fh.lp(lik.sigma2, likp.sigma2) + log(lik.sigma2);
  end
end

function lpg = lik_qgp_lpg(lik)
%LIK_QGP_LPG  Evaluate gradient of the log prior with respect
%                  to the parameters.
%
%  Description
%    LPG = LIK_QGP_LPG(LIK) takes a QGP likelihood
%    function structure LIK and returns LPG = d log (p(th))/dth,
%    where th is the vector of parameters. This subfunction is
%    needed when there are likelihood parameters.
%
%  See also
%    LIK_QGP_PAK, LIK_QGP_UNPAK, LIK_QGP_E, GP_G

  lpg = [];

  if ~isempty(lik.p.sigma2)
    likp=lik.p;

    lpgs = likp.sigma2.fh.lpg(lik.sigma2, likp.sigma2);
    lpg = lpgs(1).*lik.sigma2 + 1;
    if length(lpgs) > 1
      lpg = [lpg lpgs(2:end)];
    end
  end
end

function ll = lik_qgp_ll(lik, y, f, z)
%LIK_QGP_LL  Log likelihood
%
%  Description
%    LL = LIK_QGP_LL(LIK, Y, F, Z) takes a likelihood
%    structure LIK, observations Y and latent values F.
%    Returns the log likelihood, log p(y|f,z). This subfunction
%    is needed when using Laplace approximation or MCMC for
%    inference with non-Gaussian likelihoods. This subfunction
%    is also used in information criteria (DIC, WAIC) computations.
%
%  See also
%    LIK_QGP_LLG, LIK_QGP_LLG3, LIK_QGP_LLG2, GPLA_E

  tau=lik.quantile;
  sigma=sqrt(lik.sigma2);
  ll = sum(log(tau*(1-tau)/sigma) - (y-f)./sigma.*(tau-(y<=f)));
end

function llg = lik_qgp_llg(lik, y, f, param, z)
%LIK_QGP_LLG  Gradient of the log likelihood
%
%  Description
%    LLG = LIK_QGP_LLG(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, observations Y and latent values F. Returns
%    the gradient of the log likelihood with respect to PARAM.
%    At the moment PARAM can be 'param' or 'latent'. This subfunction
%    is needed when using Laplace approximation or MCMC for inference
%    with non-Gaussian likelihoods.
%
%  See also
%    LIK_QGP_LL, LIK_QGP_LLG2, LIK_QGP_LLG3, GPLA_E

  tau=lik.quantile;
  sigma2=sqrt(lik.sigma2);
  switch param
    case 'param'
      llg = sum(-1/(2.*sigma2) + (y-f)./(2.*sigma2^(3/2)).*(tau-(y<=f)));

      % correction for the log transformation
      llg = llg.*lik.sigma2;
    case 'latent'
      llg = (tau-(y<=f))/sqrt(sigma2);
  end
end

function llg2 = lik_qgp_llg2(lik, y, f, param, z)
%LIK_QGP_LLG2  Second gradients of the log likelihood
%
%  Description
%    LLG2 = LIK_QGP_LLG2(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, observations Y and latent values F. Returns
%    the Hessian of the log likelihood with respect to PARAM.
%    At the moment PARAM can be 'param' or 'latent'. LLG2 is
%    a vector with diagonal elements of the Hessian matrix
%    (off diagonals are zero). This subfunction is needed
%    when using Laplace approximation or EP for inference
%    with non-Gaussian likelihoods.
%
%  See also
%    LIK_QGP_LL, LIK_QGP_LLG, LIK_QGP_LLG3, GPLA_E

  tau=lik.quantile;
  sigma2=lik.sigma2;
  switch param
    case 'param'
      llg2 = sum(1/(2*sigma2^2) - 3.*(tau-(y<=f)).*(y-f)./(4.*sigma2^(5/2)));

      % correction due to the log transformation
      llg2 = llg2.*lik.sigma2;
    case 'latent'
      llg2 = zeros(size(f));
    case 'latent+param'
      llg2 = -(tau-(y<=f))./(2*sigma2^(3/2));

      % correction due to the log transformation
      llg2 = llg2.*lik.disper;
  end
end

function llg3 = lik_qgp_llg3(lik, y, f, param, z)
%LIK_QGP_LLG3  Third gradients of the log likelihood
%
%  Description
%    LLG3 = LIK_QGP_LLG3(LIK, Y, F, PARAM) takes a likelihood
%    structure LIK, observations Y and latent values F and
%    returns the third gradients of the log likelihood with
%    respect to PARAM. At the moment PARAM can be 'param' or
%    'latent'. LLG3 is a vector with third gradients. This
%    subfunction is needed when using Laplace approximation for
%    inference with non-Gaussian likelihoods.
%
%  See also
%    LIK_QGP_LL, LIK_QGP_LLG, LIK_QGP_LLG2, GPLA_E, GPLA_G

  tau=lik.quantile;
  sigma2=lik.sigma2;
  switch param
    case 'param'
      llg3 = sum(-1/sigma2^3 + 15.*(tau-(y<=f)).*(y-f)./(8.*sigma2^(7/2)));
    case 'latent'
      llg3 = 0;
    case 'latent2+param'
      llg3 = 0;

      % correction due to the log transformation
      llg3 = llg3.*lik.sigma2;
  end
end

function [logM_0, m_1, sigm2hati1] = lik_qgp_tiltedMoments(lik, y, i1, sigma2_i, myy_i, z)
%LIK_QGP_TILTEDMOMENTS  Returns the marginal moments for EP algorithm
%
%  Description
%    [M_0, M_1, M2] = LIK_QGP_TILTEDMOMENTS(LIK, Y, I, S2,
%    MYY, Z) takes a likelihood structure LIK, observations
%    Y, index I and cavity variance S2 and mean MYY. Returns
%    the zeroth moment M_0, mean M_1 and variance M_2 of the
%    posterior marginal (see Rasmussen and Williams (2006):
%    Gaussian processes for Machine Learning, page 55). This
%    subfunction is needed when using EP for inference with
%    non-Gaussian likelihoods.
%
%  See also
%    GPEP_E

  yy = y(i1);
  sigma2 = lik.sigma2;
  tau=lik.quantile;
  logM_0=zeros(size(yy));
  m_1=zeros(size(yy));
  sigm2hati1=zeros(size(yy));

  for i=1:length(i1)
    if isscalar(sigma2_i)
      sigma2_ii = sigma2_i;
    else
      sigma2_ii = sigma2_i(i);
    end

    % get a function handle of an unnormalized tilted distribution
    % (likelihood * cavity = Quantile-GP * Gaussian)
    % and useful integration limits
    [tf,minf,maxf]=init_qgp_norm(yy(i),myy_i(i),sigma2_ii,sigma2,tau);

    % Integrate with quadrature
    RTOL = 1.e-6;
    ATOL = 1.e-10;
    [m_0, m_1(i), m_2] = quad_moments(tf, minf, maxf, RTOL, ATOL);
    if isnan(m_0)
      logM_0=NaN;
      return
    end
    sigm2hati1(i) = m_2 - m_1(i).^2;

    % If the second central moment is less than cavity variance
    % integrate more precisely. Theoretically for log-concave
    % likelihood should be sigm2hati1 < sigm2_i.
    if sigm2hati1(i) >= sigma2_ii
      ATOL = ATOL.^2;
      RTOL = RTOL.^2;
      [m_0, m_1(i), m_2] = quad_moments(tf, minf, maxf, RTOL, ATOL);
      sigm2hati1(i) = m_2 - m_1(i).^2;
      if sigm2hati1(i) >= sigma2_ii
        error('lik_qgp_tilted_moments: sigm2hati1 >= sigm2_i');
      end
    end
    logM_0(i) = log(m_0);
  end
end

function [g_i] = lik_qgp_siteDeriv(lik, y, i1, sigm2_i, myy_i, z)
%LIK_QGP_SITEDERIV  Evaluate the expectation of the gradient
%                      of the log likelihood term with respect
%                      to the likelihood parameters for EP
%
%  Description [M_0, M_1, M2] =
%    LIK_QGP_SITEDERIV(LIK, Y, I, S2, MYY, Z) takes a
%    likelihood structure LIK, observations Y, index I
%    and cavity variance S2 and mean MYY. Returns E_f
%    [d log p(y_i|f_i) /d a], where a is the likelihood
%    parameter and the expectation is over the marginal posterior.
%    This term is needed when evaluating the gradients of
%    the marginal likelihood estimate Z_EP with respect to
%    the likelihood parameters (see Seeger (2008):
%    Expectation propagation for exponential families).This
%    subfunction is needed when using EP for inference with
%    non-Gaussian likelihoods and there are likelihood parameters.
%
%  See also
%    GPEP_G

  yy = y(i1);
  sigma2=lik.sigma2;
  tau=lik.quantile;

  % get a function handle of an unnormalized tilted distribution
  % (likelihood * cavity = Quantile-GP * Gaussian)
  % and useful integration limits
  [tf,minf,maxf]=init_qgp_norm(yy,myy_i,sigm2_i,sigma2,tau);
  % additionally get function handle for the derivative
  td = @deriv;

  % Integrate with quadgk
  [m_0, fhncnt] = quadgk(tf, minf, maxf);
  [g_i, fhncnt] = quadgk(@(f) td(f).*tf(f)./m_0, minf, maxf);
  g_i = g_i.*sigma2;

  function g = deriv(f)

    g = -1/(2.*sigma2) + (yy-f)./(2.*sigma2^(3/2)).*(tau-(yy<=f));

  end
end

function [lpy, Ey, Vary] = lik_qgp_predy(lik, Ef, Varf, yt, zt)
%LIK_QGP_PREDY  Returns the predictive mean, variance and density of y
%
%  Description
%    LPY = LIK_QGP_PREDY(LIK, EF, VARF YT, ZT)
%    Returns logarithm of the predictive density PY of YT, that is
%        p(yt | zt) = \int p(yt | f, zt) p(f|y) df.
%    This subfunction is needed when computing posterior predictive
%    distributions for future observations.
%
%    [LPY, EY, VARY] = LIK_QGP_PREDY(LIK, EF, VARF) takes a
%    likelihood structure LIK, posterior mean EF and posterior
%    Variance VARF of the latent variable and returns the
%    posterior predictive mean EY and variance VARY of the
%    observations related to the latent variables. This
%    subfunction is needed when computing posterior predictive
%    distributions for future observations.
%

%
%  See also
%    GPLA_PRED, GPEP_PRED, GPMC_PRED

  sigma2=lik.sigma2;
  tau=lik.quantile;

  Ey=[];
  Vary=[];

  % Evaluate the posterior predictive densities of the given observations
  lpy = zeros(length(yt),1);
  if (size(Ef,2) > 1) && (size(Ef,2) > 1) && size(yt,2) == 1
    % Approximate integral with sum of grid points when using corrected
    % marginal posterior
    for i1=1:length(yt)
      py = arrayfun(@(f) exp(lik.fh.ll(lik, yt(i1), f, [])), Ef(i1,:));
      pf = Varf(i1,:)./sum(Varf(i1,:));
      lpy(i1) = log(sum(py.*pf));
    end
  else
    for i1=1:length(yt)
      % get a function handle of the likelihood times posterior
      % (likelihood * posterior = Quantile-GP * Gaussian)
      % and useful integration limits
      [pdf,minf,maxf]=init_qgp_norm(...
        yt(i1),Ef(i1),Varf(i1),sigma2, tau);
      % integrate over the f to get posterior predictive distribution
      lpy(i1) = log(quadgk(pdf, minf, maxf));
    end
  end
end

function [df,minf,maxf] = init_qgp_norm(yy,myy_i,sigm2_i,sigma2,tau)
%INIT_QGP_NORM
%
%  Description
%    Return function handle to a function evaluating
%    Quantile-GP * Gaussian which is used for evaluating
%    (likelihood * cavity) or (likelihood * posterior) Return
%    also useful limits for integration. This is private function
%    for lik_qgp. This subfunction is needed by subfunctions
%    tiltedMoments, siteDeriv and predy.
%
%  See also
%    LIK_QGP_TILTEDMOMENTS, LIK_QGP_SITEDERIV,
%    LIK_QGP_PREDY

  sigma=sqrt(sigma2);
% avoid repetitive evaluation of constant part
  ldconst = log(tau*(1-tau)/sigma) ...
            - log(sigm2_i)/2 - log(2*pi)/2;
  % Create function handle for the function to be integrated
  df = @qgp_norm;
  % use log to avoid underflow, and derivates for faster search
  ld = @log_qgp_norm;
  ldg = @log_qgp_norm_g;
%   ldg2 = @log_qgp_norm_g2;

  % Set the limits for integration
  % Quantile-GP likelihood is log-concave so the qgp_norm
  % function is unimodal, which makes things easier
  if yy==0
    % with yy==0, the mode of the likelihood is not defined
    % use the mode of the Gaussian (cavity or posterior) as a first guess
    modef = myy_i;
  else
    % use precision weighted mean of the Gaussian approximation
    % of the Quantile-GP likelihood and Gaussian
    modef = (myy_i/sigm2_i + yy/sigma2)/(1/sigm2_i + 1/sigma2);
  end
  % find the mode of the integrand using Newton iterations
  % few iterations is enough, since the first guess in the right direction
  niter=8;       % number of Newton iterations

  minf=modef-6*sigm2_i;
  while ldg(minf) < 0
    minf=minf-2*sigm2_i;
  end
  maxf=modef+6*sigm2_i;
  while ldg(maxf) > 0
    maxf=maxf+2*sigm2_i;
  end
  for ni=1:niter
%     h=ldg2(modef);
    modef=0.5*(minf+maxf);
    if ldg(modef) < 0
      maxf=modef;
    else
      minf=modef;
    end
  end
  % integrand limits based on Gaussian approximation at mode
  minf=modef-6*sqrt(sigm2_i);
  maxf=modef+6*sqrt(sigm2_i);
  modeld=ld(modef);
  iter=0;
  % check that density at end points is low enough
  lddiff=20; % min difference in log-density between mode and end-points
  minld=ld(minf);
  step=1;
  while minld>(modeld-lddiff)
    minf=minf-step*sqrt(sigm2_i);
    minld=ld(minf);
    iter=iter+1;
    step=step*2;
    if iter>100
      error(['lik_qgp -> init_qgp_norm: ' ...
             'integration interval minimun not found ' ...
             'even after looking hard!'])
    end
  end
  maxld=ld(maxf);
  iter=0;
  step=1;
  while maxld>(modeld-lddiff)
    maxf=maxf+step*sqrt(sigm2_i);
    maxld=ld(maxf);
    iter=iter+1;
    step=step*2;
    if iter>100
      error(['lik_qgp -> init_qgp_norm: ' ...
             'integration interval maximun not found ' ...
             'even after looking hard!'])
    end
  end

  function integrand = qgp_norm(f)
  % Quantile-GP * Gaussian
    integrand = exp(ldconst ...
                    -(yy-f)./sqrt(sigma2).*(tau-(yy<=f)) ...
                    -0.5*(f-myy_i).^2./sigm2_i);
  end

  function log_int = log_qgp_norm(f)
  % log(Quantile-GP * Gaussian)
  % log_qgp_norm is used to avoid underflow when searching
  % integration interval
    log_int = ldconst...
              -(yy-f)./sqrt(sigma2).*(tau-(yy<=f)) ...
              -0.5*(f-myy_i).^2./sigm2_i;
  end

  function g = log_qgp_norm_g(f)
  % d/df log(Quantile-GP * Gaussian)
  % derivative of log_qgp_norm
    g = (tau-(yy<=f))/sqrt(sigma2) ...
        + (myy_i - f)./sigm2_i;
  end

end

function mu = lik_qgp_invlink(lik, f, z)
%LIK_QGP_INVLINK  Returns values of inverse link function
%
%  Description
%    MU = LIK_QGP_INVLINK(LIK, F) takes a likelihood structure LIK and
%    latent values F and returns the values MU of inverse link function.
%    This subfunction is needed when using function gp_predprctmu.
%
%     See also
%     LIK_QGP_LL, LIK_QGP_PREDY

  mu = f;
end

function reclik = lik_qgp_recappend(reclik, ri, lik)
%RECAPPEND  Append the parameters to the record
%
%  Description
%    RECLIK = LIK_QGP_RECAPPEND(RECLIK, RI, LIK) takes a
%    likelihood record structure RECLIK, record index RI and
%    likelihood structure LIK with the current MCMC samples of
%    the parameters. Returns RECLIK which contains all the old
%    samples and the current samples from LIK.  This subfunction
%    is needed when using MCMC sampling (gp_mc).
%
%  See also
%    GP_MC

  if nargin == 2
    % Initialize the record
    reclik.type = 'Quantile-GP';

    % Initialize parameter
    reclik.sigma2 = [];

    % Set the function handles
    reclik.fh.pak = @lik_qgp_pak;
    reclik.fh.unpak = @lik_qgp_unpak;
    reclik.fh.lp = @lik_qgp_lp;
    reclik.fh.lpg = @lik_qgp_lpg;
    reclik.fh.ll = @lik_qgp_ll;
    reclik.fh.llg = @lik_qgp_llg;
    reclik.fh.llg2 = @lik_qgp_llg2;
    reclik.fh.llg3 = @lik_qgp_llg3;
    reclik.fh.tiltedMoments = @lik_qgp_tiltedMoments;
    reclik.fh.predy = @lik_qgp_predy;
    reclik.fh.invlink = @lik_qgp_invlink;
    reclik.fh.recappend = @lik_qgp_recappend;
    reclik.p=[];
    reclik.p.sigma2=[];
    if ~isempty(ri.p.sigma2)
      reclik.p.sigma2 = ri.p.sigma2;
    end
  else

    % Append to the record
    reclik.sigma2(ri,:)=lik.sigma2;
    if ~isempty(lik.p.sigma2)
      reclik.p.sigma2 = lik.p.sigma2.fh.recappend(reclik.p.sigma2, ri, lik.p.sigma2);
    end
    reclik.quantile = lik.quantile;
  end
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function roi = roiLoad(pth, format, referenceMR);
% Load a mrVista2 ROI.
%
% roi = roiLoad(<pth='dialog'>, <format='Inplane'>, <referenceMR>);
%
% Possible formats are 'Inplane' (mrVista Inplane ROI),
% 'Volume' (mrVista Volume/Gray ROI), and 'mat' (saved MATLAB file,
% mrVista2 format).
%
% ras, 02/24/07.
if notDefined('pth')
	pth = mrvSelectFile('r', 'mat', [], 'Select an ROI to Load');
end

if notDefined('format'),	format = 'Inplane';		end

if exist(pth, 'file') | exist([pth '.mat'], 'file')
	roi = load(pth);
else
	error(sprintf('ROI path %s not found. ', pth))
end

% if we successfully loaded an ROI file, check for
% the format: first check if it's a mrVista 1 ROI:
fields = fieldnames(roi);
if ismember(fields, 'ROI')
	% it's a mrVista 1 ROI: use mrVista-specific tools
	% to parse it
	mrGlobals2;
	roi = roi.ROI; % subsume ROI struct, only thing of interest
	if isequal(lower(roi.viewType), 'inplane')
		reference = 'Inplane';
		referenceMR = fullfile(HOMEDIR, 'Inplane', 'anat.mat');
		voxelSize = mrSESSION(1).inplanes.voxelSize;
	elseif ismember(lower(roi.viewType), {'volume', 'gray'})
		reference = 'I|P|R';
		if notDefined('referenceMR')
			referenceMR = mrLoadHeader(getVAnatomyPath);
		else
			referenceMR = mrParse(referenceMR);
		end
		voxelSize = referenceMR.voxelSize(1:3);
	else
		error('ROI view type not supported yet.')

	end

	template = roiCreate(reference, roi.coords, 'color', roi.color, ...
					'name', roi.name, 'voxelSize', voxelSize, ...
					'referenceMR', referenceMR);
	for f = setdiff( fieldnames(template), fieldnames(roi) )'
		roi.(f{1}) = template.(f{1});
	end

	for f = setdiff( fieldnames(roi), fieldnames(template) )'
		roi = rmfield(roi, f{1});
	end

	return
end

return


text:
function tokens = opt_getToken(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

author: PatternRecognition

repo: OpenBMI


text:
classdef TestFisheyeEstimateNewCameraMatrixForUndistortRectify
    %TestFisheyeEstimateNewCameraMatrixForUndistortRectify

    methods (Static)
        function test_1
            %TODO
        end

        function test_error_argnum
            try
                cv.fisheyeEstimateNewCameraMatrixForUndistortRectify();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end


text:
classdef SMEA < ALGORITHM
% <multi> <real/integer>
% Self-organizing multiobjective evolutionary algorithm
% D    ---     --- Number of neurons in each dimension of the latent space
% tau0 --- 0.7 --- Initial learning rate
% H    ---   5 --- Size of neighborhood mating pools

%------------------------------- Reference --------------------------------
% H. Zhang, A. Zhou, S. Song, Q. Zhang, X. Gao, and J. Zhang, A self-
% organizing multiobjective evolutionary algorithm, IEEE Transactions on
% Evolutionary Computation, 2016, 20(5): 792-806.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            [D,tau0,H] = Algorithm.ParameterSet(repmat(ceil(Problem.N.^(1/(Problem.M-1))),1,Problem.M-1),0.7,5);
            Problem.N  = prod(D);
            sigma0     = sqrt(sum(D.^2)/(Problem.M-1))/2;

            %% Generate random population
            Population = Problem.Initialization();
            FrontNo    = NDSort(Population.objs,inf);

            %% Initialize the SOM
            % Training set
            S = Population.decs;
            % Weight vector of each neuron
            W = S;
            % Position of each neuron
            D = arrayfun(@(S)1:S,D,'UniformOutput',false);
            eval(sprintf('[%s]=ndgrid(D{:});',sprintf('c%d,',1:length(D))))
            eval(sprintf('Z=[%s];',sprintf('c%d(:),',1:length(D))))
            % Distance between each two neurons in latent space
            LDis = pdist2(Z,Z);
            % H nearest neurons of each neuron in latent space
            [~,B] = sort(LDis,2);
            B     = B(:,2:min(H+1,end));

            %% Optimization
            while Algorithm.NotTerminated(Population)
                % Update SOM
                for s = 1 : size(S,1)
                    sigma  = sigma0*(1-(Problem.FE+s)/Problem.maxFE);
                    tau    = tau0*(1-(Problem.FE+s)/Problem.maxFE);
                    [~,u1] = min(pdist2(S(s,:),W));
                    U      = LDis(u1,:) < sigma;
                    W(U,:) = W(U,:) + tau.*repmat(exp(-LDis(u1,U))',1,size(W,2)).*(repmat(S(s,:),sum(U),1)-W(U,:));
                end

                % Associate each solution with a neuron
                A  = 1 : Problem.N;
                U  = 1 : Problem.N;
                XU = zeros(1,Problem.N);
                for i = 1 : Problem.N
                    x        = randi(length(A));
                    [~,u]    = min(pdist2(Population(A(x)).dec,W(U,:)));
                    XU(U(u)) = A(x);
                    A(x)     = [];
                    U(u)     = [];
                end

                % Evolution
                A = Population.decs;
                for u = 1 : Problem.N
                    drawnow('limitrate');
                    if rand < 0.9
                        Q = XU(B(u,:));
                    else
                        Q = 1 : Problem.N;
                    end
                    Q = Q(randperm(end,2));
                    y = OperatorDE(Problem,Population(u),Population(Q(1)),Population(Q(2)));
                    [Population,FrontNo] = Select(Population,FrontNo,y);
                end

                % Update the training set
                S = setdiff(Population.decs,A,'rows');
            end
        end
    end
end


text:
function fd1d_heat_explicit_test03 ( )

%*****************************************************************************80
%
%% FD1D_HEAT_EXPLICIT_TEST03 does a simple test problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD1D_HEAT_EXPLICIT_TEST03:\n' );
  fprintf ( 1, '  Compute an approximate solution to the time-dependent\n' );
  fprintf ( 1, '  one dimensional heat equation:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    dH/dt - K * d2H/dx2 = f(x,t)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Run a simple test case.\n' );
%
%  Heat coefficient.
%
  k = k_test03 ( );
%
%  X_NUM is the number of equally spaced nodes to use between 0 and 1.
%
  x_num = 21;
  x_min = -5.0;
  x_max = +5.0;
  dx = ( x_max - x_min ) / ( x_num - 1 );
  x = linspace ( x_min, x_max, x_num );
%
%  T_NUM is the number of equally spaced time points between 0 and 10.0.
%
  t_num = 81;
  t_min = 0.0;
  t_max = 4.0;
  dt = ( t_max - t_min ) / ( t_num - 1 );
  t = linspace ( t_min, t_max, t_num );
%
%  Get the CFL coefficient.
%
  cfl = fd1d_heat_explicit_cfl ( k, t_num, t_min, t_max, x_num, x_min, x_max );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of X nodes = %d\n', x_num );
  fprintf ( 1, '  X interval is [%f,%f]\n', x_min, x_max );
  fprintf ( 1, '  X spacing is %f\n', dx );
  fprintf ( 1, '  Number of T values = %d\n', t_num );
  fprintf ( 1, '  T interval is [%f,%f]\n', t_min, t_max );
  fprintf ( 1, '  T spacing is %f\n', dt );
  fprintf ( 1, '  Constant K = %g\n', k );
  fprintf ( 1, '  CFL coefficient = %g\n', cfl );
%
%  Running the code produces an array H of temperatures H(t,x),
%  and vectors x and t.
%

  hmat = zeros ( x_num, t_num );

  for j = 1 : t_num
    if ( j == 1 )
      h = ic_test03 ( x_num, x, t(j) );
      h = bc_test03 ( x_num, x, t(j), h );
    else
      h = fd1d_heat_explicit ( x_num, x, t(j-1), dt, cfl, @rhs_test03, @bc_test03, h );
    end
    hmat(1:x_num,j) = h(1:x_num);
  end
%
%  Plot the data.
%
  figure ( 3 )
  [ tmat, xmat ] = meshgrid ( t, x );
  mesh ( tmat, xmat, hmat );
  title ( 'H(X,T) for TEST03 computed by FD1D\_HEAT\_EXPLICIT' );
  xlabel ( '<-- Time -->' );
  ylabel ( '<-- X -->' );
  zlabel ( '<--  H(X,T) -->' );
%
%  Write the data to files.
%
  r8mat_write ( 'h_test03.txt', x_num, t_num, hmat );
  r8vec_write ( 't_test03.txt', t_num, t );
  r8vec_write ( 'x_test03.txt', x_num, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  H(X,T) written to "h_test03.txt"\n' );
  fprintf ( 1, '  T values written to "t_test03.txt"\n' );
  fprintf ( 1, '  X values written to "x_test3.txt"\n' );
  return
end


text:
%PNTSET Set Dirichlet point constraints.
%
%   [ F, INDROW, AMAT, T ] = PNTSET( PROB, F, AMAT, ISYMM, SET_NULL )
%   Sets Dirichlet point constraints conditions in the right hand side
%   load vector F and global matrix AMAT with the information in the finite
%   element problem struct PROB.
%
%       Input       Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       prob        struct                 Finite element problem struct
%       f           (neq,1)                Right hand side/load vector
%       amat        (n_a,n_a)              System matrix (sparse or triplet format)
%       isymm       scalar/{0}             Symmetrize BCs if applicable.
%       set_null    scalar/{0}             Set zeros in f vector.
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       f           (neq,1)                Modified right hand side/load vector
%       indrow      (neq,1)                Index to rows (dofs) in rhs which were set
%       amat        (n_a,n_a)              Modified system matrix
%       t           scalar                 Time spent in function

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics