text:
% Add Manopt to the path and make all manopt components available.

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Jan. 3, 2013.
% Contributors:
% Change log:
%   Aug.  7, 2013 (NB): Changed to work without the import command
%                      (new structure of the toolbox).
%   Aug.  8, 2013 (NB): Changed to use addpath_recursive, home brewed.
%   Aug. 22, 2013 (NB): Using genpath instead of homecooked
%                       addpath_recursive.

addpath(pwd);

% Recursively add Manopt directories to the Matlab path.
cd manopt;
addpath(genpath(pwd));
cd ..;

author: andrewssobral

repo: lrslibrary


text:
function plotHistConv(model, sample, rxnNames, nSubSamples)
% Plots convergence of sample histograms
%
% USAGE:
%
%    plotHistConv(model, sample, rxnNames, nSubSamples)
%
% INPUTS:
%    model:          COBRA model structure
%    sample:         Sampled fluxes
%    rxnNames:       List of reactions to plot
%    nSubSamples:    Number of sub samples
%
% EXAMPLE:
%
%    example 1:
%    rxnNames = {'R1', 'R2'}
%    plotHistConv(model, sample, rxnNames, nSubSamples)
%
% .. Author: - Markus Herrgard 8/14/06

nSkip = 10;
nBin = 20;

[nRxns,nSamples] = size(sample);

[isInModel,rxnInd] = ismember(rxnNames,model.rxns);
rxnInd = rxnInd(isInModel);

nPlotRxn = sum(isInModel);
rxnNames = rxnNames(isInModel);

nCol = ceil(sqrt(nPlotRxn));
nRow = ceil(nPlotRxn/nCol);

subSampleSize = floor(nSamples/nSubSamples);

