text:
%kvdiskstr 'Generates a disk structuring element'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vdiskstr.pane file
%
% Parameters:
% OutputFile: o 'Output Struc File', required: 'output resulting struc file definition'
% InputFile: trigger 'Trigger Input ', optional: 'trigger input'
% Integer: d 'Disk Diameter ', default: 5: 'diameter of the disk to be created'
% Integer: off 'Nonplanar offset', default: 0: 'value to be added to nonplanar elements'
% MultiChoice: t 'Select Metric', default: 1: 'choose euclidean metric'
% String: t 'Euclidean ', default: 'euclidean': 'choose city-block metric'
% String: t 'City-Block', default: 'city-block': 'choose chessboard metric'
% String: t 'Chessboard', default: 'chessboard': 'lets user choose which kind of disk is to be created'
%    Choices are:
%
% Example: o = kvdiskstr(trigger, {'o','';'trigger','';'d',5;'off',0;'t',1;'t','euclidean';'t','city-block';'t','chessboard'})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vdiskstr - Generates a disk structuring element
%
%  DESCRIPTION
% .I vdiskstr
% generates special structuring elements which are disks in one of three well-known
% metrics (2-D case):
%
%  _
% Euclidean : D(p,q) = sqrt( (px-qx)^2 + (py-qy)^2 )
%
%  _
% City-Block: D(p,q) = |px-qx| + |py-qy|
%
%  _
% Chessboard: D(p,q) = max( |px-qx| , |py-qy| )
%
% Disks in each one of the above metrics are respectively circles,
% diamonds and squares. For a concise review of this subject, see Gonzales &
% Wintz, "Digital Image Processing", chapter 2.
%
% The disks generated may be planar or nonplanar, 2-dimensional or 3-dimensional.
% The default is euclidean planar 2-D.
%
% The offset makes sense only for the nonplanar case, in which case the value
% specified is added to the element. If no offset or 0 is specified, the element
% values will be adjusted so as to have zero in the borders. For example, a
% 5x5 chessboard nonplanar 2-D element will look like this:
%                           0 0 0 0 0
%                           0 1 1 1 0
%                           0 1 2 1 0
%                           0 1 1 1 0
%                           0 0 0 0 0
%
% The disk-like structuring elements  are very useful in many morphological operations
% such as  distance transforms, granulometry, skeleton algorithms etc.
%
%
%
%  EXAMPLES
%
% vdiskstr -r 3  -t "city-block"
% Creates a 2-D disk-like structuring element of radius 3 according to the
% city-block metric.
%
%  "SEE ALSO"
% lvdiskstr(3)
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993-1997 Junior Barrera, Roberto Lotufo.  All rights reserved.
%

