text:
function VO = spm_smoothto8bit(V,fwhm)
% 3 dimensional convolution of an image to 8bit data in memory
% FORMAT VO = spm_smoothto8bit(V,fwhm)
% V     - mapped image to be smoothed
% fwhm  - FWHM of Guassian filter width in mm
% VO    - smoothed volume in a form that can be used by the
%         spm_*_vol.mex* functions.
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id$

if nargin>1 & fwhm>0,
    VO = smoothto8bit(V,fwhm);
else,
    VO = V;
end;
return;
%_______________________________________________________________________

%_______________________________________________________________________
function VO = smoothto8bit(V,fwhm)
% 3 dimensional convolution of an image to 8bit data in memory
% FORMAT VO = smoothto8bit(V,fwhm)
% V     - mapped image to be smoothed
% fwhm  - FWHM of Guassian filter width in mm
% VO    - smoothed volume in a form that can be used by the
%         spm_*_vol.mex* functions.
%_______________________________________________________________________

vx   = sqrt(sum(V.mat(1:3,1:3).^2));
s    = (fwhm./vx./sqrt(8*log(2)) + eps).^2;
r    = cell(1,3);
for i=1:3,
    r{i}.s = ceil(3.5*sqrt(s(i)));
    x      = -r{i}.s:r{i}.s;
    r{i}.k = exp(-0.5 * (x.*x)/s(i))/sqrt(2*pi*s(i));
    r{i}.k = r{i}.k/sum(r{i}.k);
end;

buff = zeros([V.dim(1:2) r{3}.s*2+1]);

VO        = V;
VO.dt     = [spm_type('uint8') spm_platform('bigend')];
V0.dat    = uint8(0);
V0.dat(VO.dim(1:3)) = uint8(0);
VO.pinfo  = [];