clf
for rxnID = 1:nPlotRxn
    subplot(nRow,nCol,rxnID);
    hold on
    maxx = -1e9;
    minx = 1e9;
    maxy = 0;
    for subID = 1:nSubSamples
        [n,x] = hist(sample(rxnInd(rxnID), 1:nSkip:(subSampleSize * subID))', nBin);
        plot(x, n/sum(n));
        maxx = max([max(x) maxx]);
        minx = min([min(x) minx]);
        maxy = max([max(n/sum(n)) maxy]);
    end
    axis([minx maxx 0 maxy]);
    title(rxnNames{rxnID});
end

author: opencobra

repo: cobratoolbox


text:
function A = TriArea3D(X1,X2,X3)

% Xi n*3 are the three points of triangles where n is the number of
% triangle

x1 = X1(:,1);  y1 = X1(:,2);   z1 = X1(:,3);
x2 = X2(:,1);  y2 = X2(:,2);   z2 = X2(:,3);
x3 = X3(:,1);  y3 = X3(:,2);   z3 = X3(:,3);
A = 1/2*(((x1-x3).*(y2-y1) - (x1-x2).*(y3-y1)).^2 + ...
    ((y1-y3).*(z2-z1) - (y1-y2).*(z3-z1)).^2 + ...
    ((z1-z3).*(x2-x1) - (z1-z2).*(x3-x1)).^2).^(1/2);

A = abs(A);

end

repo: ifem


text:
function theta = vector_separation_nd ( dim_num, v1, v2, theta )

%*****************************************************************************80
%
%% VECTOR_SEPARATION_ND finds the angular separation between vectors in ND.
%
%  Discussion:
%
%    Any two vectors lie in a plane, and are separated by a plane angle.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the vectors.
%
%    Input, real V1(DIM_NUM), V2(DIM_NUM), the two vectors.
%
%    Output, real THETA, the angle between the two vectors.
%
  v1_norm = sqrt ( sum ( v1(1:dim_num).^2 ) );

  v2_norm = sqrt ( sum ( v2(1:dim_num).^2 ) );

  cos_theta = ( v1(1:dim_num) * v2(1:dim_num)' ) / ( v1_norm * v2_norm );

  theta = r8_acos ( cos_theta );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function mri = ctf_avw2mri(avw)

% mri = ctf_avw2mri(avw)
%
% The purpose of this function is to convert an Analyze volume into a CTF
% .mri volume.  It currently requires that the Analyze volume is
% 256x256x256, 1x1x1 mm voxels.  The avw_read function can
% handle various Analyze orientations, so this function assumes that the
% avw.img volume has axial unflipped orientation (it always is when
% returned by avw_read, regardless of the format in the .img file). The
% returned mri struct in the matlab workspace can be saved using
% ctf_write_mri.
%
% This function depends on the avw* functions availabe in the mri_toolbox
% at http://eeg.sf.net/.  In particular, the avw_read and associated
% functions for reading an Analyze volume into a matlab avw struct, used as
% input to this function.
%
%      <>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> %
%      <                                                      > %
%      <                    DISCLAIMER:                       > %
%      <                                                      > %
%      < THIS PROGRAM IS INTENDED FOR RESEARCH PURPOSES ONLY. > %
%      < THIS PROGRAM IS IN NO WAY INTENDED FOR CLINICAL OR   > %
%      <                    OFFICIAL USE.                     > %
%      <                                                      > %
%      <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<> %
%

% $Revision: 1.1 $ $Date: 2009-01-30 03:49:26 $

% Copyright (C) 2003  Darren L. Weber
%
% 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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

% History:  08/2003, Darren.Weber_at_radiology.ucsf.edu
%                    - adapted from an appendex to CTF document
%                    MRIConverter.pdf, which is copied at the end of this
%                    function.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mri = ctf_mri_make;

mri.file = [avw.fileprefix,'.mri'];

ver = '[$Revision: 1.1 $]';
fprintf('\nCTF_AVW2MRI [v%s]\n',ver(12:16));  tic;

% these checks for the volume dims/pixel size could be replaced with an
% interpolation function that could take any Analyze volume and convert it
% to 256x256x256, 1x1x1 mm
if avw.hdr.dime.dim(2) ~= 256,
    error('avw.hdr.dime.dim(2) ~= 256');
end
if avw.hdr.dime.dim(3) ~= 256,
    error('avw.hdr.dime.dim(3) ~= 256');
end
if avw.hdr.dime.dim(4) ~= 256,
    error('avw.hdr.dime.dim(4) ~= 256');
end

if avw.hdr.dime.pixdim(2) ~= 1,
    error('avw.hdr.dime.dim(2) ~= 256');
end
if avw.hdr.dime.pixdim(3) ~= 1,
    error('avw.hdr.dime.dim(3) ~= 256');
end
if avw.hdr.dime.pixdim(4) ~= 1,
    error('avw.hdr.dime.dim(4) ~= 256');
end

% mri.hdr.dataSize = 1 or 2 (bytes), 8 or 16 bits
if avw.hdr.dime.bitpix == 8,
    mri.hdr.dataSize = 1;
    mri.hdr.clippingRange = 255;
else
    mri.hdr.dataSize = 2;
    mri.hdr.clippingRange = 65536;
end

% This next step should always work correctly, given that avw_read always
% converts any Analyze orientation into axial unflipped in the matlab
% workspace variable avw.img; the axial unflipped orientation is
% (+X left,  +Y anterior,  +Z superior); the ctf orientation is
% (+X right, +Y posterior, +Z inferior), the complete opposite.
temp = flipdim(avw.img,1);
temp = flipdim(temp,2);
temp = flipdim(temp,3);
mri.img = temp;

t=toc; fprintf('...done (%5.2f sec).\n\n',t);

return

author: PatternRecognition

repo: OpenBMI


text:
% Test file for trigtech/real.m

function pass = test_real(pref)

% Get preferences.
if ( nargin < 1 )
    pref = trigtech.techPref();
end

testclass = trigtech();

% Test a scalar-valued function.
f = testclass.make(@(x) exp(20i*pi*x) + 1i*sin(100*pi*x), [], pref);
g = testclass.make(@(x) cos(20*pi*x), [], pref);
h = real(f);
g = prolong(g, length(h));
pass(1) = norm(h.coeffs - g.coeffs, inf) < 10*vscale(h)*eps;

% Test an array-valued function.
f = testclass.make(@(x) [exp(20i*pi*x) + 1i*sin(100*pi*x), -exp(10i*pi*x)], [], pref);
g = testclass.make(@(x) [cos(20*pi*x), -real(exp(10i*pi*x))], [], pref);
h = real(f);
n = max(length(g),length(h));
g = prolong(g,n); h = prolong(h,n);
pass(2) = norm(h.coeffs - g.coeffs, inf) < 10*max(vscale(h)*eps);

% Test a real function.
f = 1i*testclass.make(@(x) cos(30*pi*x), [], pref);
g = real(f);
pass(3) = numel(g.coeffs) == 1 && g.coeffs == 0;

% Test an array-valued real function.
f = 1i*testclass.make(@(x) [cos(99*pi*x), sin(99*pi*x), exp(cos(pi*x))], [], pref);
g = real(f);
pass(4) = all(size(g.coeffs) == [1, 3]) && all(g.coeffs == 0);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% patchExtrude
% Below is a demonstration of the features of the |patchExtrude| function

%%
clear; close all; clc;

%% Syntax
% |[Et,Vt]=patchExtrude(F,V,zRange);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function addCOBRABinaryPathToSystemPath()
% Sets up the path to system executables shipped by the toolbox.
%
% USAGE:
%     addCOBRABinaryPathToSystemPath()
%

global CBTDIR

arch = computer('arch');

%Set thefolders
binaryDir = [CBTDIR filesep 'binary' filesep arch filesep 'bin'];
libDir = [CBTDIR filesep 'binary' filesep arch filesep 'lib'];

%get all Subfolders
bindirs = rdir([binaryDir filesep '**' filesep '*']);
libdirs = rdir([libDir filesep '**' filesep '*']);
binnames = {bindirs.name};
libnames = {libdirs.name};
bindirs = [{binaryDir},binnames([bindirs.isdir])]; %Add all subdirectories and the top level
libdirs = [{libDir},libnames([libdirs.isdir])];

if ispc
    setenv('Path', [getenv('Path') ';' strjoin(union(bindirs,libdirs),';')  ]);
else
    setenv('PATH', [getenv('PATH') ':' strjoin(bindirs,':')]);
    setenv('LD_LIBRARY_PATH', [getenv('LD_LIBRARY_PATH') ':' strjoin(libdirs,':')]);
end

end


text:
function varargout = drawRay(ray, varargin)
%DRAWRAY Draw a ray on the current axis.
%
%   drawRay(RAY)
%   With RAY having the syntax: [x0 y0 dx dy], draws the ray starting from
%   point (x0 y0) and going to direction (dx dy), clipped with the current
%   window axis.
%
%   drawRay(RAY, PARAMS, VALUE)
%   Can specify param-pair values.
%
%   H = drawRay(...)
%   Returns handle on line object
%
%   See also
%     rays2d, drawLine, clipRay
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2003-10-31
% Copyright 2003-2022 INRA - TPV URPOI - BIA IMASTE

% extract handle of axis to draw in
if isAxisHandle(ray)
    ax = ray;
    ray = varargin{1};
    varargin(1) = [];
else
    ax = gca;
end

% get bounding box limits
box = axis(ax);

% compute clipped shapes
[clipped, isInside] = clipRay(ray, box);

% allocate memory for handle
h = -ones(size(ray, 1), 1);

% draw visible rays
h(isInside) = drawEdge(ax, clipped(isInside, :), varargin{:});

% process output
if nargout > 0
    varargout = {h};
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [ prob, ier ] = mdbeta ( x, p, q )

%*****************************************************************************80
%
%% MDBETA evaluates the incomplete beta function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2008
%
%  Author:
%
%    Original FORTRAN77 version by Oliver Ludwig;
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Oliver Ludwig,
%    Algorithm 179:
%    Incomplete Beta Ratio,
%    Communications of the ACM,
%    Volume 6, Number 6, June 1963, page 314.
%
%  Parameters:
%
%    Input, real X, the value to which function is to be
%    integrated.  X must be in the range [0,1] inclusive.
%
%    Input, real P, the first parameter.  P must be greater
%    than 0.0.
%
%    Input, real Q, the second parameter.  Q must be greater
%    than 0.0.
%
%    Output, real PROB.  The probability that a random variable
%    from a Beta distribution having parameters P and Q will be less than
%    or equal to X.
%
%    Output, integer IER, error parameter.
%    0, normal exit.
%    1, X is not in the range [0,1] inclusive.
%    2, P or Q is less than or equal to 0.
%
%  Local parameters:
%
%    Local, real ALEPS, the logarithm of EPS1.
%
%    Local, real EPS, the machine precision.
%
%    Local, real EPS1, the smallest representable number.
%
  aleps = - 179.6016;
  eps = 2.2E-16;
  eps1 = 1.0E-78;
%
%  Check ranges of the arguments.
%
  prob = 0.0;
  y = x;

  if ( x < 0.0 || 1.0 < x )
    ier = 1;
    return
  end

  if ( p <= 0.0 || q <= 0.0 )
    ier = 2;
    return
  end

  ier = 0;

  if ( x <= 0.5 )
    interval = 0;
  else
    interval = 1;
    temp = p;
    p = q;
    q = temp;
    y = 1.0 - y;
  end

  if ( x == 0.0 || x == 1.0 )

    prob = 0.0;

    if ( interval ~= 0 )
      prob = 1.0 - prob;
      temp = p;
      p = q;
      q = temp;
    end

    return
  end

  ib = q;
  temp = ib;
  ps = q - ib;

  if ( q == temp )
    ps = 1.0;
  end

  dp = p;
  dq = q;
  px = dp * log ( y );
  pq = alogam ( dp + dq );
  p1 = alogam ( dp );
  c = alogam ( dq );
  d4 = log ( dp );
  xb = px + alogam ( ps + dp ) - alogam ( ps ) - d4 - p1;
%
%  Scaling
%
  ib = floor ( xb / aleps );
  infsum = 0.0;
%
%  First term of a decreasing series will underflow.
%
  if ( ib == 0 )

    infsum = exp ( xb );
    cnt = infsum * dp;
%
%  CNT will equal dexp ( temp ) * ( 1.d0 - ps ) * i * p * y**i / factorial ( i ).
%
    wh = 0.0;

    while ( 1 )

      wh = wh + 1.0;
      cnt = cnt * ( wh - ps ) * y / wh;
      xb = cnt / ( dp + wh );
      infsum = infsum + xb;

      if ( xb / eps < infsum )
        break
      end

    end

  end

  finsum = 0.0;

  if ( dq <= 1.0 )

    prob = finsum + infsum;

    if ( interval ~= 0 )
      prob = 1.0 - prob;
      temp = p;
      p = q;
      q = temp;
    end

    return
  end

  xb = px + dq * log ( 1.0 - y ) + pq - p1 - log ( dq ) - c;
%
%  Scaling.
%
  ib = floor ( xb / aleps );

  if ( ib < 0 )
    ib = 0;
  end

  c = 1.0 / ( 1.0 - y );
  cnt = exp ( xb - ib * aleps );
  ps = dq;
  wh = dq;

  while ( 1 )

    wh = wh - 1.0;

    if ( wh <= 0.0 )

      prob = finsum + infsum;

      if ( interval ~= 0 )
        prob = 1.0 - prob;
        temp = p;
        p = q;
        q = temp;
      end

      break

    end

    px = ( ps * c ) / ( dp + wh );

    if ( px <= 1.0 )

      if ( cnt / eps <= finsum || cnt <= eps1 / px )

        prob = finsum + infsum;

        if ( interval ~= 0 )
          prob = 1.0 - prob;
          temp = p;
          p = q;
          q = temp;
        end

        break

      end

    end

    cnt = cnt * px;
%
%  Rescale.
%
    if ( 1.0 < cnt )
      ib = ib - 1;
      cnt = cnt * eps1;
    end

    ps = wh;

    if ( ib == 0 )
      finsum = finsum + cnt;
    end

  end

  return
end


text:
function aedat = RemoveInvalidEvents(aedat)

%{
This function removes any events which are flagged as invalid, and then
removes the 'valid' fields. It's possible that all events get removed. In
this case the data type field itself needs to be removed, but this should
be handled by a subsequent call to NumEventsByType.
%}

dbstop if error

if ~isfield(aedat, 'data')
    return
end
% Special

if isfield(aedat.data, 'special') && isfield(aedat.data.special, 'valid')
    aedat.data.special.timeStamp    = aedat.data.special.timeStamp  (aedat.data.special.valid);
    aedat.data.special.address      = aedat.data.special.address    (aedat.data.special.valid);
    aedat.data.special = rmfield(aedat.data.special, 'valid');
end

% Polarity

if isfield(aedat.data, 'polarity') && isfield(aedat.data.polarity, 'valid')
    aedat.data.polarity.timeStamp   = aedat.data.polarity.timeStamp (aedat.data.polarity.valid);
    aedat.data.polarity.x           = aedat.data.polarity.x         (aedat.data.polarity.valid);
    aedat.data.polarity.y           = aedat.data.polarity.y         (aedat.data.polarity.valid);
    aedat.data.polarity.polarity    = aedat.data.polarity.polarity  (aedat.data.polarity.valid);
    aedat.data.polarity = rmfield(aedat.data.polarity, 'valid');
end
% Frames

if isfield(aedat.data, 'frame') && isfield(aedat.data.frame, 'valid')
    aedat.data.frame.timeStampStart    = aedat.data.frame.timeStampStart (aedat.data.frame.valid);
    aedat.data.frame.timeStampEnd      = aedat.data.frame.timeStampEnd   (aedat.data.frame.valid);
    aedat.data.frame.samples           = aedat.data.frame.samples        (aedat.data.frame.valid);
    aedat.data.frame.xLength           = aedat.data.frame.xLength        (aedat.data.frame.valid);
    aedat.data.frame.yLength           = aedat.data.frame.yLength        (aedat.data.frame.valid);
    aedat.data.frame.xPosition         = aedat.data.frame.xPosition      (aedat.data.frame.valid);
    aedat.data.frame.yPosition         = aedat.data.frame.yPosition      (aedat.data.frame.valid);
    aedat.data.frame.roiId             = aedat.data.frame.roiId          (aedat.data.frame.valid);
    aedat.data.frame.colorChannels     = aedat.data.frame.colorChannels  (aedat.data.frame.valid);
    aedat.data.frame.colorFilter       = aedat.data.frame.colorFilter    (aedat.data.frame.valid);
    aedat.data.frame = rmfield(aedat.data.frame, 'valid');
end
% Imu6

if isfield(aedat.data, 'imu6') && isfield(aedat.data.imu6, 'valid')
    aedat.data.imu6.timeStamp    = aedat.data.imu6.timeStamp    (aedat.data.imu6.valid);
    aedat.data.imu6.accelX       = aedat.data.imu6.accelX       (aedat.data.imu6.valid);
    aedat.data.imu6.accelY       = aedat.data.imu6.accelY       (aedat.data.imu6.valid);
    aedat.data.imu6.accelZ       = aedat.data.imu6.accelZ       (aedat.data.imu6.valid);
    aedat.data.imu6.gyroX        = aedat.data.imu6.gyroX        (aedat.data.imu6.valid);
    aedat.data.imu6.gyroY        = aedat.data.imu6.gyroY        (aedat.data.imu6.valid);
    aedat.data.imu6.gyroZ        = aedat.data.imu6.gyroZ        (aedat.data.imu6.valid);
    aedat.data.imu6.temperature  = aedat.data.imu6.temperature  (aedat.data.imu6.valid);
    aedat.data.imu6 = rmfield(aedat.data.imu6, 'valid');
end

if isfield(aedat.data, 'sample') && isfield(aedat.data.sample, 'valid')
    aedat.data.sample.timeStamp     = aedat.data.sample.timeStamp   (aedat.data.sample.valid);
    aedat.data.sample.sampleType	= aedat.data.sample.sampleType  (aedat.data.sample.valid);
    aedat.data.sample.sample		= aedat.data.sample.sample      (aedat.data.sample.valid);
    aedat.data.sample = rmfield(aedat.data.sample, 'valid');
end

if isfield(aedat.data, 'ear') && isfield(aedat.data.ear, 'valid')
    aedat.data.ear.timeStamp	= aedat.data.ear.timeStamp  (aedat.data.ear.valid);
    aedat.data.ear.position     = aedat.data.ear.position   (aedat.data.ear.valid);
    aedat.data.ear.channel		= aedat.data.ear.channel	(aedat.data.ear.valid);
    aedat.data.ear.neuron		= aedat.data.ear.neuron		(aedat.data.ear.valid);
    aedat.data.ear.filter		= aedat.data.ear.filter		(aedat.data.ear.valid);
    aedat.data.ear = rmfield(aedat.data.ear, 'valid');
end

if isfield(aedat.data, 'point1D') && isfield(aedat.data.point1D, 'valid')
    aedat.data.point1D.type         = aedat.data.point1D.type         (aedat.data.point1D.valid);
    aedat.data.point1D.timeStamp    = aedat.data.point1D.timeStamp    (aedat.data.point1D.valid);
    aedat.data.point1D.x            = aedat.data.point1D.x            (aedat.data.point1D.valid);
    aedat.data.point1D = rmfield(aedat.data.point1D, 'valid');
end

if isfield(aedat.data, 'point2D') && isfield(aedat.data.point2D, 'valid')
    aedat.data.point2D.type         = aedat.data.point2D.type       (aedat.data.point2D.valid);
    aedat.data.point2D.timeStamp    = aedat.data.point2D.timeStamp  (aedat.data.point2D.valid);
    aedat.data.point2D.x            = aedat.data.point2D.x          (aedat.data.point2D.valid);
    aedat.data.point2D.y            = aedat.data.point2D.y          (aedat.data.point2D.valid);
    aedat.data.point2D = rmfield(aedat.data.point2D, 'valid');
end

if isfield(aedat.data, 'point3D') && isfield(aedat.data.point3D, 'valid')
    aedat.data.point3D.type         = aedat.data.point3D.type       (aedat.data.point3D.valid);
    aedat.data.point3D.timeStamp    = aedat.data.point3D.timeStamp  (aedat.data.point3D.valid);
    aedat.data.point3D.x            = aedat.data.point3D.x          (aedat.data.point3D.valid);
    aedat.data.point3D.y            = aedat.data.point3D.y          (aedat.data.point3D.valid);
    aedat.data.point3D.z            = aedat.data.point3D.z          (aedat.data.point3D.valid);
    aedat.data.point3D = rmfield(aedat.data.point3D, 'valid');
end

author: panpanfei


text:
function Data=FilterJumps(Dates,Data,Name)
Chgs=diff(Data);
DatesChgs=Dates(2:end);
ChgsNoJumps=Chgs;
DatesChgsNoJumps=DatesChgs;

Reject=find(abs(Chgs)>7*median(abs(Chgs)));  % reject outliers
ChgsNoJumps(Reject)=[];
DatesChgsNoJumps(Reject)=[];

ChgsJumps=0*DatesChgs;
ChgsJumps(Reject)=Chgs(Reject);

PlotSeries(DatesChgs,Chgs,DatesChgsNoJumps,ChgsNoJumps,ChgsJumps,Name)

Data=cumsum([ChgsNoJumps]);               % reconstruct series

author: Sable

repo: mcbench-benchmarks


text:
% sgwt_check_connected : Check connectedness of graph
%
% function r=sgwt_check_connected(A)
%
% returns 1 if graph is connected, 0 otherwise
% Uses boost graph library breadth first search
%
% Inputs :
% A - adjacency matrix
%
% Outputs :
% r - result
%

% This file is part of the SGWT toolbox (Spectral Graph Wavelet Transform toolbox)
% Copyright (C) 2010, David K. Hammond.
%
% The SGWT toolbox 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.
%
% The SGWT toolbox 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 the SGWT toolbox.  If not, see <http://www.gnu.org/licenses/>.

function r=sgwt_check_connected(A)
  d=bfs(A,1);
  r=~any(d==-1);

repo: gspbox


text:
function params_struct = get_param_val(params_struct, direct_params, ...
                                       field_name, default_val)
% used by all params_* files to get the value of a particular parameter
%
% @authors:     Ahmad Humayun
% @contact:     ahumayun@cc.gatech.edu
% @affiliation: Georgia Institute of Technology
% @date:        Fall 2013 - Summer 2014

    if ~isfield(params_struct, field_name)
        params_struct.(field_name) = default_val;
    end

    % if user gave the parameter directly, then override everything else
    if any(strcmpi(direct_params, field_name))
        % the first specified value is picked
        param_pos = find(strcmpi(direct_params, field_name), 1, 'first');
        params_struct.(field_name) = direct_params{param_pos + 1};
    end
end

author: Cloud-CV

repo: object-proposals


text:
%kkprstats 'Print statistics to Cantata variables'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kprstats.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'First Input data object'
% String: wmin_var ' ', default: 'width_min': 'Name of minimum width variable'
% String: wmax_var ' ', default: 'width_max': 'Name of maximum width variable'
% String: hmin_var ' ', default: 'height_min': 'Name of minimum height variable'
% String: hmaxvar ' ', default: 'height_max': 'Name of maximum height variable'
% String: dmin_var ' ', default: 'depth_min': 'Name of minimum depth variable'
% String: dmax_var ' ', default: 'depth_max': 'Name of maximum depth variable'
% String: tmin_var ' ', default: 'time_min': 'Name of minimum time variable'
% String: tmax_var ' ', default: 'time_max': 'Name of maximum time variable'
% String: emin_var ' ', default: 'elem_min': 'Name of minimum element variable'
% String: emax_var ' ', default: 'elem_max': 'Name of maximum element variable'
% String: minval_var ' ', default: 'val_min': 'Name of minimum value variable'
% String: maxval_var ' ', default: 'val_max': 'Name of maximum value variable'
% String: mean_var ' ', default: 'mean': 'Name for mean variable'
% String: wsize_var ' ', default: 'width_size': 'Name of total width variable (KSTATS_WSIZE)'
% String: var_var ' ', default: 'variance': 'Name for variance variable'
% String: hsize_var ' ', default: 'height_size': 'Name of total height variable (KSTATS_HSIZE)'
% String: sd_var ' ', default: 'std_dev': 'Name for standard deviation variable'
% String: dsize_var ' ', default: 'depth_size': 'Name of total depth variable'
% String: rms_var ' ', default: 'rms': 'Name for rms variable'
% String: tsize_var ' ', default: 'time_size': 'Name of total time variable'
% String: psum_var ' ', default: 'psum': 'Name of positive sum variable'
% String: esize_var ' ', default: 'elem_size': 'Name of total element variable'
% String: nsum_var ' ', default: 'nsum': 'Name for negative sum variable'
% String: sum_var ' ', default: 'sum': 'Name for sum variable'
% String: skew_var ' ', default: 'skew': 'Name for skew variable'
% String: kur_var ' ', default: 'kurtosis': 'Name for kurtosis variable'
% String: pts_var ' ', default: 'total_pts': 'Name for total points variable'
% String: ppts_var ' ', default: 'pos_pts': 'Name for positive points variable'
% String: npts_var ' ', default: 'neg_pts': 'Name for negative points variable'
% String: zpts_var ' ', default: 'zero_pts': 'Name for zero-valued points variable'
%
% Example: kkprstats(i, {'i','';'wmin_var','width_min';'wmax_var','width_max';'hmin_var','height_min';'hmaxvar','height_max';'dmin_var','depth_min';'dmax_var','depth_max';'tmin_var','time_min';'tmax_var','time_max';'emin_var','elem_min';'emax_var','elem_max';'minval_var','val_min';'maxval_var','val_max';'mean_var','mean';'wsize_var','width_size';'var_var','variance';'hsize_var','height_size';'sd_var','std_dev';'dsize_var','depth_size';'rms_var','rms';'tsize_var','time_size';'psum_var','psum';'esize_var','elem_size';'nsum_var','nsum';'sum_var','sum';'skew_var','skew';'kur_var','kurtosis';'pts_var','total_pts';'ppts_var','pos_pts';'npts_var','neg_pts';'zpts_var','zero_pts'})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kprstats - Print statistics to Cantata variables
%
%  DESCRIPTION
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkprstats(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkprstats(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'wmin_var', 'width_min';'wmax_var', 'width_max';'hmin_var', 'height_min';'hmaxvar', 'height_max';'dmin_var', 'depth_min';'dmax_var', 'depth_max';'tmin_var', 'time_min';'tmax_var', 'time_max';'emin_var', 'elem_min';'emax_var', 'elem_max';'minval_var', 'val_min';'maxval_var', 'val_max';'mean_var', 'mean';'wsize_var', 'width_size';'var_var', 'variance';'hsize_var', 'height_size';'sd_var', 'std_dev';'dsize_var', 'depth_size';'rms_var', 'rms';'tsize_var', 'time_size';'psum_var', 'psum';'esize_var', 'elem_size';'nsum_var', 'nsum';'sum_var', 'sum';'skew_var', 'skew';'kur_var', 'kurtosis';'pts_var', 'total_pts';'ppts_var', 'pos_pts';'npts_var', 'neg_pts';'zpts_var', 'zero_pts'};
maxval={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
minval={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
istoggle=[0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String','String'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(0);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
callKhoros([w 'kprstats"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function data = boxesData( varargin )
% Get ground truth data for object proposal bounding box evaluation.
%
% Used to get dataset information for evaluation of bounding box object
% proposals using boxesEval.m. This requires separate download of the
% appropriate dataset first. Currently only set up for PASCAL VOC 2007 but
% it would be fairly simple to extend to other datasets.
%
% As a first step, it is necessary to download the PASCAL VOC 2007 dataset.
% Go to http://pascallin.ecs.soton.ac.uk/challenges/VOC/voc2007/ and get:
%  VOCdevkit_08-Jun-2007, VOCtrainval_06-Nov-2007, VOCtest_06-Nov-2007
% After downloading these files extract them to dataDir='boxes/VOCdevkit'.
% Once complete this function will process the data and return information
% in a convenient format for boxesEval.m.
%
% USAGE
%  data = boxesData( opts )
%
% INPUTS
%  opts       - parameters (struct or name/value pairs)
%   .resDir     - ['boxes/'] location for results and evaluation
%   .dataDir    - ['boxes/VOCdevkit/'] dir containing PASCAL VOC 2007
%   .split      - ['val'] data split for evaluation
%
% OUTPUTS
%  data       - dataset split information
%   .split      - data split for evaluation
%   .n          - number of images
%   .ids        - list of string ids
%   .imgs       - list of image filenames
%   .gt         - ground truth for each image
%
% EXAMPLE
%
% See also edgeBoxesDemo, edgeBoxes, boxesEval
%
% Structured Edge Detection Toolbox      Version 3.0
% Copyright 2014 P. Dollar and L. Zitnick.  [pdollar-at-microsoft.com]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the MSR-LA Full Rights License [see license.txt]

% get default parameters (unimportant parameters are undocumented)
dfs={ 'resDir','boxes/', 'dataDir','boxes/VOCdevkit/', 'split','val' };
o=getPrmDflt(varargin,dfs,1);

% locations of PASCAL VOC dataset
if(~exist(o.dataDir,'dir')), error('dataset not found, see help'); end
imDir=[o.dataDir 'VOC2007/JPEGImages/'];
gtDir=[o.dataDir 'VOC2007/Annotations/'];
imList=[o.dataDir 'VOC2007/ImageSets/Main/' o.split '.txt'];
addpath(genpath([o.dataDir 'VOCcode']));

% check if data already exists, if yes load and return
dataNm=[o.resDir '/GroundTruth' '-' o.split '.mat'];
if(exist(dataNm,'file')), data=load(dataNm); data=data.data; return; end

% get list of image ids
if(~exist(imList,'file')), error('ids file not found'); end
f=fopen(imList); ids=textscan(f,'%s %*s'); ids=ids{1}; fclose(f);

% generate list of image and gt filenames then load gt
n=length(ids); imgs=cell(n,1); gt=cell(n,1);
for i=1:n, imgs{i}=[imDir ids{i} '.jpg']; end
for i=1:n, gt{i}=[gtDir ids{i} '.xml']; end
if(~exist(imgs{1},'file')), error('images not found'); end
if(~exist(gt{1},'file')), error('annotations not found'); end
parfor i=1:n, [~,gt{i}]=bbGt('bbLoad',gt{i},'format',1); end

% create output structure and cache to disk
data=struct('split',o.split,'n',n,'ids',{ids},'imgs',{imgs},'gt',{gt});
if(~exist(o.resDir,'dir')), mkdir(resDir); end; save(dataNm,'data');

end


text:
function [indices, pot] = find_most_prob_entry(pot)
% function [indices, pot] = find_most_prob_entry(pot)
% function [indices, pot] = find_most_prob_entry(pot)
% Find the indices of the argmax, and set all other enties to 0.

%indices = argmax(pot.T);
[m i] = max(pot.T(:));
indices = ind2subv(pot.sizes, i);
pot.T = 0*myones(pot.sizes);
pot.T(i) = m;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function GUOdemo(Action);

persistent GUOdemoFig1 GUOdemoFig2 GUO1 GUO2 GUO1TextPos GUO2SmallPos GUO1Ready;

% function GUOdemo(Action);
%
% Demonstrates the graphicuserobject class.
% The Action argument selects the button callback code
% and should not be supplied from the command line.
%
% Copyright (c) SINUS Messtechnik GmbH 2002
% www.sinusmess.de - Sound & Vibration Instrumentation
%                  - PCB Services
%                  - Electronic Design & Production

if nargin == 0
   % No Action parameter supplied:  create figures, GUOs and other controls

   % Create figures
   GUO1Ready = 0;
   GUOdemoFig1 = figure('Position', [50 300 600 400], ...
                        'ResizeFcn', 'GUOdemo(''ResizeFig1'')');
   GUOdemoFig2 = figure('Position', [350 50 600 400]);

   % Create GUO1:  graphicuserobject with visible frame (default units=normalized)
   GUO1 = graphicuserobject(GUOdemoFig1, 'Position', [1/12 0.7 7/12 0.2], ...
                                         'Visible', 'on');
   GUO1 = uicontrol(GUO1, 'Tag', 'GUO1TextNormalized', ...
                          'Style', 'text', ...
                          'Position', [0.1 0.6 0.8 0.3], ...
                          'String', 'This text and this GUO have units=normalized');
   GUO1TextPos = [10 10 250 20];
   GUO1 = uicontrol(GUO1, 'Tag', 'GUO1TextPixels', ...
                          'Style', 'text', ...
                          'Units', 'pixels', ...
                          'Position', GUO1TextPos, ...
                          'String', 'This text and the other GUO have units=pixels');

   % Create child GUOs within GUO1
   GUOPB1N = PlaybackButtons('Tag', 'PBBinGUO1N', ...
                             'Position', [0.45 0.35 0.25 0.3], ...
                             'Visible', 'on');
   GUOPB1N = SetPBBCallbacks(GUOPB1N);  % SetPBBCallbacks is a subfunction of GUOdemo
   GUO1 = addchildguo(GUO1, GUOPB1N);   % Add GUOPB1N within GUO1
   GUOPB1P = PlaybackButtons('Tag', 'PBBinGUO1P', ...
                             'Units', 'pixels', ...
                             'Position', [260 10 80 20], ...
                             'Visible', 'on');
   GUOPB1P = SetPBBCallbacks(GUOPB1P);  % SetPBBCallbacks is a subfunction of GUOdemo
   GUO1 = addchildguo(GUO1, GUOPB1P);   % Add GUOPB1P within GUO1
   GUO1Ready = 1;  % Avoid "Child not found" error from initial resize (see below)

   % Create GUO2:  graphicuserobject with visible frame (units=pixels)
   GUO2SmallPos = [100 60 350 200];
   GUO2 = graphicuserobject(GUOdemoFig1, 'Units', 'Pixels', ...
                                         'Position', GUO2SmallPos, ...
                                         'Visible', 'on');
   % Pushbuttons are normalized relative to graphicuserobject (default units=normalized)
   GUO2 = uicontrol(GUO2, 'Style', 'pushbutton', ...
                          'Position', [0.1 0.8 0.15 0.1], ...
                          'String', 'Big', ...
                          'Callback', 'GUOdemo(''Big'')');
   GUO2 = uicontrol(GUO2, 'Style', 'pushbutton', ...
                          'Position', [0.1 0.7 0.15 0.1], ...
                          'String', 'Small', ...
                          'Callback', 'GUOdemo(''Small'')');
   GUO2 = uicontrol(GUO2, 'Style', 'pushbutton', ...
                          'Position', [0.1 0.6 0.15 0.1], ...
                          'String', 'Edit', ...
                          'Callback', 'GUOdemo(''Edit'')');
   GUO2 = uicontrol(GUO2, 'Style', 'pushbutton', ...
                          'Position', [0.1 0.5 0.15 0.1], ...
                          'String', 'Plot', ...
                          'Callback', 'GUOdemo(''Plot'')');
   GUO2 = uicontrol(GUO2, 'Style', 'pushbutton', ...
                          'Position', [0.1 0.4 0.15 0.1], ...
                          'String', 'Figure 1', ...
                          'Callback', 'GUOdemo(''Fig1'')');
   GUO2 = uicontrol(GUO2, 'Style', 'pushbutton', ...
                          'Position', [0.1 0.3 0.15 0.1], ...
                          'String', 'Figure 2', ...
                          'Callback', 'GUOdemo(''Fig2'')');
   GUO2 = uicontrol(GUO2, 'Style', 'pushbutton', ...
                          'Position', [0.1 0.2 0.15 0.1], ...
                          'String', 'Close', ...
                          'Callback', 'GUOdemo(''Close'')');
   % Axes are normalized relative to graphicuserobject (default units=normalized)
   GUO2 = axes(GUO2, 'Tag', 'PlotAxes', ...
                     'Position', [0.45 0.5 0.25 0.4], ...
                     'XLim', [0 100], ...
                     'YLim', [0 100]);
   axes(GUO2, '');  % Make GUO2 frame the current axes for rectangle
   rectangle('Position', [0.075 0.175 0.2 0.75], ...
             'Curvature', [0.1]);

   % Create child GUOs within GUO2
   GUOPB2N = PlaybackButtons('Tag', 'PBBinGUO2N', ...
                             'Position', [0.45 0.2 0.25 0.14], ...
                             'Visible', 'on');
   GUOPB2N = SetPBBCallbacks(GUOPB2N);  % SetPBBCallbacks is a subfunction of GUOdemo
   GUO2 = addchildguo(GUO2, GUOPB2N);   % Add GUOPB2N within GUO2
   GUOPB2P = PlaybackButtons('Tag', 'PBBinGUO2P', ...
                             'Units', 'pixels', ...
                             'Position', [260 10 80 20], ...
                             'Visible', 'on');
   GUOPB2P = SetPBBCallbacks(GUOPB2P);  % SetPBBCallbacks is a subfunction of GUOdemo
   GUO2 = addchildguo(GUO2, GUOPB2P);   % Add GUOPB2P within GUO2

   % Create and initialise text at fixed position relative to GUO2
   GUO2 = uicontrol(GUO2, 'Tag', 'FixedText', ...
                          'Style', 'text', ...
                          'Units', 'pixels', ...
                          'Position', [10 10 130 20]);
   GUOdemo('SetText');  % Update displayed position of GUO2 (pixels)

   % Demonstrate Display function
   GUO1
   GUO2

   % Create pushbuttons on figure GUOdemoFig1 (not in GUO)
   figure(GUOdemoFig1);
   uicontrol('Style', 'pushbutton', ...
             'Position', [100 25 100 20], ...
             'String', 'Hide GUO1', ...
             'Callback', 'GUOdemo(''HideGUO1'')');
   uicontrol('Style', 'pushbutton', ...
             'Position', [200 25 100 20], ...
             'String', 'Show GUO1', ...
             'Callback', 'GUOdemo(''ShowGUO1'')');
   uicontrol('Style', 'pushbutton', ...
             'Position', [300 25 100 20], ...
             'String', 'Hide Controls 1', ...
             'Callback', 'GUOdemo(''HideControls1'')');
   uicontrol('Style', 'pushbutton', ...
             'Position', [400 25 100 20], ...
             'String', 'Show Controls 1', ...
             'Callback', 'GUOdemo(''ShowControls1'')');
   uicontrol('Style', 'pushbutton', ...
             'Position', [100  5 100 20], ...
             'String', 'Hide GUO2', ...
             'Callback', 'GUOdemo(''HideGUO2'')');
   uicontrol('Style', 'pushbutton', ...
             'Position', [200  5 100 20], ...
             'String', 'Show GUO2', ...
             'Callback', 'GUOdemo(''ShowGUO2'')');
   uicontrol('Style', 'pushbutton', ...
             'Position', [300  5 100 20], ...
             'String', 'Hide Controls 2', ...
             'Callback', 'GUOdemo(''HideControls2'')');
   uicontrol('Style', 'pushbutton', ...
             'Position', [400  5 100 20], ...
             'String', 'Show Controls 2', ...
             'Callback', 'GUOdemo(''ShowControls2'')');

else
   % Action argument supplied (selects callback function)

   switch Action
   case 'ResizeFig1'  % ResizeFcn for GUOdemoFig1
      if GUO1Ready    % GUO1Ready avoids "Child not found" error from initial resize.
         GUO1 = resizefcn(GUO1);  % resizefcn necessary because GUO1 is normalized
      end                         % but has children with Units=pixels
   case 'Big'         % "Big" pushbutton in GUO2 resizes GUO2 and normalized children
      GUO2 = set(GUO2, 'Position', [50 60 500 300]);
      GUOdemo('SetText');  % Update displayed position of GUO2 (pixels)
   case 'Small'       % "Small" pushbutton in GUO2 resizes GUO2 and normalized children
      GUO2 = set(GUO2, 'Position', GUO2SmallPos);
      GUOdemo('SetText');  % Update displayed position of GUO2 (pixels)
   case 'Edit'        % "Edit" pushbutton in GUO2 demonstrates access to child uicontrol
      try             % "EditText":  creation, getting child properties and deletion
         MsgText = getchild(GUO2, 'EditText', 'String');
      catch
         MsgText = [];
         GUO2 = uicontrol(GUO2, 'Tag', 'EditText', ...
                                'Style', 'edit', ...
                                'Position', [0.2 0.85 0.6 0.1], ...
                                'HorizontalAlignment', 'left');
      end
      if length(MsgText) > 0
         GUO2 = deletechild(GUO2, 'EditText');
         msgbox(['EditText contained: ' MsgText], 'GUOdemo');
      end
   case 'Plot'        % "Plot" pushbutton in GUO2 demonstrates access to child axes
      % Demonstrate indexing by number instead of Tag...
      axes(GUO2, 8);  % ... but "axes(GUO2, 'PlotAxes')" would normally be better
      plot([0:100], [0:100], 'r');
      text(50, 50, 'XYZ');  % Create text relative to PlotAxes
   case 'Fig1'        % "Fig1" pushbutton in GUO2 moves GUO2 to GUOdemoFig1
      GUO2 = set(GUO2, 'Parent', GUOdemoFig1);
   case 'Fig2'        % "Fig2" pushbutton in GUO2 moves GUO2 to GUOdemoFig2
      GUO2 = set(GUO2, 'Parent', GUOdemoFig2);
   case 'Close'       % "Close" pushbutton in GUO2 ends demonstration
      GUO2 = delete(GUO2);  % Demonstrate GUO deletion
      close(GUOdemoFig1);
      close(GUOdemoFig2);
   case 'SetText'     % Update displayed position of GUO2 (pixels)
      PosText = ['Pos:  ' num2str(get(GUO2, 'Position'))];
      GUO2 = setchild(GUO2, 'FixedText', 'String', PosText);
   case {'PBBnormalized', 'PBBpixels'}  % Callback for all pushbuttons in
      Tag = get(gcbo, 'Tag');           % PlaybackButttons GUOs - see SetPBBCallbacks
      msgbox(['"' Tag '" button in child GUO with units = ' Action(4:end)], 'GUOdemo');

   % The remaining Actions are for the pushbuttons on figure GUOdemoFig1,
   % which demonstrate the hiding and showing of GUOs and their children.
   case 'HideGUO1'       % Hide GUO1
      GUO1 = hide(GUO1);
   case 'ShowGUO1'       % Show GUO1
      GUO1 = show(GUO1);
   case 'HideControls1'  % Hide normalized child GUO & control in GUO1
      GUO1 = guoeval(GUO1, 'PBBinGUO1N', 'hide');
      GUO1 = setchild(GUO1, 'GUO1TextNormalized', 'Visible', 'off');
   case 'ShowControls1'  % Show normalized child GUO & control in GUO1
      GUO1 = guoeval(GUO1, 'PBBinGUO1N', 'show');
      GUO1 = setchild(GUO1, 'GUO1TextNormalized', 'Visible', 'on');
   case 'HideGUO2'       % Hide GUO2
      GUO2 = hide(GUO2);
   case 'ShowGUO2'       % Show GUO2
      GUO2 = show(GUO2);
   case 'HideControls2'  % Hide normalized child GUO & control in GUO2
      GUO2 = guoeval(GUO2, 'PBBinGUO2N', 'hide');
      GUO2 = setchild(GUO2, 'PlotAxes', 'Visible', 'off');
      % Line and text are children of PlotAxes
      set(getchild(GUO2, 'PlotAxes', 'Children'), 'Visible', 'off');
   case 'ShowControls2'  % Show normalized child GUO & control in GUO2
      GUO2 = guoeval(GUO2, 'PBBinGUO2N', 'show');
      GUO2 = setchild(GUO2, 'PlotAxes', 'Visible', 'on');
      % Line and text are children of PlotAxes
      set(getchild(GUO2, 'PlotAxes', 'Children'), 'Visible', 'on');

   otherwise  % This should never happen...
      msgbox(['Unknown Action: ' Action], 'GUOdemo', 'error');
   end
end

function  GUOPBB = SetPBBCallbacks(GUOPBB);

% Subfunction of GUOdemo:  sets the same callback function for each of
% the four pushbuttons within GUOPBB (which is a PlaybackButtons GUO).

PBBunits = get(GUOPBB, 'Units');
for Tag = {'First', 'Previous', 'Next', 'Last'}
   GUOPBB = setchild(GUOPBB, Tag{:}, 'Callback', ['GUOdemo(''PBB' PBBunits ''')']);
end


text:
clear

scrsz = get(0,'ScreenSize');
figure1 = figure('Position',[20 50 3*scrsz(3)/4 0.9*scrsz(4)]);

set(figure1,'Units','Inches');
pos = get(figure1,'Position');
set(figure1,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])

% Create axes
axes1 = axes('Parent',figure1,'FontSize',40,'FontName','Helvetica');

line_width = 6;
hold on;

load('../../matlab_version/experiments_300VW/results/cfss_errors.mat');
[error_x, error_y] = cummErrorCurve(cfss_error_66_cat_1);
plot(error_x, error_y, 'DisplayName', 'CFSS', 'LineWidth',line_width);

load('../../matlab_version/experiments_300VW/results/cfan_errors.mat');
[error_x, error_y] = cummErrorCurve(cfan_error_66_cat_1);
plot(error_x, error_y, 'DisplayName', 'CFAN', 'LineWidth',line_width);
%
load('../../matlab_version/experiments_300VW/results/iccr_errors.mat');
[error_x, error_y] = cummErrorCurve(iccr_error_66_cat_1);
plot(error_x, error_y, 'DisplayName', 'iCCR', 'LineWidth',line_width);

load('results/300VW_OpenFace.mat');
[error_x, error_y] = cummErrorCurve(clnf_error_66_cat_1);
plot(error_x, error_y, 'DisplayName', 'OpenFace - CLNF', 'LineWidth',line_width);

% Make sure CE-CLM is drawn on top
[error_x, error_y] = cummErrorCurve(ceclm_error_66_cat_1);
plot(error_x, error_y, 'r', 'LineWidth',line_width, 'DisplayName', 'OpenFace - CE-CLM');

% Make it looks nice and print to a pdf
set(gca,'xtick',[0.01:0.01:0.08])
xlim([0.01,0.08]);
xlabel('IOD normalized MAE','FontName','Helvetica');
ylabel('Proportion of images','FontName','Helvetica');
grid on
ax=legend('show', 'Location', 'SouthEast');
ax.FontSize = 30;

print -dpdf results/300VWres_66_cat1.pdf

%%
clear

scrsz = get(0,'ScreenSize');
figure1 = figure('Position',[20 50 3*scrsz(3)/4 0.9*scrsz(4)]);

set(figure1,'Units','Inches');
pos = get(figure1,'Position');
set(figure1,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])

% Create axes
axes1 = axes('Parent',figure1,'FontSize',40,'FontName','Helvetica');

line_width = 6;
hold on;

load('../../matlab_version/experiments_300VW/results/cfss_errors.mat');
[error_x, error_y] = cummErrorCurve(cfss_error_66_cat_2);
plot(error_x, error_y, 'DisplayName', 'CFSS', 'LineWidth',line_width);

load('../../matlab_version/experiments_300VW/results/cfan_errors.mat');
[error_x, error_y] = cummErrorCurve(cfan_error_66_cat_2);
plot(error_x, error_y, 'DisplayName', 'CFAN', 'LineWidth',line_width);
%
load('../../matlab_version/experiments_300VW/results/iccr_errors.mat');
[error_x, error_y] = cummErrorCurve(iccr_error_66_cat_2);
plot(error_x, error_y, 'DisplayName', 'iCCR', 'LineWidth',line_width);

load('results/300VW_OpenFace.mat');
[error_x, error_y] = cummErrorCurve(clnf_error_66_cat_2);
plot(error_x, error_y, 'DisplayName', 'OpenFace - CLNF', 'LineWidth',line_width);

% Make sure CE-CLM is drawn on top
[error_x, error_y] = cummErrorCurve(ceclm_error_66_cat_2);
plot(error_x, error_y, 'r', 'LineWidth',line_width, 'DisplayName', 'OpenFace - CE-CLM');

% Make it looks nice and print to a pdf
set(gca,'xtick',[0.01:0.01:0.08])
xlim([0.01,0.08]);
xlabel('IOD normalized MAE','FontName','Helvetica');
ylabel('Proportion of images','FontName','Helvetica');
grid on
ax=legend('show', 'Location', 'SouthEast');
ax.FontSize = 30;

print -dpdf results/300VWres_66_cat2.pdf

%%
clear

scrsz = get(0,'ScreenSize');
figure1 = figure('Position',[20 50 3*scrsz(3)/4 0.9*scrsz(4)]);

set(figure1,'Units','Inches');
pos = get(figure1,'Position');
set(figure1,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])

% Create axes
axes1 = axes('Parent',figure1,'FontSize',40,'FontName','Helvetica');

line_width = 6;
hold on;

load('../../matlab_version/experiments_300VW/results/cfss_errors.mat');
[error_x, error_y] = cummErrorCurve(cfss_error_66_cat_3);
plot(error_x, error_y, 'DisplayName', 'CFSS', 'LineWidth',line_width);

load('../../matlab_version/experiments_300VW/results/cfan_errors.mat');
[error_x, error_y] = cummErrorCurve(cfan_error_66_cat_3);
plot(error_x, error_y, 'DisplayName', 'CFAN', 'LineWidth',line_width);
%
load('../../matlab_version/experiments_300VW/results/iccr_errors.mat');
[error_x, error_y] = cummErrorCurve(iccr_error_66_cat_3);
plot(error_x, error_y, 'DisplayName', 'iCCR', 'LineWidth',line_width);

load('results/300VW_OpenFace.mat');
[error_x, error_y] = cummErrorCurve(clnf_error_66_cat_3);
plot(error_x, error_y, 'DisplayName', 'OpenFace - CLNF', 'LineWidth',line_width);

% Make sure CE-CLM is drawn on top
[error_x, error_y] = cummErrorCurve(ceclm_error_66_cat_3);
plot(error_x, error_y, 'r', 'LineWidth',line_width, 'DisplayName', 'OpenFace - CE-CLM');

% Make it looks nice and print to a pdf
set(gca,'xtick',[0.01:0.01:0.08])
xlim([0.01,0.08]);
xlabel('IOD normalized MAE','FontName','Helvetica');
ylabel('Proportion of images','FontName','Helvetica');
grid on
ax=legend('show', 'Location', 'SouthEast');
ax.FontSize = 30;

print -dpdf results/300VWres_66_cat3.pdf

author: TadasBaltrusaitis

repo: OpenFace


text:
% This is a test script which demonstrates the usage of the "PhdTrackInitiatorX" class.
% =========================================================================>

% Load the ground truth data
load('multi-target-tracking-1.mat');
GroundTruthStateSequence = GroundTruthStateSequence(~cellfun('isempty',GroundTruthStateSequence));

% Plot settings
ShowPlots = 0;              % Set to 0 to hide plots
numTrueTracks = 3;

% Model parameter shortcuts
lambdaV = 10; % Expected number of clutter measurements over entire surveillance region
V = 5000^2;     % Volume of surveillance region (10x10 2D-grid)
V_bounds = [-2500 2500 -2500 2500]; % [x_min x_max y_min y_max]
P_D = 0.9;    % Probability of detection
timestep_duration = duration(0,0,1);

%% Models
transition_model = ConstantVelocityX('VelocityErrVariance', 5^2,...
                                     'NumDims', 2,...
                                     'TimestepDuration', timestep_duration);
% measurement_model = LinearGaussianX('NumMeasDims', 2,...
%                                     'NumStateDims', 4,...
%                                     'MeasurementErrVariance', 10^2,...
%                                     'Mapping', [1 3]);
% clutter_model = PoissonRateUniformPositionX('ClutterRate',lambdaV,...
%                                             'Limits',[V_bounds(1:2);...
%                                                       V_bounds(3:4)]);
measurement_model = RangeBearing2CartesianX('NumStateDims', 4,...
                                    'MeasurementErrVariance', [(pi/180)^2, 10^2],...
                                    'Mapping', [1 3]);
clutter_model = PoissonRateUniformPositionX('ClutterRate',lambdaV,...
                                            'Limits',[[-pi, pi];...
                                                      [0, 2500]]);
detection_model = ConstantDetectionProbabilityX('DetectionProbability',P_D);

birth_model = DistributionBasedBirthModelX('Distribution', UniformDistributionX([V_bounds(1:2); ...
                                                                                [-10 10 ];...
                                                                                V_bounds(3:4);...
                                                                                [-10 10 ]]),...
                                           'BirthIntensity', 0.0001);

% Compile the State-Space model
model = StateSpaceModelX(transition_model, measurement_model,...
                        'Clutter',clutter_model,...
                        'Detection', detection_model,...
                        'Birth', birth_model);

%% Simulate measurements, based on ground-truth
meas_simulator = MultiTargetMeasurementSimulatorX('Model',model);
DataList = meas_simulator.simulate(GroundTruthStateSequence);
N = numel(DataList);

%% Base Filter
obs_covar= measurement_model.covar();
PriorState = GaussianStateX(zeros(4,1), transition_model.covar() + blkdiag(obs_covar(1,1), 0, obs_covar(2,2),0));
base_filter = UnscentedKalmanFilterX('Model', model, 'StatePrior', PriorState);
% base_filter = ParticleFilterX('Model', model, 'StatePrior', PriorState);

%% Data Associator
config.ClutterModel = clutter_model;
config.DetectionModel = detection_model;
config.Clusterer = NaiveClustererX();
config.Gater = EllipsoidalGaterX(2,'GateLevel',10)';
config.DetectionProbability = P_D;
jpdaf = JointIntegratedProbabilisticDataAssocX(config);

%% Track Initiator

% Initiate Data Associator
config_phd.Model = model;
[priorParticles, priorWeights] = model.Birth.random(50000);
config_phd.StatePrior = ParticleStateX(priorParticles,10*priorWeights);
config_phd.BirthScheme = {'Expansion', 5000};
config_phd.SurvivalProbability = 0.99;

% Instantiate PHD filter
myphd = SMC_PHDFilterX(config_phd);

% Initiate Tag Generator
tag_gen = UuidTagGeneratorX();

% Prepare initiator parameters
config_ti.TagGenerator = tag_gen;
config_ti.InitFilter = base_filter;
config_ti.PhdFilter = myphd;
config_ti.ConfirmThreshold = 0.8;

% Create the track initiator
myti = PhdTrackInitiatorX(config_ti);

%% Track Deleter
config_td = struct('Fieldname', 'ExistenceProbability',...
                   'ReferenceValue', 0.1,...
                   'ReferenceOperand', 'lt');
mytd = FieldBasedDeleterX(config_td);

%% Metric Generators
ospa = OSPAX('CutOffThreshold',100,'Order',2);
gospa = GOSPAX('CutOffThreshold',100,'Order',2);

%% START OF SIMULATION
%  ===================>

% Create figure windows
if(ShowPlots)
    img = imread('maze.png');

    % set the range of the axes
    % The image will be stretched to this.
    min_x = 0;
    max_x = 10;
    min_y = 0;
    max_y = 10;

    % make data to plot - just a line.
    x = min_x:max_x;
    y = (6/8)*x;

    figure('units','normalized','outerposition',[0 0 .5 1])
    ax(1) = gca;
end

TrackList = [];
gospa_vals= zeros(N,4);
timestamp_km1 = DataList(1).Timestamp;
for k=2:N
    fprintf('Iteration = %d/%d\n================>\n',k,N);

    %% Extract DataList at time k
    MeasurementList = DataList(k);
    if(MeasurementList.NumMeasurements == 0)
        continue;
    end
    timestamp_km1 = DataList(k-1).Timestamp;
    timestamp_k = MeasurementList.Timestamp;
    dt = timestamp_k - timestamp_km1;
    transition_model.TimestepDuration = dt;
    fprintf('Timestamp = %s\n================>\n',timestamp_k);

    %% Process JPDAF
    jpdaf.MeasurementList = MeasurementList;
    jpdaf.TrackList = TrackList;
    jpdaf.predictTracks();
    jpdaf.associate(TrackList, MeasurementList);
    jpdaf.updateTracks();

    %% Perform Track initiation
    [TrackList] = myti.initiateTracks(jpdaf.TrackList, MeasurementList, jpdaf.AssocWeightsMatrix);

    %% Perform Track deletion
    TrackList = mytd.deleteTracks(TrackList);

    %% Evaluate metric
    [gospa_vals(k,1), gospa_vals(k,2), gospa_vals(k,3), gospa_vals(k,4)]= ...
        gospa.evaluate(GroundTruthStateSequence{k},jpdaf.TrackList);
    %% Plot update step results
    if(ShowPlots)

        cla(ax(1));
        %imagesc(ax(1),[min_x max_x], [min_y max_y], flipud(img));
        hold on;
        if(exist('data_plot','var'))
            delete(data_plot);
        end
        data_inv = measurement_model.finv(MeasurementList.Vectors);
        data_plot = plot(ax(1), data_inv(1,:), data_inv(3,:),'k*','MarkerSize', 10);

        % Plot confirmed tracks
        for j=1:numel(TrackList)
            means = [TrackList{j}.Trajectory.Mean];
            h2 = plot(ax(1), means(1,:),means(3,:),'-','LineWidth',1);
            h2 = plotgaussellipse(TrackList{j}.Filter.StatePosterior.Mean([1 3]),...
                                  TrackList{j}.Filter.StatePosterior.Covar([1 3],[1 3]),...
                                  'Color','r',...
                                  'Axis',ax(1));
        end

        % set the y-axis back to normal.
        set(ax(1),'ydir','normal');
        str = sprintf('Robot positions (Update)');
        title(ax(1),str)
        xlabel('X position (m)')
        ylabel('Y position (m)')
        axis(ax(1),V_bounds)
        pause(0.01)
    end
end
figure
title("GOSPA");
subplot(2,3,[1 2 3]), plot(1:k,gospa_vals(1:k,1)), title("GOSPA Metric");
subplot(2,3,4), plot(1:k,gospa_vals(1:k,2)), title("GOSPA Localisation");
subplot(2,3,5), plot(1:k,gospa_vals(1:k,3)), title("GOSPA Missed");
subplot(2,3,6), plot(1:k,gospa_vals(1:k,4)), title("GOSPA False");

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function [qt] = gal2qt(gal)
% Convert volume from US liquid gallons to US liquid quarts.
% Chad Greene 2012
qt = gal*4;


text:
function newnonl( prob, ncones, arg ) %#ok
error( nargchk( 2, 3, nargin ) );
global cvx___
cones = cvx___.cones;
if nargin == 3,
    ncones = struct( 'type', ncones, 'indices', arg );
end
for k = 1 : length( ncones ),
    ncone = ncones(k);
    if any( cvx___.reserved( ncone.indices( : ) ) ),
        error( 'Variables placed in nonlinearities must be free.' );
    else
        cvx___.reserved( ncone.indices ) = 1;
    end
    if isequal( ncone.type, 'nonnegative' ),
        ncone.indices = ncone.indices(:)';
    end
    if isempty( cones ),
        cones = ncone;
    else
        match = find( strcmp( { cones.type }, ncone.type ) );
        if ~isempty( match ),
            nlsiz = size( ncone.indices, 1 );
            match = match( cellfun( 'size', { cones(match).indices }, 1 ) == nlsiz );
            if isempty( match ),
                cones = [ cones, ncone ];
            else
                match = match(1);
                cones(match).indices = [ cones(match).indices, ncone.indices ];
            end
        else
            cones = [ cones, ncone ];
        end
    end
end
cvx___.cones = cones;

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function [Dsubj, nf, r, thdemean, M] = surf_registration_stats(dirname, arg2, arg3)

%
% surf_registration_stats.m
%
% Original Author: Laurence Wastiaux
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:13 $
%    $Revision: 1.3 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

LoadTHMat=1;
D=[];
th=90;
if (nargin==3)
    D=arg2;
    th=arg3;
    LoadTHMat=0;
end
if (nargin==2)
    D=arg2;
    th=0.90;
    LoadTHMat=0;
end
files=dir(dirname);
nf=files;
length(files)
nv = 163842;
%D=zeros(length(files), nv);

%% 1- Load the thickness matrix (N_vertices x N_subjects) %%
if(LoadTHMat==1);
    for i_sub=5:300
        subject_name=files(i_sub).name;
        disp(subject_name)
        SubjectDir=strcat(dirname, '/', subject_name);
        SurfDir=strcat(SubjectDir, '/surf');

        % a- Ressample the thickness curv file to sphere %
        trglh_file=sprintf('~/lh.thickness.average7.mgh');
        trgrh_file=sprintf('~/rh.thickness.average7.mgh');

        if(exist(trglh_file))
            cmdrm=sprintf('rm %s', trglh_file);
            u=unix(cmdrm);
        end
        if(exist(trglh_file))
            cmdrm=sprintf('rm %s', trglh_file);
            u=unix(cmdrm);
        end

        lhthickness_file=strcat(SurfDir, '/lh.thickness');
        rhthickness_file=strcat(SurfDir, '/rh.thickness');
        if (isempty(SurfDir) | ~exist(lhthickness_file) | ~exist(rhthickness_file) )
            fprintf('thickness files do not exist\n');
            i_sub=i_sub+1;
        else
            %cmd_lh=sprintf('mri_surf2surf --srcsubject %s --srcsurfval thickness --src_type curv --trgsubject average7 --trgsurfval %s --trg_type mgh --hemi lh --nsmooth-in 100 --noreshape',subject_name, trglh_file);
            %cmd_rh=sprintf('mri_surf2surf --srcsubject %s --srcsurfval thickness --src_type curv --trgsubject average7 --trgsurfval %s --trg_type mgh --hemi rh --nsmooth-in 100 --noreshape',subject_name, trgrh_file);
            cmd_lh=sprintf('~/work/dev/mris_average_curvature/mris_average_curvature -a 100 -o average7 thickness lh sphere.reg %s %s', subject_name, trglh_file);
            u1=unix(cmd_lh);
            %u2=unix(cmd_rh);

            % b- Load the ressampled surface %
            if (exist(trglh_file)==0)
                fprintf('could not create the ressampled surface\n');
                i_sub=i_sub+1;
            else
                %[vol]=load_mgh(trglh_file);
                vol=read_curv(trglh_file);
                sz=size(vol)
                if(sz(1)~= nv)
                    msgerr=sprintf('Wrong ressampled surface');
                    error(msgerr)
                end
                D= [D ; vol'];
                nf(size(D,1)).name=subject_name;
            end
        end
    end
end
Dsubj=D;
%save('/space/okapi/3/data/laurence/surf_registration/thickness_avg100.mat', 'Dsubj');
%save('/space/okapi/3/data/laurence/surf_registration/Subj_Names.mat', 'nf');
sd=size(D);

nsubjs=sd(1);
nns=1:nsubjs;

%% 2- compute stats %%

% demean
thdemean = Dsubj - repmat(mean(Dsubj,2),[1 nv]);

% Covariance
M = thdemean*thdemean';

% % Jackknife
for jksubj = 1:nsubjs
  fprintf('%d *\n',jksubj);
  indjk = find(nns ~= jksubj);

  % SVD
  [u s v] = fast_svd(thdemean(indjk,:),M(indjk,indjk));
  fcvs(:,jksubj) = cumsum((diag(s.^2))/(sum(diag(s.^2))));

  % keep enough to explain 85% of the data. The fcvs does not
  % change that much regardless of which data set is jked
  nkeep(jksubj) = max(find(fcvs(:,jksubj)<th));
  v = v(:,1:nkeep(jksubj));

  % Extract the jkth subjects thickness and compute var
  thsubj = thdemean(jksubj,:);
  %thsubj = randn(size(thsubj));
  d0(jksubj) = var(thsubj);

  % Project out the first nkeep EVs and recompute var
  thsubjproj = thsubj - (thsubj*v)*v';
  d1(jksubj) = var(thsubjproj);

  fprintf('  nkeep = %3d, d0=%g, d1=%g, r=%g, rwhite = %g\n',...
	  nkeep(jksubj),d0(jksubj),d1(jksubj),d1(jksubj)/d0(jksubj),1-nkeep(jksubj)/size(thdemean,1));

end

% Compute percent residual variance
r = 100*d1./d0;
%save('/space/okapi/3/data/laurence/ADF/surf_registration/residual_var80.mat', 'r');
% Plot the histogram
% indr = find(r<25);
% rr = r(indr);
 [h x] = hist(r);
 p = h/sum(h);
 ph = plot(x,p,'r-');
% set(ph,'Linewidth',3);
% xlabel('Percent Residual Variance');
% title('Distribution of Percent Residual Variance in Thickness')
% set(gcf,'color',[1 1 1]);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef MMF5 < PROBLEM
% <multi> <real> <multimodal>
% Multi-modal multi-objective test function

%------------------------------- Reference --------------------------------
% C. Yue, B. Qu, and J. Liang, A multi-objective particle swarm optimizer
% using ring topology for solving multimodal multiobjective Problems, IEEE
% Transactions on Evolutionary Computation, 2018, 22(5): 805-817.
%------------------------------- 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".
%--------------------------------------------------------------------------

    properties
        POS;    % Pareto optimal set for IGDX calculation
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            obj.D = 2;
            obj.lower    = [1,-1];
            obj.upper    = [3,3];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            temp = X(:,2) <= 1;
            y    = zeros(size(X,1),1);
            y(temp)  = X(temp,2) - sin(6*pi*abs(X(temp,1)-2)+pi);
            y(~temp) = X(~temp,2) - 2 - sin(6*pi*abs(X(~temp,1)-2)+pi);
            PopObj(:,1) = abs(X(:,1)-2);
            PopObj(:,2) = 1 - sqrt(PopObj(:,1)) + 2*y.^2;
        end
        %% Generate Pareto optimal solutions
        function R = GetOptimum(obj,N)
            % Generate points in Pareto optimal set
            obj.POS(:,1) = linspace(1,3,N/2)';
            obj.POS(:,2) = sin(6*pi*abs(obj.POS(:,1)-2)+pi);
            obj.POS = [obj.POS;obj.POS(:,1),obj.POS(:,2)+2];
            % Generate points on Pareto front
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - sqrt(R(:,1));
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R(:,1) = linspace(0,1,100)';
            R(:,2) = 1 - sqrt(R(:,1));
        end
        %% Calculate the metric value
        function score = CalMetric(obj,metName,Population)
            switch metName
                case 'IGDX'
                    score = feval(metName,Population,obj.POS);
                otherwise
                    score = feval(metName,Population,obj.optimum);
            end
        end
        %% Display a population in the objective space
        function DrawObj(obj,Population)
            PopDec = Population.decs;
            temp1  = PopDec(:,1)<=2;
            temp2  = PopDec(:,2)<=1;
            Draw(Population(temp1&temp2).objs,'o','MarkerSize',6,'Marker','o','Markerfacecolor',[1 .5 .5],'Markeredgecolor',[1 .2 .2],{'\it f\rm_1','\it f\rm_2',[]});
            Draw(Population(temp1&~temp2).objs+0.05,'o','MarkerSize',6,'Marker','o','Markerfacecolor',[.5 .5 1],'Markeredgecolor',[.2 .2 1]);
            Draw(Population(~temp1&temp2).objs+0.1,'o','MarkerSize',6,'Marker','o','Markerfacecolor',[.5 1 .5],'Markeredgecolor',[.2 1 .2]);
            Draw(Population(~temp1&~temp2).objs+0.15,'o','MarkerSize',6,'Marker','o','Markerfacecolor',[1 .5 1],'Markeredgecolor',[1 .2 1]);
            Draw(obj.PF,'-','LineWidth',1,'Color',[1 .2 .2]);
            Draw(obj.PF+0.05,'-','LineWidth',1,'Color',[.2 .2 1]);
            Draw(obj.PF+0.1,'-','LineWidth',1,'Color',[.2 1 .2]);
            Draw(obj.PF+0.15,'-','LineWidth',1,'Color',[1 .2 1]);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% sq_dist - a function to compute a matrix of all pairwise squared distances
% between two sets of vectors, stored in the columns of the two matrices, a
% (of size D by n) and b (of size D by m). If only a single argument is given
% or the second matrix is empty, the missing matrix is taken to be identical
% to the first.
%
% Special functionality: If an optional third matrix argument Q is given, it
% must be of size n by m, and in this case a vector of the traces of the
% product of Q' and the coordinatewise squared distances is returned.
%
% NOTE: The program code is written in the C language for efficiency and is
% contained in the file sq_dist.c, and should be compiled using matlabs mex
% facility. However, this file also contains a (less efficient) matlab
% implementation, supplied only as a help to people unfamiliar with mex. If
% the C code has been properly compiled and is avaiable, it automatically
% takes precendence over the matlab code in this file.
%
% Usage: C = sq_dist(a, b)
%    or: C = sq_dist(a)  or equiv.: C = sq_dist(a, [])
%    or: c = sq_dist(a, b, Q)
% where the b matrix may be empty.
%
% where a is of size D by n, b is of size D by m (or empty), C and Q are of
% size n by m and c is of size D by 1.
%
% Copyright (c) 2003, 2004, 2005 and 2006 Carl Edward Rasmussen. 2006-03-09.

function C = kafbox_sq_dist(a, b, Q);

if nargin < 1 | nargin > 3 | nargout > 1
  error('Wrong number of arguments.');
end

if nargin == 1 | isempty(b)                   % input arguments are taken to be
  b = a;                                   % identical if b is missing or empty
end

[D, n] = size(a);
[d, m] = size(b);
if d ~= D
  error('Error: column lengths must agree.');
end

if nargin < 3
  C = zeros(n,m);
  for d = 1:D
    C = C + (repmat(b(d,:), n, 1) - repmat(a(d,:)', 1, m)).^2;
  end
  % C = repmat(sum(a.*a)',1,m)+repmat(sum(b.*b),n,1)-2*a'*b could be used to
  % replace the 3 lines above; it would be faster, but numerically less stable.
else
  if [n m] == size(Q)
    C = zeros(D,1);
    for d = 1:D
      C(d) = sum(sum((repmat(b(d,:), n, 1) - repmat(a(d,:)', 1, m)).^2.*Q));
    end
  else
    error('Third argument has wrong size.');
  end
end

repo: kafbox


text:
function opts = SRSolverParams(varargin)
%SRSOLVERPARAMS Parameters for numerical optizimation.
%   SRSOLVERPARAMS returns a parameter structure used for numerical
%   optimization and computations in different super-resolution algorithms.
%
%   The parameter structure consist of the following parameters:
%       - maxFunEvals:  Maximum number of objective function evaluations.
%       - maxIter:      Maximum number of iterations for non-linear
%                       optimization.
%       - maxIrlsIter:  Maximum number of iterations for iteratively
%                       re-weighted least squares algorithms.
%       - tolX:         Measure of the absolute precision required for the
%                       super-resolved pixel values (termination
%                       tolerance).
%       - tolF:         Measure of the absolute precision required for
%                       objective functions in non-linear optimization
%                       (termination tolerance).
%       - gradCheck:    Set to 1 to check the analytic derived gradients
%                       for debug purposes.
%       - verbose:      Level for debug messages. Set to -1 for no
%                       messages, to 0 for only warning messages and 1 for
%                       error messages.

    % Set default value to parameter structure
    opts = struct('maxFunEvals',         50, ...    % Maximim number of function evals
                  'maxIter',             50, ...    % Maximum number of iterations
                  'maxIrlsIter',         15, ...    % Maximum number of iteration for iteratively re-weighted least squares
                  'tolX',                1e-3, ...  % Tolerance criterion for image pixels
                  'tolF',                1e-3, ...  % Tolerance criterion for objective function value
                  'verbose',             false, ... % Use debug outputs
                  'gradCheck',           false);    % Check objective function gradient

    % Update with user-defined parameters
    for k = 1:2:(nargin - 1)

        param = varargin{k};
        value = varargin{k+1};
        if isfield(opts, param)
            opts = setfield(opts, param, value);
        else
            error( sprintf('Invalid parameter %s', param) );
        end

    end

author: thomas-koehler

repo: SupER


text:
function [err, commonSize, numElements] = statsizechk(nparams,varargin)
%STATSIZECHK Check for compatible array sizes.
%   [ERR,COMMONSIZE,NUMELEMENTS] = STATSIZECHK(NPARAMS,A,B,...,M,N,...) or
%   [ERR,COMMONSIZE,NUMELEMENTS] = STATSIZECHK(NPARAMS,A,B,...,[M,N,...])
%   in effect computes size( A + B + ... + zeros(M,N,...) ), and catches
%   any size mismatches.  NPARAMS is the number of array input arguments.

%   Copyright 1993-2004 The MathWorks, Inc.
%   $Revision: 1.1.6.2 $  $Date: 2004/01/24 09:36:35 $
%
%   Mex file.

% try
%     tmp = 0;
%     for argnum = 1:nparams
%         tmp = tmp + varargin{argnum};
%     end
%     if nargin > nparams+1
%         tmp = tmp + zeros(varargin{nparams+1:end});
%     end
%     err = 0;
%     commonSize = size(tmp);
%     numElements = numel(tmp);
%
% catch
%     err = 1;
%     commonSize = [];
%     numElements = 0;
% end

author: zouchuhang

repo: LayoutNet


text:
% SP_EVALUATE_ELEMENT_LIST_PARAM: compute the basis functions, in the parametric domain, in a given list of elements.
%
%     sp = sp_evaluate_element_list_param (space, msh_elems, 'option1', value1, ...)
%
% INPUTS:
%
%    space:   object defining the space of discrete functions (see sp_scalar)
%    msh_elems: msh structure containing the information of quadrature or
%               visualization points, for a given list of elements
%               (see msh_cartesian/msh_evaluate_element_list)
%   'option', value: additional optional parameters, currently available options are:
%
%              Name     |   Default value |  Meaning
%           ------------+-----------------+----------------------------------
%            value      |      true       |  compute shape_functions
%            gradient   |      false      |  compute shape_function_gradients
%            hessian    |      false      |  compute shape_function_hessians
%
% OUTPUT:
%
%    sp: struct representing the discrete function space, with the following fields:
%              (see the article for a detailed description)
%
%    FIELD_NAME      (SIZE)                                        DESCRIPTION
%    ncomp           (scalar)                                      number of components of the functions of the space (actually, 1)
%    ndof            (scalar)                                      total number of degrees of freedom
%    ndof_dir        (1 x ndim vector)                             degrees of freedom along each direction
%    nsh_max         (scalar)                                      maximum number of shape functions per element
%    nsh             (1 x msh_elems.nel vector)                    actual number of shape functions per each element
%    connectivity    (nsh_max x msh_elems.nel vector)              indices of basis functions that do not vanish in each element
%    shape_functions (msh_elems.nqn x nsh_max x msh_elems.nel)     basis functions evaluated at each quadrature node in each element
%    shape_function_gradients
%         (ndim x msh_elems.nqn x nsh_max x msh_elems.nel)         basis function gradients evaluated at each quadrature node in each element
%    shape_function_hessians
%         (ndim x ndim x msh_elems.nqn x nsh_max x msh_elems.nel)  basis function hessians evaluated at each quadrature node in each element
%
% Copyright (C) 2009, 2010, 2011 Carlo de Falco
% Copyright (C) 2011, 2015, 2019 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function sp = sp_evaluate_element_list_param (space, msh, varargin)

value = true;
gradient = false;
hessian = false;
if (~isempty (varargin))
  if (~rem (length (varargin), 2) == 0)
    error ('sp_evaluate_element_list_param: options must be passed in the [option, value] format');
  end
  for ii=1:2:length(varargin)-1
    if (strcmpi (varargin {ii}, 'value'))
      value = varargin {ii+1};
    elseif (strcmpi (varargin {ii}, 'gradient'))
      gradient = varargin {ii+1};
    elseif (strcmpi (varargin {ii}, 'hessian'))
      hessian = varargin {ii+1};
    else
      warning ('Ignoring unknown option %s', varargin {ii});
    end
  end
end

sp_univ = space.sp_univ;

elem_list = msh.elem_list;

elem_ind = cell (msh.ndim, 1);
[elem_ind{:}] = ind2sub (msh.nel_dir, elem_list);
elem_ind = cell2mat (elem_ind);

nsh = 1;
for idim = 1:msh.ndim
  nsh = nsh .* sp_univ(idim).nsh(elem_ind(idim,:));
end

for idim = 1:msh.ndim
  csize = ones (1, msh.ndim);
  csize(idim) = sp_univ(idim).nsh_max;
  crep = [sp_univ.nsh_max];
  crep(idim) = 1;
  conn{idim} = reshape (sp_univ(idim).connectivity(:,elem_ind(idim,:)), [csize, msh.nel]);
  conn{idim} = repmat (conn{idim}, [crep, 1]);
end

connectivity = zeros (space.nsh_max, msh.nel);
indices = ones (size (conn{1}));
for idim = 1:msh.ndim
  indices = indices & conn{idim} ~= 0;
end
% for idim = 1:msh.ndim
%   conn{idim} = conn{idim}(indices);
% end
connectivity(indices) = sub2ind ([space.ndof_dir, 1], conn{:}); % The extra 1 makes things work in any dimension
connectivity = reshape (connectivity, space.nsh_max, msh.nel);

clear conn csize crep indices

sp = struct('nsh_max', space.nsh_max, 'nsh', nsh, 'ndof', space.ndof,  ...
            'ndof_dir', space.ndof_dir, 'connectivity', connectivity, ...
            'ncomp', 1, 'degree', space.degree);

if (value || gradient || hessian)
shp = cell(1,msh.ndim); shg = cell(1,msh.ndim); shh = cell(1,msh.ndim);
  for idim = 1:msh.ndim
    ssize = ones (1, 2*msh.ndim);
    ssize([idim, msh.ndim+idim]) = [msh.nqn_dir(idim), sp_univ(idim).nsh_max];
    srep = [msh.nqn_dir, sp_univ.nsh_max];
    srep([idim, msh.ndim+idim]) = 1;
    shp{idim} = reshape (sp_univ(idim).shape_functions(:,:,elem_ind(idim,:)), [ssize, msh.nel]);
    shp{idim} = repmat (shp{idim}, [srep, 1]);
    shp{idim} = reshape (shp{idim}, msh.nqn, space.nsh_max, msh.nel);
    shg{idim} = reshape (sp_univ(idim).shape_function_gradients(:,:,elem_ind(idim,:)), [ssize, msh.nel]);
    shg{idim} = repmat (shg{idim}, [srep, 1]);
    shg{idim} = reshape (shg{idim}, msh.nqn, space.nsh_max, msh.nel);
    shh{idim} = reshape (sp_univ(idim).shape_function_hessians(:,:,elem_ind(idim,:)), [ssize, msh.nel]);
    shh{idim} = repmat (shh{idim}, [srep, 1]);
    shh{idim} = reshape (shh{idim}, msh.nqn, space.nsh_max, msh.nel);
  end

  if (value)
    sp.shape_functions = 1;
    for idim = 1:msh.ndim
      sp.shape_functions = sp.shape_functions .* shp{idim};
    end
  end

  if (gradient)
    for idim = 1:msh.ndim
      shape_fun_grad = shg{idim};
      for jdim = setdiff (1:msh.ndim, idim)
        shape_fun_grad = shape_fun_grad .* shp{jdim};
      end
      sp.shape_function_gradients(idim,:,:,:) = shape_fun_grad;
    end
    sp.shape_function_gradients = reshape (sp.shape_function_gradients, ...
                                  msh.ndim, msh.nqn, sp.nsh_max, msh.nel);
  end

  if (hessian && (isfield (msh, 'geo_map_der2') || msh.nel == 0))
    for idim = 1:msh.ndim
      shape_fun_hess = shh{idim};
      for jdim = setdiff (1:msh.ndim, idim)
        shape_fun_hess = shape_fun_hess .* shp{jdim};
      end
      sp.shape_function_hessians(idim,idim,:,:,:) = shape_fun_hess;

      for jdim = setdiff (1:msh.ndim, idim)
        shape_fun_hess = shg{idim} .* shg{jdim};
        for kdim = setdiff (1:msh.ndim, [idim, jdim])
          shape_fun_hess = shape_fun_hess .* shp{kdim};
        end
        sp.shape_function_hessians(idim,jdim,:,:,:) = shape_fun_hess;
      end
    end
  end

  clear shp shg shh

  if (strcmpi (space.space_type, 'NURBS'))
    sp = bsp_2_nrb__ (sp, msh, space.weights);
  end
end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
%GETSIZE_IN Get size_in field of mapping
%
%    SIZE_IN = GETSIZE_IN(W)

author: marianux

repo: ecg-kit

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


text:
function [A, C, Q, R, initx, initV, LL] = ...
    learn_kalman(data, A, C, Q, R, initx, initV, max_iter, diagQ, diagR, ARmode, constr_fun, varargin)
% LEARN_KALMAN Find the ML parameters of a stochastic Linear Dynamical System using EM.
%
% [A, C, Q, R, INITX, INITV, LL] = LEARN_KALMAN(DATA, A0, C0, Q0, R0, INITX0, INITV0) fits
% the parameters which are defined as follows
%   x(t+1) = A*x(t) + w(t),  w ~ N(0, Q),  x(0) ~ N(init_x, init_V)
%   y(t)   = C*x(t) + v(t),  v ~ N(0, R)
% A0 is the initial value, A is the final value, etc.
% DATA(:,t,l) is the observation vector at time t for sequence l. If the sequences are of
% different lengths, you can pass in a cell array, so DATA{l} is an O*T matrix.
% LL is the "learning curve": a vector of the log lik. values at each iteration.
% LL might go positive, since prob. densities can exceed 1, although this probably
% indicates that something has gone wrong e.g., a variance has collapsed to 0.
%
% There are several optional arguments, that should be passed in the following order.
% LEARN_KALMAN(DATA, A0, C0, Q0, R0, INITX0, INITV0, MAX_ITER, DIAGQ, DIAGR, ARmode)
% MAX_ITER specifies the maximum number of EM iterations (default 10).
% DIAGQ=1 specifies that the Q matrix should be diagonal. (Default 0).
% DIAGR=1 specifies that the R matrix should also be diagonal. (Default 0).
% ARMODE=1 specifies that C=I, R=0. i.e., a Gauss-Markov process. (Default 0).
% This problem has a global MLE. Hence the initial parameter values are not important.
%
% LEARN_KALMAN(DATA, A0, C0, Q0, R0, INITX0, INITV0, MAX_ITER, DIAGQ, DIAGR, F, P1, P2, ...)
% calls [A,C,Q,R,initx,initV] = f(A,C,Q,R,initx,initV,P1,P2,...) after every M step. f can be
% used to enforce any constraints on the params.
%
% For details, see
% - Ghahramani and Hinton, "Parameter Estimation for LDS", U. Toronto tech. report, 1996
% - Digalakis, Rohlicek and Ostendorf, "ML Estimation of a stochastic linear system with the EM
%      algorithm and its application to speech recognition",
%       IEEE Trans. Speech and Audio Proc., 1(4):431--442, 1993.

%    learn_kalman(data, A, C, Q, R, initx, initV, max_iter, diagQ, diagR, ARmode, constr_fun, varargin)
if nargin < 8, max_iter = 10; end
if nargin < 9, diagQ = 0; end
if nargin < 10, diagR = 0; end
if nargin < 11, ARmode = 0; end
if nargin < 12, constr_fun = []; end
verbose = 1;
thresh = 1e-4;

if ~iscell(data)
  N = size(data, 3);
  data = num2cell(data, [1 2]); % each elt of the 3rd dim gets its own cell
else
  N = length(data);
end

N = length(data);
ss = size(A, 1);
os = size(C,1);

alpha = zeros(os, os);
Tsum = 0;
for ex = 1:N
  %y = data(:,:,ex);
  y = data{ex};
  T = length(y);
  Tsum = Tsum + T;
  alpha_temp = zeros(os, os);
  for t=1:T
    alpha_temp = alpha_temp + y(:,t)*y(:,t)';
  end
  alpha = alpha + alpha_temp;
end

previous_loglik = -inf;
loglik = 0;
converged = 0;
num_iter = 1;
LL = [];

% Convert to inline function as needed.
if ~isempty(constr_fun)
  constr_fun = fcnchk(constr_fun,length(varargin));
end

while ~converged & (num_iter <= max_iter)

  %%% E step

  delta = zeros(os, ss);
  gamma = zeros(ss, ss);
  gamma1 = zeros(ss, ss);
  gamma2 = zeros(ss, ss);
  beta = zeros(ss, ss);
  P1sum = zeros(ss, ss);
  x1sum = zeros(ss, 1);
  loglik = 0;

  for ex = 1:N
    y = data{ex};
    T = length(y);
    [beta_t, gamma_t, delta_t, gamma1_t, gamma2_t, x1, V1, loglik_t] = ...
	Estep(y, A, C, Q, R, initx, initV, ARmode);
    beta = beta + beta_t;
    gamma = gamma + gamma_t;
    delta = delta + delta_t;
    gamma1 = gamma1 + gamma1_t;
    gamma2 = gamma2 + gamma2_t;
    P1sum = P1sum + V1 + x1*x1';
    x1sum = x1sum + x1;
    %fprintf(1, 'example %d, ll/T %5.3f\n', ex, loglik_t/T);
    loglik = loglik + loglik_t;
  end
  LL = [LL loglik];
  if verbose, fprintf(1, 'iteration %d, loglik = %f\n', num_iter, loglik); end
  %fprintf(1, 'iteration %d, loglik/NT = %f\n', num_iter, loglik/Tsum);
  num_iter =  num_iter + 1;

  %%% M step

  % Tsum =  N*T
  % Tsum1 = N*(T-1);
  Tsum1 = Tsum - N;
  A = beta * inv(gamma1);
  %A = (gamma1' \ beta')';
  Q = (gamma2 - A*beta') / Tsum1;
  if diagQ
    Q = diag(diag(Q));
  end
  if ~ARmode
    C = delta * inv(gamma);
    %C = (gamma' \ delta')';
    R = (alpha - C*delta') / Tsum;
    if diagR
      R = diag(diag(R));
    end
  end
  initx = x1sum / N;
  initV = P1sum/N - initx*initx';

  if ~isempty(constr_fun)
    [A,C,Q,R,initx,initV] = feval(constr_fun, A, C, Q, R, initx, initV, varargin{:});
  end

  converged = em_converged(loglik, previous_loglik, thresh);
  previous_loglik = loglik;
end

%%%%%%%%%

function [beta, gamma, delta, gamma1, gamma2, x1, V1, loglik] = ...
    Estep(y, A, C, Q, R, initx, initV, ARmode)
%
% Compute the (expected) sufficient statistics for a single Kalman filter sequence.
%

[os T] = size(y);
ss = length(A);

if ARmode
  xsmooth = y;
  Vsmooth = zeros(ss, ss, T); % no uncertainty about the hidden states
  VVsmooth = zeros(ss, ss, T);
  loglik = 0;
else
  [xsmooth, Vsmooth, VVsmooth, loglik] = kalman_smoother(y, A, C, Q, R, initx, initV);
end

delta = zeros(os, ss);
gamma = zeros(ss, ss);
beta = zeros(ss, ss);
for t=1:T
  delta = delta + y(:,t)*xsmooth(:,t)';
  gamma = gamma + xsmooth(:,t)*xsmooth(:,t)' + Vsmooth(:,:,t);
  if t>1 beta = beta + xsmooth(:,t)*xsmooth(:,t-1)' + VVsmooth(:,:,t); end
end
gamma1 = gamma - xsmooth(:,T)*xsmooth(:,T)' - Vsmooth(:,:,T);
gamma2 = gamma - xsmooth(:,1)*xsmooth(:,1)' - Vsmooth(:,:,1);

x1 = xsmooth(:,1);
V1 = Vsmooth(:,:,1);


text:
function [ImTrans,tau,ImTrans_color] = align_c(I2,I1,tau,iter)

c=size(I2,3);
ImTrans = I2;
if nargin<3
  tau = zeros(6,1);
end
if nargin<4
  iter = 5;
end

numLevel = fix(log(size(I2,1)*size(I2,2))/log(2)/2-3)-1;

% I1=mean(Imref,3);
% I2=mean(Im,3);
% I2=luminance_transfer(I1,I2);
[ImTrans,ImTrans_color,tau] = regMGNC_c(I1,I2,tau,numLevel,iter);

end

author: csjcai

repo: RealSR

save_path: github-repos/MATLAB/csjcai-RealSR