function varargout = kvdiskstr(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kvdiskstr(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'o', '__output';'trigger', '__input';'d', 5;'off', 0;'t', 1;'t', 'euclidean';'t', 'city-block';'t', 'chessboard'};
maxval={0,1,2,0,0,0,0,0};
minval={0,1,2,0,0,0,0,0};
istoggle=[0,1,0,0,0,1,1,1];
was_set=istoggle * 0;
paramtype={'OutputFile','InputFile','Integer','Integer','MultiChoice','String','String','String'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

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

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% op_zeropad.m
% Jamie Near, McGill University 2014.
%
% USAGE:
% out=op_zeropad(in,zpFactor);
%
% DESCRIPTION:
% Apply zeropadding (a.k.a. zero-filling) to MRS data.
%
% INPUTS:
% in         = input data in matlab structure format.
% zpFactor   = the factor by which the number of points in the fid will be
%             increased.  ie.  if zpFactor =2, then the number of zeros
%             added to the end of the fid will be equal to the number of
%             points in the original spectrum.
%
% OUTPUTS:
% out        = Output dataset following zeropadding.

function out=op_zeropad(in,zpFactor);

if in.flags.zeropadded
    cont=input('WARNING:  Zero padding has already been performed!  Continue anyway?  (y or n)','s');
    if cont=='y'
        %continue;
    else
        error('STOPPING');
    end
end

%calculate how many zeros to add
zp=ceil((in.sz(1)*zpFactor)-in.sz(1));

%Add zeros using MATLAB array zeropadding function;
fids=padarray(in.fids,zp,'post');

%Calculate Specs using fft
specs=fftshift(ifft(fids,[],in.dims.t),in.dims.t);

%recalculate the sz vector
sz=size(fids);

%Now re-calculate t and ppm arrays using the calculated parameters:
f=[(-in.spectralwidth/2)+(in.spectralwidth/(2*sz(1))):...
    in.spectralwidth/(sz(1)):...
    (in.spectralwidth/2)-(in.spectralwidth/(2*sz(1)))];

ppm=-f/(in.Bo*42.577);
ppm=ppm+4.65;

t=[0:in.dwelltime:(sz(1)-1)*in.dwelltime];

%FILLING IN DATA STRUCTURE
out=in;
out.fids=fids;
out.specs=specs;
out.sz=sz;
out.ppm=ppm;
out.t=t;
out.n=sz(1);

%FILLING IN THE FLAGS
out.flags=in.flags;
out.flags.writtentostruct=1;
out.flags.zeropadded=1;

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function str = showmydoc(item, indent)

% function str = showmydoc(item, indent)
% Display help text for a cfg_choice and all of its options, without
% recursive calls to child nodes.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: showmydoc.m 2085 2008-09-12 10:26:59Z volkmar $

rev = '$Rev: 2085 $'; %#ok
str1 = showdoc(item.cfg_item, indent);
citems = subsref(item, substruct('.','val'));
str2{1} = sprintf('This branch contains %d items:', numel(citems));
% Display short listing of branch items first
str3 = cell(1,numel(citems));
for k = 1:numel(citems)
    str3{k} = sprintf('* %s', subsref(citems{k}, substruct('.','name')));
end;
str = [str1(:); {''}; str2(:); str3(:)];

author: spm


text:
function X = sigmrnd(P)
%     X = double(1./(1+exp(-P)))+1*randn(size(P));
    X = double(1./(1+exp(-P)) > rand(size(P)));
end

author: rasmusbergpalm

repo: DeepLearnToolbox


text:
function p = gam_cdf(x,a,b)
%GAM_CDF Cumulative of Gamma probability density function (cdf).
%
%   P = GAM_CDF(X,A,B) Returns the gamma cdf with
%   shape A and inverse scale B, at the values in X.
%
%   The size of X is the common size of the input arguments. A
%   scalar input functions as a constant matrix of the same size as
%   the other inputs.
%
%   The parameterization is as in Gelman et al. (2004). Bayesian Data
%   Analysis (second edition)

% Copyright (c) 1998-2004 Aki Vehtari

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

if nargin < 3,
  error('Requires three input arguments.');
end

p=gammainc(b*x,a);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%%% Projection onto constrained force space + Constraint Stabilization

% Course: Robotic Manipulation and Mobility
% Advisor: Dr. V. Krovi
%
% Homework Number: MIDTERM
%
% Names: Sourish Chakravarty
% 	Hrishi Lalit Shah

function [dX]= ROBO_midterm_f3(t,X)

global l1 lc1 m1 j1 tau1 l2 lc2 m2 j2 tau2 g Lambda2 omega zeta

th1=X(1);
x2=X(2);
y2=X(3);
th2=X(4);

th1d= X(5);
x2d=X(6);
y2d=X(7);
th2d= X(8);
%%%%%%%%%%%%%%%%%%%%%%% CREATING IMPORTANT MATRICES IN THE GOVERNING EQUATION
s1=sin(th1);
s2=sin(th2);
c1=cos(th1);
c2=cos(th2);

%%%% Element - 1
M1= [j1+m1*(lc1^2)];
E1= [1,-1];
G1= [m1*lc1*c1*g];
% B11 = [-l1*s1, l1*c1];
% B12 = [-m1*lc1*c1];

%%%% Element - 2
M2= [m2, 0, 0;
     0, m2, 0;
     0, 0, j2];
E2= [ 0; 0; 1];
G2= [0; m2*g; 0];
% B21= [-1, 0;
%      0, -1;
%      -lc2*s2, lc2*c2];
% B22= [0, -m2, 0]';

%%%% Constraints
C= [x2-l1*c1-lc2*c2;
    y2-l1*s1-lc2*s2];
A= [l1*s1, 1, 0, lc2*s2;
    -l1*c1, 0, 1, -lc2*c2];% Jacobian of constraint matrix
Ad = [l1*c1*th1d, 0, 0, lc2*c2*th2d;
     l1*s1*th1d, 0, 0, lc2*s2*th2d]; % Derivative of matrix A
S = [1, 0;
    -l1*s1, -lc2*s2;
    l1*c1, lc2*c2;
    0, 1]; % Null space of A or Feasible Motion
Sd = [0, 0;
    -l1*c1*th1d, -lc2*c2*th2d;
    -l1*s1*th1d, -lc2*s2*th2d;
    0, 0]; % Derivative of S matrix

%%%%%%%%%%%%%%%%%%%%%%% PROJECTION ONTO CONSTRAINED FORCE SPACE USING
%%%%%%%%%%%%%%%%%%%%%%% CONSTRAINT STABILIZATION
M=[M1, zeros(1,3);
    zeros(3,1), M2];
E=[E1;
    zeros(3,1),E2];
G= [G1;G2];
T =[tau1; tau2];
V=zeros(4,1);

Qd = [th1d, x2d, y2d, th2d]';

t1= inv([M, A';
        A, zeros(2,2)]); % since m=2;
%     size(Ad)
%     size(2*zeta*A)
%     size((Ad+2*zeta*A)*Qd)
%     size((omega^2)*C)
t2= [E*T-V-G; -(Ad+2*zeta*A)*Qd-(omega^2)*C];%%% Applying constraint stabilization
t3= t1*t2;

dX(1:4,1)= [th1d,x2d,y2d,th2d]';
dX(5:8,1)= t3(1:4);
lam = t3(5:end);
Lambda2 = [Lambda2;lam];%%% Passing instaneous values of Lagrange multipliers back to main function
return

author: Sable

repo: mcbench-benchmarks


text:
function neqn = p24_neqn ( )

%*****************************************************************************80
%
%% P24_NEQN returns the number of equations for problem P24.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NEQN, the number of equations.
%
  neqn = 2;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = VBA_spm_detrend(x,p)
% Polynomial detrending over columns
% FORMAT y = spm_detrend(x,p)
% x   - data matrix
% p   - order of polynomial [default: 0]
%
% y   - detrended data matrix
%__________________________________________________________________________
%
% spm_detrend removes linear and nonlinear trends from column-wise data
% matrices.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_detrend.m 5219 2013-01-29 17:07:07Z spm $

% defaults
%--------------------------------------------------------------------------
[m,n] = size(x);
if ~m || ~n
    y = [];
    return
end
if nargin == 1
    p = 0;
end

% centre columns
%--------------------------------------------------------------------------
if ~p
    y = x - ones(m,1)*mean(x);
    return
end

% polynomial adjustment
%--------------------------------------------------------------------------
G     = zeros(m,p+1);
for i = 0:p
    d = (1:m).^i;
    G(:,i+1) = d(:);
end
y     = x - G*(pinv(full(G))*x);

author: MBB-team

repo: VBA-toolbox


text:
function out=mp_euler(precision)

if nargin==0
 out=mpEuler(mp(0));
else
 out=mpEuler(mp(precision));
end

author: opencobra

repo: cobratoolbox


text:
function [newVar newMean state] = hlp_adaptExpWinMovVar(varargin)
% estimate exponential window moving variance (and mean). The
% moving variance is calculated in place based on new data values and
% previous data.

g = arg_define(varargin, ...
        arg_norep({'values'},[],[],'data values. can be scalar, vector or matrix'), ...
        arg_nogui({'instate','State'},[],[],'state'), ...
        arg_sub({'adaptOpts'},{},@hlp_scaleLimits,'Adaptation options'), ...
        arg({'reset'},false,[],'Reset adaptation state') ...
        );

state = g.instate;

if g.reset || ~isfield(g,'instate') || isempty(state)
    % initialize state
    state.lastVar      = 1;
    state.lastMean     = 0;
    state.numRunsSoFar = 0;
end

% adapt limits using exponential window moving average
[state.lastVar state.lastMean] = hlp_expWinMovVar(g.values,state.lastVar, state.lastMean, ...
                                     state.numRunsSoFar,rmfield(g.adaptOpts,'arg_direct'));
state.numRunsSoFar = state.numRunsSoFar + 1;

newVar  = state.lastVar;
newMean = state.lastMean;

repo: REST


text:
function [err,elemErr] = getL2errorRT0(node,elem,exactSigma,sigmah,d)
%% GETL2ERRORRT0 L2 norm of RT0 element.
%
%  The input exactSigma can be a function bundle or vector array.
%  When exactSigma is a vector array, it should be size NT*2, and
%  exactSigma = \nabla u (or exactSigma= K\nabla u), vector array of size
%  NT by 2.
%
%  [err,elemErr] = getL2errorRT0(node,elem,exactSigma,sigmah).
%  err gives the L2 norm between exact flux and approximate one from RT0,
%  elemErr, an NT by 1 array, gives the square of L2 norm between exact
%  flux and approximate one from RT0 on each element.
%
%  [err,elemErr] = getL2errorRT0(node,elem,exactSigma,sigmah,d).
%  If the coefficient 'd' is inlcuded, then the input should be in the way
%  sigmah     = K\nabla u_h  (including K),
%  exactSigma =  \nabla u    (function bundle, don't including K),    or
%  exactSigma = K\nabla u    (vector array, including K).
%  err=||K exactSigma - sigma_h||_{K^(-1)} = |||u-u_h||| is the energy norm.
%  elemErr, an NT by 1 array, gives the square of L2 norm between exact
%  flux and approximate one from RT0 on each element
%
% Example
%
%     maxIt = 5;
%     node = [1,0; 1,1; 0,1; -1,1; -1,0; -1,-1; 0,-1; 0,0]; % nodes
%     elem = [1,2,8; 3,8,2; 8,3,5; 4,5,3; 7,8,6; 5,6,8];    % elements
%     bdFlag = setboundary(node,elem,'Dirichlet');
%     pde = mixBCdata;
%     err = zeros(maxIt,1);
%     h = zeros(maxIt,1);
%     for k = 1:maxIt
%         [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
%         [u,sigma] = PoissonRT0(node,elem,bdFlag,pde);
%         err(k) = getL2errorRT0(node,elem,pde.Du,sigma);
%         h(k) = 1./(sqrt(size(node,1))-1);
%     end
%     r1 = showrateh(h,err,2);
%     legend('||\sigma - \sigma_h||',['h^{' num2str(r1) '}'],...
%            'LOCATION','Best');
%
% See also getHdiverrorRT0, getL2error3RT0.
%
% Created by Ming Wang at Jan 17, 2011, M-lint modified at May 15, 2011.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if exist('d','var') && ~isempty(d)
    if isreal(d)
        K = d;                  % d is an array
    else
    center = (node(elem(:,1),:) + node(elem(:,2),:) + ...
              node(elem(:,3),:))/3;
    K = d(center);              % d is a function
    end
else
    K = [];
end

%% Construct Data Structure
[elem2dof,~,elem2edgeSign] = dofedge(elem);
NT = size(elem,1);% Ndof = max(elem2dof(:)); %N = size(node,1);
[Dlambda,area] = gradbasis(node,elem);
locEdge = [2 3; 3 1; 1 2];

%% Compute square of the L2 error element-wise
[lambda,w] = quadpts(3); % quadrature order is 3
nQuad = size(lambda,1);
err = zeros(NT,1);
rotMat = [0 -1; 1 0]; % rotation matrix for computing rotLambda.
for p = 1:nQuad
    pxy = lambda(p,1)*node(elem(:,1),:) ...
        + lambda(p,2)*node(elem(:,2),:) ...
        + lambda(p,3)*node(elem(:,3),:);
    if ~isnumeric(exactSigma)
        sigmap = exactSigma(pxy);
        if(nargin >=5)&& ~isempty(K) % multiply coeff. if fun bundle flux.
            sigmap = repmat(K,1,2).*sigmap;
        end
    else
        sigmap = exactSigma;         %  K*\nabla u_h.
    end
    sigmahp = zeros(NT,2);
    for k = 1:3 % for each basis
        i = locEdge(k,1); j = locEdge(k,2);
        % phi_k = lambda_iRot_j - lambda_jRot_i;
        sigmahp = sigmahp + repmat(elem2edgeSign(:,k).*sigmah(elem2dof(:,k)),1,2).*...
                   (lambda(p,i)*Dlambda(:,:,j)*rotMat-lambda(p,j)*Dlambda(:,:,i)*rotMat);
    end
    err = err + w(p)*sum((sigmap - sigmahp).^2,2);
end
err = err.*area;               % ||sigma - K\nabla u_h||^2
if(nargin >=5)&& ~isempty(K)   % ||sigma - K\nabla u_h||^2_{K^(-1)}
    err = err./K;
end
elemErr = err;           % ||sigma - K\nabla u_h||^2_{K^(-1)}
% modify the error
err(isnan(err)) = 0;
err = sqrt(sum(err));

repo: ifem


text:
function [E] = gdf(V1, V2)
%
% GDF   Ground distance between two vectors
%    [E] = GDF(F1, F2) is the ground distance between two feature vectors.
%
%    Example:
%    -------
%        v1 = [100, 40, 22];
%        v2 = [50, 100, 80];
%        ...
%        [e] = gdf(v1, v2);
%        ...
%
%    This file and its content belong to Ulas Yilmaz.
%    You are welcome to use it for non-commercial purposes, such as
%    student projects, research and personal interest. However,
%    you are not allowed to use it for commercial purposes, without
%    an explicit written and signed license agreement with Ulas Yilmaz.
%    Berlin University of Technology, Germany 2006.
%    http://www.cv.tu-berlin.de/~ulas/RaRF
%

E = norm(V1 - V2, 2);

end


text:
%IMPEXP_DOLFIN Import/export grid in FEniCS/Dolfin XML format.
%
%   [ GRID ] = IMPEXP_DOLFIN( FILE_NAME, MODE, DATA, USE_MESHFCN, FID )
%   Import or export of a FEniCS/Dolfin (ASCII .xml) grid and data
%   format. FILE_NAME is a string specifying the (root) file name to
%   process. MODE can either be a string indicating import (no
%   boundary reconstruction), import_bdr (boundary reconstruction with
%   gridbdr), or export. For export, DATA can be either a full fea
%   struct or just the grid struct. USE_MESHFCN is a boolean flag to
%   write subdomain numbers in a mesh function instead of the
%   (default) domain field. A GRID struct is output when importing
%   FID is an optional log file handle for message output
%   (negative for gui output or empty for no output).
%
%   Example:
%
%      1) Export and re-import of a unit square with a hole grid.
%
%      grid1 = quad2tri(holegrid());
%      impexp_hdf5( 'featool-fenics-mesh.xml', 'export', grid1 )
%      grid2 = impexp_hdf5( 'featool-fenics-mesh.h5', 'import' )
%      subplot(1,2,1), plotgrid(grid1), title('grid1')
%      subplot(1,2,2), plotgrid(grid2), title('grid2')
%      is_ok = gridcheck(grid2) == 0
%
%   See also FENICS, FENICS_IMPORT, IMPEXP_HDF5

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function [X, nh, hrf] = designMatrix(stim,params,tSeries,reshapeFlag);
%
% Create a design matrix appropriate for the hrf option
% specified in the selected event-related params
% (see eventParamsDefault for more info on the params struct).
%
% Usage:
% [X, nh, hrf] = designMatrix(stim, [params], [tSeries or nFrames], ...
%                                   [reshapeFlag]);
%
%
% If the params.glmHRF option is 0, the design matrix X will
% be appropriate for a deconvolution GLM in which time
% courses relative to trial onset for each condition are
% estimated as part of the GLM. If this flag is positive (1-3),
% a design matrix X will be returned in which there is a single
% predictor function for each condition, using an assumed form
% of the hemodynamic response function (HRF).
%
% Entering the optional tSeries argument is needed if the
% params.glmHRF option is 1 -- estimate HRF from mean
% response to all stim. If entered, the design matrix
% will also be clipped to the # of frames in the tSeries.
% You can enter the # of frames directly as the third
% argument instead of the whole tSeries if you are using
% a different HRF option, but want to clip to the # of frames.
%
% The optional reshape flag, if set to 1 [default 0], will
% cause the otherwise-2D matrix to be set as a 3D matrix
% with the third dimension being different scans. This is
% appropriate for much of the new GLM code.
%
% Also returns nh, the # of time points to use in the
% hemodynamic response window for estimating a GLM.
% (see glm, applyGlm); and hrf, the response function
% used (empty if deconvolving).
%
%
% ras, 04/05
% ras, 08/05 -- imported into mrVista 2.0.
if notDefined('params'),      params = eventParamsDefault;    end
if notDefined('reshapeFlag'), reshapeFlag = 0;                end
if notDefined('tSeries'),     tSeries = [];                   end

tr = params.framePeriod;

% figure out whether an entire tSeries
% was passed, or just a # of frames:
if length(tSeries)==1
    % nFrames is specified, rather than tSeries
    nFrames = tSeries;
    tSeries = [];
elseif ~isempty(tSeries)
   nFrames = size(tSeries,1);
   if nFrames==1
       % need it as a column vector
       tSeries = tSeries';
       nFrames = size(tSeries,1);
   end
end

if isempty(tSeries)
	% default is max frames specified in stim struct
    nFrames = stim.onsetFrames(end);
end

% decide whether we're deconvolving (getting
% estimated time courses for each condition)
% or fitting an HIRF (getting only a single beta value
% for each condition) based on the selected event-related
% hrf parameter. Get a corresponding stimulus matrix:
if params.glmHRF==0
    framesPerScan = max(stim.framesPerRun);

    % make a delta-function matrix for onsets of
    % different conditions:
    s = delta_function_from_parfile(stim.parfiles,tr,framesPerScan);

    % create Toeplitz matrix for estimating deconvolved responses
    % (see papers on 'selective-averaging', e.g. by Randy Buckner
    % or Douglas Greve for the Freesurfer code)
    fw = unique(round(params.timeWindow/tr)); % time window
    hrf = [];
    nConds =  size(s,2);
    [X nh] = eventDeconvolutionMatrix(s,fw);

    % reshape to 2D (will undo later if needed)
    X = reshape(permute(X,[1 3 2]),[size(X,1)*size(X,3) size(X,2)]);
else
    % apply HRF -- first, construct the impulse response function:
    switch params.glmHRF
        case 1,
            % estimate hrf as mean time course to all stim
            if isempty(tSeries)
                errmsg = 'Need tSeries to estimate mean trial response (glmHRF option 1)';
                error(errmsg);
            end

            disp('Estimating mean trial response for GLM')
            tSeries = squeeze(mean(tSeries,2)); % mean across voxels
            anal = eventTimeCourses(tSeries(:),stim,params);
            postStim = find(anal.timeWindow>=0);
            hrf = anal.meanTcs(postStim,anal.condNums>0);
            hrf = hrf(:,params.snrConds);
            if size(hrf,2) > 1
                hrf = nanmeanDims(hrf,2);
            end
            hrf = hrf ./ sum(hrf); % normalize to integral 1
        case 2,  % boynton gamma function
            hrf = 'boynton';
        case 3,  % spm difference-of-gammas
            hrf = 'spm';
        case 4, % dale & buckner HRF
            hrf = 'dale';
    end

    % apply HRF -- return a 2D matrix covering whole time course
    [X, hrf] = eventConvolveHRF(stim,params,hrf,nFrames);

    % we're only returning one predictor per condition: the
    % nh variable should note this:
    nh = 1;
end

if reshapeFlag==1
    % return a 3D matrix w/ runs as the 3rd dimension

    % figure out max # frames per run to
    % use for reshaping
    for run = unique(stim.run)
        ind = find(stim.run==run);
        onsets{run} = stim.onsetFrames(ind);
        framesPerRun(run) = onsets{run}(end)-onsets{run}(1)+1;
    end
    maxFrames = max(framesPerRun);
    nRuns = length(framesPerRun);
    nPredictors = size(X,2);

    % init 3D X matrix
    oldX = X;
    X = zeros(maxFrames,nPredictors,nRuns);

    % reshape (allow for different-length runs)
    for run = 1:nRuns
        rng = onsets{run}(1):onsets{run}(end);
        X(1:framesPerRun(run),:,run) = oldX(rng,:);
    end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef specimenSymmetry < symmetry
% defines a specimen symmetry
%
% usually specimen symmetry is either
% triclinic or orthorhombic
%

properties
  axes = [xvector,yvector,zvector];
  plotOptions = {}
end

  methods

    function s = specimenSymmetry(varargin)
      % defines a specimen symmetry
      %
      % usually specimen symmetry is either triclinic or orthorhombic
      %

      axes = getClass(varargin,'vector3d',[xvector,yvector,zvector]);

      if nargin == 0

        id = 1;
        rot = rotation.id;

      elseif isa(varargin{1},'quaternion') % define the symmetry just by rotations

        rot = varargin{1};

        if check_option(varargin,'pointId')
          id = get_option(varargin,'pointId');
        else
          id = symmetry.rot2pointId(rot,axes);
        end

      else

        id = symmetry.extractPointId(varargin{:});
        rot = symmetry.calcQuat(id,axes);

      end

      s = s@symmetry(id,rot);
      s.axes = axes;

    end

    function display(s)
      disp(' ');
      disp([inputname(1) ' = ' char(s.lattice) ' ' doclink(s) ' ' docmethods(inputname(1))]);
      disp(' ');
    end

  end

  methods (Static = true)

    function cs = loadobj(s)
      % called by Matlab when an object is loaded from an .mat file
      % this overloaded method ensures compatibility with older MTEX
      % versions

       % maybe there is nothing to do
      if isa(s,'specimenSymmetry'), cs = s; return; end

      if isfield(s,'rot')
        rot = s.rot;
      else
        rot = rotation(s.a,s.b,s.c,s.d,s.i);
      end

      if isfield(s,'axes')
        axes = s.axes;
      else
        axes = [];
      end

      if isfield(s,'id')
        id = {'pointId',s.id};
      else
        id = {};
      end

      cs = specimenSymmetry(rot,id{:},axes);

      if isfield(s,'opt'), cs.opt = s.opt; end
      if isfield(s,'plotOptions'), cs.plotOptions = s.plotOptions; end

    end

  end

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
classdef mme_gen3p_opf < mp.mme_gen3p

%   MATPOWER
%   Copyright (c) 2022, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

%     properties
%     end

    methods
        function x0 = interior_x0(obj, mm, nm, dm, x0)
        end
    end     %% methods
end         %% classdef

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function [h0, h1] = dfilters(fname, type)
% DFILTERS	Generate directional 2D filters
%
%	[h0, h1] = dfilters(fname, type)
%
% Input:
%	fname:	Filter name.  Available 'fname' are:
%		'haar':		the "Haar" filters
%		'vk':		McClellan transformed of the filter from
%				    the VK book
%		'ko':		orthogonal filter in the Kovacevic's paper
%		'kos':		smooth 'ko' filter
%		'lax':		17 x 17 by Lu, Antoniou and Xu
%		'sk':		9 x 9 by Shah and Kalker
%		'cd':		7 and 9 McClellan transformed by
%				                       Cohen and Daubechies
%		'pkva':		ladder filters by Phong et al.
%		'oqf_362':	regular 3 x 6 filter
%       'dvmlp'  :  regular linear phase biorthogonal filter with 3 dvm
%		'sinc':		ideal filter (*NO perfect recontruction*)
%       'dmaxflat': diamond maxflat filters obtained from a three stage ladder
%
%	     type:	'd' or 'r' for decomposition or reconstruction filters
%
% Output:
%	h0, h1:	diamond filter pair (lowpass and highpass)

% To test those filters (for the PR condition for the FIR case), verify that:
% conv2(h0, modulate2(h1, 'b')) + conv2(modulate2(h0, 'b'), h1) = 2
% (replace + with - for even size filters)
%
% To test for orthogonal filter
% conv2(h, reverse2(h)) + modulate2(conv2(h, reverse2(h)), 'b') = 2

% The diamond-shaped filter pair
switch fname
    case {'haar'}
	if lower(type(1)) == 'd'
	    h0 = [1, 1] / sqrt(2);
	    h1 = [-1, 1] / sqrt(2);
	else
	    h0 = [1, 1] / sqrt(2);
	    h1 = [1, -1] / sqrt(2);
	end

    case 'vk'	% in Vetterli and Kovacevic book
	if lower(type(1)) == 'd'
	    h0 = [1, 2, 1] / 4;
	    h1 = [-1, -2, 6, -2, -1] / 4;
	else
	    h0 = [-1, 2, 6, 2, -1] / 4;
	    h1 = [-1, 2, -1] / 4;
	end

	% McClellan transfrom to obtain 2D diamond filters
	t = [0, 1, 0; 1, 0, 1; 0, 1, 0] / 4;	% diamond kernel
	h0 = ftrans2(h0, t);
	h1 = ftrans2(h1, t);

    case 'ko'	% orthogonal filters in Kovacevic's thesis
	a0 = 2; a1 = 0.5; a2 = 1;

	h0 = [0,      -a1, -a0*a1, 0;
	      -a2, -a0*a2, -a0,    1;
	      0, a0*a1*a2, -a1*a2, 0];

	% h1 = qmf2(h0);
	h1 = [0, -a1*a2, -a0*a1*a2, 0;
	      1,     a0, -a0*a2, a2;
	      0, -a0*a1,  a1, 0];

	% Normalize filter sum and norm;
	norm = sqrt(2) / sum(h0(:));

	h0 = h0 * norm;
	h1 = h1 * norm;

	if type == 'r'
	    % Reverse filters for reconstruction
	    h0 = h0(end:-1:1, end:-1:1);
	    h1 = h1(end:-1:1, end:-1:1);
	end

    case 'kos'	% Smooth orthogonal filters in Kovacevic's thesis
	a0 = -sqrt(3); a1 = -sqrt(3); a2 = 2+sqrt(3);

	h0 = [0,      -a1, -a0*a1, 0;
	      -a2, -a0*a2, -a0,    1;
	      0, a0*a1*a2, -a1*a2, 0];

	% h1 = qmf2(h0);
	h1 = [0, -a1*a2, -a0*a1*a2, 0;
	      1,     a0, -a0*a2, a2;
	      0, -a0*a1,  a1, 0];

	% Normalize filter sum and norm;
	norm = sqrt(2) / sum(h0(:));

	h0 = h0 * norm;
	h1 = h1 * norm;

	if type == 'r'
	    % Reverse filters for reconstruction
	    h0 = h0(end:-1:1, end:-1:1);
	    h1 = h1(end:-1:1, end:-1:1);
	end

    case 'lax'	% by Lu, Antoniou and Xu
	h = [-1.2972901e-5,  1.2316237e-4, -7.5212207e-5,  6.3686104e-5, ...
	      9.4800610e-5, -7.5862919e-5,  2.9586164e-4, -1.8430337e-4; ...

	      1.2355540e-4, -1.2780882e-4, -1.9663685e-5, -4.5956538e-5, ...
	     -6.5195193e-4, -2.4722942e-4, -2.1538331e-5, -7.0882131e-4; ...

	     -7.5319075e-5, -1.9350810e-5, -7.1947086e-4,  1.2295412e-3, ...
	      5.7411214e-4,  4.4705422e-4,  1.9623554e-3,  3.3596717e-4; ...

	      6.3400249e-5, -2.4947178e-4,  4.4905711e-4, -4.1053629e-3, ...
	     -2.8588307e-3,  4.3782726e-3, -3.1690509e-3, -3.4371484e-3; ...

	      9.6404973e-5, -4.6116254e-5,  1.2371871e-3, -1.1675575e-2, ...
	      1.6173911e-2, -4.1197559e-3,  4.4911165e-3,  1.1635130e-2; ...

	     -7.6955555e-5, -6.5618379e-4,  5.7752252e-4,  1.6211426e-2, ...
	      2.1310378e-2, -2.8712621e-3, -4.8422645e-2, -5.9246338e-3; ...

	      2.9802986e-4, -2.1365364e-5,  1.9701350e-3,  4.5047673e-3, ...
	     -4.8489158e-2, -3.1809526e-3, -2.9406153e-2,  1.8993868e-1; ...

	     -1.8556637e-4, -7.1279432e-4,  3.3839195e-4,  1.1662001e-2, ...
	     -5.9398223e-3, -3.4467920e-3,  1.9006499e-1,  5.7235228e-1];

	h0 = sqrt(2) * [h, h(:, end-1:-1:1); ...
			h(end-1:-1:1, :), h(end-1:-1:1, end-1:-1:1)];

	h1 = modulate2(h0, 'b');

    case 'sk'	% by Shah and Kalker
	h = [ 0.621729,    0.161889,  -0.0126949, -0.00542504, 0.00124838; ...
	      0.161889,   -0.0353769, -0.0162751, -0.00499353, 0; ...
	     -0.0126949,  -0.0162751,  0.00749029, 0, 0; ...
	     -0.00542504,  0.00499353, 0, 0, 0; ...
	      0.00124838, 0, 0, 0, 0];

	h0 = sqrt(2) * [h(end:-1:2, end:-1:2), h(end:-1:2, :); ...
			h(:, end:-1:2), h];

	h1 = modulate2(h0, 'b');

    case 'dvmlp'
    q= sqrt(2); b=.02; b1 = b*b;
    h  = [b/q 0 -2*q*b 0 3*q*b 0 -2*q*b 0 b/q;
          0 -1/(16*q) 0 9/(16*q) 1/q  9/(16*q) 0 -1/(16*q) 0;
          b/q 0 -2*q*b 0 3*q*b 0 -2*q*b 0 b/q];
    g0 = [-b1/q   0   4*b1*q   0    -14*q*b1    0   28*q*b1    0       -35*q*b1  0    28*q*b1    0     -14*q*b1   0    4*b1*q  0     -b1/q;
           0     b/(8*q)   0  -13*b/(8*q) b/q  33*b/(8*q) -2*q*b -21*b/(8*q) 3*q*b -21*b/(8*q)   -2*q*b 33*b/(8*q)  b/q -13*b/(8*q) 0  b/(8*q) 0;
          -q*b1  0  -1/(256*q) + 8*q*b1  0   9/(128*q) - 28*q*b1   -1/(q*16)   -63/(256*q) + 56*q*b1   9/(16*q)   87/(64*q)-70*q*b1      9/(16*q)   -63/(256*q) + 56*q*b1    -1/(q*16)   9/(128*q) - 28*q*b1   0   -1/(256*q) + 8*q*b1   0   -q*b1;
          0     b/(8*q)   0  -13*b/(8*q) b/q  33*b/(8*q) -2*q*b -21*b/(8*q) 3*q*b -21*b/(8*q)   -2*q*b 33*b/(8*q)  b/q -13*b/(8*q) 0  b/(8*q) 0;
          -b1/q   0   4*b1*q   0    -14*q*b1    0   28*q*b1    0       -35*q*b1  0    28*q*b1    0     -14*q*b1   0    4*b1*q  0     -b1/q];
    h1 = modulate2(g0,'b' );
    h0 = h;
    if lower(type(1)) == 'r'
        h1 = modulate2(h ,'b');
        h0 = g0;
    end

    case {'cd', '7-9'}	% by Cohen and Daubechies
	% 1D prototype filters: the '7-9' pair
	h0 = [0.026748757411, -0.016864118443, -0.078223266529, ...
	      0.266864118443, 0.602949018236, 0.266864118443, ...
	      -0.078223266529, -0.016864118443, 0.026748757411];
	g0 = [-0.045635881557, -0.028771763114, 0.295635881557, ...
	      0.557543526229, 0.295635881557, -0.028771763114, ...
	      -0.045635881557];

	if lower(type(1)) == 'd'
	    h1 = modulate2(g0, 'c');
	else
	    h1 = modulate2(h0, 'c');
	    h0 = g0;
	end

	% Use McClellan to obtain 2D filters
	t = [0, 1, 0; 1, 0, 1; 0, 1, 0] / 4;	% diamond kernel
	h0 = sqrt(2) * ftrans2(h0, t);
	h1 = sqrt(2) * ftrans2(h1, t);

    case {'pkva', 'ldtest'}	% Filters from the ladder structure
	% Allpass filter for the ladder structure network
	beta = ldfilter(fname);

	% Analysis filters
	[h0, h1] = ld2quin(beta);

	% Normalize norm
	h0 = sqrt(2) * h0;
	h1 = sqrt(2) * h1;

	% Synthesis filters
	if lower(type(1)) == 'r'
	    f0 = modulate2(h1, 'b');
	    f1 = modulate2(h0, 'b');

	    h0 = f0;
	    h1 = f1;
	end

 case {'pkva-half4'}	% Filters from the ladder structure
	% Allpass filter for the ladder structure network
	beta = ldfilterhalf( 4);

	% Analysis filters
	[h0, h1] = ld2quin(beta);

	% Normalize norm
	h0 = sqrt(2) * h0;
	h1 = sqrt(2) * h1;

	% Synthesis filters
	if lower(type(1)) == 'r'
	    f0 = modulate2(h1, 'b');
	    f1 = modulate2(h0, 'b');

	    h0 = f0;
	    h1 = f1;
	end

    case {'pkva-half6'}	% Filters from the ladder structure
	% Allpass filter for the ladder structure network
	beta = ldfilterhalf( 6);

	% Analysis filters
	[h0, h1] = ld2quin(beta);

	% Normalize norm
	h0 = sqrt(2) * h0;
	h1 = sqrt(2) * h1;

	% Synthesis filters
	if lower(type(1)) == 'r'
	    f0 = modulate2(h1, 'b');
	    f1 = modulate2(h0, 'b');

	    h0 = f0;
	    h1 = f1;
	end

    case {'pkva-half8'}	% Filters from the ladder structure
	% Allpass filter for the ladder structure network
	beta = ldfilterhalf( 8);

	% Analysis filters
	[h0, h1] = ld2quin(beta);

	% Normalize norm
	h0 = sqrt(2) * h0;
	h1 = sqrt(2) * h1;

	% Synthesis filters
	if lower(type(1)) == 'r'
	    f0 = modulate2(h1, 'b');
	    f1 = modulate2(h0, 'b');

	    h0 = f0;
	    h1 = f1;
	end

    case 'sinc'	% The "sinc" case, NO Perfect Reconstruction
	% Ideal low and high pass filters
	flength = 30;

	h0 = fir1(flength, 0.5);
	h1 = modulate2(h0, 'c');

	% Use McClellan to obtain 2D filters
	t = [0, 1, 0; 1, 0, 1; 0, 1, 0] / 4;	% diamond kernel
	h0 = sqrt(2) * ftrans2(h0, t);
	h1 = sqrt(2) * ftrans2(h1, t);

    case {'oqf_362'}	% Some "home-made" filters!
        h0 = sqrt(2) / 64 * ...
             [	sqrt(15),	-3,	0; ...
                 0,		5,	-sqrt(15); ...
                 -2*sqrt(15),	30,	0; ...
                 0,		30,	2*sqrt(15); ...
                 sqrt(15),	5, 	0; ...
                 0, 		-3,	-sqrt(15)]';

         h1 = -reverse2(modulate2(h0, 'b'));

	if type == 'r'
	    % Reverse filters for reconstruction
	    h0 = h0(end:-1:1, end:-1:1);
	    h1 = h1(end:-1:1, end:-1:1);
	end

    case {'test'}	% Only for the shape, not for PR
	h0 = [0, 1, 0; 1, 4, 1; 0, 1, 0];
	h1 = [0, -1, 0; -1, 4, -1; 0, -1, 0];

    case {'testDVM'}	% Only for directional vanishing moment
	h0 = [1, 1; 1, 1] / sqrt(2);
	h1 = [-1, 1; 1, -1] / sqrt(2);

    case 'qmf'	% by Lu, Antoniou and Xu
	% ideal response
    % window
    m=2;
    n=2;
    w=[];
    w1d=kaiser(4*m+1,2.6);
    for n1=-m:m
        for n2 = -n:n
            w(n1+m+1,n2+n+1) = w1d(2*m+n1+n2+1)*w1d(2*m+n1-n2+1) ;
        end
    end
    h=[];
    for n1=-m:m
        for n2 = -n:n
            h(n1+m+1,n2+n+1) =  .5*sinc((n1+n2)/2)*.5*sinc((n1-n2)/2);
        end
    end
    c=sum(sum(h));
    h = sqrt(2)*h/c;
    h0=h.*w;
    h1 = modulate2(h0,'b');

    %h0 = modulate2(h,'r');
    %h1 = modulate2(h,'b');

     case 'qmf2'	% by Lu, Antoniou and Xu
	% ideal response
    % window

    h=[-.001104 .002494 -0.001744 0.004895 -0.000048 -.000311;
    0.008918 -0.002844 -0.025197 -0.017135 0.003905 -0.000081;
    -0.007587 -0.065904 00.100431 -0.055878 0.007023 0.001504;
    0.001725 0.184162 0.632115 0.099414 -0.027006 -0.001110;
    -0.017935 -0.000491 0.191397 -0.001787 -0.010587 0.002060;
    .001353 0.005635 -0.001231 -0.009052 -0.002668 0.000596];
    h0 = h./sum(sum(h));
    h1 = modulate2(h0,'b');

    %h0 = modulate2(h,'r');
    %h1 = modulate2(h,'b');

   case 'dmaxflat4'
       M1 = 1/sqrt(2); M2=M1;
       k1=1-sqrt(2);k3=k1;
       k2 = M1;
       h  = [.25*k2*k3 .5*k2 1+.5*k2*k3]*M1; h = [h fliplr(h(1:end-1))];
       g  = [-.125*k1*k2*k3 0.25*k1*k2 (-0.5*k1-0.5*k3-0.375*k1*k2*k3) 1+ .5*k1*k2]*M2;
       g = [g fliplr(g(1:end-1))];
       B  = dmaxflat(4,0);
       h0 = mctrans(h,B);
       g0 = mctrans(g,B);
       h0 = sqrt(2)*(h0./sum(h0(:)));
       g0 = sqrt(2)*(g0./sum(g0(:)));

       h1 = modulate2(g0,'b' );
    if lower(type(1)) == 'r'
       h1 = modulate2(h0 ,'b');
       h0 = g0;
    end
   case 'dmaxflat5'
       M1 = 1/sqrt(2); M2=M1;
       k1=1-sqrt(2);k3=k1;
       k2 = M1;
       h  = [.25*k2*k3 .5*k2 1+.5*k2*k3]*M1; h = [h fliplr(h(1:end-1))];
       g  = [-.125*k1*k2*k3 0.25*k1*k2 (-0.5*k1-0.5*k3-0.375*k1*k2*k3) 1+ .5*k1*k2]*M2;
       g = [g fliplr(g(1:end-1))];
       B  = dmaxflat(5,0);
       h0 = mctrans(h,B);
       g0 = mctrans(g,B);
       h0 = sqrt(2)*(h0./sum(h0(:)));
       g0 = sqrt(2)*(g0./sum(g0(:)));

       h1 = modulate2(g0,'b' );
    if lower(type(1)) == 'r'
       h1 = modulate2(h0 ,'b');
       h0 = g0;
    end

   case 'dmaxflat6'
       M1 = 1/sqrt(2); M2=M1;
       k1=1-sqrt(2);k3=k1;
       k2 = M1;
       h  = [.25*k2*k3 .5*k2 1+.5*k2*k3]*M1; h = [h fliplr(h(1:end-1))];
       g  = [-.125*k1*k2*k3 0.25*k1*k2 (-0.5*k1-0.5*k3-0.375*k1*k2*k3) 1+ .5*k1*k2]*M2;
       g = [g fliplr(g(1:end-1))];
       B  = dmaxflat(6,0);
       h0 = mctrans(h,B);
       g0 = mctrans(g,B);
       h0 = sqrt(2)*(h0./sum(h0(:)));
       g0 = sqrt(2)*(g0./sum(g0(:)));

       h1 = modulate2(g0,'b' );
    if lower(type(1)) == 'r'
       h1 = modulate2(h0 ,'b');
       h0 = g0;
    end

   case 'dmaxflat7'
       M1 = 1/sqrt(2); M2=M1;
       k1=1-sqrt(2);k3=k1;
       k2 = M1;
       h  = [.25*k2*k3 .5*k2 1+.5*k2*k3]*M1; h = [h fliplr(h(1:end-1))];
       g  = [-.125*k1*k2*k3 0.25*k1*k2 (-0.5*k1-0.5*k3-0.375*k1*k2*k3) 1+ .5*k1*k2]*M2;
       g = [g fliplr(g(1:end-1))];
       B  = dmaxflat(7,0);
       h0 = mctrans(h,B);
       g0 = mctrans(g,B);
       h0 = sqrt(2)*(h0./sum(h0(:)));
       g0 = sqrt(2)*(g0./sum(g0(:)));

       h1 = modulate2(g0,'b' );
    if lower(type(1)) == 'r'
       h1 = modulate2(h0 ,'b');
       h0 = g0;
    end

    otherwise
	% Assume the "degenerated" case: 1D wavelet filters
	[h0, h1] = wfilters(fname, type);
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [hPa] = MPa2hPa(MPa)
% Convert pressure from megapascals to hectopascals.
% Chad Greene 2012
hPa = MPa*10000.0;


text:
%BAGGINGC Bootstrapping and aggregation of classifiers
%
%    W = BAGGINGC (A,CLASSF,N,ACLASSF,T)
%
% INPUT
%   A         Training dataset.
%   CLASSF    The base classifier (default: nmc)
%   N         Number of base classifiers to train (default: 100)
%   ACLASSF   Aggregating classifier (default: meanc), [] for no aggregation.
%   T         Tuning set on which ACLASSF is trained (default: [], meaning use A)
%
% OUTPUT
%    W        A combined classifier (if ACLASSF given) or a stacked
%             classifier (if ACLASSF []).
%
% DESCRIPTION
% Computation of a stabilised version of a classifier by bootstrapping and
% aggregation ('bagging'). In total N bootstrap versions of the dataset A
% are generated and used for training of the untrained classifier CLASSF.
% Aggregation is done using the combining classifier specified in CCLASSF.
% If ACLASSF is a trainable classifier it is trained by the tuning dataset
% T, if given; else A is used for training. The default aggregating classifier
% ACLASSF is MEANC. Default base classifier CLASSF is NMC.
%
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS, NMC, MEANC

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: baggingc.m,v 1.3 2010/06/01 08:47:05 duin Exp $

function w = baggingc (a,clasf,n,rule,t)
		if (nargin < 5),
		prwarning(2,'no tuning set supplied, using training set for tuning (risk of overfit)');
		t = [];
	end
	if (nargin < 4)
		prwarning(2,'aggregating classifier not specified, assuming meanc');
		rule = meanc;
	end
	if (nargin < 3) | isempty(n),
		prwarning(2,'number of repetitions not specified, assuming 100');
		n = 100;
	end
	if (nargin < 2) | isempty(clasf),
		prwarning(2,'base classifier not specified, assuming nmc');
		clasf = nmc;
	end
	if ((nargin < 1) | isempty(a))
		w = prmapping('baggingc',{clasf,n,rule});
		return
	end

	iscomdset(a,t); % test compatibility training and tuning set

	% Concatenate N classifiers on bootstrap samples (100%) taken
	% from the training set.

	w = [];
	for i = 1:n
		w = [w gendat(a)*clasf];
	end

	% If no aggregating classifier is given, just return the N classifiers...

	if (~isempty(rule))

		% ... otherwise, train the aggregating classifier on the train or
		% tuning set.

		if (isempty(t))
			w = traincc(a,w,rule);
		else
			w = traincc(t,w,rule);
		end
	end

	w = setcost(w,a);

	return

author: marianux

repo: ecg-kit

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


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

SCALE = 10; % centimeters instead of millimeters
MAXTIMESTAMPERR = .01;

success = false;
msg = '';

%% set parameters

[inmoviefile,expfile,...
  expdir,moviefilestr,trxfilestr,perframedirstr,...
  arenatype,arenacenterx,arenacentery,...
  arenaradius,arenawidth,arenaheight,...
  dosoftlink,dotransposeimage] = myparse(varargin,...
  'inmoviefile','','expfile','',...
  'expdir','','moviefilestr','movie.ufmf','trxfilestr','trx.mat','perframedirstr','perframe',...
  'arenatype','None','arenacenterx',0,'arenacentery',0,...
  'arenaradius',123,'arenawidth',123,'arenaheight',123,...
  'dosoftlink',false,...
  'dotransposeimage',false);

%% make sure that MAGATAnalyzer code is on the path

res = which('Experiment');
if isempty(res),
  msg = 'MAGATAnalyzer-Matlab-Analysis code must be on the MATLAB path';
  return;
end
res = which('deg2rad');
if isempty(res),
  msg = '"MAGATAnalyzer-Matlab-Analysis/utility functions" code must be on the MATLAB path';
  return;
end

%% check for files

if isempty(expfile),
  msg = 'Experiment file is empty';
  return;
end

if isempty(inmoviefile),
  msg = 'Video file is empty';
  return;
end

if ~exist(expfile,'file'),
  msg = sprintf('Experiment file %s does not exist',expfile);
  return;
end

if ~exist(inmoviefile,'file'),
  msg = sprintf('Video file %s does not exist',inmoviefile);
  return;
end

% output file locations
moviefile = fullfile(expdir,moviefilestr);
trxfile = fullfile(expdir,trxfilestr);
perframedir = fullfile(expdir,perframedirstr);

%% load in experiment data

try
  expdata = load(expfile);
catch ME,
  msg = getReport(ME);
  return;
end

%% compute trajectories

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

nflies = numel(expdata.experiment_1.track);
trx = [];
perframedata = struct;
timestamps = expdata.experiment_1.elapsedTime';
dt = diff(timestamps);

allframeidx = cell(1,nflies);

for i = 1:nflies,

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

  trk = struct;

  % sometimes frames are skipped
  [errt,frameidx] = min(dist2(timestamps',[expdata.experiment_1.track(i).pt.et]'),[],1);
  allframeidx{i} = frameidx;
  maxerrt = max(errt);
  if maxerrt > MAXTIMESTAMPERR,
    msg = sprintf('Could not match timestamps for larva %d',i);
    return;
  end
  trk.firstframe = min(frameidx);
  trk.endframe = max(frameidx);
  trk.nframes = trk.endframe - trk.firstframe+1;
  % make this start counting at 1
  frameidx = frameidx - trk.firstframe + 1;

  if trk.nframes > numel(expdata.experiment_1.track(i).pt),
    fprintf('%d frames missing for larva %d\n',trk.nframes - numel(expdata.experiment_1.track(i).pt),i);
  end

  % centroid
  loc = double(cat(2,expdata.experiment_1.track(i).pt.loc));
  trk.x_mm = nan(1,trk.nframes);
  trk.x_mm(frameidx) = loc(1,:);
  trk.y_mm = nan(1,trk.nframes);
  trk.y_mm(frameidx) = 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};
  area = double(cat(2,expdata.experiment_1.track(i).pt.area));

  % 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 = nan(1,trk.nframes);
  trk.x(frameidx) = loc_px(XIND,:)+1;
  trk.y = nan(1,trk.nframes);
  trk.y(frameidx) = loc_px(YIND,:)+1;

  % area
  trk.area=nan(1,trk.nframes);
  trk.area(frameidx)=area;
  % midpoint
  xmid_cm = mid(XIND,:);
  ymid_cm = mid(YIND,:);
  perframedata.xmid_mm{i} = nan(1,trk.nframes);
  perframedata.xmid_mm{i}(frameidx) = xmid_cm*SCALE;
  perframedata.ymid_mm{i} = nan(1,trk.nframes);
  perframedata.ymid_mm{i}(frameidx) = ymid_cm*SCALE;
  mid_px = cat(1,expdata.experiment_1.camcalinfo.r2cX(mid(1,:),mid(2,:)),...
    expdata.experiment_1.camcalinfo.r2cY(mid(1,:),mid(2,:)));
  trk.xmid = nan(1,trk.nframes);
  trk.xmid(frameidx) = mid_px(XIND,:)+1;
  trk.ymid = nan(1,trk.nframes);
  trk.ymid(frameidx) = mid_px(YIND,:)+1;

  % head
  xhead_cm = head(XIND,:);
  yhead_cm = head(YIND,:);
  perframedata.xhead_mm{i} = nan(1,trk.nframes);
  perframedata.xhead_mm{i}(frameidx) = xhead_cm*SCALE;
  perframedata.yhead_mm{i} = nan(1,trk.nframes);
  perframedata.yhead_mm{i}(frameidx) = yhead_cm*SCALE;
  head_px = cat(1,expdata.experiment_1.camcalinfo.r2cX(head(1,:),head(2,:)),...
    expdata.experiment_1.camcalinfo.r2cY(head(1,:),head(2,:)));
  trk.xhead = nan(1,trk.nframes);
  trk.xhead(frameidx) = head_px(XIND,:)+1;
  trk.yhead = nan(1,trk.nframes);
  trk.yhead(frameidx) = head_px(YIND,:)+1;

  % tail
  xtail_cm = tail(XIND,:);
  ytail_cm = tail(YIND,:);
  perframedata.xtail_mm{i} = nan(1,trk.nframes);
  perframedata.xtail_mm{i}(frameidx) = xtail_cm*SCALE;
  perframedata.ytail_mm{i} = nan(1,trk.nframes);
  perframedata.ytail_mm{i}(frameidx) = ytail_cm*SCALE;
  tail_px = cat(1,expdata.experiment_1.camcalinfo.r2cX(tail(1,:),tail(2,:)),...
    expdata.experiment_1.camcalinfo.r2cY(tail(1,:),tail(2,:)));
  trk.xtail = nan(1,trk.nframes);
  trk.xtail(frameidx) = tail_px(XIND,:)+1;
  trk.ytail = nan(1,trk.nframes);
  trk.ytail(frameidx) = 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(:,frameidx(isspine)) = cell2mat(cellfun(@(x) double(x(1,:)'),spine(isspine),'UniformOutput',false));
  yspine_cm(:,frameidx(isspine)) = cell2mat(cellfun(@(x) double(x(2,:)'),spine(isspine),'UniformOutput',false));
  % spines are backwards
  xspine_cm = xspine_cm(end:-1:1,:);
  yspine_cm = yspine_cm(end:-1:1,:);
  if dotransposeimage,
    trk.xspine = expdata.experiment_1.camcalinfo.r2cY(xspine_cm,yspine_cm)+1;
    trk.yspine = expdata.experiment_1.camcalinfo.r2cX(xspine_cm,yspine_cm)+1;
    trk.xspine_mm = xspine_cm*SCALE;
    trk.yspine_mm = yspine_cm*SCALE;
  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;
    trk.xspine_mm = yspine_cm*SCALE;
    trk.yspine_mm = xspine_cm*SCALE;
  end

  % contour
  xcontour_cm = cellfun(@(x) double(x(1,:)),contour,'UniformOutput',false);
  ycontour_cm = cellfun(@(x) double(x(2,:)),contour,'UniformOutput',false);
  trk.xcontour = cell(1,trk.nframes);
  trk.ycontour = cell(1,trk.nframes);
  perframedata.xcontour_mm{i} = cell(1,trk.nframes);
  perframedata.ycontour_mm{i} = cell(1,trk.nframes);
  if dotransposeimage,
    trk.xcontour(frameidx) = cellfun(@(x,y) expdata.experiment_1.camcalinfo.r2cY(x,y)+1, xcontour_cm,ycontour_cm,'UniformOutput',false);
    trk.ycontour(frameidx) = cellfun(@(x,y) expdata.experiment_1.camcalinfo.r2cX(x,y)+1, xcontour_cm,ycontour_cm,'UniformOutput',false);
    perframedata.xcontour_mm{i}(frameidx) = cellfun(@(x) x*SCALE,xcontour_cm,'UniformOutput',false);
    perframedata.ycontour_mm{i}(frameidx) = cellfun(@(x) x*SCALE,ycontour_cm,'UniformOutput',false);
  else
    trk.ycontour(frameidx) = cellfun(@(x,y) expdata.experiment_1.camcalinfo.r2cX(x,y)+1, xcontour_cm,ycontour_cm,'UniformOutput',false);
    trk.xcontour(frameidx) = cellfun(@(x,y) expdata.experiment_1.camcalinfo.r2cY(x,y)+1, xcontour_cm,ycontour_cm,'UniformOutput',false);
    perframedata.xcontour_mm{i}(frameidx) = cellfun(@(x) x*SCALE,ycontour_cm,'UniformOutput',false);
    perframedata.ycontour_mm{i}(frameidx) = cellfun(@(x) x*SCALE,xcontour_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-indexing correction because the
  % difference will cancel out
  thetahead = atan2(head_px(YIND,:)-tail_px(YIND,:),head_px(XIND,:)-tail_px(XIND,:));
  thetaflip = modrange(thetahead+modrange(theta-thetahead,-pi/2,pi/2),-pi,pi);
  trk.a = nan(1,trk.nframes);
  trk.a(frameidx) = a/2;
  trk.b = nan(1,trk.nframes);
  trk.b(frameidx) = b/2;
  trk.theta = nan(1,trk.nframes);
  trk.theta(frameidx) = thetaflip;

  % convert a, b, theta to mm
  x = loc_px(XIND,:);
  y = loc_px(XIND,:);
  if dotransposeimage,
    x1 = x + a.*cos(thetaflip);
    x2 = x - a.*cos(thetaflip);
    y1 = y + a.*sin(thetaflip);
    y2 = y - a.*sin(thetaflip);
  else
    y1 = x + a.*cos(thetaflip);
    y2 = x - a.*cos(thetaflip);
    x1 = x + a.*sin(thetaflip);
    x2 = x - 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 = nan(1,trk.nframes);
  trk.theta_mm(frameidx) = atan2(y2_cm-y1_cm,x2_cm-x1_cm);
  trk.a_mm = nan(1,trk.nframes);
  trk.a_mm(frameidx) = sqrt((y2_cm-y1_cm).^2+(x2_cm-x1_cm).^2)/4;
  if dotransposeimage,
    x1 = x + b.*cos(thetaflip+pi/2);
    x2 = x - b.*cos(thetaflip+pi/2);
    y1 = y + b.*sin(thetaflip+pi/2);
    y2 = y - b.*sin(thetaflip+pi/2);
  else
    y1 = x + b.*cos(thetaflip+pi/2);
    y2 = x - b.*cos(thetaflip+pi/2);
    x1 = y + b.*sin(thetaflip+pi/2);
    x2 = y - 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 = nan(1,trk.nframes);
  trk.b_mm(frameidx) = 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;

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

  % convert are to square mm from pixels
  trk.area_mm=trk.area/(trk.pxpermm*trk.pxpermm);
  % frame rate
  trk.dt = dt(trk.firstframe:trk.endframe-1);
  %trk.dt = diff([expdata.experiment_1.track(i).pt.et]);

  trk.off = 1-trk.firstframe;
  trk.id = expdata.experiment_1.track(i).trackNum;
  trx = structappend(trx,trk);

end

%% arena parameters

switch lower(arenatype),

  case 'circle',
    arenacenterx_mm = expdata.experiment_1.camcalinfo.c2rX(arenacentery,arenacenterx)*SCALE;
    arenacentery_mm = expdata.experiment_1.camcalinfo.c2rY(arenacentery,arenacenterx)*SCALE;
    arenaradius_mm = arenaradius/trx(1).pxpermm;

    for i = 1:numel(trx),
      trx(i).arena = struct;
      trx(i).arena.arena_radius_mm = arenaradius_mm;
      trx(i).arena.arena_center_mm_x = arenacenterx_mm;
      trx(i).arena.arena_center_mm_y = arenacentery_mm;
    end

  case 'rectangle',

    tl = [arenacenterx - arenawidth/2,arenacentery - arenaheight/2];
    tr = [arenacenterx + arenawidth/2,arenacentery - arenaheight/2];
    bl = [arenacenterx - arenawidth/2,arenacentery + arenaheight/2];
    br = [arenacenterx + arenawidth/2,arenacentery + arenaheight/2];

    tl_mm = [expdata.experiment_1.camcalinfo.c2rX(tl(2),tl(1)),...
      expdata.experiment_1.camcalinfo.c2rY(tl(2),tl(1))]*SCALE;
    tr_mm = [expdata.experiment_1.camcalinfo.c2rX(tr(2),tr(1)),...
      expdata.experiment_1.camcalinfo.c2rY(tr(2),tr(1))]*SCALE;
    bl_mm = [expdata.experiment_1.camcalinfo.c2rX(bl(2),bl(1)),...
      expdata.experiment_1.camcalinfo.c2rY(bl(2),bl(1))]*SCALE;
    br_mm = [expdata.experiment_1.camcalinfo.c2rX(br(2),br(1)),...
      expdata.experiment_1.camcalinfo.c2rY(br(2),br(1))]*SCALE;

    for i = 1:numel(trx),
      trx(i).arena = struct;
      trx(i).arena.tl = tl_mm;
      trx(i).arena.tr = tr_mm;
      trx(i).arena.bl = bl_mm;
      trx(i).arena.br = br_mm;
    end

end

%% create the experiment directory

if ~exist(expdir,'dir'),
  [success1,msg1] = mkdir(expdir);
  if ~success1,
    msg = msg1;
    return;
  end
end

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

%% copy/soft-link movie

if strcmp(fullfile(inmoviefile),fullfile(moviefile)),
  fprintf('Input and out movie files are the same, not copying/linking.\n');
else

  if dosoftlink,
    if isunix,
      cmd = sprintf('ln -s %s %s',inmoviefile,moviefile);
      unix(cmd);
      % test to make sure it worked
      [status,result] = unix(sprintf('readlink %s',moviefile));
      result = strtrim(result);
      if status ~= 0 || ~strcmp(result,inmoviefile),
        warndlg(sprintf('Failed to make soft link, copying %s to %s instead',inmoviefile,moviefile));
        dosoftlink = false;
      end
    elseif ispc,
      cmd = sprintf('mkshortcut.vbs /target:"%s" /shortcut:"%s"',inmoviefile,moviefile);
      fprintf('Making a Windows shortcut file at "%s" with target "%s"\n',inmoviefile,moviefile);
      system(cmd);
      % test to make sure that worked
      [equalmoviefile,didfind] = GetPCShortcutFileActualPath(moviefile);
      if ~didfind || ~strcmp(equalmoviefile,inmoviefile),
        warndlg(sprintf('Failed to make shortcut, copying %s to %s instead',inmoviefile,moviefile));
        dosoftlink = false;
      end
    else
      warndlg(sprintf('Unknown OS, not soft-linking movie file %s',inmoviefile));
      dosoftlink = false;
    end
  end

  if ~dosoftlink,
    [success1,msg1] = copyfile(inmoviefile,moviefile);
    if ~success1,
      msg = msg1;
      success = false;
      return;
    end
  end

end

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

%% per-frame data

% save the perframedata we computed above
fns = fieldnames(perframedata);
% everything is currently in mm
units = parseunits('mm'); %#ok<NASGU>
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 = 'ma_vel';
% fielddict.speed = 'ma_speed';
% fielddict.vnorm = 'ma_vnorm';
% fielddict.theta = 'ma_theta';
% fielddict.adjspeed = 'ma_adjspeed';
% fielddict.lrdtheta = 'ma_lrdtheta';
% fielddict.pathLength = 'ma_pathLength';
% fielddict.covRatio = 'ma_covRatio';
% fielddict.covTheta = 'ma_covTheta';
% fielddict.covMinor = 'ma_covMinor';
% fielddict.covMajor = 'ma_covMajor';
% fielddict.dcovRatio = 'ma_dcovRatio';
% fielddict.sarea = 'ma_aarea';
%
% 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,
%     tmp = expdata.experiment_1.track(i).getDerivedQuantity(fnin);
%     data{i} = nan(size(tmp,1),trx(i).nframes);
%     data{i}(:,allframeidx{i}) = tmp(:,idx{i});
%   end
%   save(matfilename,'data','units');
% end

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

success = true;

author: kristinbranson

repo: JAABA


text:
function region = p23_region ( )

%*****************************************************************************80
%
%% P23_REGION returns the name of the integration region for problem 23.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string REGION, the name of the integration region.
%
  region = 'simplex';

  return
end


text:
function [cumRet, cumProdRet, dailyRet, dailyPortfolio] ...
    = ubah(fid, data, varargins, opts)
% This file is an entry for the market strategy.
%
% function [cumRet, cumProdRet, dailyRet, dailyPortfolio] ...
%     = ubah(fid, data, varargins, opts)
% cumRet: cumulative wealth achived at the end of a period.
% cumProdRet: cumulative wealth achieved till the end each period.
% dailyRet: daily return achieved by a strategy.
% dailyPortfolio: daily portfolios
%
% data: market sequence vectors
% fid: handle for write log file
% varargins: variable parameters
% opts: option parameter for behvaioral control
%
% Example: [cumRet, cumProdRet, dailyRet, dailyPortfolio] ...
%          = market_start(fid, data, {0}, opts);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Bin LI, Steven C.H. Hoi
% Contributors:
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extract the parameters
tc = varargins{1};      % transaction cost fee rate

% Run the market algorithm
[cumRet, cumProdRet, dailyRet, dailyPortfolio] ...
    = ubah_run(fid, data, tc, opts);

end
%%%%%%%%%%%%%%End%%%%%%%%%%%%%%%%%%%%%%

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function c8_atan_test ( )

%*****************************************************************************80
%
%% C8_ATAN_TEST tests C8_ATAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 February 2015
%
%  Author:
%
%    John Burkardt
%
  seed = 123456678;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'C8_ATAN_TEST\n' );
  fprintf ( 1, '  C8_ATAN computes the inverse tangent of a C8.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '       C1=C8_UNIFORM_01          C2 = C8_ATAN(C1)           C3 = C8_TAN(C2)\n' );
  fprintf ( 1, '     ---------------------     ---------------------     ---------------------\n' );
  fprintf ( 1, '\n' );

  for test = 1 : 10

    [ c1, seed ] = c8_uniform_01 ( seed );

    c2 = c8_atan ( c1 );

    c3 = c8_tan ( c2 );

    fprintf ( 1, '  (%12f  %12f)  (%12f  %12f)  (%12f  %12f)\n', ...
      real ( c1 ), imag ( c1 ), real ( c2 ), imag ( c2 ), real ( c3 ), imag ( c3 ) );

  end

  return
end


text:
classdef EventRate
%EventRate Event Rate class constructor.
%
%    EventRate is a class that has been developed around plotting earthquake
%    counts - i.e. the rate of events per unit time. It has evolved to compute
%    other metrics such the hourly mean event rate, median event rate, mean
%    magnitude and cumulative magnitude, which are important metrics for an AVO
%    swarm tracking system.
%
%    EventRate can import information from:
%    (1) a Catalog object.
%    (2) a Datascope database written in the "swarms1.0" schema, defined at AVO.
%        This is the format used by the swarm tracking system (Thompson &
%        West, 2010).
%
%    ER = EventRate(Catalog_OBJECT, 'binsize', BINSIZE) creates an eventrate object
%    from a Catalog object using non-overlapping bins of BINSIZE days.
%
%    ER = EventRate(Catalog_OBJECT, 'binsize', BINSIZE, 'stepsize', STEPSIZE) creates an eventrate object
%    using overlapping bins. If omitted STEPSIZE==BINSIZE.
%
%%   EXAMPLES:
%
%       First create a catalog object from the demo database:
%           dbpath = demodb('avo')
%           catalogObject = readEvents('datascope', 'dbpath', dbpath, ...
%                  'dbeval', ...
%                  'deg2km(distance(lat, lon, 60.4853, -152.7431))<15.0' ...
%                  );
%
%       (1) Create an eventrate object using a binsize of 1 day:
%           erobj = catalogObject.eventrate('binsize', 1);
%
%       (2) Create an eventrate object using a binsize of 1 hour:
%           erobj = catalogObject.eventrate('binsize', 1/24);
%
%       (3) Create an eventrate object using a binsize of 1 hour but a stepsize of 5 minutes:
%           erobj = catalogObject.eventrate('binsize', 1/24, 'stepsize', 5/1440);
%
%%   PROPERTIES
%
%    For a list of all properties type properties(EventRate)
%
%    time                % (array) time of the center of each bin as a DATENUM
%
%    METRICS:
%        counts 		     % (array) number of events in each bin
%        mean_rate           % (array) number of events per hour in each bin
%        median_rate	     % (array) reciprocal of the median time interval between events. Represented as an hourly rate.
%        cum_mag		     % (array) total sum of energy in each bin, represented as a magnitude.
%        mean_mag		     % (array) mean magnitude of events in each bin
%        median_mag          % (array) median magnitude of events in each bin
%        min_mag             % (array) smallest magnitude in each bin
%        max_mag             % (array) largest magnitude in each bin
%
%    SUMMARY DATA:
%        numbins             % (scalar) number of bins used for grouping
%                                events
%        total_counts        % (scalar) sum of counts
%        total_mag           % (scalar) total sum of energy of all catalogObjects, represented as a magnitude
%
%    METADATA:
%        etype               % event type/classification.
%        snum                % (scalar) start date/time in DATENUM format
%        enum                % (scalar) end date/time in DATENUM format
%        binsize             % (scalar) bin size in days
%        stepsize            % (scalar) step size in days
%        region              % (4-element vector) [minlon maxlon minlat maxlat]
%        minmag              % (scalar) magnitudes smaller than this were eliminated
%        dbroot              % path to the original data on disk
%        archiveformat       % indicates if the source is a flat file, or
%                              'daily' or 'monthly' volumes
%        auth                % auth of the events
%
%%   METHODS
%
%    For a list of all methods type methods EventRate
%
%
%%   See also Catalog, Catalog_lite
%
%% AUTHOR: Glenn Thompson

% $Date: 2014-05-06 14:52:40 -0800 (Tue, 06 May 2014) $
% $Revision: 404 $

% I don't think these parts work anymore
%       (4) Create a vector of eventrate objects subclassified using event types 'r', 'e', 'l', 'h', 't':
%               erobj = eventrate(catalogObject, 1, 'etypes', 'relht');
%           To plot counts on separate figures:
%               erobj.plot()
%           To plot counts and energy panels, each event type as a separate figure:
%               erobj.plot('metric', {'counts';'energy'});
%           To plot counts and energy panels on separate figures, each event type as panels:
%               erobj.plot('metric', {'counts';'energy'}, 'plotmode', 'panels');
%           To plot counts and energy panels on separate figures, each event type stacked:
%               erobj.plot('metric', {'counts';'energy'}, 'plotmode', 'stacked');
%
%       (5) A full example:
%               catalogObject = catalog(fullfile(MVO_DATA, 'mbwh_catalog'), 'seisan', 'snum', datenum(1996,10,1), 'enum', datenum(2004,3,1), 'region', 'Montserrat')
%               erobj = eventrate(catalogObject, 365/12, 'stepsize', 1, 'etypes', 'thlr');
%               erobj.plot('metric', {'counts';'energy'}, 'plotmode', 'stacked');
%

%% PROPERTIES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    properties(GetAccess = 'public', SetAccess = 'public')
        time = [];          % (array) in datenum format
        counts = []; 		% (array) number of events in each bin
		mean_rate = [];      % (array) number of events per hour in each bin
		median_rate = [];	% (array) reciprocal of the median time interval between events. Represented as an hourly rate.
		cum_mag = [];		% (array) total sum of energy in each bin, represented as a magnitude.
		mean_mag = [];		% (array)
        median_mag = [];     % (array)
        energy = [];
        total_counts = [];   % (scalar) sum of counts
		total_mag = [];      % (scalar) total sum of energy of all catalogObjects, represented as a magnitude
        numbins = [];        % (scalar)
        min_mag = [];
        max_mag = [];
        etype = '*';
        snum = 0;
        enum = now;
        binsize = 1;
        stepsize = 1;
        misc_fields = {};
        misc_values = {};
    end

 %% PUBLIC METHODS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	methods
        %% CONSTRUCTOR
        function self = EventRate(time, counts, energy, median_energy, ...
                smallest_energy, biggest_energy, median_time_interval, total_counts, ...
                snum, enum, etypes, binsize, stepsize, numbins)
            self.time = time;
            self.counts = counts;
            self.median_rate = 1 ./ (median_time_interval * 24);
            self.median_rate(counts<10) = 0;
            self.median_rate = max([self.counts / (24 * binsize); self.median_rate]);
            self.median_mag = magnitude.eng2mag(median_energy);
            self.energy = energy;
            self.total_counts = total_counts;
            self.numbins = numbins;
            self.min_mag = magnitude.eng2mag(smallest_energy);
            self.max_mag = magnitude.eng2mag(biggest_energy);
            self.etype = etypes;
            self.snum = snum;
            self.enum = enum;
            self.binsize = binsize;
            self.stepsize = stepsize;
            if (enum-snum) < binsize
                error('binsize cannot be bigger than data time range');
            end
        end

        %% ----------------------------------------------
        %% GETTERS
        function cum_mag = get.cum_mag(erobj)
            cum_mag = magnitude.eng2mag(erobj.energy);
        end
        function mean_mag = get.mean_mag(erobj)
            mean_mag = magnitude.eng2mag(erobj.energy./erobj.counts);
        end
        function mean_rate = get.mean_rate(erobj)
            mean_rate = erobj.counts / (24 * erobj.binsize);
        end
        function total_mag = get.total_mag(erobj)
            total_mag = magnitude.eng2mag(sum(erobj.energy));
        end

    end % methods

    methods(Static)
        cookbook()
    end

end

author: geoscience-community-codes

repo: GISMO


text:
classdef FuncRegClass < matlab.mixin.Copyable

    properties
        userfuncdir
        userfuncfiles
        entries
        type
        typeNameIdx
        config
        logger
    end

    methods

        % ----------------------------------------------------------------------------------
        function obj = FuncRegClass(varargin)
            %
            % Syntax:
            %    obj = FuncRegClass(type);
            %    obj = FuncRegClass(obj2);
            %
            % Description:
            %    Generates a list of FuncRegEntryClass objects of type
            %    'run', 'subj', or 'group' for any hmr* user files it finds
            %    in the folder <Rootdir>/FuncRegistry/UserFunctions.
            %    <Rootdir> is relative to the location of this file,
            %    FuncRegClass.m
            %
            % Example:
            %    Generate array of registry entries for all run-level users
            %    functions
            %
            %    freg = FuncRegClass('run')
            %
            %         ===> FuncRegClass with properties:
            %
            %         userfuncdir: 'c:/jdubb/workspaces/homer3/FuncRegistry/UserFunctions/'
            %       userfuncfiles: {1x10 cell}
            %             entries: [1x10 FuncRegEntryClass]
            %                type: 'run'
            %
            %    freg.entries(1)
            %
            %         ===> FuncRegEntryClass with properties:
            %
            %                name: 'hmrR_BandpassFilt'
            %              uiname: 'Bandpass_Filter'
            %        usageoptions: {'Bandpass_Filter'  'dod = hmrR_BandpassFilt( dod, t, hpf, lpf )'  'hmrR_BandpassFilt dod (dod,t hpf %0.3f 0.010 lpf %0.3f 0.500'}
            %              params: {2x2 cell}
            %                help: [1x1 FuncHelpClass]
            %
            global logger

            logger = InitLogger(logger);
            obj.logger = logger;

            obj.entries = FuncRegEntryClass().empty();

            if nargin==0
                return;
            end

            obj.typeNameIdx = 1;
            if isa(varargin{1}, 'FuncRegClass')
                obj.Copy(varargin{1});
                return;
            end

            obj.type = varargin{1};
            if nargin > 1
                obj.typeNameIdx = varargin{2};
            end

            % Get the parameter items from config file relevant to this class
            obj.userfuncdir = FindUserFuncDir(obj);
            obj.userfuncfiles = [];
            obj.Load();
        end

        % ----------------------------------------------------------------------------------
        function Load(obj)
            heading = sprintf('Registry loading %s-level user functions:', [upper(obj.type(1)), obj.type(2:end)]);
            fprintf('%s\n', heading);

            files = cell(length(obj.userfuncdir),1);
            for jj = 1:length(obj.userfuncdir)
                expr = sprintf('%shmr%s_*.m', obj.userfuncdir{jj}, upper(obj.type(obj.typeNameIdx)));
                files{jj} = dir(expr);
                for iF = 1:length(files{jj})
                    files{jj}(iF).folder = obj.userfuncdir{jj};
                end
                fprintf('User functions folder %s\n', obj.userfuncdir{jj});
            end

            h = waitbar(0, heading);

            kk = 1;
            for jj = 1:length(files)
                N = length(files{jj});
                for ii = 1:N
                    if strfind(files{jj}(ii).name, '_result.m')
                        continue;
                    end
                    progressmsg = sprintf('Parsing %s', files{jj}(ii).name);
                    obj.logger.Write('%s\n', progressmsg);
                    waitbar_improved(kk/N, h, sprintf(progressmsg));
                    obj.AddEntry(files{jj}(ii).name, jj);
                    kk = kk+1;
                end
            end

            close(h)
            fprintf('\n');
        end

        % -------------------------------------------------------------------------------
        function err = AddEntry(obj, funcname, jj)
            err = 0;
            if nargin<3
                jj = 1;
            end
            idx = obj.GetIdx(funcname);
            if ~isempty(idx)
                q = menu('Function entry already exists in registry. Do you want to reload it?', 'YES','NO');
                if q==2
                    return;
                end
            end
            tmp = FuncRegEntryClass(funcname);
            if tmp.IsValid()
                idx = length(obj.entries)+1;
                obj.entries(idx) = FuncRegEntryClass(tmp);
                obj.userfuncfiles{idx} = [obj.userfuncdir{jj}, funcname];
            else
                fprintf('  Discarding %s. Function is not valid registry entry.\n', tmp.name)
            end
        end

        % -------------------------------------------------------------------------------
        function err = ReloadEntry(obj, funcname)
            err = 0;
            idx = obj.GetIdx(funcname);
            obj.entries(idx) = FuncRegEntryClass(funcname);
        end

        % ----------------------------------------------------------------------------------
        function Copy(obj, obj2)
            obj.userfuncdir = obj2.userfuncdir;
            obj.userfuncfiles = obj2.userfuncfiles;
            obj.type = obj2.type;
            obj.config = obj2.config;
            for ii=1:length(obj2.entries)
                obj.entries(ii) = FuncRegEntryClass(obj2.entries(ii));
            end
        end

        % ----------------------------------------------------------------------------------
        function b = IsEmpty(obj)
            b = true;
            if isempty(obj)
                return
            end
            if isempty(obj.userfuncdir)
                return;
            end
            if isempty(obj.userfuncfiles)
                return;
            end
            if isempty(obj.entries)
                return;
            end
            b = false;
        end

        % ------------------------------------------------------
        function idx = GetIdx(obj, key)
            idx = [];
            if ~exist('key','var') || isempty(key)
                key = 1;
            end
            if ischar(key)
                k = [];
                for ii=1:length(obj.entries)
                    if strcmp(obj.entries(ii).GetName(), key)
                        k = ii;
                        break;
                    end
                end
                if isempty(k)
                    return
                end
                key = k(1);
            end
            if ~iswholenum(key) || key<1 || key>length(obj.entries)
                return;
            end
            idx = key;
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Methods to retrieve all or multiple entries
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % ----------------------------------------------------------------------------------
        function usagestrs = GetUsageStrs(obj, addnewline)
            if ~exist('addnewline','var')
                addnewline = 0;
            end
            if addnewline
                nl = sprintf('\n');
            else
                nl = '';
            end
            usagestrs = cell(2*length(obj.entries), 1);
            kk = 1;
            for ii=1:length(obj.entries)
                for jj=1:obj.entries(ii).GetOptionsNum()
                    usagestrs{kk} = sprintf('@ %s%s', obj.entries(ii).GetUsageStr(jj), nl);
                    kk=kk+1;
                end
            end
            % Remove extra empty cells
            usagestrs = usagestrs(~cellfun('isempty',usagestrs));
        end

        % ----------------------------------------------------------------------------------
        function names = GetFuncNames(obj)
            names = cell(length(obj.entries), 1);
            for ii=1:length(obj.entries)
                names{ii} = obj.entries(ii).GetName();
            end
        end

        % ----------------------------------------------------------------------------------
        function fcallsstr = GetFuncCallsEncoded(obj, procStream)
            fcallsstr = {};
            if ~isa(fcall,'FuncCallClass')
                return
            end
            if isempty(procStream)
                return;
            end
            fcalls = procStream.fcalls;
            kk=1;
            for iFcall=1:length(fcalls)
                idx = obj.GetIdx(fcall(iFcall).GetName());
                if isempty(idx)
                    continue;
                end
                str = obj.entries(idx).GetFuncCallsEncoded(fcall(iFcall));
                if ~isempty(str)
                    fcallsstr{kk} = str;
                    kk=kk+1;
                end
            end
        end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Methods to retrieve individual entries
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    methods

        % ----------------------------------------------------------------------------------
        function [name, idx] = GetFuncName(obj, key)
            idx = obj.GetIdx(key);
            if isempty(idx)
                return;
            end
            name = obj.entries(idx).GetName();
        end

        % ----------------------------------------------------------------------------------
        function helpstr = GetFuncHelp(obj, key)
            helpstr = '';
            idx = obj.GetIdx(key);
            if isempty(idx)
                return;
            end
            helpstr = obj.entries(idx).GetHelpStr();
        end

        % ----------------------------------------------------------------------------------
        function fcallstr = GetFuncCallStrDecoded(obj, key, usagename)
            fcallstr = '';
            idx = obj.GetIdx(key);
            if isempty(idx)
                return;
            end
            fcallstr = obj.entries(idx).GetFuncCallStrDecoded(usagename);
        end

        % ----------------------------------------------------------------------------------
        function fcall = GetFuncCallDecoded(obj, key, usagename)
            fcall = FuncCallClass().empty();
            idx = obj.GetIdx(key);
            if isempty(idx)
                return;
            end
            fcall = obj.entries(idx).GetFuncCallDecoded(usagename);
        end

        % ----------------------------------------------------------------------------------
        function fcall = FindClosestMatch(obj, fcall0)
            fcall = FuncCallClass().empty();
            if isempty(obj)
                return;
            end
            if nargin<2
                return;
            end
            for ii=1:length(obj.entries)
                fcall = obj.entries(ii).FindClosestMatch(fcall0);
                if ~isempty(fcall)
                    break
                end
            end
        end

        % ----------------------------------------------------------------------------------
        function paramtxt = GetParamText(obj, key)
            paramtxt = '';
            idx = obj.GetIdx(key);
            if isempty(idx)
                return;
            end
            paramtxt = obj.entries(idx).GetParamText();
        end

        % ----------------------------------------------------------------------------------
        function usagename = GetUsageName(obj, fcall)
            usagename = '';
            if ~isa(fcall,'FuncCallClass')
                return
            end
            idx = obj.GetIdx(fcall.GetName());
            if isempty(idx)
                return;
            end
            usagename = obj.entries(idx).GetUsageName(fcall);
        end

        % ----------------------------------------------------------------------------------
        function usagenames = GetUsageNames(obj, funcname)
            usagenames = {};
            idx = obj.GetIdx(funcname);
            if isempty(idx)
                return;
            end
            usagenames = obj.entries(idx).GetUsageNames();
        end

        % ----------------------------------------------------------------------------------
        function n = GetNumUsages(obj, funcname)
            n = length(obj.GetUsageNames(funcname));
        end

        % ----------------------------------------------------------------------------------
        function usagestr = GetUsageStrDecorated(obj, funcname, usageparam)
            usagestr = '';
            if ~exist('funcname','var')
                return;
            end
            if ~exist('usageparam','var')
                usageparam='';
            end
            idx = obj.GetIdx(funcname);
            if isempty(idx)
                return;
            end
            usagestr = obj.entries(idx).GetUsageStrDecorated(usageparam);
        end

        % ----------------------------------------------------------------------------------
        function lastdt = DateLastModified(obj)
            try
                lastdt = datetime(1970,01,01);
                for ii = 1:length(obj.entries)
                    for jj = 1:length(obj.userfuncdir)
                        file = dir([obj.userfuncdir{jj}, '/', obj.entries(ii).GetName(), '.m']);
                        if ~isempty(file)
                            if  datetime(file.date,'locale','system') > lastdt
                                lastdt =  datetime(file.date,'locale','system');
                            end
                            break;
                        end
                    end
                end
            catch
                lastdt = -1;
            end
        end

    end
end

author: BUNPC


text:
function colors = distinguishable_colors(n_colors, bg, func)
% DISTINGUISHABLE_COLORS: pick colors that are maximally perceptually distinct
%
% When plotting a set of lines, you may want to distinguish them by color.
% By default, Matlab chooses a small set of colors and cycles among them,
% and so if you have more than a few lines there will be confusion about
% which line is which. To fix this problem, one would want to be able to
% pick a much larger set of distinct colors, where the number of colors
% equals or exceeds the number of lines you want to plot. Because our
% ability to distinguish among colors has limits, one should choose these
% colors to be "maximally perceptually distinguishable."
%
% This function generates a set of colors which are distinguishable
% by reference to the "Lab" color space, which more closely matches
% human color perception than RGB. Given an initial large list of possible
% colors, it iteratively chooses the entry in the list that is farthest (in
% Lab space) from all previously-chosen entries. While this "greedy"
% algorithm does not yield a global maximum, it is simple and efficient.
% Moreover, the sequence of colors is consistent no matter how many you
% request, which facilitates the users' ability to learn the color order
% and avoids major changes in the appearance of plots when adding or
% removing lines.
%
% Syntax:
%   colors = distinguishable_colors(n_colors)
% Specify the number of colors you want as a scalar, n_colors. This will
% generate an n_colors-by-3 matrix, each row representing an RGB
% color triple. If you don't precisely know how many you will need in
% advance, there is no harm (other than execution time) in specifying
% slightly more than you think you will need.
%
%   colors = distinguishable_colors(n_colors,bg)
% This syntax allows you to specify the background color, to make sure that
% your colors are also distinguishable from the background. Default value
% is white. bg may be specified as an RGB triple or as one of the standard
% "ColorSpec" strings. You can even specify multiple colors:
%     bg = {'w','k'}
% or
%     bg = [1 1 1; 0 0 0]
% will only produce colors that are distinguishable from both white and
% black.
%
%   colors = distinguishable_colors(n_colors,bg,rgb2labfunc)
% By default, distinguishable_colors uses the image processing toolbox's
% color conversion functions makecform and applycform. Alternatively, you
% can supply your own color conversion function.
%
% Example:
%   c = distinguishable_colors(25);
%   figure
%   image(reshape(c,[1 size(c)]))
%
% Example using the file exchange's 'colorspace':
%   func = @(x) colorspace('RGB->Lab',x);
%   c = distinguishable_colors(25,'w',func);

% Copyright 2010-2011 by Timothy E. Holy

% Parse the inputs
if (nargin < 2)
    bg = [1 1 1];  % default white background
else
    if iscell(bg)
        % User specified a list of colors as a cell aray
        bgc = bg;
        for i = 1:length(bgc)
            bgc{i} = parsecolor(bgc{i});
        end
        bg = cat(1,bgc{:});
    else
        % User specified a numeric array of colors (n-by-3)
        bg = parsecolor(bg);
    end
end

% Generate a sizable number of RGB triples. This represents our space of
% possible choices. By starting in RGB space, we ensure that all of the
% colors can be generated by the monitor.
n_grid = 30;  % number of grid divisions along each axis in RGB space
x = linspace(0,1,n_grid);
[R,G,B] = ndgrid(x,x,x);
rgb = [R(:) G(:) B(:)];
if (n_colors > size(rgb,1)/3)
    error('You can''t readily distinguish that many colors');
end

% Convert to Lab color space, which more closely represents human
% perception
if (nargin > 2)
    lab = func(rgb);
    bglab = func(bg);
else
    C = makecform('srgb2lab');
    lab = applycform(rgb,C);
    bglab = applycform(bg,C);
end

% If the user specified multiple background colors, compute distances
% from the candidate colors to the background colors
mindist2 = inf(size(rgb,1),1);
for i = 1:size(bglab,1)-1
    dX = bsxfun(@minus,lab,bglab(i,:)); % displacement all colors from bg
    dist2 = sum(dX.^2,2);  % square distance
    mindist2 = min(dist2,mindist2);  % dist2 to closest previously-chosen color
end

% Iteratively pick the color that maximizes the distance to the nearest
% already-picked color
colors = zeros(n_colors,3);
lastlab = bglab(end,:);   % initialize by making the "previous" color equal to background
for i = 1:n_colors
    dX = bsxfun(@minus,lab,lastlab); % displacement of last from all colors on list
    dist2 = sum(dX.^2,2);  % square distance
    mindist2 = min(dist2,mindist2);  % dist2 to closest previously-chosen color
    [~,index] = max(mindist2);  % find the entry farthest from all previously-chosen colors
    colors(i,:) = rgb(index,:);  % save for output
    lastlab = lab(index,:);  % prepare for next iteration
end

%-------------------------------------------------------------------------
function c = parsecolor(s)
if ischar(s)
    c = colorstr2rgb(s);
elseif isnumeric(s) && size(s,2) == 3
    c = s;
else
    error('MATLAB:InvalidColorSpec','Color specification cannot be parsed.');
end

%-------------------------------------------------------------------------
function c = colorstr2rgb(c)
% Convert a color string to an RGB value.
% This is cribbed from Matlab's whitebg function.
% Why don't they make this a stand-alone function?
rgbspec = [1 0 0;0 1 0;0 0 1;1 1 1;0 1 1;1 0 1;1 1 0;0 0 0];
cspec = 'rgbwcmyk';
k = find(cspec==c(1));
if isempty(k)
    error('MATLAB:InvalidColorString','Unknown color string.');
end
if k~=3 || length(c)==1,
    c = rgbspec(k,:);
elseif length(c)>2,
    if strcmpi(c(1:3),'bla')
        c = [0 0 0];
    elseif strcmpi(c(1:3),'blu')
        c = [0 0 1];
    else
        error('MATLAB:UnknownColorString', 'Unknown color string.');
    end
end


text:
function FinishFortify(s,~,troops)

if get(troops, 'UserData') ~= 0
    uiwait(msgbox('There still are troops on the move!'))
else
    uiresume(get(s, 'Parent'))
end


text:
function this = MrDimInfo_get_add_remove(this)
% Unit test for MrDimInfo get, add and remove method
%
%   Y = MrUnitTest()
%   run(Y, 'MrDimInfo_get_add_remove')
%
% This is a method of class MrUnitTest.
%
% IN
%
% OUT
%
% EXAMPLE
%   MrDimInfo_methods
%
%   See also MrUnitTest

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

% construct MrDimInfo object from sampling points
dimInfo = this.make_dimInfo_reference(0);

% define expected solution
expSolution = dimInfo.copyobj;

% get, remove and add dims

% get dimInfo along x
dimInfoX = dimInfo.get_dims('x');

% remove x from dimInfo
dimInfo.remove_dims('x');

% add x back to dimInfo
dimInfo.add_dims(1, 'samplingPoints', dimInfoX.samplingPoints{1}, ...
    'dimLabels', dimInfoX.dimLabels{1}, 'units', dimInfoX.units{1});

% define actual solution
actSolution = dimInfo;

% verify whether expected and actual solution are identical
% Note: convert to struct, since the PublicPropertyComparator (to allow
% nans to be treated as equal) does not compare properties of objects that
% overload subsref

warning('off', 'MATLAB:structOnObject');
this.verifyEqual(struct(actSolution), struct(expSolution), 'absTol', 10e-7);
warning('on', 'MATLAB:structOnObject');

end

author: translationalneuromodeling

repo: tapas