for i=1:V.dim(3)+r{3}.s,
    if i<=V.dim(3),
        img      = spm_slice_vol(V,spm_matrix([0 0 i]),V.dim(1:2),0);
        msk      = find(~isfinite(img));
        img(msk) = 0;
        buff(:,:,rem(i-1,r{3}.s*2+1)+1) = ...
            conv2(conv2(img,r{1}.k,'same'),r{2}.k','same');
    else,
        buff(:,:,rem(i-1,r{3}.s*2+1)+1) = 0;
    end;

    if i>r{3}.s,
        kern    = zeros(size(r{3}.k'));
        kern(rem((i:(i+r{3}.s*2))',r{3}.s*2+1)+1) = r{3}.k';
        img     = reshape(buff,[prod(V.dim(1:2)) r{3}.s*2+1])*kern;
        img     = reshape(img,V.dim(1:2));
        ii      = i-r{3}.s;
        mx      = max(img(:));
        mn      = min(img(:));
        if mx==mn, mx=mn+eps; end;
        VO.pinfo(1:2,ii) = [(mx-mn)/255 mn]';
        VO.dat(:,:,ii)   = uint8(round((img-mn)*(255/(mx-mn))));
    end;
end;


text:
function [R,P,K]=pfrac(G)
G=tf(G);
[R,P,K]=residue(G.num{1},G.den{1});
for i=1:length(R),
   if imag(P(i))>eps
      a=real(R(i)); b=imag(R(i));
      R(i)=-2*sqrt(a^2+b^2);
      R(i+1)=atan2(-a,b);
   end
end

author: Sable

repo: mcbench-benchmarks


text:
% eegplugin_brainmovie() - EEGLAB plugin for importing ASC INStep files.
%
% Usage:
%   >> eegplugin_brainmovie(menu);
%   >> eegplugin_brainmovie(menu, trystrs, catchstrs);
%
% Inputs:
%   menu       - [float]  EEGLAB menu handle
%   trystrs    - [struct] "try" strings for menu callbacks. See notes on EEGLab plugins.
%                (http://www.sccn.ucsd.edu/eeglab/contrib.html)
%   catchstrs  - [struct] "catch" strings for menu callbacks. See notes on EEGLab plugins.
%                (http://www.sccn.ucsd.edu/eeglab/contrib.html)
%
% Authors: Arnaud Delorme, SCCN, INC, UCSD, 2001-

% Copyright (C) 2001 Arnaud Delorme
%
% 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

function vers = eegplugin_brainmovie(fig, trystrs, catchstrs)

    vers = 'brainmovie0.1b';

    if nargin < 3
        error('eegplugin_brainmovie requires 3 arguments');
    end;

    % add folder to path
    % ------------------
    if ~exist('brainmovie')
        p = which('eegplugin_brainmovie.m');
        p = p(1:findstr(p,'eegplugin_brainmovie.m')-1);
        addpath( p );
    end;

author: PatternRecognition

repo: OpenBMI


text:
function safeSave(fname, varargin)

% SAFESAVE Safe save
% FORMAT
% DESC Saves given variables in such a way that existing files are
%   only overwritten with a complete save, thus making sure that a
%   complete copy exists even if the call is interrupted
% ARG fname : name of the file to save to (should include extension)
% ARG ... : the variables to save (either as names (strings) or as
%   variables
%
% SEEALSO save
%
% COPYRIGHT : Antti Honkela, 2010

% NDLUTIL

  saveargs = cell(size(varargin));
  for k=1:length(varargin),
    if isempty(inputname(k+1)),
      saveargs{k} = varargin{k};
    else
      saveargs{k} = inputname(k+1);
    end
  end
  argstring = sprintf(', ''%s''', saveargs{:});

  % First time saving to this file
  if ~exist(fname, 'file'),
    callstr = sprintf('save(''%s''%s)', fname, argstring);
    evalin('caller', callstr);
    return;
  end

  % If file exists, first write to a different name
  [pathstr, name, ext] = fileparts(fname);
  tmpfname = fullfile(pathstr,[name '_tmp' ext]);
  callstr = sprintf('save(''%s''%s)', tmpfname, argstring);
  evalin('caller', callstr);

  % Then move to overwrite the original
  movefile(tmpfname, fname);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


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

%%
clear; close all; clc;

%% Syntax
% |[T]=writeBasicSmesh(smeshStruct);|

%% 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 [errmsg,x,values] = dfevaluate(fitNames,x,fun,wantBounds,confLevel,plotFun,dum)
%DFEVALUATE Evaluate fits for DFITTOOL

%   $Revision: 1.1.6.3 $  $Date: 2004/01/24 09:35:32 $
%   Copyright 1993-2004 The MathWorks, Inc.

% If the function is th empty string, clear the plot (if there is one) and
% clear any saved data.
if isempty(fitNames)
    dfevaluateplot(false,false); % closes the plot window
    dfgetset('evaluateResults', []);
    return
end

nfits = length(fitNames);
try
    x = sprintf('[ %s ]', x); % allow an unbracketed list of numbers to work
    x = eval(x);
    confLevel = eval(confLevel) ./ 100;
catch
    x = [];
    values = zeros(0, nfits*(1+2*wantBounds));
    errmsg = sprintf('Invalid MATLAB expression: %s',lasterr);
    h = [];
    return
end
errmsg = '';

x = x(:);
n = length(x);

% % This is enforced by the evaluate panel.
% switch fun
% case {'pdf' 'hazrate' 'condmean'}
%     % No bounds allowed for pdf, hazrate, or conditional mean.
%     wantBounds = false;
% otherwise % {'cdf' 'icdf' 'survivor' 'cumhazard' 'probplot'}
%     % bounds are allowed
% end

% Output table will have first column for data, then for each fit, one
% column for function, two columns for bounds (if requested).
values = repmat(NaN, n, nfits*(1+2*wantBounds));

fitdb = getfitdb;
for i = 1:nfits
    fit = find(fitdb, 'name', fitNames{i});
    % Cannot compute bounds for kernel smooths and certain parametric fits.
    getBounds = wantBounds && ...
        (~isequal(fit.fittype, 'smooth') && fit.distspec.hasconfbounds);

    % Evaluate the requested function for this fit.
    if getBounds
        [y,ylo,yup] = eval(fit,x,fun,confLevel);
        values(:,3*i-2) = y;
        values(:,3*i-1) = ylo;
        values(:,3*i) = yup;
    else
        y = eval(fit,x,fun);
        if wantBounds
            values(:,3*i-2) = y;
        else
            values(:,i) = y;
        end
    end
end

% Save the results for SaveToWorkSpace.  The plot function can be called
% directly from java, so it uses those saved results as well.
dfgetset('evaluateResults', [x,values]);

% Save information about the fits that we've evaluated for the plot function.
dfgetset('evaluateFun', fun);
dfgetset('evaluateInfo', struct('fitNames',{fitNames}, 'wantBounds',wantBounds));

dfevaluateplot(plotFun,dum);

author: zouchuhang

repo: LayoutNet


text:
function [RV,IM,J,IMF] = remove_unreferenced(V,F)
  % REMOVE_UNREFERENCED Removes any rows in V that are not referenced in R.
  % **Warning:** This implementation is O(#V) regardless of #F
  %
  % [RV,IM,J,IMF] = remove_unreferenced(V,F)
  %
  % Inputs:
  %   V  #V by dim list of "vertex positions"
  %   F  #F by anything list of indices into V (will be treated as F(:))
  % Outputs:
  %  RV  #unique(F) by dim vertex positions, order such that if the jth vertex is
  %    some face in F, and the kth vertex is not then j comes before k
  %  IM  #V by 1 list of indices such that: RF = IM(F) and RT = IM(T)
  %    and V(find(IM~=-1),:) = RV
  %  J  #RV by 1 list, such that RV = V(J,:)
  %
  % Examples:
  %   % Tet mesh in (V,T,F)
  %   [RV,I] = remove_unreferenced(V,[T(:);F(:)]);
  %   T = I(T);
  %   F = I(F);
  %   ... % do some computation on RV
  %   % replace back into V
  %   V(find(IM<=size(SV,1)),:) = V
  %

  if isempty(F)
    RV = zeros(0,size(V,2));
    IM = (1:size(V,1))';
    J = zeros(1,0);
    IMF = zeros(1,0);
    return;
  end
  % get list of unique vertex indices that occur in faces
  %U = unique(F(:));
  % Slightly faster unique if we don't have infs or nans
  sF = sort(F(:));
  I = [true;diff(sF)~=0];
  U = sF(I);
  % get list of vertices that do not occur in faces
  %NU = (1:size(V,1))';
  %NU = NU(~ismember(NU,U));
  n = size(V,1);
  NU = find(0==sparse(U,1,1,n,1));
  assert((size(U,1) + size(NU,1)) == n);
  % allocate space for an indexmap so that IM[i] gives new index of vertex i
  IM = zeros(n,1);
  % reindex vertices that occur in faces to be first
  IM(U) = 1:size(U,1);
  % reindex vertices that do not occur in faces to come after those that do
  IM(NU) = size(U,1) + (1:size(NU,1));
  % reorder vertices
  RV(IM,:) = V;
  % Remove unreferenced
  RV = RV(1:max(IM(F(:))),:);
  J(IM) = 1:n;
  IMF = IM(F);
  J = J(1:max(IMF(:)));

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function f = isrowvector(m)

% function f = isrowvector(m)
%
% <m> is a matrix
%
% return whether <m> is 1 x n where n >= 0.
% specifically:
%   f = isvector(m) & size(m,1)==1;
%
% example:
% isrowvector([1 2])
% isrowvector([1])
% isrowvector(zeros(1,0))
% ~isrowvector([])

f = isvector(m) & size(m,1)==1;

author: cvnlab

repo: GLMsingle

save_path: github-repos/MATLAB/cvnlab-GLMsingle


text:
function evaluateAP(predidxs)
% implementation of AP measure,
% as defined in [Pishchulin et al., arXiv'15]

fprintf('evaluateAP()\n');

thresh = 0.5;

tableTex = cell(length(predidxs)+1,1);

% load ground truth
p = getExpParams(-1);
load([p.gtDir '/annolist_dataset_v12'],'annolist');
load([p.gtDir '/mpii_human_pose_v1_u12'],'RELEASE');
% load groups of people
load('./groups_v12','groups');
bTrain = false;
% obtain testing multi-person images and rectangles
[imgidxs_multi_test,rectidxs_multi_test] = getMultiPersonGroups(groups,RELEASE, bTrain);
% all testing images
annolist_test = annolist(RELEASE.img_train == 0);
% multi-person images
annolist_test_multi = annolist_test(imgidxs_multi_test);
% multi-person rectangles
for imgidx = 1:length(annolist_test_multi)
    annolist_test_multi(imgidx).annorect = annolist_test_multi(imgidx).annorect(rectidxs_multi_test{imgidx});
end

apAll = zeros(15,length(predidxs));
markerAll = repmat({'-'},length(predidxs),1);

tableDir = [fileparts(p.predFilename) '/latex']; if (~exist(tableDir,'dir')), mkdir(tableDir); end

for i = 1:length(predidxs);

    p = getExpParams(predidxs(i));
    % load predictions
    load([fileparts(p.predFilename) '/pred_keypoints_mpii_multi' ],'pred');
    assert(length(annolist_test_multi) == length(pred));
    % assign predicted poses to GT poses
    [scoresAll, labelsAll, nGTall] = assignGTmulti(pred,annolist_test_multi,thresh(end));
    % compute average precision (AP) per part
    ap = zeros(size(nGTall,1)+1,1);
    for j = 1:size(nGTall,1)
      scores = []; labels = [];
      for imgidx = 1:length(annolist_test_multi)
        scores = [scores; scoresAll{j}{imgidx}];
        labels = [labels; labelsAll{j}{imgidx}];
      end
      % compute precision/recall
      [precision,recall] = getRPC(scores,labels,sum(nGTall(j,:)));
      % compute AP
      ap(j) = VOCap(recall,precision)*100;
    end
    % compute mean AP
    ap(end) = mean(ap(1:end-1));
    columnNames = p.partNames;
    save([fileparts(p.predFilename) '/apAll'],'ap','columnNames');
    % plot results
    [row, header] = genTableAP(ap,p.name);
    tableTex{1} = header;
    tableTex{i+1} = row;
    apAll(:,i) = ap;
end

fid = fopen([tableDir '/ap.tex'],'wt');assert(fid ~= -1);
for i=1:length(tableTex),fprintf(fid,'%s',tableTex{i}); end; fclose(fid);

end

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
classdef tris

methods(Static)

function x = forward_row(L,b)
    % Solvex L x = b problem by forward substitution row wise
    % b gets overwritten in the process
    % GVL4: algorithm 3.1.1
    n = length(b);
    % start from last row
    for i=1:n
        b(i) = (b(i) - L(i, 1:i-1)* b(1:i-1)) / L(i,i);
    end
    x = b;
end % function

function x = back_row(U,b)
    % Solvex U x = b problem by back substitution row wise
    % b gets overwritten in the process
    % GVL4: algorithm 3.1.2
    n = length(b);
    % start from last row
    for i=n:-1:1
        b(i) = (b(i) - U(i, i+1:n) * b(i+1:n)) / U(i,i);
    end
    x = b;
end % function

function x = forward_col(L,b)
    % Solvex L x = b problem by forward substitution column wise
    % b gets overwritten in the process
    % GVL4: algorithm 3.1.3
    n = length(b);
    % start from last row
    for j=1:n
        b(j) = b(j) / L(j, j);
        b(j+1:n) = b(j+1:n) - L(j+1:n, j)*b(j);
    end
    x = b;
end % function

function x =  back_col(U, b)
    % Solves U x = b problem by back substitution column wise
    % b gets overwritten during the process
    % GVL4: algorithm 3.1.4
    n = length(b);
    % we iterate backwards from last column onwards
    for j=n:-1:1
        % Scale with the diagonal entry in this column
        b(j) = b(j)/U(j,j);
        % subtract rest of the column
        b(1:j-1) = b(1:j-1) - U(1:j-1,j)*b(j);
    end
    x = b;
end % function

end % methods

end % classdef

author: indigits

repo: sparse-plex


text:
%isiscalar - Test if parameter is a scalar (integer) satisfying an optional list of tests.
%
%  USAGE
%
%    test = isiscalar(x,test1,test2,...)
%
%    x              parameter to test
%    test1...       optional list of additional tests
%
%  EXAMPLES
%
%    % Test if x is a scalar (double)
%    isiscalar(x)
%
%    % Test if x is a strictly positive scalar (double)
%    isiscalar(x,'>0')
%
%    % Test if x is a scalar (double) included in [2,3]
%    isiscalar(x,'>=2','<=3')
%
%  NOTE
%
%    The tests ignore NaN, e.g. isiscalar(nan), isiscalar(nan,'>0') and isiscalar(nan,'<=0')
%    all return 1.
%
%  SEE ALSO
%
%    See also isdmatrix, isdvector, isdscalar, isimatrix, isivector, isstring_FMAT,
%    islscalar, islvector, islmatrix.
%

% Copyright (C) 2010 by Michaël Zugaro
%
% 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.

function test = isiscalar(x,varargin)

% Check number of parameters
if nargin < 1,
  error('Incorrect number of parameters (type ''help isiscalar'' for details).');
end

% Test: double, scalar
test = isa(x,'double') & isscalar(x);
if ~test, return; end

% Test: integers?
test = test & round(x)==x;

% Optional tests
for i = 1:length(varargin),
	try
		if ~eval(['x' varargin{i} ';']), test = false; return; end
	catch err
		error(['Incorrect test ''' varargin{i} ''' (type ''help isiscalar'' for details).']);
	end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function e = end(X,k,n)
%END Last index of indexing expression for ttensor.
%
%   The expression X(end,:,:) will call END(X,1,3) to determine
%   the value of the first index.
%
%   See also TTENSOR, TTENSOR/SUBSREF, END.
%
%MATLAB Tensor Toolbox.
%Copyright 2010, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by Brett Bader and Tamara Kolda.
% http://csmr.ca.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2010) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in tensor_toolbox/LICENSE.txt
% $Id: end.m,v 1.7 2010/03/19 23:46:31 tgkolda Exp $

% Note that this only works with {} because () is not supported by
% subsref.
e = ndims(X);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [uniques,numUnique] = countUnique(x,option)
%COUNT_UNIQUE  Determines unique values, and counts occurrences
%   [uniques,numUnique] = count_unique(x)
%
%   This function determines unique values of an array, and also counts the
%   number of instances of those values.
%
%   This uses the MATLAB builtin function accumarray, and is faster than
%   MATLAB's unique function for intermediate to large sizes of arrays for integer values.
%   Unlike 'unique' it cannot be used to determine if rows are unique or
%   operate on cell arrays.
%
%   If float values are passed, it uses MATLAB's logic builtin unique function to
%   determine unique values, and then to count instances.
%
%   Descriptions of Input Variables:
%   x:  Input vector or matrix, N-D.  Must be a type acceptable to
%       accumarray, numeric, logical, char, scalar, or cell array of
%       strings.
%   option: Acceptable values currently only 'float'.  If 'float' is
%           specified, the input x vector will be treated as containing
%           decimal values, regardless of whether it is a float array type.
%
%   Descriptions of Output Variables:
%   uniques:    sorted unique values
%   numUnique:  number of instances of each unique value
%
%   Example(s):
%   >> [uniques] = countUnique(largeArray);
%   >> [uniques,numUnique] = countUnique(largeArray);
%
%   See also: unique, accumarray

% Author: Anthony Kendall
% Contact: anthony [dot] kendall [at] gmail [dot] com, http://www.mathworks.com/matlabcentral/fileexchange/23333
% Created: 2009-03-17

if isempty(x)
    uniques = x;
    numUnique = x;
    return;
elseif isscalar(x)
    uniques = x;
    numUnique = 1;
    return;
end

testFloat = false;
if nargin == 2 && strcmpi(option,'float')
    testFloat = true;
end

nOut = nargout;
if testFloat
    if nOut < 2
        [uniques] = float_cell_unique(x,nOut);
    else
        [uniques,numUnique] = float_cell_unique(x,nOut);
    end
else
    try %this will fail if the array is float or cell
        if nOut < 2
            [uniques] = int_log_unique(x,nOut);
        else
            [uniques,numUnique] = int_log_unique(x,nOut);
        end
    catch %default to standard approach
        if nOut < 2
            [uniques] = float_cell_unique(x,nOut);
        else
            [uniques,numUnique] = float_cell_unique(x,nOut);
        end
    end
end

end

function [uniques,numUnique] = int_log_unique(x,nOut)
%Check to see if accumarray is appropriate for this function
maxVal = max(x(:));
if maxVal / numel(x) > 1000
    throw(MException('countUnique','Accumarray is inefficient for arrays when ind values are >> than the number of elements'))
end
%First, determine the offset for negative values
minVal = min(x(:));
if minVal < 1
    %Now, offset to get the index
    index = x(:) - minVal + 1;

    %Get the number of duplicates with accumarray
    numUnique = accumarray(index,1);

    %Get the sum of those duplicate values
    sumDups = accumarray(index,x(:));
else
    %Get the number of duplicates with accumarray
    numUnique = accumarray(x(:),1);

    %Get the sum of those duplicate values
    sumDups = accumarray(x(:),x(:));
end

%Find numUnique > 0
test = (numUnique > 0);

%Determine the unique values
uniques = sumDups(test) ./ (numUnique(test));

if nOut == 2
    %Trim the numUnique array
    numUnique = numUnique(test);
end
end

function [uniques,numUnique] = float_cell_unique(x,nOut)

if ~iscell(x)
    %First, sort the input vector
    x = sort(x(:));
    numelX = numel(x);

    %Check to see if the array type needs to be converted to double
    currClass = class(x);
    isdouble = strcmp(currClass,'double');

    if ~isdouble
        x = double(x);
    end

    %Check to see if there are any NaNs or Infs, sort returns these either at
    %the beginning or end of an array
    if isnan(x(1)) || isinf(x(1)) || isnan(x(numelX)) || isinf(x(numelX))
        %Check to see if the array contains nans or infs
        xnan = isnan(x);
        xinf = isinf(x);
        testRep = xnan | xinf;

        %Remove all of these from the array
        x = x(~testRep);
    end

    %Determine break locations of unique values
    uniqueLocs = [true;diff(x) ~= 0];
else
    isdouble = true; %just to avoid conversion on finish

    %Sort the rows of the cell array
    x = sort(x(:));

    %Determine unique location values
    uniqueLocs = [true;~strcmp(x(1:end-1),x(2:end)) ~= 0] ;
end

%Determine the unique values
uniques = x(uniqueLocs);

if ~isdouble
    x = feval(currClass,x);
end

%Count the number of duplicate values
if nOut == 2
    numUnique = diff([find(uniqueLocs);length(x)+1]);
end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function val=logAXPlusY(logA,logX,logY)
%%LOGAXPLUSY  Find the element-by-element natural logarithm of a.*x+y, for
%             non-negative, vector or scalar a,x, and y, given the
%             logarithms of a, x, and y, trying to avoid overflow and
%             underflow errors in the product and sum. This lets one
%             evaluate log(a*x+y) even if e^(log(A), e^(log(X)) and
%             e^(log(Y)) all would overflow.
%
%INPUTS: logA, logX, logY The real natural logarithms of the non-zero
%                         a, x, and y that one want in the evaluation of
%                         log(a.*x+y). They can be scalars or vectors/
%                         matrices. The logarithms are element-by element,
%                         not matrix.
%
%OUTPUTS: val The value log(a.*x+y) for each corresponding element in a,x,
%             and y.
%
%The algorithm just uses some identities of logarithms and exponentials to
%minimize overflow and underflow problems. Large enough logarithmic values
%can still eventually have underflow/ overflow problems, but at values
%significantly larger than when a, x and y individually would overflow. The
%algorithm is still subject to the same loss of precision one would expect
%to see when evaluating log(x+1) and x is close to zero (e.g x=eps(1)).
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

logm=min(logA+logX,logY);
logM=max(logA+logX,logY);

%Underflow problems can occur cause 1+exp(logm-logM)=1. However, this seems
%inevitable under any method chosen.
val=logM+log(1+exp(logm-logM));
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function pass = sampleTest(op, f, pref)
%SAMPLETEST   Test an evaluation of input OP against a TRIGTECH approximation.
%   SAMPLETEST(OP, F) evaluates both the function OP and its TRIGTECH
%   representation F at one or more points within [-1,1). The difference of
%   these values is computed, and if this is sufficiently small (relative to
%   VSCALE(F)) the test passes and returns TRUE. If the difference is large,
%   it returns FALSE.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Get the interpolation points:
n = length(f);

% Set a tolerance:
tol =sqrt(max(pref.chebfuneps, eps));
tol = tol.*max(vscale(f));

% pseudo random sample points
xeval = [-0.357998918959666; 0.036785641195074];

% Evaluate the TRIGTECH:
vFun = feval(f, xeval);

% Evaluate the op:
vOp = feval(op, xeval);

% If the TRIGTECH evaluation differs from the op evaluation, SAMPLETEST failed:
err = abs(vOp - vFun); % Relative (to vscale) error.
if ( all(max(err) <= tol) )
    pass = true;  % :)
else
    pass = false; % :(
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function fem1d_heat_implicit_test01 ( )

%*****************************************************************************80
%
%% FEM1D_HEAT_IMPLICIT_TEST01 runs a simple test.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM1D_HEAT_IMPLICIT_TEST01:\n' );
  fprintf ( 1, '  Using the finite element method,\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' );
%
%  Set the nodes.
%
  x_num = 21;
  x_min = 0.0;
  x_max = 1.0;
  dx = ( x_max - x_min ) / ( x_num - 1 );
  x = ( linspace ( x_min, x_max, x_num ) )';
%
%  Set the times.
%
  t_num = 201;
  t_min = 0.0;
  t_max = 80.0;
  dt = ( t_max - t_min ) / ( t_num - 1 );
  t = ( linspace ( t_min, t_max, t_num ) )';
%
%  Set finite element information.
%
  element_num = x_num - 1;
  element_node(1,1:element_num) = 1 : x_num - 1;
  element_node(2,1:element_num) = 2 : x_num;
  quad_num = 3;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of X nodes = %d\n', x_num );
  fprintf ( 1, '  X interval = [ %f, %f ]\n', x_min, x_max );
  fprintf ( 1, '  X step size = %f\n', dx );
  fprintf ( 1, '  Number of T steps = %d\n', t_num );
  fprintf ( 1, '  T interval = [ %f, %f ]\n', t_min, t_max );
  fprintf ( 1, '  T step size = %f\n', dt );
  fprintf ( 1, '  Number of elements = %d\n', element_num );
  fprintf ( 1, '  Number of quadrature points = %d\n', quad_num );

  u_mat = zeros ( x_num, t_num );

  for j = 1 : t_num

    if ( j == 1 )

      u(1:x_num,1) = ic_test01 ( x_num, x, t(j) );
      u(1:x_num,1) = bc_test01 ( x_num, x, t(j), u(1:x_num,1) );

    else

      u(1:x_num,1) = fem1d_heat_implicit ( x_num, x, t(j), dt, @k_test01, ...
        @rhs_test01, @bc_test01, element_num, element_node, quad_num, u );

    end

    u(1:x_num,1) = bc_test01 ( x_num, x, t(j), u(1:x_num,1) );
    u_mat(1:x_num,j) = u(1:x_num,1);

  end
%
%  Make a product grid of T and X for plotting.
%
  [ t_mat, x_mat ] = meshgrid ( t, x );
%
%  Make a mesh plot of the solution.
%
  mesh ( t_mat, x_mat, u_mat )
  xlabel ( '<--Time-->' )
  ylabel ( '<--X-->' )
  zlabel ( '<--H(X,T)-->' )
  title ( 'H(X,T) for TEST01, by FEM1D\_HEAT\_IMPLICIT' )
%
%  Write the data to files.
%
  r8mat_write ( 'h_test01.txt', x_num, t_num, u_mat );
  r8vec_write ( 't_test01.txt', t_num, t );
  r8vec_write ( 'x_test01.txt', x_num, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  H(X,T) written to "h_test01.txt"\n' );
  fprintf ( 1, '  T values written to "t_test01.txt"\n' );
  fprintf ( 1, '  X values written to "x_test01.txt"\n' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [vertex,face,d,c] = read_ply(filename)

% read_ply - read data from PLY file.
%
%   [vertex,face] = read_ply(filename);
%
%   'vertex' is a 'nb.vert x 3' array specifying the position of the vertices.
%   'face' is a 'nb.face x 3' array specifying the connectivity of the mesh.
%
%
%   Copyright (c) 2003 Gabriel Peyr?

[d,c] = plyread(filename);
vi = d.face.vertex_indices;
nf = length(vi);
% http://stackoverflow.com/a/6210539/148668
% This is slow...
%lengths = cellfun(@(x)numel(x),vi);
lengths = cellfun('length',vi);
maxLength=max(lengths);
if all(maxLength == lengths)
  face = cell2mat(vi)+1;
else
  face = cell2mat( ...
    cellfun( ...
      @(x)cat(2,x,-ones(1,maxLength-length(x))), ...
      vi, ...
      'UniformOutput',false))+1;
end
vertex = [d.vertex.x, d.vertex.y, d.vertex.z];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Elements,varargout] = plyread(Path,Str)
%PLYREAD   Read a PLY 3D data file.
%   [DATA,COMMENTS] = PLYREAD(FILENAME) reads a version 1.0 PLY file
%   FILENAME and returns a structure DATA.  The fields in this structure
%   are defined by the PLY header; each element type is a field and each
%   element property is a subfield.  If the file contains any comments,
%   they are returned in a cell string array COMMENTS.
%
%   [TRI,PTS] = PLYREAD(FILENAME,'tri') or
%   [TRI,PTS,DATA,COMMENTS] = PLYREAD(FILENAME,'tri') converts vertex
%   and face data into triangular connectivity and vertex arrays.  The
%   mesh can then be displayed using the TRISURF command.
%
%   Note: This function is slow for large mesh files (+50K faces),
%   especially when reading data with list type properties.
%
%   Example:
%   [Tri,Pts] = PLYREAD('cow.ply','tri');
%   trisurf(Tri,Pts(:,1),Pts(:,2),Pts(:,3));
%   colormap(gray); axis equal;
%
%   See also: PLYWRITE

% Pascal Getreuer 2004
% Markus Leuthold 2017-09-21: make reader compatible with InstantMeshes output

[fid,Msg] = fopen(Path,'rt');	% open file in read text mode

if fid == -1, error(Msg); end

Buf = fscanf(fid,'%s',1);
if ~strcmp(Buf,'ply')
   fclose(fid);
   error('Not a PLY file.');
end

%%% read header %%%

Position = ftell(fid);
Format = '';
NumComments = 0;
Comments = {};				% for storing any file comments
NumElements = 0;
NumProperties = 0;
Elements = [];				% structure for holding the element data
ElementCount = [];		% number of each type of element in file
PropertyTypes = [];		% corresponding structure recording property types
ElementNames = {};		% list of element names in the order they are stored in the file
PropertyNames = [];		% structure of lists of property names

while 1
   Buf = fgetl(fid);   								% read one line from file
   BufRem = Buf;
   Token = {};
   Count = 0;

   while ~isempty(BufRem)								% split line into tokens
      [tmp,BufRem] = strtok(BufRem);

      if ~isempty(tmp)
         Count = Count + 1;							% count tokens
         Token{Count} = tmp;
      end
   end

   if Count 		% parse line
      switch lower(Token{1})
      case 'format'		% read data format
         if Count >= 2
            Format = lower(Token{2});

            if Count == 3 & ~strcmp(Token{3},'1.0')
               fclose(fid);
               error('Only PLY format version 1.0 supported.');
            end
         end
      case 'comment'		% read file comment
         NumComments = NumComments + 1;
         Comments{NumComments} = '';
         for i = 2:Count
            Comments{NumComments} = [Comments{NumComments},Token{i},' '];
         end
      case 'element'		% element name
         if Count >= 3
            if isfield(Elements,Token{2})
               fclose(fid);
               error(['Duplicate element name, ''',Token{2},'''.']);
            end

            NumElements = NumElements + 1;
            NumProperties = 0;
   	      Elements = setfield(Elements,Token{2},[]);
            PropertyTypes = setfield(PropertyTypes,Token{2},[]);
            ElementNames{NumElements} = Token{2};
            PropertyNames = setfield(PropertyNames,Token{2},{});
            CurElement = Token{2};
            ElementCount(NumElements) = str2double(Token{3});

            if isnan(ElementCount(NumElements))
               fclose(fid);
               error(['Bad element definition: ',Buf]);
            end
         else
            error(['Bad element definition: ',Buf]);
         end
      case 'property'	% element property
         if ~isempty(CurElement) & Count >= 3
            NumProperties = NumProperties + 1;
            eval(['tmp=isfield(Elements.',CurElement,',Token{Count});'],...
               'fclose(fid);error([''Error reading property: '',Buf])');

            if tmp
               error(['Duplicate property name, ''',CurElement,'.',Token{2},'''.']);
            end

            % add property subfield to Elements
            eval(['Elements.',CurElement,'.',Token{Count},'=[];'], ...
               'fclose(fid);error([''Error reading property: '',Buf])');
            % add property subfield to PropertyTypes and save type
            eval(['PropertyTypes.',CurElement,'.',Token{Count},'={Token{2:Count-1}};'], ...
               'fclose(fid);error([''Error reading property: '',Buf])');
            % record property name order
            eval(['PropertyNames.',CurElement,'{NumProperties}=Token{Count};'], ...
               'fclose(fid);error([''Error reading property: '',Buf])');
         else
            fclose(fid);

            if isempty(CurElement)
               error(['Property definition without element definition: ',Buf]);
            else
               error(['Bad property definition: ',Buf]);
            end
         end
      case 'end_header'	% end of header, break from while loop
         break;
      end
   end
end

%%% set reading for specified data format %%%

if isempty(Format)
	warning('Data format unspecified, assuming ASCII.');
   Format = 'ascii';
end

switch Format
case 'ascii'
   Format = 0;
case 'binary_little_endian'
   Format = 1;
case 'binary_big_endian'
   Format = 2;
otherwise
   fclose(fid);
   error(['Data format ''',Format,''' not supported.']);
end

if ~Format
   Buf = fscanf(fid,'%f');		% read the rest of the file as ASCII data
   BufOff = 1;
else
   % reopen the file in read binary mode
   fclose(fid);

   if Format == 1
      fid = fopen(Path,'r','ieee-le.l64');		% little endian
   else
      fid = fopen(Path,'r','ieee-be.l64');		% big endian
   end

   % find the end of the header again (using ftell on the old handle doesn't give the correct position)
   BufSize = 8192;
   Buf = [blanks(10),char(fread(fid,BufSize,'uchar')')];
   i = [];
   tmp = -11;

   while isempty(i)
   	i = findstr(Buf,['end_header',13,10]);			% look for end_header + CR/LF
   	i = [i,findstr(Buf,['end_header',10])];		% look for end_header + LF

      if isempty(i)
         tmp = tmp + BufSize;
         Buf = [Buf(BufSize+1:BufSize+10),char(fread(fid,BufSize,'uchar')')];
      end
   end

   % seek to just after the line feed
   fseek(fid,i + tmp + 11 + (Buf(i + 10) == 13),-1);
end

%%% read element data %%%

% PLY and MATLAB data types (for fread)
PlyTypeNames = {'char','uchar','short','ushort','int','uint','float','double', ...
    'char8','uchar8','short16','ushort16','int32','uint32','float32','double64','uint8'};
MatlabTypeNames = {'schar','uchar','int16','uint16','int32','uint32','single','double','uint8'};
SizeOf = [1,1,2,2,4,4,4,8];	% size in bytes of each type

for i = 1:NumElements
   % get current element property information
   eval(['CurPropertyNames=PropertyNames.',ElementNames{i},';']);
   eval(['CurPropertyTypes=PropertyTypes.',ElementNames{i},';']);
   NumProperties = size(CurPropertyNames,2);

%   fprintf('Reading %s...\n',ElementNames{i});

   if ~Format	%%% read ASCII data %%%
      for j = 1:NumProperties
         Token = getfield(CurPropertyTypes,CurPropertyNames{j});

         if strcmpi(Token{1},'list')
            Type(j) = 1;
         else
            Type(j) = 0;
			end
      end

      % parse buffer
      if ~any(Type)
         % no list types
         Data = reshape(Buf(BufOff:BufOff+ElementCount(i)*NumProperties-1),NumProperties,ElementCount(i))';
         BufOff = BufOff + ElementCount(i)*NumProperties;
      else
         ListData = cell(NumProperties,1);

         for k = 1:NumProperties
            ListData{k} = cell(ElementCount(i),1);
         end

         % list type
		   for j = 1:ElementCount(i)
   	      for k = 1:NumProperties
      	      if ~Type(k)
         	      Data(j,k) = Buf(BufOff);
            	   BufOff = BufOff + 1;
	            else
   	            tmp = Buf(BufOff);
      	         ListData{k}{j} = Buf(BufOff+(1:tmp))';
         	      BufOff = BufOff + tmp + 1;
            	end
            end
         end
      end
   else		%%% read binary data %%%
      % translate PLY data type names to MATLAB data type names
      ListFlag = 0;		% = 1 if there is a list type
      SameFlag = 1;     % = 1 if all types are the same

      for j = 1:NumProperties
         Token = getfield(CurPropertyTypes,CurPropertyNames{j});

         if ~strcmp(Token{1},'list')			% non-list type
	         tmp = rem(strmatch(Token{1},PlyTypeNames,'exact')-1,8)+1;

            if ~isempty(tmp)
               TypeSize(j) = SizeOf(tmp);
               Type{j} = MatlabTypeNames{tmp};
               TypeSize2(j) = 0;
               Type2{j} = '';

               SameFlag = SameFlag & strcmp(Type{1},Type{j});
	         else
   	         fclose(fid);
               error(['Unknown property data type, ''',Token{1},''', in ', ...
                     ElementNames{i},'.',CurPropertyNames{j},'.']);
         	end
         else											% list type
            if length(Token) == 3
               ListFlag = 1;
               SameFlag = 0;
               tmp = rem(strmatch(Token{2},PlyTypeNames,'exact')-1,8)+1;
               tmp2 = rem(strmatch(Token{3},PlyTypeNames,'exact')-1,8)+1;

               if ~isempty(tmp) & ~isempty(tmp2)
                  TypeSize(j) = SizeOf(tmp);
                  Type{j} = MatlabTypeNames{tmp};
                  TypeSize2(j) = SizeOf(tmp2);
                  Type2{j} = MatlabTypeNames{tmp2};
	       else
   	   	 fclose(fid);
               	 error(['Unknown property data type, ''list ',Token{2},' ',Token{3},''', in ', ...
                   ElementNames{i},'.',CurPropertyNames{j},'.']);
               end
            else
               fclose(fid);
               error(['Invalid list syntax in ',ElementNames{i},'.',CurPropertyNames{j},'.']);
            end
         end
      end

      % read file
      if ~ListFlag
         if SameFlag
            % no list types, all the same type (fast)
            Data = fread(fid,[NumProperties,ElementCount(i)],Type{1})';
         else
            % no list types, mixed type
            Data = zeros(ElementCount(i),NumProperties);

         	for j = 1:ElementCount(i)
        			for k = 1:NumProperties
               	Data(j,k) = fread(fid,1,Type{k});
              	end
         	end
         end
      else
         ListData = cell(NumProperties,1);

         for k = 1:NumProperties
            ListData{k} = cell(ElementCount(i),1);
         end

         if NumProperties == 1
            BufSize = 512;
            SkipNum = 4;
            j = 0;

            % list type, one property (fast if lists are usually the same length)
            while j < ElementCount(i)
               Position = ftell(fid);
               % read in BufSize count values, assuming all counts = SkipNum
               [Buf,BufSize] = fread(fid,BufSize,Type{1},SkipNum*TypeSize2(1));
               Miss = find(Buf ~= SkipNum);					% find first count that is not SkipNum
               fseek(fid,Position + TypeSize(1),-1); 		% seek back to after first count

               if isempty(Miss)									% all counts are SkipNum
                  Buf = fread(fid,[SkipNum,BufSize],[int2str(SkipNum),'*',Type2{1}],TypeSize(1))';
                  fseek(fid,-TypeSize(1),0); 				% undo last skip

                  for k = 1:BufSize
                     ListData{1}{j+k} = Buf(k,:);
                  end

                  j = j + BufSize;
                  BufSize = floor(1.5*BufSize);
               else
                  if Miss(1) > 1									% some counts are SkipNum
                     [Buf2,cc] = fread(fid,[SkipNum,Miss(1)-1],[int2str(SkipNum),'*',Type2{1}],TypeSize(1));
                     Buf2 = Buf2';

                     for k = 1:Miss(1)-1
                        ListData{1}{j+k} = Buf2(k,:);
                     end

                     j = j + k;
                  end
                  % Alec: check if done and rewind one step
                  if j >= ElementCount(i)
                    fseek(fid,-1,0);
                    break;
                  end

                  % read in the list with the missed count
                  SkipNum = Buf(Miss(1));
                  j = j + 1;
                  ListData{1}{j} = fread(fid,[1,SkipNum],Type2{1});
                  BufSize = ceil(0.6*BufSize);
               end
            end
         else
            % list type(s), multiple properties (slow)
            Data = zeros(ElementCount(i),NumProperties);

            for j = 1:ElementCount(i)
         		for k = 1:NumProperties
            		if isempty(Type2{k})
               		Data(j,k) = fread(fid,1,Type{k});
            		else
               		tmp = fread(fid,1,Type{k});
               		ListData{k}{j} = fread(fid,[1,tmp],Type2{k});
		            end
      		   end
      		end
         end
      end
   end

   % put data into Elements structure
   for k = 1:NumProperties
   	if (~Format & ~Type(k)) | (Format & isempty(Type2{k}))
      	eval(['Elements.',ElementNames{i},'.',CurPropertyNames{k},'=Data(:,k);']);
      else
      	eval(['Elements.',ElementNames{i},'.',CurPropertyNames{k},'=ListData{k};']);
		end
   end
end

clear Data ListData;
fclose(fid);

if (nargin > 1 & strcmpi(Str,'Tri')) | nargout > 2
   % find vertex element field
   Name = {'vertex','Vertex','point','Point','pts','Pts'};
   Names = [];

   for i = 1:length(Name)
      if any(strcmp(ElementNames,Name{i}))
         Names = getfield(PropertyNames,Name{i});
         Name = Name{i};
         break;
      end
   end

   if any(strcmp(Names,'x')) & any(strcmp(Names,'y')) & any(strcmp(Names,'z'))
      eval(['varargout{1}=[Elements.',Name,'.x,Elements.',Name,'.y,Elements.',Name,'.z];']);
   else
      varargout{1} = zeros(1,3);
	end

   varargout{2} = Elements;
   varargout{3} = Comments;
   Elements = [];

   % find face element field
   Name = {'face','Face','poly','Poly','tri','Tri'};
   Names = [];

   for i = 1:length(Name)
      if any(strcmp(ElementNames,Name{i}))
         Names = getfield(PropertyNames,Name{i});
         Name = Name{i};
         break;
      end
   end

   if ~isempty(Names)
      % find vertex indices property subfield
	   PropertyName = {'vertex_indices','vertex_indexes','vertex_index','indices','indexes'};

   	for i = 1:length(PropertyName)
      	if any(strcmp(Names,PropertyName{i}))
         	PropertyName = PropertyName{i};
	         break;
   	   end
      end

      if ~iscell(PropertyName)
         % convert face index lists to triangular connectivity
         eval(['FaceIndices=varargout{2}.',Name,'.',PropertyName,';']);
  			N = length(FaceIndices);
   		Elements = zeros(N*2,3);
   		Extra = 0;

			for k = 1:N
   			Elements(k,:) = FaceIndices{k}(1:3);

   			for j = 4:length(FaceIndices{k})
      			Extra = Extra + 1;
	      		Elements(N + Extra,:) = [Elements(k,[1,j-1]),FaceIndices{k}(j)];
   			end
         end
         Elements = Elements(1:N+Extra,:) + 1;
      end
   end
else
   varargout{1} = Comments;
end

author: alecjacobson

repo: gptoolbox


text:
function res = checkfisx(fis)
%CHECKFISX Checks the fuzzy inference system properties for legal values.
%
%   See checkfis for syntax and explanation.
%
%   It is completely based on MATLAB's checkfis
%   with some modifications, so it's compatible
%   with extendent fuzzy rule structure.
%
%   Compare also code of this function
%   with code of the original checkfis.

%   Per Konstantin A. Sidelnikov, 2009.

for i = 1 : getfisx(fis, 'numrules')
    if isempty(fis.rule(i).weight)
        error('weight of rule %d is empty.', i);
    end
    if isempty(fis.rule(i).antecedent)
        error('antecedent of rule %d is empty.', i);
    end
    if isempty(fis.rule(i).consequent)
        error('consequent of rule %d is empty.', i);
    end
    if isempty(fis.rule(i).connection)
        error('connection of rule %d is empty.', i);
    end
end

for i = 1 : getfisx(fis, 'numinputs')
    if isempty(fis.input(i).name)
        error('name of input %d is empty.', i);
    end
    if isempty(fis.input(i).range)
        error('range of input %d is empty.', i);
    end
end

for i = 1 : getfisx(fis, 'numoutputs')
    if isempty(fis.output(i).name)
        error('name of output %d is empty.', i);
    end
    if isempty(fis.output(i).range)
        error('range of output %d is empty.', i);
    end
end

res = 1;


text:
       clear
       load exampleRobots

author: petercorke

repo: robotics-toolbox-matlab


text:
function Y = sinh(X)
    % Symbolic hyperbolic sine.

    % Convert inputs to SymExpression
    % X = SymExpression(X);

    % construct the operation string
    sstr = ['Sinh[' X.s ']'];

    % create a new object with the evaluated string
    Y = SymExpression(sstr);
end

author: ayonga

repo: frost-dev

save_path: github-repos/MATLAB/ayonga-frost-dev


text:
function r8_hyper_2f1_test ( )

%*****************************************************************************80
%
%% TEST0425 tests R8_HYPER_2F1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 October 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_HYPER_2F1_TEST:\n' );
  fprintf ( 1, '  R8_HYPER_2F1 evaluates the hypergeometric 2F1 function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      A       B       C       X      ' );
  fprintf ( 1, ' 2F1                       2F1                     DIFF\n' );
  fprintf ( 1, '                                     ' );
  fprintf ( 1, '(tabulated)               (computed)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, a, b, c, x, fx ] = hyper_2f1_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r8_hyper_2f1 ( a, b, c, x );

    fprintf ( 1, '  %6.2f  %6.2f  %6.2f  %6.2f  %24.16e  %24.16e  %10.4e\n', ...
      a, b, c, x, fx, fx2, abs ( fx - fx2 ) );

  end

  return
end


text:
function rundivcrl(vxyz,cordname,tn)
% rundivcrl(vxyz,cordname,tn)
% This program tests the vector operations of divergence
% and curl in a general curvilinear coordinate system. The
% analysis process employs a symbolic vector vxyz expressed
% as a function of cartesian coordinates x,y,z. Then a
% curvilinear coordinate system defined by function cordname
% is chosen along with a numerical vector tn to specify a
% particular set of curvilinear coordinate values for which
% the divergence and curl of vxyz are computed using both
% cartesian and curvilinear coordinates. Agreement of values
% obtained by the two methods indicates that quantities such
% as the base vectors, metric tensors, Christoffel symbols
% and covariant derivatives are evaluated correctly.
% vxyz      - a symbolic vector typified by the form
%             [x^2*y*z+1; x*y^2*z+2; x*y*z^2+3]
% cordname  - a function specifying the cartesian vector X=
%             [x;y;z] as a function of curvilinear coordinate
%             coordinate variables. See function cone as an
%             example.
% tn        - a numerical vector having three curvilinear
%             coordinate values for which numerical values
%             of divergence and curl are evaluated.
%
% Among the functions provided to define cordinate systems
% are: cylin, sphr, cone, parab, oblate, elipcyl, toroid, and
% elipsod.
%
% In addition the function cordname, which is input to specify
% the coordinate system, other functions used by the program are
% crldivxyz, metric, ndiverge, curls, covardifn, and curv2cart.
% The program call list is indicated by rundivcrl(vxyz,cordname).
% If no input parameters are given, these are read interactively.

%             by Howard Wilson, July, 2007
syms x y z real
disp(' ')
disp('TEST OF DIVERGENCE AND CURL COMPUTATION IN')
disp(' CARTESIAN AND IN CURVILINEAR COORDINATES')
disp(' ')
if nargin==0
  disp('Input the cartesian components of a vector')
  disp('in symbolic form such as: x*y, y*z, z*x')
  str=input('> ?  ','s');
  vxyz=sym(['[',str,']']);
  disp(' ')
  disp('Input a function name for the coordinate system')
  disp('(such as sphr for spherical coordinates)')
  str=input('> ?  ','s');
  cordname=str2func(str);
  disp(' ')
  disp('Input curvilinear coordinate values to specify')
  disp('a position for the divergence and curl to be')
  disp('computed ( for example, typical values for')
  disp('function sphr could be 2.5,pi/3,pi/4 )')
  str=input('> ?  ','s');
  tn=eval(['[',str,']']);
elseif isstr(cordname)
  cordname=str2func(cordname);
end
vxyz=vxyz(:);
% Obtain the divergence and curl of the input
% vector expressed in cartesian form
[vdivxyz,vcrlxyz]=crldivxyz(vxyz(1),vxyz(2),vxyz(3));

% Get base vectors, metric tensor,and Christoffel
% symbols
[X,t,bco,bcn,gco,gcn,cs1,cs2]=metric(cordname);

fprintf('Coordinate function:')
type([func2str(cordname),'.m']), disp(' ')
disp('Test vector V(x,y,z):'), disp(vxyz(:))
disp('Numerical curvilinear coordinate values:')
disp(tn(:).')

% Express the cartesian vector in terms of
% curvilinear coordinates
vxyzc=subs(vxyz,{x,y,z},{X(1),X(2),X(3)});

% Transform vcxyz to curvilinear coordinates. Get
% both covariant and contravariant components.
vcovar=simple(bco.'*vxyzc); vcontr=simple(bcn.'*vxyzc);
%%vcovar=simple(bco'*vxyzc); vcontr=simple(bcn'*vxyzc);

% Obtain the divergence and the curl
vdivcn=ndiverge(vcontr,t,bco,tn);
% pcontr=curls(vcovar,names,bcovar,cs2,tn)
% gives contravariant components in numerical form
vcrlcn=curls(vcovar,t,bco,cs2,tn);
bcon=subs(bco,{t(1),t(2),t(3)},{tn(1),tn(2),tn(3)});
detbcn=double(det(bcon));
bcnn=subs(bcn,{t(1),t(2),t(3)},{tn(1),tn(2),tn(3)});
vcrlcnxyz=curv2cart(vcrlcn,1,bcnn,bcon,t);
vcrlcnxyz=double(vcrlcnxyz);

% Make numerical evaluations using the test vector
% in cartesian form. First evaluate the x,y,z values
% corresponding to the chosen numerical values of the
% curvilinear coordinates.
Xn=subs(X,{t(1),t(2),t(3)},{tn(1),tn(2),tn(3)});

% Then evaluate the vector vxyz, the divergence and
% the curl at position Xn
vxyzn=subs(vxyz,{x,y,z},{Xn(1),Xn(2),Xn(3)});

% The divergence of vxyz at position Xn
vdivxyzn=subs(vdivxyz,{x,y,z},{Xn(1),Xn(2),Xn(3)});

% The curl of vxyz at position Xn
vcrlxyzn=subs(vcrlxyz,{x,y,z},{Xn(1),Xn(2),Xn(3)});

% Print results for the chosen vector and
% coordinate position
vcrlxyzn=vcrlxyzn(:)'; vcrlcnxyz=vcrlcnxyz(:)';
disp('Cartesian coordinate position for test')
disp('vector evaluation:'), disp(Xn(:).')
disp('Jacobian determinant value')
disp(detbcn)
disp('Divergence using cartesian coordinates:')
disp(vdivxyzn)
disp('Divergence using curvilinear coordinates:')
disp(vdivcn)
disp('Curl components using cartesian coordinates:')
disp(vcrlxyzn)
disp('Curl components computed in curvilinear coordinates')
disp('and then transformed back to cartesian components:')
disp(vcrlcnxyz)
disp('All Done')

author: Sable

repo: mcbench-benchmarks


text:
function w = i4vec_direct_product2 ( factor_index, factor_order, ...
  factor_value, factor_num, point_num, w )

%*****************************************************************************80
%
%% I4VEC_DIRECT_PRODUCT2 creates a direct product of I4VEC's.
%
%  Discussion:
%
%    To explain what is going on here, suppose we had to construct
%    a multidimensional quadrature rule as the product of K rules
%    for 1D quadrature.
%
%    The product rule will be represented as a list of points and weights.
%
%    The J-th item in the product rule will be associated with
%      item J1 of 1D rule 1,
%      item J2 of 1D rule 2,
%      ...,
%      item JK of 1D rule K.
%
%    In particular,
%      X(J) = ( X(1,J1), X(2,J2), ..., X(K,JK))
%    and
%      W(J) = W(1,J1) * W(2,J2) * ... * W(K,JK)
%
%    So we can construct the quadrature rule if we can properly
%    distribute the information in the 1D quadrature rules.
%
%    This routine carries out that task for the weights W.
%
%    Another way to do this would be to compute, one by one, the
%    set of all possible indices (J1,J2,...,JK), and then index
%    the appropriate information.  An advantage of the method shown
%    here is that you can process the K-th set of information and
%    then discard it.
%
%  Example:
%
%    Rule 1:
%      Order = 4
%      W(1:4) = ( 2, 3, 5, 7 )
%
%    Rule 2:
%      Order = 3
%      W(1:3) = ( 11, 13, 17 )
%
%    Rule 3:
%      Order = 2
%      W(1:2) = ( 19, 23 )
%
%    Product Rule:
%      Order = 24
%      W(1:24) =
%        ( 2 * 11 * 19 )
%        ( 3 * 11 * 19 )
%        ( 4 * 11 * 19 )
%        ( 7 * 11 * 19 )
%        ( 2 * 13 * 19 )
%        ( 3 * 13 * 19 )
%        ( 5 * 13 * 19 )
%        ( 7 * 13 * 19 )
%        ( 2 * 17 * 19 )
%        ( 3 * 17 * 19 )
%        ( 5 * 17 * 19 )
%        ( 7 * 17 * 19 )
%        ( 2 * 11 * 23 )
%        ( 3 * 11 * 23 )
%        ( 5 * 11 * 23 )
%        ( 7 * 11 * 23 )
%        ( 2 * 13 * 23 )
%        ( 3 * 13 * 23 )
%        ( 5 * 13 * 23 )
%        ( 7 * 13 * 23 )
%        ( 2 * 17 * 23 )
%        ( 3 * 17 * 23 )
%        ( 5 * 17 * 23 )
%        ( 7 * 17 * 23 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer FACTOR_INDEX, the index of the factor being processed.
%    The first factor processed must be factor 1.
%
%    Input, integer FACTOR_ORDER, the order of the factor.
%
%    Input, integer FACTOR_VALUE(FACTOR_ORDER), the factor values for
%    factor FACTOR_INDEX.
%
%    Input, integer FACTOR_NUM, the number of factors.
%
%    Input, integer POINT_NUM, the number of elements in the direct product.
%
%    Input/output, integer W(POINT_NUM), the elements of the
%    direct product, updated by the latest factor.
%
%  Local Parameters:
%
%    Local, integer START, the first location of a block of values to set.
%
%    Local, integer CONTIG, the number of consecutive values to set.
%
%    Local, integer SKIP, the distance from the current value of START
%    to the next location of a block of values to set.
%
%    Local, integer REP, the number of blocks of values to set.
%
  persistent contig;
  persistent rep;
  persistent skip;

  if ( factor_index == 1 )
    contig = 1;
    skip = 1;
    rep = point_num;
    w(1:point_num) = 1;
  end

  rep = rep / factor_order;
  skip = skip * factor_order;

  for j = 1 : factor_order

    start = 1 + ( j - 1 ) * contig;

    for k = 1 : rep
      w(start:start+contig-1) = w(start:start+contig-1) * factor_value(j);
      start = start + skip;
    end

  end

  contig = contig * factor_order;

  return
end


text:
function [first, last] = whatcoef2QL(level, colorF, o, L)
%------------------------------------------------------------------------------
%
% This function determines where to find a 2D gridfunction in the
% 1D storage vector (heap).
% The output is the first and last index within the heap.
% They are supposed to be uniquely determined by the integer level L and
% character o describing the type of the 2D gridfunction.
%
% This function is a two-dimensional lifting scheme utility.
%
% first = integer
%
% last  = integer
%
% level = integer designated as the level of the 2D gridfunction
%
% o = character, should be either 'a' or 'd',
%     describing the type of the 2D gridfunction:
%     'a' relates to approximation (coefficients) and
%     'd' relates to detail (coefficients)
%
% L = 2D integer array of bookkeeping, see function storeR.
%
% See also: storeR, retrieveQ etc.
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: May 22, 2002.
% (c) 1999-2002 Stichting CWI, Amsterdam
%------------------------------------------------------------------------------
if isempty(L)
  first = 1;
  last  = 0;
else
  [nL, mL] = size(L);
  if mL ~= 6
    disp([' whatcoef2QL - ERROR at type ' o ' level ' num2str(level)]);
    error(' whatcoef2QL - books do not fit ')
  else
%   Choose here according to level
    lo1 = lower(o(1));
    if mod(level, 2) == 0
%     At even level the grids are rectangular
      jL = -1;
      j = 1;
      foundit = 0;
      while j <= nL && ~foundit
         if L(j, 1) == level
           switch lo1
             case 'a' , if L(j, 3) == 0 && L(j, 2) == 0
                          jL = j;
                          foundit = 1;
                        end
             case 'd' , if L(j, 3) == 1 && L(j, 2) == 0
                          jL = j;
                          foundit = 1;
                        end
             otherwise
               disp([' whatcoef2QL - ERROR at type ' o ...
                     ' level ' num2str(level)]);
               error(' whatcoef2QL - unknown type of coefficients ')
           end
         end
         j = j + 1;
      end
      if jL == -1
        disp([' whatcoef2QL - ERROR at type ' lo1 ' level ' num2str(level)]);
        error(' whatcoef2QL - no such coefficients in heap')
      else
        nF = L(jL, 4);
        if nF < 1
          disp([' whatcoef2QL - ERROR at type ' lo1 ' level ' num2str(level)]);
          error(' whatcoef2QL - unvalid 1st dimension of target')
        end
        mF = L(jL, 5);
        if mF < 1
          disp([' whatcoef2QL - ERROR at type ' lo1 ' level ' num2str(level)]);
          error(' whatcoef2QL - unvalid 2nd dimension of target')
        end
        heaptr = L(jL, 6);
        if heaptr < 2
          disp([' whatcoef2QL - ERROR at type ' lo1 ' level ' num2str(level)]);
          error(' whatcoef2QL - bookkeeping error ')
        else
          first = heaptr-nF*mF;
          if first < 1
            disp([' whatcoef2QL - ERROR at type ' lo1 ' level ' num2str(level)]);
            error(' whatcoef2QL - heap not that large, dimensions?')
          else
            last  = heaptr-1;
          end
        end
      end
    else
%     At odd level the grids are quincunx-shaped.
      switch colorF
        case '00' , icolorF = 1;
        case '11' , icolorF = 2;
        case '01' , icolorF = 3;
        case '10' , icolorF = 4;
        otherwise
          disp([' whatcoef2QL - ERROR at type ' lo1 ' color ' colorF ...
                ' with level ' num2str(level)]);
          error(' whatcoef2QL - unknown type of coefficients ')
      end
      fix2even = 999;
      FL = [];
      for j=1:nL
        if L(j, 2) == icolorF
           FL = [FL; [ (L(j,1)+fix2even) 0 L(j,3:6)]];
        end
      end
      if isempty(FL)
        first = 1;
        last  = 0;
      else
        [first, last] = whatcoef2QL((level+fix2even), colorF, lo1, FL);
%       Here colorF does not matter as (level+fix2even) is even.
      end
    end
  end
end
%------------------------------------------------------------------------------


text:
function planC = importVarianFDF(dirT,optFile)
% function scanS = importVarianFDF(dirT,optFile)
% Input: directory containing FDF files and CERRoptions
% file or optS structure
% Output: CERR plan (planC)
% Example
% dirT = 'C:\Projects\FDF reader\Test_data_Joel\gems_trans_pre-CA.dat';
% optFile = 'CERROptions.m';
% planC = importVarianFDF(dirT,optFile);
%
%  Written DK 02/16/2010
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

if exist('dirT')
    pathT=dirT;
else
    pathT = uigetdir(pwd','Select Single Data FDF directory:');
    if pathT == 0
        return
    end
end

contents = dir(fullfile(pathT,'*.fdf'));

hWait = waitbar(0,'Reading Slice');
for i=1:length(contents)
    waitbar(i/length(contents),hWait,['Reading Slice ',num2str(i)]);
    [outheader, data] = readVarianFDF(pathT, contents(i).name);
    header(i)=outheader;
    scanArray(:,:,i) = data;
end
close(hWait)

% Check the format of optFile.
if exist('optFile')
    if ~isstruct(optFile)
        optS = opts4Exe(optFile);
    else
        optS = optFile;
    end
else
    optS = opts4Exe([getCERRPath,'CERROptions.json']);
end

planC = initializeCERR;
indexS = planC{end};
planC{indexS.CERROptions} = optS;
planC{indexS.indexS}     = indexS;

scanInfo = initializeScanInfo;

[ZSlice, zOrder]   = sort([header.zValue]);
header(1:end) = header(zOrder);
scanArray(:,:,1:end) = scanArray(:,:,zOrder);

for i=1:length(header)
    scanInfo(i).grid1Units=header(i).grid1Units;
    scanInfo(i).grid2Units=header(i).grid2Units;
    scanInfo(i).sizeOfDimension1=header(i).sizeOfDimension1;
    scanInfo(i).sizeOfDimension2=header(i).sizeOfDimension2;
    scanInfo(i).xOffset=header(i).xOffset;
    scanInfo(i).yOffset=header(i).yOffset;
    scanInfo(i).zValue=header(i).zValue;
    scanInfo(i).CTOffset=1000;
    scanInfo(i).sliceThickness=[];
    scanInfo(i).imageType = header(i).imageType;
    scanInfo(i).patientName = 'FDF^Scan';
end

minScanArray = min(scanArray(:));
maxScanArray = max(scanArray(:));
scanArray = (scanArray-minScanArray)/maxScanArray*4095;

planC{indexS.scan}(1).scanInfo = scanInfo;
planC{indexS.scan}.scanArray = scanArray;
planC{indexS.scan}.scanType = 'MRI';
planC{indexS.scan}.scanUID = createUID('scan');
planC = setUniformizedData(planC,optS);

if nargout == 1
    return
else
    save_planC(planC,optS);
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
classdef SymFunction < SymExpression
    % This class wraps a symbolic expression with additional information
    % for the convenience of compiling and exporting the symbolic
    % expression to a C/C++ source files.
    %
    % Copyright (c) 2016, AMBER Lab
    % All right reserved.
    %
    % Redistribution and use in source and binary forms, with or without
    % modification, are permitted only in compliance with the BSD 3-Clause
    % license, see
    % http://www.opensource.org/licenses/bsd-license.php

    properties (SetAccess = protected, GetAccess = public)
        % An identification name of the function
        %
        % @type char
        Name

        % The symbolic representation of the dependent variables
        %
        % @type SymVariable
        Vars

        % The symbolic representation of the constant parameters
        %
        % @type SymVariable
        Params

    end

    properties (SetAccess = protected, GetAccess = public)
        % The status flags of the SymFunction
        %
        % These flags stores the status of assignment, export and
        % compilation of the associated symbolic function.
        %
        % @type struct
        Status
    end

    properties (SetAccess=private, GetAccess=public)
        % File names of the functions.
        %
        % Each field of ''Funcs'' specifies the name of a function that
        % used for a certain computation of the function.
        %
        % Required fields of funcs:
        %   Func: a string of the function that computes the
        %   function value @type char
        %   Jac: a string of the function that computes the
        %   function Jacobian @type char
        %   JacStruct: a string of the function that computes the
        %   sparsity structure of function Jacobian @type char
        %   Hess: a string of the function that computes the
        %   function Hessian @type char
        %   HessStruct: a string of the function that computes the
        %   sparsity structure of function Hessian @type char
        % @type struct
        Funcs
    end

    methods

        function obj = SymFunction(name, expr, vars, params)
            % The class constructor function.
            %
            % Parameters:
            %  name: the name of the function. @type char
            %  expr: the symbolic expression of the function
            %  @type char
            %  vars: the symbolic representation of variables in the
            %  expression @type cell
            %  params: the symbolic representation of parameters in the
            %  expression @type cell

            obj = obj@SymExpression(expr);

            % validate name string
            assert(isvarname(name), 'The name must be a char vector.');
            obj.Name = name;

            obj.Funcs = struct(...
                'Jac', ['J_' obj.Name],...
                'JacStruct', ['Js_' obj.Name],...
                'Hess', ['H_' obj.Name],...
                'HessStruct', ['Hs_' obj.Name]);

            if nargin > 2
                if ~isempty(vars)
                    if ~iscell(vars), vars = {vars}; end
                    assert(all(cellfun(@(x)isa(x,'SymVariable'), vars)),...
                        'SymFunction:invalidVariables', ...
                        'The dependent variables must be valid SymVariable objects.');
                    obj.Vars = vars;
                else
                    obj.Vars = {};
                end
                % vars = cellfun(@(x)flatten(x), vars,'UniformOutput',false);
                % obj.Vars = tovector([vars{:}]);
            else
                obj.Vars = {};
            end

            if nargin > 3
                if ~isempty(params)
                    if ~iscell(params), params = {params}; end
                    assert(all(cellfun(@(x)isa(x,'SymVariable'), params)),...
                        'SymFunction:invalidVariables', ...
                        'The dependent variables must be valid SymVariable objects.');
                    obj.Params = params;
                else
                    obj.Params = {};
                end
                % params = cellfun(@(x)flatten(x), params,'UniformOutput',false);
                % obj.Params = tovector([params{:}]);
            else
                obj.Params = {};
            end

            obj.Status = struct();
            obj.Status.FunctionExported = false;
            obj.Status.JacobianExported = false;
            obj.Status.HessianExported  = false;
            obj.Status.FunctionSaved    = false;
            obj.Status.FunctionLoaded   = false;
        end

    end

    % methods defined in external files
    methods

        f = export(obj, export_path, varargin);

        [J,Js] = exportJacobian(obj, export_path, varargin);

        [H,Hs] = exportHessian(obj, export_path, varargin);

        file = save(obj, export_path, varargin);

        obj = load(obj, file_path);
    end
end


function [ R,t,s ] = POS( xp,x )
%POS Estimate scaled orthographic projection parameters from 2D-3D
%correspondences
%   xp - 2*n matrix of 2D feature point positions
%   x  - 3*n matrix of 3D feature point positions
%   R, s, t - rotation, scale, translate parameters
%
% The algorithm is loosely based on the first iteration of the POSIT
% algorithm, although I enforce a valid rotation matrix via SVD. This
% function only provides a linear estimate and the reprojection error can
% be improved by subsequently refining the pose estimates using a
% non-linear optimiser. This is achieved by RefineOrthoCam.
%
% Note that SOP has 6 DOF: 3 for rotation, 2 for 2D translation, 1 for
% scale so 3 2D-3D correspondences should be enough. However, our linear
% solution does not enforce any constraints on the 8 estimated parameters
% and so a minimum of 4 points are required.

npts = size(xp,2);

% Build linear system of equations in 8 unknowns of projection matrix
A = zeros(2*npts,8);

A(1:2:2*npts-1,1:3)=x';
A(1:2:2*npts-1,4)=1;

A(2:2:2*npts,5:7)=x';
A(2:2:2*npts,8)=1;

b = reshape(xp,2*npts,1);

% Solve linear system
k = (A\b);

% Extract params from recovered vector
R1 = k(1:3);
R2 = k(5:7);
sTx = k(4);
sTy = k(8);
s = (norm(R1)+norm(R2))/2;
r1 = R1./norm(R1);
r2 = R2./norm(R2);
r3 = cross(r1,r2);
R = [r1'; r2'; r3'];
% Set R to closest orthogonal matrix to estimated rotation matrix
[U,S,V]=svd([r1'; r2'; r3']);
R = U*V';
% Determinant of R must = 1
if (det(R)<0)
    U(3,:)=-U(3,:);
    R = U*V';
end
% Re-estimate scale and translation based on R

% Centroid of rotated 3D points
% Rx = R(1:2,:)*x;
% c3d = mean(Rx,2);
% c2d = mean(xp,2);
% c2d-c3d
%
% a(1,:) = Rx(1,:)-c3d(1);
% a(2,:) = Rx(2,:)-c3d(2);
% b(1,:) = x(1,:)-c2d(1);
% b(2,:) = x(2,:)-c2d(2);
% s2 = a(:)'\b(:);
% t2 = (c2d-c3d)./s;

% Remove scale from translations
t(1,1)=sTx/s;
t(2,1)=sTy/s;

end


text:
function [outexpdir] = ConvertLarvaeSamuel2Trx(inexpdir,rootoutputdir,varargin)

%% set parameters

% default parameters
trxfilestr = 'trx.mat';
moviefilestr = 'movie.mmf';
perframedirstr = 'perframe';
DEBUG = false;
makesoftlink = true;
dotransposeimage = false;

[trxfilestr,moviefilestr,perframedirstr,DEBUG,makesoftlink,dotransposeimage] = ...
  myparse(varargin,...
  'trxfilestr',trxfilestr,...
  'moviefilestr',moviefilestr,...
  'perframedirstr',perframedirstr,...
  'debug',DEBUG,...
  'makesoftlink',makesoftlink,...
  'dotranspose',dotransposeimage);

%% parse experiment directory name
[rootdatadir,inexperiment_name] = fileparts(inexpdir); %#ok<ASGLU>
timestamp = inexperiment_name;

% find experiment file
files = dir(fullfile(inexpdir,['*experiment_',timestamp,'*.mat']));
if isempty(files),
  error('Could not find experiment file in directory %s',inexpdir);
end
for i = 1:numel(files),
  metadata = regexp(files(i).name,['^(?<line_name>.+)@(?<effector>.+)_experiment_',timestamp,'\.mat$'],'names','once');
  if ~isempty(metadata),
    experimentfile = fullfile(inexpdir,files(i).name);
    break;
  end
end
if isempty(metadata),
  error('Could not find experiment file in directory %s',inexpdir);
end
metadata.exp_datetime = timestamp;

%% find movie file
moviefilepattern = sprintf('%s@%s@%s*.mmf',metadata.exp_datetime,metadata.line_name,metadata.effector);
files = dir(fullfile(inexpdir,moviefilepattern));
if isempty(files),
  error('Could not find movie file in %s',inexpdir);
end
inmoviefile = fullfile(inexpdir,files(1).name);

%% create output directory

% make sure the root output directory exists
if ~exist(rootoutputdir,'dir'),
  mkdir(rootoutputdir);
end

% create the experiment directory
outexperiment_name = sprintf('%s_%s_%s',metadata.line_name,metadata.effector,metadata.exp_datetime);
outexpdir = fullfile(rootoutputdir,outexperiment_name);
if ~exist(outexpdir,'dir'),
  mkdir(outexpdir);
end

%% load in experiment data
expdata = load(experimentfile);

%% compute trajectories
nflies = numel(expdata.experiment_1.track);
trx = [];
perframedata = struct;
SCALE = 10; % centimeters instead of millimeters

if dotransposeimage,
  XIND = 1;
  YIND = 2;
else
  XIND = 2;
  YIND = 1;
end

for i = 1:nflies,

  fprintf('Larva %d...\n',i);

  trk = struct;
  trk.nframes = expdata.experiment_1.track(i).npts;

  % centroid
  loc = double(cat(2,expdata.experiment_1.track(i).pt.loc));
  trk.x_mm = loc(1,:);
  trk.y_mm = loc(2,:);

  % grab other data about the larva
  mid = double(cat(2,expdata.experiment_1.track(i).pt.mid));
  head = double(cat(2,expdata.experiment_1.track(i).pt.head));
  tail = double(cat(2,expdata.experiment_1.track(i).pt.tail));
  spine = {expdata.experiment_1.track(i).pt.spine};
  contour = {expdata.experiment_1.track(i).pt.contour};

  % approximate pxpermm
  trk.pxpermm = 1/expdata.experiment_1.camcalinfo.realUnitsPerPixel;

  % convert to pixels, as everything else will be for plotting purposes
  loc_px = cat(1,expdata.experiment_1.camcalinfo.r2cX(loc(1,:),loc(2,:)),...
    expdata.experiment_1.camcalinfo.r2cY(loc(1,:),loc(2,:)));
  trk.x = loc_px(XIND,:)+1;
  trk.y = loc_px(YIND,:)+1;

  % midpoint
  xmid_cm = mid(1,:);
  ymid_cm = mid(2,:);
  mid_px = cat(1,expdata.experiment_1.camcalinfo.r2cX(mid(1,:),mid(2,:)),...
    expdata.experiment_1.camcalinfo.r2cY(mid(1,:),mid(2,:)));
  trk.xmid = mid_px(XIND,:)+1;
  trk.ymid = mid_px(YIND,:)+1;

  % head
  xhead_cm = head(1,:);
  yhead_cm = head(2,:);
  head_px = cat(1,expdata.experiment_1.camcalinfo.r2cX(head(1,:),head(2,:)),...
    expdata.experiment_1.camcalinfo.r2cY(head(1,:),head(2,:)));
  trk.xhead = head_px(XIND,:)+1;
  trk.yhead = head_px(YIND,:)+1;

  % tail
  xtail_cm = tail(1,:);
  ytail_cm = tail(2,:);
  tail_px = cat(1,expdata.experiment_1.camcalinfo.r2cX(tail(1,:),tail(2,:)),...
    expdata.experiment_1.camcalinfo.r2cY(tail(1,:),tail(2,:)));
  trk.xtail = tail_px(XIND,:)+1;
  trk.ytail = tail_px(YIND,:)+1;

  % spine
  nspinepts = max(cellfun(@(x) size(x,2),spine));
  isspine = ~cellfun(@isempty,spine);
  xspine_cm = nan(nspinepts,trk.nframes);
  yspine_cm = nan(nspinepts,trk.nframes);
  xspine_cm(:,isspine) = cell2mat(cellfun(@(x) double(x(1,:)'),spine(isspine),'UniformOutput',false));
  yspine_cm(:,isspine) = cell2mat(cellfun(@(x) double(x(2,:)'),spine(isspine),'UniformOutput',false));
  if dotransposeimage,
    trk.xspine = expdata.experiment_1.camcalinfo.r2cX(xspine_cm,yspine_cm)+1;
    trk.yspine = expdata.experiment_1.camcalinfo.r2cY(xspine_cm,yspine_cm)+1;
  else
    trk.yspine = expdata.experiment_1.camcalinfo.r2cX(xspine_cm,yspine_cm)+1;
    trk.xspine = expdata.experiment_1.camcalinfo.r2cY(xspine_cm,yspine_cm)+1;
  end

  % contour
  xcontour_cm = cellfun(@(x) double(x(1,:)),contour,'UniformOutput',false);
  ycontour_cm = cellfun(@(x) double(x(2,:)),contour,'UniformOutput',false);
  if dotransposeimage,
    trk.xcontour = cellfun(@(x,y) expdata.experiment_1.camcalinfo.r2cX(x,y)+1, xcontour_cm,ycontour_cm,'UniformOutput',false);
    trk.ycontour = cellfun(@(x,y) expdata.experiment_1.camcalinfo.r2cY(x,y)+1, xcontour_cm,ycontour_cm,'UniformOutput',false);
  else
    trk.ycontour = cellfun(@(x,y) expdata.experiment_1.camcalinfo.r2cX(x,y)+1, xcontour_cm,ycontour_cm,'UniformOutput',false);
    trk.xcontour = cellfun(@(x,y) expdata.experiment_1.camcalinfo.r2cY(x,y)+1, xcontour_cm,ycontour_cm,'UniformOutput',false);
  end

  % covariance matrix is already in pixels; use it to get major, minor,
  % orientation
  S = cellfun(@(x) [x(1),x(2);x(2),x(3)], {expdata.experiment_1.track(i).pt.cov},'UniformOutput',false);
  S = cat(3,S{:});
  if ~dotransposeimage,
    S = [S(2,2,:),S(2,1,:);S(1,2,:),S(1,1,:)];
  end
  [a,b,theta] = cov2ell(S);
  % note that we don't worry about the zero-indexin correction because the
  % difference will cancel out
  thetahead = atan2(trk.yhead-trk.ytail,trk.xhead-trk.xtail);
  thetaflip = modrange(thetahead+modrange(theta-thetahead,-pi/2,pi/2),-pi,pi);
  trk.a = a/2;
  trk.b = b/2;
  trk.theta = thetaflip;

  % convert a, b, theta to mm
  if dotransposeimage,
    x1 = trk.x-1 + a.*cos(thetaflip);
    x2 = trk.x-1 - a.*cos(thetaflip);
    y1 = trk.y-1 + a.*sin(thetaflip);
    y2 = trk.y-1 - a.*sin(thetaflip);
  else
    y1 = trk.x-1 + a.*cos(thetaflip);
    y2 = trk.x-1 - a.*cos(thetaflip);
    x1 = trk.y-1 + a.*sin(thetaflip);
    x2 = trk.y-1 - a.*sin(thetaflip);
  end
  x1_cm = expdata.experiment_1.camcalinfo.c2rX(x1,y1);
  y1_cm = expdata.experiment_1.camcalinfo.c2rY(x1,y1);
  x2_cm = expdata.experiment_1.camcalinfo.c2rX(x2,y2);
  y2_cm = expdata.experiment_1.camcalinfo.c2rY(x2,y2);
  trk.theta_mm = atan2(y2_cm-y1_cm,x2_cm-x1_cm);
  trk.a_mm = sqrt((y2_cm-y1_cm).^2+(x2_cm-x1_cm).^2)/4;
  if dotransposeimage,
    x1 = trk.x-1 + b.*cos(thetaflip+pi/2);
    x2 = trk.x-1 - b.*cos(thetaflip+pi/2);
    y1 = trk.y-1 + b.*sin(thetaflip+pi/2);
    y2 = trk.y-1 - b.*sin(thetaflip+pi/2);
  else
    y1 = trk.x-1 + b.*cos(thetaflip+pi/2);
    y2 = trk.x-1 - b.*cos(thetaflip+pi/2);
    x1 = trk.y-1 + b.*sin(thetaflip+pi/2);
    x2 = trk.y-1 - b.*sin(thetaflip+pi/2);
  end
  x1_cm = expdata.experiment_1.camcalinfo.c2rX(x1,y1);
  y1_cm = expdata.experiment_1.camcalinfo.c2rY(x1,y1);
  x2_cm = expdata.experiment_1.camcalinfo.c2rX(x2,y2);
  y2_cm = expdata.experiment_1.camcalinfo.c2rY(x2,y2);
  trk.b_mm = sqrt((y2_cm-y1_cm).^2+(x2_cm-x1_cm).^2)/4;

  % all this was actually in centimeters, so convert to millimeters
  trk.x_mm = trk.x_mm*SCALE;
  trk.y_mm = trk.y_mm*SCALE;
  trk.a_mm = trk.a_mm*SCALE;
  trk.b_mm = trk.b_mm*SCALE;
  trk.pxpermm = trk.pxpermm/SCALE;

  % frame rate
  trk.dt = diff([expdata.experiment_1.track(i).pt.et]);

  trk.firstframe = double(expdata.experiment_1.track(i).startFrame+1);
  trk.endframe = trk.firstframe+trk.nframes-1;
  trk.off = 1-trk.firstframe;
  trk.id = expdata.experiment_1.track(i).trackNum;
  trx = structappend(trx,trk);

  % TODO: compute some useful per-frame data

  % speed of each body point
  perframedata.headspeed{i} = sqrt(diff(yhead_cm).^2+diff(xhead_cm).^2)*SCALE./trk.dt;
  perframedata.midspeed{i} = sqrt(diff(ymid_cm).^2+diff(xmid_cm).^2)*SCALE./trk.dt;
  perframedata.tailspeed{i} = sqrt(diff(ytail_cm).^2+diff(xtail_cm).^2)*SCALE./trk.dt;

  % head direction
  perframedata.headdir{i} = atan2(yhead_cm-ymid_cm,xhead_cm-xmid_cm);
  perframedata.dheaddir{i} = diff(perframedata.headdir{i})./trk.dt;

  % tail to mid direction
  perframedata.taildir{i} = atan2(ymid_cm-ytail_cm,xmid_cm-xtail_cm);
  perframedata.dtaildir{i} = diff(perframedata.taildir{i})./trk.dt;

  % bend
  perframedata.bend{i} = modrange(perframedata.headdir{i}-perframedata.taildir{i},-pi,pi);

  % change in bend angle
  perframedata.dbend{i} = modrange(diff(perframedata.bend{i}),-pi,pi);

  % skeleton length
  perframedata.spinelength{i} = sum(sqrt(diff(yspine_cm).^2+diff(xspine_cm).^2))*SCALE;

  % change in skeleton length
  perframedata.dspinelength{i} = diff(perframedata.spinelength{i}) ./ trk.dt;

  % velocity relative to head direction
  perframedata.headvelforward_headdir{i} = (diff(xhead_cm).*cos(perframedata.headdir{i}(1:end-1)) + ...
    diff(yhead_cm).*sin(perframedata.headdir{i}(1:end-1))) * SCALE ./ trk.dt;
  perframedata.headvelsideways_headdir{i} = (diff(xhead_cm).*cos(perframedata.headdir{i}(1:end-1)+pi/2) + ...
    diff(yhead_cm).*sin(perframedata.headdir{i}(1:end-1))+pi/2) * SCALE ./ trk.dt;
  perframedata.midvelforward_headdir{i} = (diff(xmid_cm).*cos(perframedata.headdir{i}(1:end-1)) + ...
    diff(ymid_cm).*sin(perframedata.headdir{i}(1:end-1))) * SCALE ./ trk.dt;
  perframedata.midvelsideways_headdir{i} = (diff(xmid_cm).*cos(perframedata.headdir{i}(1:end-1)+pi/2) + ...
    diff(ymid_cm).*sin(perframedata.headdir{i}(1:end-1))+pi/2) * SCALE ./ trk.dt;
  perframedata.tailvelforward_headdir{i} = (diff(xtail_cm).*cos(perframedata.headdir{i}(1:end-1)) + ...
    diff(ytail_cm).*sin(perframedata.headdir{i}(1:end-1))) * SCALE ./ trk.dt;
  perframedata.tailvelsideways_headdir{i} = (diff(xtail_cm).*cos(perframedata.headdir{i}(1:end-1)+pi/2) + ...
    diff(ytail_cm).*sin(perframedata.headdir{i}(1:end-1))+pi/2) * SCALE ./ trk.dt;
  perframedata.centroidvelforward_headdir{i} = (diff(trk.x_mm).*cos(perframedata.headdir{i}(1:end-1)) + ...
    diff(trk.y_mm).*sin(perframedata.headdir{i}(1:end-1))) ./ trk.dt;
  perframedata.centroidvelsideways_headdir{i} = (diff(trk.x_mm).*cos(perframedata.headdir{i}(1:end-1)+pi/2) + ...
    diff(trk.y_mm).*sin(perframedata.headdir{i}(1:end-1))+pi/2) ./ trk.dt;

%   % velocity relative to tail dir
%   perframedata.headvelforward_taildir{i} = (diff(xhead_mm).*cos(perframedata.taildir{i}(1:end-1)) + ...
%     diff(yhead_mm).*sin(perframedata.taildir{i}(1:end-1))) * SCALE ./ trk.dt;
%   perframedata.headvelsideways_taildir{i} = (diff(xhead_mm).*cos(perframedata.taildir{i}(1:end-1)+pi/2) + ...
%     diff(yhead_mm).*sin(perframedata.taildir{i}(1:end-1))+pi/2) * SCALE ./ trk.dt;
%   perframedata.midvelforward_taildir{i} = (diff(xmid_mm).*cos(perframedata.taildir{i}(1:end-1)) + ...
%     diff(ymid_mm).*sin(perframedata.taildir{i}(1:end-1))) * SCALE ./ trk.dt;
%   perframedata.midvelsideways_taildir{i} = (diff(xmid_mm).*cos(perframedata.taildir{i}(1:end-1)+pi/2) + ...
%     diff(ymid_mm).*sin(perframedata.taildir{i}(1:end-1))+pi/2) * SCALE ./ trk.dt;
%   perframedata.tailvelforward_taildir{i} = (diff(xtail_mm).*cos(perframedata.taildir{i}(1:end-1)) + ...
%     diff(ytail_mm).*sin(perframedata.taildir{i}(1:end-1))) * SCALE ./ trk.dt;
%   perframedata.tailvelsideways_taildir{i} = (diff(xtail_mm).*cos(perframedata.taildir{i}(1:end-1)+pi/2) + ...
%     diff(ytail_mm).*sin(perframedata.taildir{i}(1:end-1))+pi/2) * SCALE ./ trk.dt;
%   perframedata.centroidvelforward_taildir{i} = (diff(trk.x_mm).*cos(perframedata.taildir{i}(1:end-1)) + ...
%     diff(trk.y_mm).*sin(perframedata.taildir{i}(1:end-1))) ./ trk.dt;
%   perframedata.centroidvelsideways_taildir{i} = (diff(trk.x_mm).*cos(perframedata.taildir{i}(1:end-1)+pi/2) + ...
%     diff(trk.y_mm).*sin(perframedata.taildir{i}(1:end-1))+pi/2) ./ trk.dt;
%
%   % velocity relative to ellipse orientation
%   perframedata.headvelforward_bodydir{i} = (diff(xhead_mm).*cos(trk.theta_mm(1:end-1)) + ...
%     diff(yhead_mm).*sin(trk.theta_mm(1:end-1))) * SCALE ./ trk.dt;
%   perframedata.headvelsideways_bodydir{i} = (diff(xhead_mm).*cos(trk.theta_mm(1:end-1)+pi/2) + ...
%     diff(yhead_mm).*sin(trk.theta_mm(1:end-1))+pi/2) * SCALE ./ trk.dt;
%   perframedata.midvelforward_bodydir{i} = (diff(xmid_mm).*cos(trk.theta_mm(1:end-1)) + ...
%     diff(ymid_mm).*sin(trk.theta_mm(1:end-1))) * SCALE ./ trk.dt;
%   perframedata.midvelsideways_bodydir{i} = (diff(xmid_mm).*cos(trk.theta_mm(1:end-1)+pi/2) + ...
%     diff(ymid_mm).*sin(trk.theta_mm(1:end-1))+pi/2) * SCALE ./ trk.dt;
%   perframedata.tailvelforward_bodydir{i} = (diff(xtail_mm).*cos(trk.theta_mm(1:end-1)) + ...
%     diff(ytail_mm).*sin(trk.theta_mm(1:end-1))) * SCALE ./ trk.dt;
%   perframedata.tailvelsideways_bodydir{i} = (diff(xtail_mm).*cos(trk.theta_mm(1:end-1)+pi/2) + ...
%     diff(ytail_mm).*sin(trk.theta_mm(1:end-1))+pi/2) * SCALE ./ trk.dt;
%   perframedata.centroidvelforward_bodydir{i} = (diff(trk.x_mm).*cos(trk.theta_mm(1:end-1)) + ...
%     diff(trk.y_mm).*sin(trk.theta_mm(1:end-1))) ./ trk.dt;
%   perframedata.centroidvelsideways_bodydir{i} = (diff(trk.x_mm).*cos(trk.theta_mm(1:end-1)+pi/2) + ...
%     diff(trk.y_mm).*sin(trk.theta_mm(1:end-1))+pi/2) ./ trk.dt;

end

%%

if DEBUG,
  colors = jet(nflies);
  readframe = get_readframe_fcn(inmoviefile);
  t = 100;
  clf;
  imagesc(readframe(t));
  axis image;
  hold on;
  for i = 1:nflies,
    firstframe = double(expdata.experiment_1.track(i).startFrame+1);
    endframe = trx(i).firstframe+trx(i).nframes-1;
    off = 1-firstframe;
    if firstframe > t || endframe < t, continue; end
    j = t+off;
    plot(trx(i).xcontour{j}+1,trx(i).ycontour{j}+1,'-','color',colors(i,:))
    plot(trx(i).xspine(:,j)+1,trx(i).yspine(:,j)+1,'-','color',colors(i,:),'linewidth',2)
    plot([trx(i).xhead(j),trx(i).xmid(j),trx(i).xtail(j)]+1,[trx(i).yhead(j),trx(i).ymid(j),trx(i).ytail(j)]+1,'.-','color',colors(i,:))
  end

end

%% create the trx file

outmatname = fullfile(outexpdir,trxfilestr);
save(outmatname,'trx');

%% per-frame data

units = parseunits('unit'); %#ok<NASGU>

perframedir = fullfile(outexpdir,perframedirstr);
if ~exist(perframedir,'dir'),
  mkdir(perframedir);
end

% save the perframedata we computed above
fns = fieldnames(perframedata);
for i = 1:numel(fns),
  fn = fns{i};
  data = perframedata.(fn); %#ok<NASGU>
  matfilename = fullfile(perframedir,[fn,'.mat']);
  fprintf('Saving data to %s...\n',matfilename);
  save(matfilename,'data','units');
end

% for some reason, the derived measurements have different indices
idx = cell(1,nflies);
for i = 1:nflies,
  idx{i} = expdata.experiment_1.track(i).getDerivedQuantity('mapptstointerped',false);
end

fielddict = struct;
fielddict.vel = 'vel_vel';
fielddict.speed = 'vel_speed';
fielddict.vnorm = 'vel_vnorm';
fielddict.theta = 'vel_theta';
fielddict.adjspeed = 'vel_adjspeed';
%fielddict.speed_diff_local = 'speed_diff_local';
fielddict.deltatheta = 'acc_deltatheta';
fielddict.ddtheta = 'acc_ddtheta';
fielddict.curv = 'acc_curv';
fielddict.lrdtheta = 'lrdtheta';
fielddict.pathLength = 'pathLength';
fielddict.covRatio = 'cov_covRatio';
fielddict.covTheta = 'cov_covTheta';
fielddict.covMinor = 'cov_covMinor';
fielddict.covMajor = 'cov_covMajor';
fielddict.dcovRatio = 'dcovRatio';
fielddict.sarea = 'sarea';

fnsin = fieldnames(fielddict);

for j = 1:numel(fnsin),
  fnin = fnsin{j};
  fnout = fielddict.(fnin);
  matfilename = fullfile(perframedir,[fnout,'.mat']);
  fprintf('%s -> %s...\n',fnin,matfilename);

  data = cell(1,nflies);

  for i = 1:nflies,
    data{i} = expdata.experiment_1.track(i).getDerivedQuantity(fnin);
  end
  save(matfilename,'data','units');
end

% also save unsmoothed area
data = cell(1,nflies);
for i = 1:nflies,
  data{i} = [expdata.experiment_1.track(i).pt.area];
end
matfilename = fullfile(perframedir,['area','.mat']);
save(matfilename,'data','units');

%% copy over the mmf file

outmoviefile = fullfile(outexpdir,moviefilestr);
if isunix && makesoftlink,
  if exist(outmoviefile,'file'),
    delete(outmoviefile);
  end
  cmd = sprintf('ln -s %s %s',inmoviefile,outmoviefile);
  unix(cmd);
else
  [success,msg] = copyfile(inmoviefile,outmoviefile);
  if ~success,
    error('Error copying file %s to %s: %s',inmoviefile,outmoviefile,msg);
  end
end

author: kristinbranson

repo: JAABA


text:
function tapas_BLT_runEnd_almost(params)
% tapas_BLT_runEnd_almost shows a screen once the task is finished but further
% measures are still needed

% Show Cue
Screen('DrawTexture', params.screen.window, params.text.endfieldmap, [], params.screen.fit, 0);
Screen('Flip', params.screen.window);

% Wait
params = wait2_escapeOption(3, params);

author: translationalneuromodeling

repo: tapas


text:
classdef MaF9 < PROBLEM
% <multi/many> <real>
% ML-DMP

%------------------------------- Reference --------------------------------
% R. Cheng, M. Li, Y. Tian, X. Zhang, S. Yang, Y. Jin, and X. Yao, A
% benchmark test suite for evolutionary many-objective optimization,
% Complex & Intelligent Systems, 2017, 3(1): 67-81.
%------------------------------- 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(Access = private)
        Points;     % Vertexes
        Polygons;   % Infeasible polygons
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            % Parameter setting
            if isempty(obj.M); obj.M = 10; end
            obj.M        = max(obj.M,3);
            obj.D        = 2;
            obj.lower    = [-10000,-10000];
            obj.upper    = [10000,10000];
            obj.encoding = ones(1,obj.D);
            % Generate vertexes
            obj.Points = [];
            [thera,rho] = cart2pol(0,1);
            [obj.Points(:,1),obj.Points(:,2)] = pol2cart(thera-(1:obj.M)*2*pi/obj.M,rho);
            % Generate invalid polygons
            head = repmat((1:obj.M)',ceil(obj.M/2-2),1);
            tail = repmat(1:ceil(obj.M/2-2),obj.M,1);
            tail = head + tail(:);
            obj.Polygons = cell(1,length(head));
            for i = 1 : length(obj.Polygons)
                obj.Polygons{i} = obj.Points(mod((head(i):tail(i))-1,obj.M)+1,:);
                obj.Polygons{i} = [obj.Polygons{i};repmat(2*Intersection(obj.Points(mod([head(i)-1,head(i),tail(i),tail(i)+1]-1,obj.M)+1,:)),size(obj.Polygons{i},1),1)-obj.Polygons{i}];
            end
        end
        %% Repair invalid solutions
        function PopDec = CalDec(obj,PopDec)
            Invalid = getInvalid(PopDec,obj.Polygons,obj.Points);
            while any(Invalid)
                PopDec(Invalid,:) = unifrnd(repmat(obj.lower,sum(Invalid),1),repmat(obj.upper,sum(Invalid),1));
                Invalid           = getInvalid(PopDec,obj.Polygons,obj.Points);
            end
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            PopObj = zeros(size(PopDec,1),size(obj.Points,1));
            for m = 1 : size(obj.Points,1)
                PopObj(:,m) = Point2Line(PopDec,obj.Points(mod(m-1:m,size(obj.Points,1))+1,:));
            end
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            [X,Y] = ndgrid(linspace(-1,1,ceil(sqrt(N))));
            ND    = inpolygon(X(:),Y(:),obj.Points(:,1),obj.Points(:,2));
            R     = obj.CalObj([X(ND),Y(ND)]);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 3
                [X,Y]    = ndgrid(linspace(-1,1,40));
                R        = obj.CalObj([X(:),Y(:)]);
                ND       = inpolygon(X(:),Y(:),obj.Points(:,1),obj.Points(:,2));
                R(~ND,:) = nan;
                R = {reshape(R(:,1),size(X)),reshape(R(:,2),size(X)),reshape(R(:,3),size(X))};
            else
                R = [];
            end
        end
        %% Display a population in the decision space
        function DrawDec(obj,Population)
            Draw(obj.Points([1:end,1],:),'-k','LineWidth',1.5,{'\it x\rm_1','\it x\rm_2',[]});
            Draw(Population.decs);
        end
    end
end

function r = Intersection(p)
    if p(1,1) == p(2,1)
        r(1) = p(1,1);
        r(2) = p(3,2)+(r(1)-p(3,1))*(p(3,2)-p(4,2))/(p(3,1)-p(4,1));
    elseif p(3,1) == p(4,1)
        r(1) = p(3,1);
        r(2) = p(1,2)+(r(1)-p(1,1))*(p(1,2)-p(2,2))/(p(1,1)-p(2,1));
    else
        k1   = (p(1,2)-p(2,2))/(p(1,1)-p(2,1));
        k2   = (p(3,2)-p(4,2))/(p(3,1)-p(4,1));
        r(1) = (k1*p(1,1)-k2*p(3,1)+p(3,2)-p(1,2))/(k1-k2);
        r(2) = p(1,2)+(r(1)-p(1,1))*k1;
    end
end

function Invalid = getInvalid(PopDec,Polygons,Points)
    Invalid = false(size(PopDec,1),1);
    for i = 1 : length(Polygons)
        Invalid = Invalid | inpolygon(PopDec(:,1),PopDec(:,2),Polygons{i}(:,1),Polygons{i}(:,2));
    end
    Invalid = Invalid & ~inpolygon(PopDec(:,1),PopDec(:,2),Points(:,1),Points(:,2));
end

function Distance = Point2Line(PopDec,Line)
    Distance = abs((Line(1,1)-PopDec(:,1)).*(Line(2,2)-PopDec(:,2))-(Line(2,1)-PopDec(:,1)).*(Line(1,2)-PopDec(:,2)))./sqrt((Line(1,1)-Line(2,1)).^2+(Line(1,2)-Line(2,2)).^2);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function make_patch_video(A,C,b,f,Y,cont,options)

defoptions = CNMFSetParms;
if nargin < 7
    options = [];
end

if isfield(options,'ind');  % indeces of components to be shown
    ind = options.ind;
else
    ind = defoptions.ind;
end

if isfield(options,'skip_frame'); % skip frames when showing the video
    skp = options.skip_frame;
else
    skp = defoptions.skip_frame;
end

if isfield(options,'sx');   % size of patch in pixels
    sx = options.sx;
else
    sx = defoptions.sx;
end

if isfield(options,'make_avi');  % save an avi video
    make_avi = options.make_avi;
else
    make_avi = defoptions.make_avi;
end

if isfield(options,'name');  % name of video
    name = options.name;
else
    name = defoptions.name;
end

if isfield(options,'show_background');
    sb = options.show_background;
else
    sb = defoptions.show_background;
end

if isfield(options,'show_contours');
    sc = options.show_contours;
else
    sc = defoptions.show_contours;
end

if isfield(options,'cmap');
    cmap = options.cmap;
else
    cmap = defoptions.cmap;
end

d1 = options.d1;
d2 = options.d2;

if isempty(cont) || nargin < 6; sc = 0; end

T = size(C,2);

if ismatrix(Y)
    Y = reshape(Y,d1,d2,T);
end
sx = min(sx,floor(min(d1,d2)/2));

cm = com(A,d1,d2);

C_rec = reshape(A*C+b*f,d1,d2,T);
C_np = reshape(b*f,d1,d2,T);
if make_avi
    vidObj = VideoWriter(name);
    set(vidObj,'FrameRate',15);
    open(vidObj);
end
fig = figure; %colormap('bone');
    set(gcf, 'PaperUnits', 'inches','Units', 'inches')
    set(gcf, 'PaperPositionMode', 'manual')
    set(gcf, 'PaperPosition',1.5*[0,0, 14, 10.5]/1.5)
    set(gcf, 'Position',1.5*[2,2, 14, 10.5]/1.5)%Yr = reshape(Y,d1,d2,T);
    colormap(cmap);
int_x = zeros(length(ind),2*sx);
int_y = zeros(length(ind),2*sx);
A_com = zeros(2*sx,2*sx,length(ind));
for i = 1:length(ind)
    int_x(i,:) = round(cm(ind(i),1)) + (-(sx-1):sx);
    if int_x(i,1)<1
        int_x(i,:) = int_x(i,:) + 1 - int_x(i,1);
    end
    if int_x(i,end)>d1
        int_x(i,:) = int_x(i,:) - (int_x(i,end)-d1);
    end
    int_y(i,:) = round(cm(ind(i),2)) + (-(sx-1):sx);
    if int_y(i,1)<1
        int_y(i,:) = int_y(i,:) + 1 - int_y(i,1);
    end
    if int_y(i,end)>d2
        int_y(i,:) = int_y(i,:) - (int_y(i,end)-d2);
    end
    A_temp = reshape(A(:,ind(i)),d1,d2);
    A_com(:,:,i) = A_temp(int_x(i,:),int_y(i,:));
end
mY = 0.7*max(Y(:));
mC = 0.7*max(C_rec(:));
mB = 0.95*max(C_np(:));

up = zeros(1,4);
for i = 1:length(ind)
    up(i) = 0.95*max(max(A_com(:,:,i)))*max(C(ind(i),:));
    ff = find(cont{ind(i)}(1,2:end)<0.1);
    cont{ind(i)}(:,ff+1) = [];
end

for t = 1:skp:T
    subplot(4,6,[1,2,7,8]); imagesc(Y(:,:,t),[0,mY]);
        title('Raw data','fontweight','bold','fontsize',16);
        axis equal; axis tight; axis off;
        if sc
            hold on;
            for i = 1:length(ind)
                plot(cont{ind(i)}(1,2:end),cont{ind(i)}(2,2:end),'w'); hold on;
                text(min(cont{ind(i)}(1,2:end))-2,min(cont{ind(i)}(2,2:end))-2,num2str(ind(i)),'color','w','fontsize',14,'fontweight','bold'); hold on;
            end
            hold off;
        end
    subplot(4,6,[3,4,9,10]); imagesc(C_rec(:,:,t) - (1-sb)*C_np(:,:,t),[0,mC]); %title(sprintf('%i',t));
        title('Denoised','fontweight','bold','fontsize',16);
        xlabel(sprintf('Timestep %i out of %i',t,T),'fontweight','bold','fontsize',16);
        axis equal; axis tight; set(gca,'XTick',[],'YTick',[]);
        if sc
            hold on;
            for i = 1:length(ind)
                plot(cont{ind(i)}(1,2:end),cont{ind(i)}(2,2:end),'w'); hold on;
            end
            hold off;
        end
    subplot(4,6,[5,6,11,12]); imagesc(Y(:,:,t)-C_rec(:,:,t),[-mC/2,mC/2]);
        title('Residual 2x','fontweight','bold','fontsize',16);
        axis equal; axis tight; axis off;
    subplot(4,6,[13,14,19,20]); imagesc(C_np(:,:,t),[0,mB]);
        title({'Background'; 'synchronized activity'},'fontweight','bold','fontsize',16);
        axis equal; axis tight; axis off;

    for i = 1:4
        subplot(4,6,14+i); imagesc(A_com(:,:,i)*C(ind(i),t),[0,up(i)/1.2]); axis square;
        if sc
            hold on;
            plot(cont{ind(i)}(1,2:end)-int_y(i,1)+1  ,cont{ind(i)}(2,2:end)-int_x(i,1)+1,'w');
            hold off;
        end
        xlabel(sprintf('Component %i',ind(i)),'fontweight','bold','fontsize',14);
        set(gca,'XTick',[],'YTick',[]);
        if i == 2
            handle = title('Representative spatial components','fontweight','bold','fontsize',16);
            set(handle,'position',get(handle,'position')./[0.43,1,1])
        end
        subplot(4,6,20+i); imagesc(Y(int_x(i,:),int_y(i,:),t),[0,up(i)+mB/2]); axis square; axis off;
        if sc
            hold on;
            plot(cont{ind(i)}(1,2:end)-int_y(i,1)+1  ,cont{ind(i)}(2,2:end)-int_x(i,1)+1,'w');
            hold off;
        end
        if i == 2
            handle = title('Raw data patches','fontweight','bold','fontsize',16);
            set(handle,'position',get(handle,'position')./[0.43,1,1])
        end
    end
    drawnow;
    if make_avi
        currFrame = getframe(fig);
        writeVideo(vidObj,currFrame);
    else
        pause(0.05);
    end
end
if make_avi
    close(vidObj);
end

end

author: flatironinstitute

repo: CaImAn-MATLAB