text:
function m=v_ccwarpf(f,n,s)
%V_CCWARPF  Warp cepstral coefficients M=(F,N,S)
% f(1) is the original sample freq, f(2) is the new sample freq
% n(1) is the original number of coefficients, n(2) is the new number
% s is a string: s(1),s(2) =l for linear, m for mel frequency, use capitals if c0 included

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: v_ccwarpf.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<3
   s='ll';
end
if length(f)<2
   f(2)=1;
end
if length(n)<2
   n(2)=n(1);
end
z=s<'a';
s=s+32*z;
if all(s=='l')
   k=1:n(2)-z(2);
   ff=((1:n(1)).'-z(1))*f(2)/f(1);
   fa=2*sin(ff*pi).*ff/pi;
   fb=ff.^2;
   ka=1-2*rem(k,2);
   kb=k.^2;
   r1=ones(n(1),1);
   c1=ones(1,n(2)-z(2));
   a=fa(:,c1).*ka(r1,:);
   b=fb(:,c1)-kb(r1,:);
   f0=find(fix(ff)==ff);
   if length(f0)
      a(f0,:)=ff(f0,c1)==k(ones(length(f0),1),:);
      b(f0,:)=1;
   end
   m=a./b;
   if z(2)
      m=[[1; 0.5*fa(2:n(1))./fb(2:n(1))] m];
   end
end


text:
% INTLAB  INTerval LABoratory.
% Version 7.1   June-24-2013   966 files, 849 m-files, 32,346 lines of Matlab code w/o comments
%                 (54,324 with comments),
%             supported by a test-case library with more than 60,000 lines of Matlab code w/o comments
%
%===> For license information and terms of use, see README.TXT.
%
%
%Directories and toolboxes
%  intval       - Interval package
%  gradient     - Automatic differentiation package
%  hessian      - Automatic Hessian package
%  taylor       - Automatic Taylor package
%  slope        - Automatic slope package
%  polynom      - Polynomial package (univariate and multivariate)
%  long         - Rudimentary long package
%  utility      - Some useful functions
%  AccSum       - Reference implementations for sum/dot routines
%  demos        - several demo routines
%
%Auxiliary files
%  intlablogo   - INTLAB logo display
%  startintlab  - Initialization of INTLAB
%  startup      - Calls startintlab: adapt to your local installation
%  intlabsetting   - Current setting of INTLAB control variables
%  INTLAB_Version_2     - Additions and changes in version 2
%  INTLAB_Version_3     - Additions and changes in version 3
%  INTLAB_Version_3.1   - Additions and changes in version 3.1
%  INTLAB_Version_4     - Additions and changes in version 4
%  INTLAB_Version_4.1   - Additions and changes in version 4.1
%  INTLAB_Version_4.1.1 - Additions and changes in version 4.1.1
%  INTLAB_Version_4.1.2 - Additions and changes in version 4.1.2
%  INTLAB_Version_5     - Additions and changes in version 5
%  INTLAB_Version_5.1   - Additions and changes in version 5.1
%  INTLAB_Version_5.2   - Additions and changes in version 5.2
%  INTLAB_Version_5.3   - Additions and changes in version 5.3
%  INTLAB_Version_5.4   - Additions and changes in version 5.4
%  INTLAB_Version_5.5   - Additions and changes in version 5.5
%  INTLAB_Version_6     - Additions and changes in version 6
%  INTLAB_Version_7     - Additions and changes in version 7
%  INTLAB_Version_7.1   - Additions and changes in version 7.1
%  FAQ          - Frequently asked questions
%  Readme       - Installation, a little tutorial and miscellaneous
%
%
%INTLAB based on
%  S.M. Rump: INTLAB - INTerval LABoratory, in "Developments in
%    Reliable Computing", ed. T. Csendes, Kluwer Academic Publishers,
%    77-104, 1999.
%
%INTLAB implementation of interval arithmetic is based on
%  S.M. Rump: Fast and Parallel Interval Arithmetic,
%    BIT 39(3), 539-560, 1999.
%and
%  S. Oishi, S.M. Rump: Fast verification of solutions of matrix equations,
%    Numerische Mathametik 90, 755-773, 2002.
%
%Real interval standard functions based on
%  S.M. Rump: Rigorous and portable standard functions,
%    BIT 41(3), 540-56, 2001.
%
%Complex interval standard functions based on
%  N.C. Boersken: Komplexe Kreis-Standardfunktionen, Freiburger
%    Intervallberichte 78/2.
%
%Accurate summation and dot product with specified _precision_
%  T. Ogita, S.M. Rump, and S. Oishi. Accurate Sum and Dot Product,
%     SIAM Journal on Scientific Computing (SISC), 26(6):1955-1988, 2005.
%  S.M. Rump, T. Ogita, and S. Oishi. Fast High Precision Summation.
%     Nonlinear Theory and Its Applications (NOLTA), IEICE, 1(1), 2010.
%  S.M. Rump: Ultimately Fast Accurate Summation, SIAM Journal on
%     Scientific Computing (SISC), 31(5):3466-3502, 2009.
%
%Accurate summation and dot product with specified _accuracy_
%  S.M. Rump, T. Ogita, and S. Oishi: Accurate Floating-point Summation I:
%     Faithful Rounding. SIAM Journal on Scientific Computing (SISC),
%     31(1): 189-224, 2008.
%  S.M. Rump, T. Ogita, and S. Oishi: Accurate Floating-point Summation II:
%     Sign, K-fold Faithful and Rounding to Nearest. SIAM Journal on Scientific
%     Computing (SISC), 31(2):1269-1302, 2008.
%  S.M. Rump: Ultimately Fast Accurate Summation, SIAM Journal on
%     Scientific Computing (SISC), 31(5):3466-3502, 2009
%
%The linear system solver including inner inclusions is completely redesigned and
%now based on
%  S.M. Rump. Accurate solution of dense linear systems, Part II: Algorithms using
%    directed rounding. Journal of Computational and Applied Mathematics (JCAM),
%    242:185–212, 2013.
%
%Least squares problems and underdetermined linear systems are based on
%  S.M. Rump. Verified Bounds for Least Squares Problems and Underdetermined Linear Systems.
%    SIAM J. Matrix Anal. Appl. (SIMAX), 33(1):130–148, 2012.
%and new componentwise error estimates on
%  S.M. Rump: Improved componentwise verified error bounds for least squares problems
%    and underdetermined linear systems, to appear.
%
%For references to verification algorithms cf. the corresponding routines. Many
%of them are discussed in
%  S.M. Rump: Verification methods: Rigorous results using floating-point arithmetic.
%    Acta Numerica, 19:287-449, 2010.
%
%Slopes based on
%  R. Krawzcyk, A. Neumaier: Interval slopes for rational functions
%    and associated centered forms, SIAM J. Numer. Anal. 22, 604-616
%    (1985)
%with improvements based on
%  S.M. Rump: Expansion and Estimation of the Range of Nonlinear Functions,
%    Math. Comp. 65(216), pp. 1503-1512, 1996.
%
%Gradients, Hessians and Taylor expansion based on standard forward mode, cf. for example
%  L.B. Rall: Automatic Differentiation: Techniques and Applications,
%    Lecture Notes in Computer Science 120, Springer, 1981.
%
%Long floating point and interval arithmetic based on standard techniques
%  with adaptation and speed up for interpretative systems.
%

% written  12/30/98     S.M. Rump
% modified 03/06/99     S.M. Rump  Version 2
% modified 11/16/99     S.M. Rump  Version 3
% modified 03/07/02     S.M. Rump  Version 3.1
% modified 12/08/02     S.M. Rump  Version 4
% modified 12/27/02     S.M. Rump  Version 4.1
% modified 01/22/03     S.M. Rump  Version 4.1.1
% modified 11/18/03     S.M. Rump  Version 4.1.2
% modified 04/04/04     S.M. Rump  Version 5
% modified 06/04/05     S.M. Rump  Version 5.1
% modified 12/20/05     S.M. Rump  Version 5.2
% modified 05/26/06     S.M. Rump  Version 5.3
% modified 05/31/07     S.M. Rump  Version 5.4
% modified 11/05/08     S.M. Rump  Version 5.5
% modified 05/08/09     S.M. Rump  Version 6
% modified 12/12/12     S.M. Rump  Version 7
% modified 06/24/13     S.M. Rump  Version 7.1
%

% Copyright (c) Siegfried M. Rump, head of the Institute for Reliable Computing,
%               Hamburg University of Technology
%
% For license information and terms of use, see README.TXT.
%
% :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
% ::  Some references to papers using INTLAB are collected on the
% ::    INTLAB homepage   http://www.ti3.tuhh.de/
% ::  If you have additional references to add, please send me
% ::  a mail ( rump [at] tuhh.de )
%

repo: OpenMAS


text:
%GEOM_REVERT_OBJECT Revert/undo geometry object operations.
%
%   [ GEOM ] = GEOM_REVERT_OBJECT( GEOM, TAG ) Reverts/undos geometry
%   object operations by replacing the geometry object TAG with the
%   objects in its children field (with uniqified tags).

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
%  INTERNAL FUNCTION
%
%  See also: UTILS.CODE.CODE2VECTOR
%

author: jmaih


text:
function jsubs = cfg2jobsubs(item, subs)

% function jsubs = cfg2jobsubs(item, subs)
% Return the subscript into the job tree for a given subscript vector into
% the val part of the cfg tree. In a cfg_choice, this is a struct reference
% to a field with the name of the tag of the corresponding child node.
%
% 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: cfg2jobsubs.m 4867 2012-08-30 13:04:51Z volkmar $

rev = '$Rev: 4867 $'; %#ok

% Only de-reference subscripts into .val{1} field
if isequal(subs(1:2), substruct('.','val','{}',{1}))
    % return possible indices into .val{1} field
    citem = subsref(item, subs(1:2));
    if numel(subs) > 2
        jsubs1 = cfg2jobsubs(citem, subs(3:end));
    else
        jsubs1 = [];
    end;
    jsubs = [substruct('.', gettag(citem)) jsubs1];
else
    cfg_message('matlabbatch:cfg2jobsubs:wrongsubs', 'Inappropriate subscript reference in item ''%s''.', item.tag);
    jsubs = struct('type',{},'subs',{});
end;

author: spm


text:
function priceOut = fcnCalCostwithDegradation(mpcModel, mpcdata, varargin) %this shows the objective funciton
%UPDATE Degradation Cost
    m   = mpcModel.mpciter;
    u = mpcModel.u;
    x = mpcModel.x;
    price = mpcdata.price(1:m,1);
    priceOut = zeros(m,1);

    A = mpcModel.battery.lifeParam(1,1);
    b = mpcModel.battery.lifeParam(1,2);

    coeff = mpcModel.battery.totalprice /(2*A*( mpcModel.battery.capacity ^b));

for i = 1:1:m
    if  u(i,1)>=0
        if u(i,2)*x(i,1)>=0   % x(1) is the cumulative kWh
            cost =  price(i)*u(i,1)  + coeff*( abs(x(i,1)+u(i,2))^b - abs(x(i,1)^b ) );
        else %0.0001*(1/u(2))^2 +
            cost =  price(i)*u(i,1)  + coeff*( abs(u(i,2))^b );
        end
    else
        if u(i,2)*x(i,1)>=0   % x(1) is the cumulative kWh
            cost =  0.8*price(i)*u(i,1)  + coeff*( abs(x(i,1)+u(i,2))^b - abs(x(i,1)^b ) );
        else %0.0001*(1/u(2))^2 +
            cost =  0.8*price(i)*u(i,1)  + coeff*( abs(u(i,2))^b );
        end
    end
    priceOut(i,1) = cost;
end

end

author: juchengquan


text:
function geometry_test0321 ( )

%*****************************************************************************80
%
%% TEST0321 tests HEXAGON_VERTICES_2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  dim_num = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0321\n' );
  fprintf ( 1, '  HEXAGON_VERTICES_2D: the vertices of the unit hexagon.\n' );

  p = hexagon_vertices_2d ( );

  r8mat_transpose_print ( dim_num, 6, p, '  Vertices:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function v = rtd(X)
% Dc RTD output
% R1 R2 R3 R4 R5 R6 R7 R8 R9 RT E1
% 1  2  3  4  5  6  7  8  9  10 11
%
Y=4; % N = 0, U = 4, and M = 1 not required for dc.
%
% The following adds to execution time, but is given for clarity.
%
R1=X(1);R2=X(2);R3=X(3);R4=X(4);R5=X(5);R6=X(6);R7=X(7);
R8=X(8);R9=X(9);RT=X(10);E1=X(11);
%
% Simplify A1 with the following substitutions
%
G1=1/R1+1/R4+1/RT;
G2=1/R5+1/R6+1/RT;
G3=1/R2+1/R3+1/R4;
G4=1/R5+1/R7;
%
A1=[G1 -1/RT -1/R4 0 -1/R1;
      -1/RT G2 -1/R5 0 0;
      -1/R4 0 G3 -1/R3 0;
      0 -1/R5 G4 0 0;
      0 0 0 0 1];
%
B2=[0; E1/R6;E1/R2;0;-E1*R9/R8];
%
V=A1\B2;v=V(Y);
% For dc circuits (no inductors or capacitors), this is as
% far as we have to go.

author: Sable

repo: mcbench-benchmarks


text:
function pass = test_sample( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e2*pref.techPrefs.chebfuneps;

% Example 1
f = ballfun(@(x,y,z)x);
X = sample(f);
Y = zeros(2,3,3);
Y(2,:,2) = [-1,0.5,0.5];
pass(1) = norm( X(:)-Y(:) ) < tol;

% Example 2
f = ballfun(@(x,y,z)z);
X = sample(f, 3, 1, 4);
r = [0;sqrt(0.5);1];
cosT = [1,0.5,-0.5,-1];
Y = reshape(r*cosT,3,1,4);
pass(2) = norm( X(:)-Y(:) ) < tol;

% Example 3
f = ballfun(@(x,y,z)z);
pass(3) = norm(ballfun(sample(f,4,4,4))-f) < tol;

if (nargout > 0)
    pass = all(pass(:));
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function d = setyname(d, n)

%tstoolbox/@description/setyname
%   Syntax:
%     * d = setyunit(d, string)
%
%   Set signal's y-name
%   e.g. d = setyunit(d, 'V')
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

narginchk(2,2);

if isa(n, 'char')
	d.yname = n;
else
	error('Wrong type of argument(s) given');
end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% SPHARM-MAT is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ============================================
% FLD.m
%
% Goal:
%   1. Calculate the optimal FLD projection
%   2. For two class problem
%
% Li Shen
% 11/13/2002 - create

function [FLD_basis, FLD_vals] = FLD(Samples,Labels)

d = size(Samples);
if (d(2)>d(1)-2)
    disp('---------------------------------');
    disp(sprintf('N=%d points, M=%d dims (IGNORE the last %d dims): a nonsigular Sb requires M<=N-c (c=2)',d,d(2)-d(1)+2));
    disp('---------------------------------');
    Samples = Samples(:,1:d(1)-2);
end

FLD_basis = [];

Sb = get_Sb(Samples,Labels);
Sw = get_Sw(Samples,Labels);

[V,D] = eig(Sb,Sw);

[eigval,ind] = sort(diag(D));

FLD_basis = V(:,ind(end));

FLD_vals = Samples*FLD_basis;

return;

%
% calculate between class scatter matrix Sb
%

function Sb = get_Sb(Samples,Labels)

c1_ind = find(Labels==1); c2_ind = find(Labels==2);
m = mean(Samples,1); m1_m = mean(Samples(c1_ind,:),1) - m ; m2_m = mean(Samples(c2_ind,:),1) - m;
Sb = length(c1_ind)*(m1_m'*m1_m) + length(c2_ind)*(m2_m'*m2_m);

return;

%
% calculate within class scatter matrix Sw
%

function Sw = get_Sw(Samples,Labels)

c1_ind = find(Labels==1); c2_ind = find(Labels==2);

% need to have at least 2 points in each class
Sw = cov(Samples(c1_ind,:))*(length(c1_ind)-1) + cov(Samples(c2_ind,:))*(length(c2_ind)-1);

return;

%
% calculate within class scatter matrix Sw (directly according to the definition)
%

function Sw = get_Sw_v2(Samples,Labels)

c1_ind = find(Labels==1); c2_ind = find(Labels==2);

% to verify the correctness of Sw
Sw = zeros(size(Samples,2));
m1 = mean(Samples(c1_ind,:),1);
xs1 = Samples(c1_ind,:) - m1(ones(1,length(c1_ind)),:);
for i = 1:size(xs1,1)
    Sw = Sw + xs1(i,:)'*xs1(i,:);
end

m2 = mean(Samples(c2_ind,:),1);
xs2 = Samples(c2_ind,:) - m2(ones(1,length(c2_ind)),:);
for i = 1:size(xs2,1)
    Sw = Sw + xs2(i,:)'*xs2(i,:);
end

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%% DEMO_febio_0019_vessel_pressure_inflate
% Below is a demonstration for:
%
% * Building geometry for a cylindrical vessel with tetrahedral elements
% * Defining the boundary conditions
% * Coding the febio structure
% * Running the model
% * Importing and visualizing the displacement results

%% Keywords
%
% * febio_spec version 3.0
% * febio, FEBio
% * vessel, cylinder
% * prescribed pressure
% * tetrahedral elements, tet4
% * tube, cylindrical
% * static, solid
% * hyperelastic, Ogden
% * displacement logfile
% * stress logfile

%%

clear; close all; clc;

%% Plot settings
fontSize=20;
faceAlpha1=0.8;
markerSize=40;
markerSize2=20;
lineWidth=3;

%% Control parameters

% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');

% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=[febioFebFileNamePart,'.txt']; %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement
febioLogFileName_stress_prin=[febioFebFileNamePart,'_stress_prin_out.txt']; %Log file name for exporting principal stress

%Specifying geometry parameters
pointSpacing=1;

radiusInner1=9;
radiusInner2=10;

radiusOuter1=10;
radiusOuter2=12;

vesselLength=50;

%Load
appliedPressure=0.1; %MPa

%Material parameter set
c1=1; %Shear-modulus-like parameter
m1=2; %Material parameter setting degree of non-linearity
k_factor=1e2; %Bulk modulus factor
k=c1*k_factor; %Bulk modulus

% FEA control settings
numTimeSteps=10; %Number of time steps desired
max_refs=25; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=12; %Optimum number of iterations
max_retries=10; %Maximum number of retires
dtmin=(1/numTimeSteps)/100; %Minimum time step size
dtmax=1/numTimeSteps; %Maximum time step size
min_residual=1e-20;

%% Creating model boundary polygons
%

nRad=round((2*pi*mean([radiusInner1 radiusInner2]))/pointSpacing); %Number of radial steps

t=linspace(0,2*pi,nRad)'; %Angles
t=t(1:end-1); %take away last which equals start
v1_Inner=[-(vesselLength/2)*ones(size(t)) radiusInner1*sin(t) radiusInner1*cos(t)]; %Circular coordinates

t=linspace(0,2*pi,nRad)'; %Angles
t=t(1:end-1); %take away last which equals start
v2_Inner=[(vesselLength/2)*ones(size(t)) radiusInner2*sin(t) radiusInner2*cos(t)]; %Circular coordinates

%% Creating model boundary surfaces

% controlStructLoft.numSteps=17;
controlStructLoft.closeLoopOpt=1;
controlStructLoft.patchType='quad';

%Meshing outer surface
[F1,V1]=polyLoftLinear(v1_Inner,v2_Inner,controlStructLoft);
F1=fliplr(F1); %Invert orientation

%%

outerRadii=V1(:,1); %x
outerRadii=outerRadii-min(outerRadii(:)); %[0 - ...]
outerRadii=outerRadii./max(outerRadii(:)); %[0 - 1]
outerRadii=radiusOuter1+(outerRadii.*(radiusOuter2-radiusOuter1)); %[0 - 1]

innerRadii=V1(:,1); %x
innerRadii=innerRadii-min(innerRadii(:)); %[0 - ...]
innerRadii=innerRadii./max(innerRadii(:)); %[0 - 1]
innerRadii=radiusInner1+(innerRadii.*(radiusInner2-radiusInner1)); %[0 - 1]

wallThickness=outerRadii-innerRadii;

%%
% Plotting model boundary polygons

cFigure; hold on;
title('Inner surface and polygons','FontSize',fontSize);
gpatch(F1,V1,outerRadii);
plotV(v1_Inner,'g.-','LineWidth',3);
plotV(v2_Inner,'g.-','LineWidth',3);
axisGeom(gca,fontSize); camlight headlight;
colormap gjet; colorbar;
gdrawnow;

%%

numSteps=ceil(max(wallThickness)./pointSpacing);
[E,V,Fp1,Fp2]=patchThick(F1,V1,1,wallThickness,numSteps);

[F,~,C_type]=element2patch(E,[],'hex8');
indBoundary=tesBoundary(F);
Fb=F(indBoundary,:);
Cb=C_type(indBoundary,:);

%%
% Plotting model boundary surfaces

cFigure;
hold on;
title('Model boundary surfaces','FontSize',fontSize);

gpatch(Fb,V,Cb);

axisGeom(gca,fontSize); camlight headlight;
colormap gjet; icolorbar;
drawnow;

%% Defining the boundary conditions
% The visualization of the model boundary shows colors for each side of the
% cube. These labels can be used to define boundary conditions.

%Define supported node set
bcSupportList=unique(Fb(ismember(Cb,[3 4]),:)); %Node set part of selected face

F_pressure=Fb(Cb==1,:);

%%
% Visualizing boundary conditions. Markers plotted on the semi-transparent
% model denote the nodes in the various boundary condition lists.

hf=cFigure;
title('Boundary conditions','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

gpatch(Fb,V,'kw','none',0.5);

hl(1)=plotV(V(bcSupportList,:),'k.','MarkerSize',markerSize);
hl(2)=gpatch(F_pressure,V,'rw','r',1);
patchNormPlot(F_pressure,V);

legend(hl,{'BC support','Pressure surface'});

axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Defining the FEBio input structure
% See also |febioStructTemplate| and |febioStruct2xml| and the FEBio user
% manual.

%Get a template with default settings
[febio_spec]=febioStructTemplate;

%febio_spec version
febio_spec.ATTR.version='3.0';

%Module section
febio_spec.Module.ATTR.type='solid';

%Control section
febio_spec.Control.analysis='STATIC';
febio_spec.Control.time_steps=numTimeSteps;
febio_spec.Control.step_size=1/numTimeSteps;
febio_spec.Control.solver.max_refs=max_refs;
febio_spec.Control.solver.max_ups=max_ups;
febio_spec.Control.solver.min_residual=min_residual;
febio_spec.Control.time_stepper.dtmin=dtmin;
febio_spec.Control.time_stepper.dtmax=dtmax;
febio_spec.Control.time_stepper.max_retries=max_retries;
febio_spec.Control.time_stepper.opt_iter=opt_iter;

%Material section
materialName1='Material1';
febio_spec.Material.material{1}.ATTR.name=materialName1;
febio_spec.Material.material{1}.ATTR.type='Ogden';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.c1=c1;
febio_spec.Material.material{1}.m1=m1;
febio_spec.Material.material{1}.c2=c1;
febio_spec.Material.material{1}.m2=-m1;
febio_spec.Material.material{1}.k=k;

%Mesh section
% -> Nodes
febio_spec.Mesh.Nodes{1}.ATTR.name='nodeSet_all'; %The node set name
febio_spec.Mesh.Nodes{1}.node.ATTR.id=(1:size(V,1))'; %The node id's
febio_spec.Mesh.Nodes{1}.node.VAL=V; %The nodel coordinates

% -> Elements
partName1='Part1';
febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of this part
febio_spec.Mesh.Elements{1}.ATTR.type='hex8'; %Element type
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E; %The element matrix

% -> Surfaces
surfaceName1='LoadedSurface';
febio_spec.Mesh.Surface{1}.ATTR.name=surfaceName1;
febio_spec.Mesh.Surface{1}.tri3.ATTR.id=(1:1:size(F_pressure,1))';
febio_spec.Mesh.Surface{1}.tri3.VAL=F_pressure;

% -> NodeSets
nodeSetName1='bcSupportList';
febio_spec.Mesh.NodeSet{1}.ATTR.name=nodeSetName1;
febio_spec.Mesh.NodeSet{1}.node.ATTR.id=bcSupportList(:);

%MeshDomains section
febio_spec.MeshDomains.SolidDomain.ATTR.name=partName1;
febio_spec.MeshDomains.SolidDomain.ATTR.mat=materialName1;

%Boundary condition section
% -> Fix boundary conditions
febio_spec.Boundary.bc{1}.ATTR.type='fix';
febio_spec.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{1}.dofs='x,y,z';

%Loads section
% -> Surface load
febio_spec.Loads.surface_load{1}.ATTR.type='pressure';
febio_spec.Loads.surface_load{1}.ATTR.surface=surfaceName1;
febio_spec.Loads.surface_load{1}.pressure.ATTR.lc=1;
febio_spec.Loads.surface_load{1}.pressure.VAL=appliedPressure;
febio_spec.Loads.surface_load{1}.symmetric_stiffness=1;

%LoadData section
% -> load_controller
febio_spec.LoadData.load_controller{1}.ATTR.id=1;
febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{1}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{1}.points.point.VAL=[0 0; 1 1];

%Output section
% -> log file
febio_spec.Output.logfile.ATTR.file=febioLogFileName;
febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.node_data{1}.VAL=1:size(V,1);

febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_stress_prin;
febio_spec.Output.logfile.element_data{1}.ATTR.data='s1;s2;s3';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';

%% Quick viewing of the FEBio input file structure
% The |febView| function can be used to view the xml structure in a MATLAB
% figure window.

%%
% |febView(febio_spec); %Viewing the febio file|
% febView(febio_spec)

%% Exporting the FEBio input file
% Exporting the febio_spec structure to an FEBio input file is done using
% the |febioStruct2xml| function.

febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode

%% Running the FEBio analysis
% To run the analysis defined by the created FEBio input file the
% |runMonitorFEBio| function is used. The input for this function is a
% structure defining job settings e.g. the FEBio input file name. The
% optional output runFlag informs the user if the analysis was run
% succesfully.

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.disp_log_on=1; %Display convergence information in the command window
febioAnalysis.runMode='external';%'internal';
febioAnalysis.t_check=0.25; %Time for checking log file (dont set too small)
febioAnalysis.maxtpi=1e99; %Max analysis time
febioAnalysis.maxLogCheckTime=10; %Max log file checking time

[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!

%% Import FEBio results

if runFlag==1 %i.e. a succesful run

    %%
    % Importing nodal displacements from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp),1,1);

    %Access data
    N_disp_mat=dataStruct.data; %Displacement
    timeVec=dataStruct.time; %Time

    %Create deformed coordinate set
    V_DEF=N_disp_mat+repmat(V,[1 1 size(N_disp_mat,3)]);

    %%
    % Importing element principal stresses from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_stress_prin),1,1);

    %Access data
    E_stress_prin_mat=dataStruct.data;
    S1_mat=E_stress_prin_mat(:,1,:);
    S2_mat=E_stress_prin_mat(:,2,:);
    S3_mat=E_stress_prin_mat(:,3,:);
    S_vm = sqrt(((S1_mat-S2_mat).^2+(S2_mat-S3_mat).^2+(S3_mat-S1_mat).^2)./2);

    %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    [~,CF_S_vm,~]=element2patch(E,S_vm(:,:,end),'hex8');
    Cb_S_vm=CF_S_vm(indBoundary,:);

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    title('Von Mises stres [MPa]','Interpreter','Latex')
    hp=gpatch(Fb,V_DEF(:,:,end),Cb_S_vm,'k',1); %Add graphics object to animate

    axisGeom(gca,fontSize);
    colormap(gjet(250)); colorbar;
    caxis([0 max(S_vm(:))]);
    axis(axisLim(V_DEF)); %Set axis limits statically
    camlight headlight;

    % Set up animation features
    animStruct.Time=timeVec; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments
        [~,CF_S_vm,~]=element2patch(E,S_vm(:,:,qt),'hex8');
        Cb_S_vm=CF_S_vm(indBoundary,:);

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp hp]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData'}; %Properties of objects to animate
        animStruct.Set{qt}={V_DEF(:,:,qt),Cb_S_vm}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;

    %%
    % Importing element principal stresses from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_stress_prin),1,1);

    %Access data
    E_stress_prin_mat=dataStruct.data;

    time_vec=dataStruct.time;

end

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function []=createCheckerBoardImage(Nrows,Ncols,squareSize,resolution,path)
%% function for writing a checker board image to be printed for camera clibration in step0.
% The function creates the image and saves it in the desired path.

% INPUT:
% * Nrows: number of rows (uneven)
% * Ncols: number of columns (even)
% * squareSize: size of each square (in meters)
% * resolution: image resolution (pixels per meter)
% * path: path where the image will be saves (including file name)

%%
squareSizePixel=squareSize*resolution;
CB = checkerboardBW(squareSizePixel,Nrows,Ncols);
figure; imshow(CB);
imwrite(CB,[path '\CB_' num2str(Nrows) '_' num2str(Nrows) '_' num2str(squareSize*1000) '.png'],'png','ResolutionUnit','meter','XResolution',resolution); %XResolution=pixels per meter

end

 %%
% MultiDIC: a MATLAB Toolbox for Multi-View 3D Digital Image Correlation
%
% License: <https://github.com/MultiDIC/MultiDIC/blob/master/LICENSE.txt>
%
% Copyright (C) 2018  Dana Solav
%
% If you use the toolbox/function for your research, please cite our paper:
% <https://engrxiv.org/fv47e>

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function p = parse_filename(filename,fields)
% Split a filename into its building constituents
% FORMAT p = parse_filename(filename,fields)
%
% Example:
%
% >> filename = '../sub-16/anat/sub-16_ses-mri_run-1_echo-2_FLASH.nii.gz';
% >> parse_filename(filename)
%
% ans =
%
%   struct with fields:
%
%     filename: 'sub-16_ses-mri_run-1_echo-2_FLASH.nii.gz'
%         type: 'FLASH'
%          ext: '.nii.gz'
%          sub: '16'
%          ses: 'mri'
%          run: '1'
%         echo: '2'
%__________________________________________________________________________

% Copyright (C) 2016-2018, Guillaume Flandin, Wellcome Centre for Human Neuroimaging
% Copyright (C) 2018--, BIDS-MATLAB developers

filename = file_utils(filename,'filename');
[parts, dummy] = regexp(filename,'(?:_)+','split','match');
p.filename = filename;
[p.type, p.ext] = strtok(parts{end},'.');
for i=1:numel(parts)-1
    [d, dummy] = regexp(parts{i},'(?:\-)+','split','match');
    if length(d)<2
%        warning([filename ' cannot be parsed. ''' d{1} ''' is not associated with a label (' d{1} '-)']);
        continue
    end
    p.(d{1}) = d{2};
end
if nargin == 2
    for i=1:numel(fields)
        if ~isfield(p,fields{i})
            p.(fields{i}) = '';
        end
    end
    try
        p = orderfields(p,['filename','ext','type',fields]);
    catch
        warning('Ignoring file "%s" not matching template.',filename);
        p = struct([]);
    end
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [MSC]=coherence_MVDR(x1,x2,L,K);

%% This program computes the coherence function between 2 signals
%% x1 and x2 with the MVDR method.
%% This algorithm is based on the paper by the same authors:
%% J. Benesty, J. Chen, and Y. Huang, "A generalized MVDR spectrum,"
%% IEEE Signal Processing letters, vol. 12, pp. 827-830, Dec. 2005.

%% x1, first signal vector of length n
%% x2, second signal vector of length n
%% L is the length of MVDR filter or window length
%% K is the resolution (the higher K, the better the resolution)

%initialization
xx1     = zeros(L,1);
xx2     = zeros(L,1);
r11     = zeros(L,1);
r22     = zeros(L,1);
r12     = zeros(L,1);
r21     = zeros(L,1);

%construction of the Fourier Matrix
F       = zeros(L,K);
l       = [0:L-1]';
f       = exp(2*pi*l*j/K);
for k = 0:K-1
    F(:,k+1) = f.^k;
end
F       = F/sqrt(L);

%number of samples, equal to the lenght of x1 and x2
n       = length(x1);

for i = 1:n
    xx1 = [x1(i);xx1(1:L-1)];
    xx2 = [x2(i);xx2(1:L-1)];
    r11 = r11 + xx1*conj(xx1(1));
    r22 = r22 + xx2*conj(xx2(1));
    r12 = r12 + xx1*conj(xx2(1));
    r21 = r21 + xx2*conj(xx1(1));
end
r11 = r11/n;
r22 = r22/n;
r12 = r12/n;
r21 = r21/n;
%
R11 = toeplitz(r11);
R22 = toeplitz(r22);
R12 = toeplitz(r12,conj(r21));
%
%for regularization
Dt1     = 0.01*r11(1)*diag(diag(ones(L)));
Dt2     = 0.01*r22(1)*diag(diag(ones(L)));
%
Ri11    = inv(R11 + Dt1);
Ri22    = inv(R22 + Dt2);
Rn12    = Ri11*R12*Ri22;
%
Si11    = real(diag(F'*Ri11*F));
Si22    = real(diag(F'*Ri22*F));
S12     = diag(F'*Rn12*F);
%
%Magnitude squared coherence function
MSC     = real(S12.*conj(S12))./(Si11.*Si22);


text:
function x=othertriangle(y);

if y(1,1)==1;
    warning('doesnt work for DIS matrices');
end

x=y*0;
if ndims(y)==2;
    x=y+rot90(flipud(triu(y)),3);
end

if ndims(y)==3;
    for n=1:size(y,3);
        y1=y(:,:,n);
        x(:,:,n)=y1+rot90(flipud(triu(y1)),3);
    end
end

if ndims(y)==4;
    for n=1:size(y,3);
        for n1=1:size(y,4);
            y1=y(:,:,n,n1);
            x(:,:,n,n1)=y1+rot90(flipud(triu(y1)),3);
        end
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
clear all; close all; clc

addpath('apm')

s = 'http://byu.apmonitor.com';
c = 'mpc';

% Connect to Arduino
tclab;

% Run time in minutes
run_time = 5.0;

% Number of cycles (1 cycle per 3 seconds)
loops = round(20*run_time);

% Temperature (degC)
T1 = ones(1,loops) * T1C(); % measured T
Tsp1 = ones(1,loops) * 30;
T2 = ones(1,loops) * T2C(); % measured T
Tsp2 = ones(1,loops) * 23;
time = zeros(1,loops);

% Changes in set point
Tsp1(50:end) = 40;
Tsp2(10:end) = 35;

% milli-volts input
Q1 = zeros(1,loops);
Q2 = zeros(1,loops);

% model predictive control initialization
mpc_init(s,c);

start_time = clock;
prev_time = start_time;

% dynamic plot (note: subplots needs to be declared here first)
figure(1)
subplot(3,1,1)
hold on, grid on
anexp1 = animatedline('LineStyle','-', 'Color', 'k', 'LineWidth', 2);
ansp1 = animatedline('LineStyle','--','Color', 'r','LineWidth', 2);
ylabel('Temperature \circC')
legend('T_1 Measured', 'T_1 Set Point', ...
    'Location', 'northwest')
subplot(3,1,2)
hold on, grid on
anexp2 = animatedline('LineStyle','-', 'Color', 'k', 'LineWidth', 2);
ansp2 = animatedline('LineStyle','--','Color', 'r','LineWidth', 2);
ylabel('Temperature \circC')
legend('T_2 Measured', 'T_2 Set Point', ...
    'Location', 'northwest')
subplot(3,1,3)
hold on, grid on
anQ1 = animatedline('LineStyle','-', 'Color', 'k', 'LineWidth', 2);
anQ2 = animatedline('LineStyle','--', 'Color', 'b', 'LineWidth', 2);
ylabel('Power Level Q (%)')
legend('Q_1', 'Q_2', 'Location', 'northwest')
xlabel('Time (sec)')

for ii = 1:loops
    % Pause Sleep time
    pause_max = 3.0;
    pause_time = pause_max - etime(clock,prev_time);
    if pause_time >= 0.0
        pause(pause_time - 0.01)
    else
        pause(0.01)
    end

    % Record time and change in time
    t = clock;
    dt = etime(t,prev_time);
    if ii>=2
        time(ii) = time(ii-1) + dt;
    end
    prev_time = t;

    % read and record from temperature controller
    T1(ii) = T1C();
    T2(ii) = T2C();

    % model predictive control
    [Q1(ii), Q2(ii)] = mpc2(T1(ii),Tsp1(ii),T2(ii),Tsp2(ii));

    % adjust power level
    h1(Q1(ii));
    h2(Q2(ii));

    % plot
    addpoints(anexp1,time(ii),T1(ii))
    addpoints(ansp1,time(ii),Tsp1(ii))
    addpoints(anexp2,time(ii),T2(ii))
    addpoints(ansp2,time(ii),Tsp2(ii))
    addpoints(anQ1,time(ii),Q1(ii))
    addpoints(anQ2,time(ii),Q2(ii))
    drawnow

    if ii==10
        apm_web(s,c);
    end
end

h1(0);
h2(0);
disp('Heaters off')
% turn off heater but keep LED on if T > 50
if (T1C() || T2C()) > 50
    led(1)
    disp(['Warning, heater temperature 1 =', num2str(T1C())])
    disp(['Warning, heater temperature 2 =', num2str(T2C())])
else
    led(0)
end

% save txt file with data
data = [time',Q1',Q2',T1',T2',Tsp1',Tsp2'];
csvwrite('data.txt',data);

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function slice3(T,varargin)
%SLICE3 Visualize a third-order tensor with slices.
%   slice3(T) visualizes the third-order tensor T by drawing its mode-1,
%   -2, and -3 slices using sliders to define their respective indices.
%   Press 'h' to show/hide the figure's controls.
%
%   slice3(T,varargin) passes the parameters varargin{:} to the plot.

%   Authors: Laurent Sorber (Laurent.Sorber@cs.kuleuven.be)
%            Marc Van Barel (Marc.VanBarel@cs.kuleuven.be)
%            Lieven De Lathauwer (Lieven.DeLathauwer@kuleuven-kulak.be)

% Check the dimensions.
if ndims(T) <= 2
    imagesc(T,varargin{:}); return;
elseif ndims(T) ~= 3
    error('slice3:T','ndims(T) should be 3.');
end

% If the data is complex, convert it to the modulus.
% Remove any NaNs and Infs.
if any(~isreal(T(:))), T = abs(T); end
T(isnan(T) | (isinf(T) & T < 0)) = min(T(~isinf(T(:))));
T(isinf(T) & T > 0) = max(T(~isinf(T(:))));

% Set up the figure.
idx = size(T); idx(2) = 1;
T = double(permute(T,[2 3 1]));
cax = newplot;
set(gcf,'Toolbar','figure');
set(datacursormode(gcf),'UpdateFcn',@datacursor);
zoom off; pan off; rotate3d off; datacursormode off;
obj{1} = uicontrol('Style','text','Position',[5 45 15 15],'String','i');
obj{2} = uicontrol('Style','slider','Position',[25 45 120 15], ...
                   'Min',1,'Max',size(T,3),'Value',size(T,3), ...
                   'SliderStep',1/(size(T,3)-1)*[1 1], ...
                   'Callback',{@redraw,1}, ...
                   'KeyPressFcn',@(obj,evt)toggle(evt.Key));
obj{3} = uicontrol('Style','text','Position',[5 25 15 15],'String','j');
obj{4} = uicontrol('Style','slider','Position',[25 25 120 15], ...
                   'Min',1,'Max',size(T,1),'Value',1, ...
                   'SliderStep',1/(size(T,1)-1)*[1 1], ...
                   'Callback',{@redraw,2}, ...
                   'KeyPressFcn',@(obj,evt)toggle(evt.Key));
obj{5} = uicontrol('Style','text','Position',[5 5 15 15],'String','k');
obj{6} = uicontrol('Style','slider','Position',[25 5 120 15], ...
                   'Min',1,'Max',size(T,2),'Value',size(T,2), ...
                   'SliderStep',1/(size(T,2)-1)*[1 1], ...
                   'Callback',{@redraw,3}, ...
                   'KeyPressFcn',@(obj,evt)toggle(evt.Key));
isVisible = true;
set(gcf,'KeyPressFcn',@(obj,evt)toggle(evt.Key));
redraw();

function toggle(key)
    if strcmpi(key,'h')
        if isVisible, val = 'off'; else val = 'on'; end
        for i = 1:length(obj), set(obj{i},'Visible',val); end
        isVisible = ~isVisible;
    end
end

function txt = datacursor(~,event_obj)
    pos = get(event_obj,'Position');
    txt = {['I: ' int2str(pos(3))], ...
           ['J: ' int2str(pos(2))], ...
           ['K: ' int2str(pos(1))], ...
           ['Value: ' num2str(T(pos(2),pos(1),pos(3)))]};
end

function redraw(hobj,~,ax)

    % Draw slices.
    if nargin >= 1, idx(ax) = round(get(hobj,'Value')); end
    slice(T,idx(3),idx(2),idx(1),varargin{:});
    shading flat;

    % Set axis properties.
    xlabel('k');
    ylabel('j');
    zlabel('i');
    xlim([1 size(T,2)]);
    ylim([1 size(T,1)]);
    zlim([1 size(T,3)]);
    set(cax,'YDir','reverse');
    set(cax,'ZDir','reverse');

    % Display grid.
    step = max(1,round(size(T)/6));
    set(cax,'XTickMode','manual','YTickMode','manual','ZTickMode','manual');
    set(cax,'XTick',[1:step(2):size(T,2)-1 size(T,2)]);
    set(cax,'YTick',[1:step(1):size(T,1)-1 size(T,1)]);
    set(cax,'ZTick',[1:step(3):size(T,3)-1 size(T,3)]);
    grid on;

end

end

author: andrewssobral

repo: lrslibrary


text:
function [V,nr,nre]=lcon2vert(A,b,Aeq,beq,TOL,checkbounds)
%An extension of Michael Kleder's con2vert function, used for finding the
%vertices of a bounded polyhedron in R^n, given its representation as a set
%of linear constraints. This wrapper extends the capabilities of con2vert to
%also handle cases where the  polyhedron has zero volume in R^n, i.e., where the
%polyhedron is defined by both equality and inequality constraints.
%
%SYNTAX:
%
%  [V,nr,nre]=lcon2vert(A,b,Aeq,beq,TOL)
%
%The rows of the N x n matrix V are a series of N vertices of the polyhedron
%in R^n, defined by the linear constraints
%
%   A*x  <= b
%   Aeq*x = beq
%
%By default, Aeq=beq=[], implying no equality constraints. The output "nr"
%lists non-redundant inequality constraints, and "nre" lists non-redundant
%equality constraints.
%
%The optional TOL argument is a tolerance used for both rank-estimation and
%for testing feasibility of the equality constraints. Default=1e-10.
%The default can also be obtained by passing TOL=[];
%
%NOTE: It is important that the region specified by the inequality system A*x<=b
%have non-zero volume in R^n. For example, A=b=[1;-1] is not legal input data,
%because the only solution to A*x<=b is x=1, which has zero volume in R^1. The
%proper way to express a zero-volume region is with the addition of
%equality constraint data, as for example Aeq=1, beq=1, A=1,b=100.
%
%EXAMPLE:
%
%The 3D region defined by x+y+z=1, x>=0, y>=0, z>=0
%is described by the following constraint data.
%
%
%     A =
%
%         0.4082   -0.8165    0.4082
%         0.4082    0.4082   -0.8165
%        -0.8165    0.4082    0.4082
%
%
%     b =
%
%         0.4082
%         0.4082
%         0.4082
%
%
%     Aeq =
%
%         0.5774    0.5774    0.5774
%
%
%     beq =
%
%         0.5774
%
%
%  >> V=lcon2vert(A,b,Aeq,beq)
%
%         V =
%
%             1.0000    0.0000    0.0000
%             0.0000    0.0000    1.0000
%            -0.0000    1.0000    0.0000
%
%

  %%initial argument parsing

  nre=[];
  nr=[];
  if nargin<5 || isempty(TOL), TOL=1e-10; end
  if nargin<6, checkbounds=true; end

  switch nargin

      case 0

           error 'At least 1 input argument required'

      case 1

         b=[]; Aeq=[]; beq=[];

      case 2

          Aeq=[]; beq=[];

      case 3

          beq=[];
          error 'Since argument Aeq specified, beq must also be specified'

  end

  b=b(:); beq=beq(:);

  if xor(isempty(A), isempty(b))
     error 'Since argument A specified, b must also be specified'
  end

  if xor(isempty(Aeq), isempty(beq))
        error 'Since argument Aeq specified, beq must also be specified'
  end

  nn=max(size(A,2)*~isempty(A),size(Aeq,2)*~isempty(Aeq));

  if ~isempty(A) && ~isempty(Aeq) && ( size(A,2)~=nn || size(Aeq,2)~=nn)

      error 'A and Aeq must have the same number of columns if both non-empty'

  end

  inequalityConstrained=~isempty(A);
  equalityConstrained=~isempty(Aeq);

 [A,b]=rownormalize(A,b);
 [Aeq,beq]=rownormalize(Aeq,beq);

  if equalityConstrained && nargout>2

        [discard,nre]=lindep([Aeq,beq].',TOL);  %#ok<ASGLU>

        if ~isempty(nre) %reduce the equality constraints

            Aeq=Aeq(nre,:);
            beq=beq(nre);

        else
            equalityConstrained=false;
        end

   end

   %%Find 1 solution to equality constraints within tolerance

   if equalityConstrained

       Neq=null(Aeq);

       x0=pinv(Aeq)*beq;

       if norm(Aeq*x0-beq)>TOL*norm(beq)  %infeasible

          nre=[]; nr=[]; %All constraints redundant for empty polytopes
          V=[];
          return;

       elseif isempty(Neq)

           if inequalityConstrained && ~all(A*x0<=b)

              nre=[]; nr=[]; %All constraints redundant for empty polytopes
              V=[];
              return;

           else %inequality constraints all satisfied, including vacuously

               V=x0(:).';
               nre=(1:nn).'; %Equality constraints determine everything.
               nr=[];%All inequality constraints are therefore redundant.
               return

           end

       end

       %rkAeq= nn - size(Neq,2);

  end

    %%
  if inequalityConstrained && equalityConstrained

   AAA=A*Neq;
   bbb=b-A*x0;

  elseif inequalityConstrained

    AAA=A;
    bbb=b;

  elseif equalityConstrained && ~inequalityConstrained

       error('Non-bounding constraints detected. (Consider box constraints on variables.)')

  end

  nnn=size(AAA,2);

  if nnn==1 %Special case

     idxu=sign(AAA)==1;
     idxl=sign(AAA)==-1;
     idx0=sign(AAA)==0;

     Q=bbb./AAA;
     U=Q;
       U(~idxu)=inf;
     L=Q;
       L(~idxl)=-inf;

     [ub,uloc]=min(U);
     [lb,lloc]=max(L);

     if ~all(bbb(idx0)>=0) || ub<lb %infeasible

         V=[]; nr=[]; nre=[];
         return

     elseif ~isfinite(ub) || ~isfinite(lb)

         error('Non-bounding constraints detected. (Consider box constraints on variables.)')

     end

     Zt=[lb;ub];

     if nargout>1
        nr=unique([lloc,uloc]); nr=nr(:);
     end

  else

          if nargout>1
           [Zt,nr]=con2vert(AAA,bbb,TOL,checkbounds);
          else
            Zt=con2vert(AAA,bbb,TOL,checkbounds);
          end

  end

  if equalityConstrained && ~isempty(Zt)

      V=bsxfun(@plus,Zt*Neq.',x0(:).');

  else

      V=Zt;

  end

  if isempty(V)
     nr=[]; nre=[];
  end

 function [V,nr] = con2vert(A,b,TOL,checkbounds)
% CON2VERT - convert a convex set of constraint inequalities into the set
%            of vertices at the intersections of those inequalities;i.e.,
%            solve the "vertex enumeration" problem. Additionally,
%            identify redundant entries in the list of inequalities.
%
% V = con2vert(A,b)
% [V,nr] = con2vert(A,b)
%
% Converts the polytope (convex polygon, polyhedron, etc.) defined by the
% system of inequalities A*x <= b into a list of vertices V. Each ROW
% of V is a vertex. For n variables:
% A = m x n matrix, where m >= n (m constraints, n variables)
% b = m x 1 vector (m constraints)
% V = p x n matrix (p vertices, n variables)
% nr = list of the rows in A which are NOT redundant constraints
%
% NOTES: (1) This program employs a primal-dual polytope method.
%        (2) In dimensions higher than 2, redundant vertices can
%            appear using this method. This program detects redundancies
%            at up to 6 digits of precision, then returns the
%            unique vertices.
%        (3) Non-bounding constraints give erroneous results; therefore,
%            the program detects non-bounding constraints and returns
%            an error. You may wish to implement large "box" constraints
%            on your variables if you need to induce bounding. For example,
%            if x is a person's height in feet, the box constraint
%            -1 <= x <= 1000 would be a reasonable choice to induce
%            boundedness, since no possible solution for x would be
%            prohibited by the bounding box.
%        (4) This program requires that the feasible region have some
%            finite extent in all dimensions. For example, the feasible
%            region cannot be a line segment in 2-D space, or a plane
%            in 3-D space.
%        (5) At least two dimensions are required.
%        (6) See companion function VERT2CON.
%        (7) ver 1.0: initial version, June 2005
%        (8) ver 1.1: enhanced redundancy checks, July 2005
%        (9) Written by Michael Kleder
%
%Modified by Matt Jacobson - March 30, 2011
%

   %%%3/4/2012 Improved boundedness test - unfortunately slower than Michael Kleder's
   if checkbounds

    [aa,bb,aaeq,bbeq]=vert2lcon(A,TOL);

    if any(bb<=0) || ~isempty(bbeq)
        error('Non-bounding constraints detected. (Consider box constraints on variables.)')
    end

    clear aa bb aaeq bbeq

   end

   dim=size(A,2);

   %%%Matt J initialization
   if strictinpoly(b,TOL)

       c=zeros(dim,1);

   else

            slackfun=@(c)b-A*c;

            %Initializer0
            c = pinv(A)*b; %02/17/2012 -replaced with pinv()
            s=slackfun(c);

            if ~approxinpoly(s,TOL) %Initializer1

                c=Initializer1(TOL,A,b,c);
                s=slackfun(c);

            end

            if  ~approxinpoly(s,TOL)  %Attempt refinement

                %disp 'It is unusually difficult to find an interior point of your polytope. This may take some time... '
                %disp ' '

                c=Initializer2(TOL,A,b,c);
                %[c,fval]=Initializer1(TOL,A,b,c,10000);
                s=slackfun(c);

            end

            if ~approxinpoly(s,TOL)
                    %error('Unable to locate a point near the interior of the feasible region.')
                    V=[];
                    nr=[];
                    return
            end

           if ~strictinpoly(s,TOL) %Added 02/17/2012 to handle initializers too close to polytope surface

                %disp 'Recursing...'

                idx=(  abs(s)<=max(s)*TOL );

                Amod=A; bmod=b;
                 Amod(idx,:)=[];
                 bmod(idx)=[];

                Aeq=A(idx,:); %pick the nearest face to c
                beq=b(idx);

                faceVertices=lcon2vert(Amod,bmod,Aeq,beq,TOL,1);
                if isempty(faceVertices)
                   disp 'Something''s wrong. Couldn''t find face vertices. Possibly polyhedron is unbounded.'
                   keyboard
                end

                c=faceVertices(1,:).';  %Take any vertex - find local recession cone vector
                s=slackfun(c);

                idx=(  abs(s)<=max(s)*TOL );

                Asub=A(idx,:); bsub=b(idx,:);

                [aa,bb,aaeq,bbeq]=vert2lcon(Asub);
                aa=[aa;aaeq;-aaeq];
                bb=[bb;bbeq;-bbeq];

                clear aaeq bbeq

                [bmin,idx]=min(bb);

                 if bmin>=-TOL
                   disp 'Something''s wrong. We should have found a recession vector (bb<0).'
                   keyboard
                 end

                Aeq2=null(aa(idx,:)).';
                beq2=Aeq2*c;  %find intersection of polytope with line through facet centroid.

                linetips = lcon2vert(A,b,Aeq2,beq2,TOL,1);

                if size(linetips,1)<2
                   disp 'Failed to identify line segment through interior.'
                   disp 'Possibly {x: Aeq*x=beq} has weak intersection with interior({x: Ax<=b}).'
                   keyboard
                end

                lineCentroid=mean(linetips);%Relies on boundedness

                clear aa bb

                c=lineCentroid(:);
                s=slackfun(c);

            end

            b = s;
   end
    %%%end Matt J initialization

    D=bsxfun(@rdivide,A,b);

    k = convhulln(D);
    nr = unique(k(:));

    G  = zeros(size(k,1),dim);
    ee=ones(size(k,2),1);
    discard=false( 1, size(k,1) );

    for ix = 1:size(k,1) %02/17/2012 - modified

        F = D(k(ix,:),:);
        if lindep(F,TOL)<dim;
            discard(ix)=1;
            continue;
        end

        G(ix,:)=F\ee;

    end

    G(discard,:)=[];

    V = bsxfun(@plus, G, c.');

    [discard,I]=unique( round(V*1e6),'rows');
    V=V(I,:);

return

function [c,fval]=Initializer1(TOL, A,b,c,maxIter)

    thresh=-10*max(eps(b));

    if nargin>4
     [c,fval]=fminsearch(@(x) max([thresh;A*x-b]), c,optimset('MaxIter',maxIter));
    else
     [c,fval]=fminsearch(@(x) max([thresh;A*x-b]), c);
    end

return

function c=Initializer2(TOL,A,b,c)
 %norm(  (I-A*pinv(A))*(s-b) )  subj. to s>=0

    maxIter=100000;

    [mm,nn]=size(A);

     Ap=pinv(A);
     Aaug=speye(mm)-A*Ap;
     Aaugt=Aaug.';

    M=Aaugt*Aaug;
    C=sum(abs(M),2);
     C(C<=0)=min(C(C>0));

    slack=b-A*c;
    slack(slack<0)=0;

        %     relto=norm(b);
        %     relto =relto + (relto==0);
        %
        %      relres=norm(A*c-b)/relto;

    IterThresh=maxIter;
    s=slack;
    ii=0;
    %for ii=1:maxIter
    while ii<=2*maxIter %HARDCODE

       ii=ii+1;
       if ii>IterThresh,
           %warning 'This is taking a lot of iterations'
           IterThresh=IterThresh+maxIter;
       end

     s=s-Aaugt*(Aaug*(s-b))./C;
     s(s<0)=0;

       c=Ap*(b-s);
       %slack=b-A*c;
       %relres=norm(slack)/relto;
       %if all(0<slack,1)||relres<1e-6||ii==maxIter, break;  end

    end

return

function [r,idx,Xsub]=lindep(X,tol)
%Extract a linearly independent set of columns of a given matrix X
%
%    [r,idx,Xsub]=lindep(X)
%
%in:
%
%  X: The given input matrix
%  tol: A rank estimation tolerance. Default=1e-10
%
%out:
%
% r: rank estimate
% idx:  Indices (into X) of linearly independent columns
% Xsub: Extracted linearly independent columns of X

   if ~nnz(X) %X has no non-zeros and hence no independent columns

       Xsub=[]; idx=[];
       return
   end

   if nargin<2, tol=1e-10; end

     [Q, R, E] = qr(X,0);

     diagr = abs(diag(R));

     %Rank estimation
     r = find(diagr >= tol*diagr(1), 1, 'last'); %rank estimation

     if nargout>1
      idx=sort(E(1:r));
        idx=idx(:);
     end

     if nargout>2
      Xsub=X(:,idx);
     end

 function [A,b]=rownormalize(A,b)
 %Modifies A,b data pair so that norm of rows of A is either 0 or 1

  if isempty(A), return; end

  normsA=sqrt(sum(A.^2,2));
  idx=normsA>0;
  A(idx,:)=bsxfun(@rdivide,A(idx,:),normsA(idx));
  b(idx)=b(idx)./normsA(idx);

 function tf=approxinpoly(s,TOL)

   smax=max(s);

   if smax<=0
      tf=false; return
   end

   tf=all(s>=-smax*TOL);

  function tf=strictinpoly(s,TOL)

   smax=max(s);

   if smax<=0
      tf=false; return
   end

   tf=all(s>=smax*TOL);

repo: paper-simulation


text:
function bitPos=findPosOfMin1Bit(n)
%%FINDPOSOFMIN1BIT Given n, which is of an integer data type, determine the
%           position of the least significant bit that is a 1.
%
%INPUTS: n An integer or a matrix of integers. These must be an integer
%          data type (i.e. not a double).
%
%OUTPUTS: bitPos The position of the first 1 in the binary representation
%                of each of the integers in n. Counting starts at 1. A
%                value of 0 means that there is no 1 in the binary
%                representation.
%
%The algorithm just keeps shifting the number (dividing by 2) until it
%the first digit is no longer 0, counting the number of shifts needed.
%
%EXAMPLE:
% n=findPosOfMin1Bit(int32([0,4,1213,2^23]))
%
%The function returns n=[0, 3, 1, 24]. Binary representations of
%the numbers (the least significant bit is to the left) are
%0    0
%4    001
%1213 10111101001
%2^23 000000000000000000000001
%
%January 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(isfloat(n)||~isreal(n))
    error('The values must be real and be an integer data type.')
end

if(ischar(n))
    n=uint8(n);
end

bitPos=ones(size(n));
numEls=numel(n);

for curEl=1:numEls
    if(n(curEl)~=0)
        while(bitand(n(curEl),1)==0)
            bitPos(curEl)=bitPos(curEl)+1;
            n(curEl)=bitshift(n(curEl),-1);
        end
    else
        bitPos(curEl)=0;
    end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [Pbest,Gbest] = GetBest(Archive,W,Z)
% Update the pbest and gbest of each particle

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    %% Select the pbest
    AObj  = Archive.objs - repmat(Z,length(Archive),1);
    NormW = sqrt(sum(W.^2,2))';
    d1    = AObj*W'./repmat(NormW,length(Archive),1);
    d2    = sqrt(repmat(sum(AObj.^2,2),1,size(W,1))-d1.^2);
    PBI   = d1 + 5*d2;
    [~,p] = min(PBI,[],1);
    Pbest = Archive(p);

    %% Select the gbest
    Gbest = Archive(randi(size(Archive,1),1,size(W,1)));
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function y = At_Subset(x,n,Omega);

y = zeros(n,1);
y(Omega) =x;

repo: REST


text:
function variance = uniform_discrete_variance ( a, b )

%*****************************************************************************80
%
%% UNIFORM_DISCRETE_VARIANCE returns the variance of the Uniform discrete PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer A, B, the parameters of the PDF.
%    A <= B.
%
%    Output, real VARIANCE, the variance of the PDF.
%
  variance = ( ( b + 1 - a )^2 - 1 ) / 12.0;

  return
end


text:
% fig_osem_mismatch.m
% examine effects of model mismatch on ML-OSEM
% Copyright 2002-6-10, Jeff Fessler, The University of Michigan

disp 'todo: NOT DONE!'

%
% generate data
%
if ~isvar('xfbp'), printm 'data'
	x = read_zubal_emis;
	[nx ny] = size(x);
	nb = 100;
	na = 10;

	arg = arg_pair('system', 12, 'nx', nx, 'ny', ny, ...
		'nb', nb, 'na', na, 'support', 'ellipse 0 0 60 60', ...
		'orbit', 180, 'orbit_start', 0, ...
		'pixel_size', 3, 'ray_spacing', 3, ...
		'obj2det_x', 2*nx, 'obj2det_y', 2*nx, ...
		'fwhm_detector', 0, ...
		'fwhm_collimator0', 0, ...
		'fwhm_slope', 0, ...
		'scale', 0);
prompt
end

%
% block object for block iterative
%
if ~isvar('Gb'), printm 'Gb'
	f.nblock = 5;
	Gb = Gblock(G, f.nblock);
prompt
end

%
% matlab iterations
%
if ~isvar('xmat'), printm 'matlab E-ML-OSEM'
	f.niter = 9;
	f.pixmax = 6;

	xinit = ones(size(xtrue));
	xmat = eml_osem(xinit(mask(:)), Gb, yi, ci, ri, f.niter, f.pixmax, ...
		[], 1.0);	% default fast precon matches aspire
	xmat = embed(xmat, mask);
	im clf, im(121, xmat, 'Matlab E-ML-OSEM iterations'), cbar horiz
prompt
end

if ~has_aspire, return, end

%
% aspire iterations
%
if ~isvar('xasp'), printm 'aspire E-ML-OSEM'

	f.init	= [f.dir 'init.fld'];
	f.out	= [f.dir 'out.fld'];
	fld_write(f.init, xinit, 'check', 0)
	if exist(f.out, 'file'), delete(f.out), end

	if ~isvar('f.sys_type')
		f.sys_type = ['2z@' f.wtr '@-'];
	end
	if ~isvar('n.n2')
		n.n2 = n.a;	% kludge
	end
	f.saver	= 'stack,1';
	f.alg	= sprintf('osemc,fast,%d,%d,1', f.nblock, n.n2);
	f.penal	= '-';
	f.method = sprintf('@%d@%s@%s', f.niter-1, f.alg, f.penal);
	f.com = sprintf(['i -chat 0 empl3 %s %s  %s %s 1 %s 1 %s -' ...
			' %s %s 0 1 %g 0 -'], ...
		f.out, f.init, f.yi, f.ci, f.ri, f.sys_type, ...
		f.method, f.saver, f.pixmax);
	os_run(f.com)

	xasp = double(fld_read(f.out));
	im(122, xasp, 'Aspire E-ML-OSEM iterations'), cbar horiz
end

if 1
	t = vcorrcoef(xasp, xmat);
	printf('corr. %g,%g', t, t-1)

	im clf, im(221, xmat, 'xhat matlab')
	im(222, xasp, 'xhat aspire')
	im(223, (xasp-xmat)/max(xmat(:)), 'aspire-matlab'), cbar

	t1 = eql_obj(xmat, G, yi(:), ci(:), ri(:), [], mask);
	t2 = eql_obj(xasp, G, yi(:), ci(:), ri(:), [], mask);

	if im
		subplot(224)
		plot(0:f.niter-1, t1-t1(1), '-o', 0:f.niter-1, t2-t1(1), '-x')
		xlabel iteration, ylabel '\Phi change', legend('mat', 'asp', 4)
		title(sprintf('E-ML-OSPS, Nsubset=%d', f.nblock))
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function hmm_out = logisticMergeHMM(hmm_marg,hmm_full,iY)
% This function takes a multiple logistic strcuture HMM and returns the
% equivalent HMM with only the marginalised component iY of the
% logisticYdim components.

Xfulldim = hmm_full.train.ndim;
Xdim = Xfulldim - hmm_full.train.logisticYdim;
Xfulldim_new=Xdim+1;
select_vec = [1:Xdim,Xdim+iY];
hmm_out=hmm_full;

% %update S and options:
% hmm_out.train.S=hmm_out.train.S(1:Xfulldim_new,1:Xfulldim_new);
% hmm_out.train.Sind=hmm_out.train.Sind(1:Xfulldim_new,1:Xfulldim_new);
% hmm_out.train.logisticYdim = 1;
% hmm_out.train.ndim = Xfulldim_new;

%update W:
for st=1:hmm_full.train.K
    hmm_out.state(st).W.Mu_W(1:Xdim,Xdim+iY) = hmm_marg.state(st).W.Mu_W(1:Xdim,Xdim+1);
    hmm_out.state(st).W.S_W(Xdim+iY,1:Xdim,1:Xdim) = hmm_marg.state(st).W.S_W(Xdim+1,1:Xdim,1:Xdim);
    hmm_out.state(st).W.iS_W(Xdim+iY,1:Xdim,1:Xdim) = hmm_marg.state(st).W.iS_W(Xdim+1,1:Xdim,1:Xdim);
end

%update alpha and sigma
for st=1:hmm_out.train.K
    %recall that for logistic setups, alpha has dimension Xdim x logisticYdim
    hmm_out.state(st).alpha.Gam_rate(1:Xdim,iY) = hmm_marg.state(st).alpha.Gam_rate;
end

%update psi:
hmm_out.psi(:,iY) = hmm_marg.psi;

end

author: OHBA-analysis

repo: HMM-MAR


text:
% DEMSPGP1D2 Do a simple 1-D regression after Snelson & Ghahramani's example.

% GP

% Fix seeds
randn('seed', 2e5);
rand('seed', 2e5);
seedVal = 2e5;
dataSetName = 'spgp1d';
experimentNo = 2;

% load data
[X, y] = mapLoadData(dataSetName, seedVal);

% Set up model
options = gpOptions('fitc');
options.numActive = 9;

% use the deterministic training conditional.
q = size(X, 2);
d = size(y, 2);

model = gpCreate(q, d, X, y, options);
model.X_u = randn(9, 1)*0.25 - 0.75;
params = gpExtractParam(model);
model = gpExpandParam(model, params);

% Optimise the model.
iters = 1000;
display = 1;
model.beta = 4/var(y);
model.kern.variance = var(y);
model.kern.inverseWidth = 1./((-min(X)+max(X))'/2).^2

model = gpOptimise(model, display, iters);

% Save the results.
capName = dataSetName;;
capName(1) = upper(capName(1));
save(['dem' capName num2str(experimentNo) '.mat'], 'model');

demSpgp1dPlot

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function varargout = fracInt(varargin)
%FRACINT  Fractional integral of an TRIGTECH.
%   Fractional integrals of TRIGTECH objects are not yet supported.

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

error('CHEBFUN:TRIGTECH:fracInt:noSupport', ...
    'Fractional integrals of TRIGTECH objects are not yet supported.')

% TODO: Implement this.

end


text:
function [] = createnewblocklibrary(CGen)
%CodeGenerator.createnewblocklibrary
%
% Creates a new Simulink library file along with its slblocks.m.
%
% Authors::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
% See also SerialLink, Link.

% Copyright (C) 1993-2012, by Peter I. Corke
% Copyright (C) 2012-2013, by Joern Malzahn
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
%
% The code generation module emerged during the work on a project funded by
% the German Research Foundation (DFG, BE1569/7-1). The authors gratefully
% acknowledge the financial support.

warning off;
CGen.rmpath;
warning on;

%% Create Simulink library file
new_system(CGen.slib,'Library', 'ErrorIfShadowed');
open_system(CGen.slib);
save_system(CGen.slib,CGen.slibpath);

%% Create slblock.m

fid = fopen(fullfile(CGen.basepath,'slblocks.m'),'w+');

fprintf(fid,'%s\n',['% This file is the slblocks file for the ', CGen.getrobfname,' Simulink library.']);
fprintf(fid,'%s\n',['% It has been automatically generated using The Robotics Toolbox for Matlab (RTB).']);
fprintf(fid,'%s\n',' ');
fprintf(fid,'%s\n\n',CGen.generatecopyrightnote);

fprintf(fid,'%s\n','function blkStruct = slblocks');
fprintf(fid,'%s\n',' ');
fprintf(fid,'%s\n', ['Browser.Library = ''',CGen.slib,''';']);
fprintf(fid,'%s\n', ['Browser.Name = ''',CGen.getrobfname,''';']);
fprintf(fid,'%s\n', ['Browser.IsFlat = 1;']);
fprintf(fid,'%s\n',' ');

fprintf(fid,'%s\n', ['blkStruct.Name = [''',CGen.getrobfname,''', sprintf(''\n''), ''Library''];']);
fprintf(fid,'%s\n', ['blkStruct.OpenFcn= ''',CGen.slib,''';']);
fprintf(fid,'%s\n', ['blkStruct.MaskDisplay = ''''',';']);
fprintf(fid,'%s\n', ['blkStruct.Bworser = Browser;']);
fprintf(fid,'%s\n',' ');

fclose(fid);

author: petercorke

repo: robotics-toolbox-matlab


text:
fprintf('** building frne mex file\n');

mex frne.c ne.c vmath.c


text:
%% MEEG dataset operations
%
% This example shows MVPA analyses performed on MEEG data.
%
% The input dataset involved a paradigm with electrical median nerve
% stimulation for durations of 2s at 20Hz.
%
%
% The code presented here can be adapted for other MEEG analyses, but
% there are a few potential caveats:
% * assignment of targets (labels of conditions) is based here on
%   stimulation periods versus pre-stimulation periods. In typical
%   analyses the targets should be based on different trial conditions, for
%   example as set a FieldTrip .trialinfo field.
% * assignment of chunks (parts of the data that are assumed to be
%   independent) is done randomly. In typical analyses they can be based on
%   different runs, with one chunk value per run.
% * the time window used for analyses is rather small. This means that in
%   particular for time-freq analysis a lot of data is missing, especially
%   for early and late timepoints in the lower frequency bands. For typical
%   analyses it may be preferred to use a wider time window.
% * the current examples do not perform baseline corrections or signal
%   normalizations, which may reduce discriminatory power.
%
% Note: running this code requires FieldTrip.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

%% get timelock data in CoSMoMVPA format

% set configuration
config=cosmo_config();
data_path=fullfile(config.tutorial_data_path,'meg_20hz');

% show dataset information
readme_fn=fullfile(data_path,'README');
cosmo_type(readme_fn);

% reset citation list
cosmo_check_external('-tic');

% load data
data_fn=fullfile(data_path,'subj102_B01_20Hz_timelock.mat');
data_tl=load(data_fn);

% convert to cosmomvpa struct
ds_tl=cosmo_meeg_dataset(data_tl);

% set the target (trial condition)
ds_tl.sa.targets=ds_tl.sa.trialinfo(:,1); % 1=pre, 2=post

% in addition give a label to each trial
index2label={'pre','post'}; % 1=pre, 2=peri/post
ds_tl.sa.labels=cellfun(@(x)index2label(x),num2cell(ds_tl.sa.targets));

% just to check everything is ok
cosmo_check_dataset(ds_tl);

%% Print some information about dataset

% print the size. Note that the number of features is the number of time
% points times the number of channels.
[ns,nf]=size(ds_tl.samples);
fprintf('There are %d samples and %d features\n', ns, nf);

fdim_labels=ds_tl.a.fdim.labels;
fdim_values=ds_tl.a.fdim.values;
ndim=numel(fdim_labels);
for dim=1:ndim
    dim_label=fdim_labels{dim};
    nvalues=numel(fdim_values{dim});
    fprintf('Dimension %d (%s) has %d values\n', dim, dim_label, nvalues);
end

%% Simple data manipulations
% select five trials from post-stim period
post_mask=cosmo_match(ds_tl.sa.labels,'post');
ds_tl_post=cosmo_slice(ds_tl,post_mask);
select_trials=2:2:18;

% the last argument indicates to slice along the first dimension.
% The value 1 is the default and can be omitted
ds_trial=cosmo_slice(ds_tl_post, select_trials, 1);

% select a single channel (multiple channels are also possible but
% makes plotting the samples directly a bit more tricky)
chan={'MEG1042'};
chan_msk=cosmo_dim_match(ds_trial,'chan', chan);

% slice again, this time along the feature dimension (dimension 2),
% to get a new dataset with just one channel
ds_trial_chan=cosmo_slice(ds_trial, chan_msk, 2);

% plot five individual trials
figure();
subplot(2,1,1);
plot(data_tl.time, ds_trial_chan.samples');
title(sprintf('%d trials of sensor %s', numel(select_trials), chan{1}));

% use original full data to select data in just one channel
chan_msk=cosmo_dim_match(ds_tl,'chan', chan);
ds_chan=cosmo_slice(ds_tl, chan_msk, 2);

ds_mean=cosmo_fx(ds_chan,@(x)mean(x,1),'targets');
% plot signal time course averaged over all samples
% note: the 'pre' data has been shifted in time to match the 'post' data
subplot(2,1,2);
[time_dim, time_index]=cosmo_dim_find(ds_mean,'time');
assert(time_dim==2); % it is a feature dimension
plot(ds_mean.a.fdim.values{time_index},ds_mean.samples);
title(sprintf('average of sensor %s', chan{1}));
legend(ds_mean.sa.labels);

% Show citation information
cosmo_check_external('-cite');

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [varargout] = p2pSimpleICP(Md, MovData, RefData, Tf0, DistThr )
if nargin == 0
    clc; close all;
    A = pcread('teapot.ply');  %pcread('E:\Dataset\Point Cloud\ply\Bunny\reconstruction\bun_zipper.ply');
    A = pcdownsample(A, 'gridAverage', 0.10 );
    MovData = double(A.Location');
    R0 = eul2rotm( deg2rad([20.0 10.0 10.0]) );
    T0 = [5 5 10.00 ]';
    RefData = Loc2Glo( MovData, R0', T0 );
    Sigma = 0.01;
    rng(100);
    MovData = MovData + Sigma * randn(3,  size(MovData, 2) );
    Md = createns(RefData');
    Tf0 = [ eye(3) mean(RefData-MovData, 2); 0 0 0 1];
    DistThr = Inf;
    [NNIdx, DD] = knnsearch( RefData', RefData', 'k', 2);
    Res = mean(DD(:, 2));
    DistThr = Inf;
    AngThr = cosd(45);
end
bVerbose = false; % true;
nDim = size(RefData, 1);
MaxIter = 30;
R = Tf0(1:nDim, 1:nDim);
T = Tf0(1:nDim, end);
for Iter = 1 : 1 : MaxIter
    %%%%%%%%%%% establish correspondence.
    AftData = Loc2Glo( MovData, R', T );
    [NNIdx, DD] = knnsearch( Md, AftData' );
    idx = find( DD' < DistThr );
    MovIdx = idx;
    RefIdx = NNIdx(idx);
    %%%%%%%%%%% obtain rotation and translation via SVD.
    [ dR, dT ] = RegFun(RefData(:, RefIdx), AftData(:, MovIdx) );
    R = dR * R;
    T = dR * T + dT;
    %%%%%%%%%%% check convergence condition.
    Err = max( norm(dR - eye(nDim)), norm(dT) );  % CalRT_Diff(TotalTf(end).Tf, TotalTf(end-1).Tf );
%     str = sprintf( 'Iter = %02d, Err = %f\n', Iter, Err );
%     disp(str);
    if Err(1) <= 1e-4
        break;
    end
end
if nargout == 1
    varargout{1} = [R T];
end
if nargout == 2
    varargout{1} = R;
    varargout{2} = T;
end
IS_SHOW = 0;
if IS_SHOW
    figure;
    hold on;
    grid on;
    axis equal;
    AftData = Loc2Glo( MovData, R', T );
    if nDim == 3
        view(3);
        showPointCloud(RefData', 'g');
        showPointCloud(MovData', 'r');
        showPointCloud(AftData', 'b');
    else
        plot(RefData(1, :), RefData(2, :), 'g.');
        plot(MovData(1, :), MovData(2, :), 'r.');
        plot(AftData(1, :), AftData(2, :), 'bo', 'markersize', 3);
    end

    if nargin == 0
        ErrR = RotationDiff( R, R0 );
        ErrT = norm(T - T0 );
        title(sprintf( 'Point-to-Point ICP, ErrR = %.4fdegs ErrT = %.4f', rad2deg(ErrR), ErrT ));
    else
        title('Point-to-Point ICP');
    end
end
bTest = 1;
end

author: DrGabor

repo: LiDAR

save_path: github-repos/MATLAB/DrGabor-LiDAR


text:
function [xPred, PPred, F, B] = ExtendedKalmanFilterX_PredictState(x,P,f,Q,u,b,Qu)
% EKALMANFILTERX_PREDICTSTATE Perform the discrete-time KF state prediction
% step, under the assumption of additive process noise.
%
% Parameters
% ----------
% x: column vector
%   The (xDim x 1) state estimate at the previous time-step.
% P: matrix
%   The (xDim x xDim) state covariance matrix at the previous
%   time-step.
% f: function handle
%   A (non-linear) state transition function.
% Q: matrix
%   The (xDim x xDim) process noise covariance matrix.
% u: column vector, optional
%   An optional (xDim x 1) control input.
%   If omitted, no control input is used.
% b: function handle, optional
%   A (non-linear) control gain function.
%   (Optional, Default = 1 if u provided, 0 otherwise)
% O: matrix, optional
%   An optional (xDim x xDim) control noise covariance
%   matrix. If omitted, Q is assumed to be 0.
%
% Returns
% -------
% xPred: column vector
%   The (xDim x 1) predicted state estimate.
% PPred: matrix
%   The (xDim x xDim) predicted state covariance matrix.
% F: matrix
%   The computed Jacobian transition matrix
% H: matrix
%   The computed (yDim x yDim) Jacobian measurement matrix
% B: matrix, optional
%   The computed Jacobian control gain matrix
%
% October 2017 Lyudmil Vladimirov, University of Liverpool.

    switch(nargin)
        case(4)
            u  = 0;
            b  = 0;
            Qu = 0;
        case(5)
            b  = 1;
            Qu = 0;
        case(6)
            Qu = 0;
    end

    % Prediction for state vector and covariance:
    [xPred,F] = ExtendedKalmanFilterX_computeJac(f,x);    %nonlinear update and linearization at current state
    PPred = F*P*F' + Q;                 %partial update

    % Compute Control Input (if applicable)
    [controlInputWithGain,B] = ExtendedKalmanFilterX_computeJac(b,u);

    % Add control input
    xPred = xPred + controlInputWithGain;
    PPred = PPred + B*Qu*B';
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function [assignment,cost] = munkres(costMat)
% MUNKRES   Munkres (Hungarian) Algorithm for Linear Assignment Problem.
%
% [ASSIGN,COST] = munkres(COSTMAT) returns the optimal column indices,
% ASSIGN assigned to each row and the minimum COST based on the assignment
% problem represented by the COSTMAT, where the (i,j)th element represents the cost to assign the jth
% job to the ith worker.
%

% This is vectorized implementation of the algorithm. It is the fastest
% among all Matlab implementations of the algorithm.

% Examples
% Example 1: a 5 x 5 example
%{
[assignment,cost] = munkres(magic(5));
disp(assignment); % 3 2 1 5 4
disp(cost); %15
%}
% Example 2: 400 x 400 random data
%{
n=400;
A=rand(n);
tic
[a,b]=munkres(A);
toc                 % about 2 seconds
%}
% Example 3: rectangular assignment with inf costs
%{
A=rand(10,7);
A(A>0.7)=Inf;
[a,b]=munkres(A);
%}
% Reference:
% "Munkres' Assignment Algorithm, Modified for Rectangular Matrices",
% http://csclab.murraystate.edu/bob.pilgrim/445/munkres.html

% version 2.2 by Yi Cao at Cranfield University on 1st March 2010

assignment = zeros(1,size(costMat,1));
cost = 0;

costMat(costMat~=costMat)=Inf;
validMat = costMat<Inf;
validCol = any(validMat,1);
validRow = any(validMat,2);

nRows = sum(validRow);
nCols = sum(validCol);
n = max(nRows,nCols);
if ~n
    return
end

maxv=10*max(costMat(validMat));

dMat = zeros(n) + maxv;
dMat(1:nRows,1:nCols) = costMat(validRow,validCol);

%*************************************************
% Munkres' Assignment Algorithm starts here
%*************************************************

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   STEP 1: Subtract the row minimum from each row.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
minR = min(dMat,[],2);
minC = min(bsxfun(@minus, dMat, minR));

%**************************************************************************
%   STEP 2: Find a zero of dMat. If there are no starred zeros in its
%           column or row start the zero. Repeat for each zero
%**************************************************************************
zP = dMat == bsxfun(@plus, minC, minR);

starZ = zeros(n,1);
while any(zP(:))
    [r,c]=find(zP,1);
    starZ(r)=c;
    zP(r,:)=false;
    zP(:,c)=false;
end

while 1
%**************************************************************************
%   STEP 3: Cover each column with a starred zero. If all the columns are
%           covered then the matching is maximum
%**************************************************************************
    if all(starZ>0)
        break
    end
    coverColumn = false(1,n);
    coverColumn(starZ(starZ>0))=true;
    coverRow = false(n,1);
    primeZ = zeros(n,1);
    [rIdx, cIdx] = find(dMat(~coverRow,~coverColumn)==bsxfun(@plus,minR(~coverRow),minC(~coverColumn)));
    while 1
        %**************************************************************************
        %   STEP 4: Find a noncovered zero and prime it.  If there is no starred
        %           zero in the row containing this primed zero, Go to Step 5.
        %           Otherwise, cover this row and uncover the column containing
        %           the starred zero. Continue in this manner until there are no
        %           uncovered zeros left. Save the smallest uncovered value and
        %           Go to Step 6.
        %**************************************************************************
        cR = find(~coverRow);
        cC = find(~coverColumn);
        rIdx = cR(rIdx);
        cIdx = cC(cIdx);
        Step = 6;
        while ~isempty(cIdx)
            uZr = rIdx(1);
            uZc = cIdx(1);
            primeZ(uZr) = uZc;
            stz = starZ(uZr);
            if ~stz
                Step = 5;
                break;
            end
            coverRow(uZr) = true;
            coverColumn(stz) = false;
            z = rIdx==uZr;
            rIdx(z) = [];
            cIdx(z) = [];
            cR = find(~coverRow);
            z = dMat(~coverRow,stz) == minR(~coverRow) + minC(stz);
            rIdx = [rIdx(:);cR(z)];
            cIdx = [cIdx(:);stz(ones(sum(z),1))];
        end
        if Step == 6
            % *************************************************************************
            % STEP 6: Add the minimum uncovered value to every element of each covered
            %         row, and subtract it from every element of each uncovered column.
            %         Return to Step 4 without altering any stars, primes, or covered lines.
            %**************************************************************************
            [minval,rIdx,cIdx]=outerplus(dMat(~coverRow,~coverColumn),minR(~coverRow),minC(~coverColumn));
            minC(~coverColumn) = minC(~coverColumn) + minval;
            minR(coverRow) = minR(coverRow) - minval;
        else
            break
        end
    end
    %**************************************************************************
    % STEP 5:
    %  Construct a series of alternating primed and starred zeros as
    %  follows:
    %  Let Z0 represent the uncovered primed zero found in Step 4.
    %  Let Z1 denote the starred zero in the column of Z0 (if any).
    %  Let Z2 denote the primed zero in the row of Z1 (there will always
    %  be one).  Continue until the series terminates at a primed zero
    %  that has no starred zero in its column.  Unstar each starred
    %  zero of the series, star each primed zero of the series, erase
    %  all primes and uncover every line in the matrix.  Return to Step 3.
    %**************************************************************************
    rowZ1 = find(starZ==uZc);
    starZ(uZr)=uZc;
    while rowZ1>0
        starZ(rowZ1)=0;
        uZc = primeZ(rowZ1);
        uZr = rowZ1;
        rowZ1 = find(starZ==uZc);
        starZ(uZr)=uZc;
    end
end

% Cost of assignment
rowIdx = find(validRow);
colIdx = find(validCol);
starZ = starZ(1:nRows);
vIdx = starZ <= nCols;
assignment(rowIdx(vIdx)) = colIdx(starZ(vIdx));
cost = trace(costMat(assignment>0,assignment(assignment>0)));

function [minval,rIdx,cIdx]=outerplus(M,x,y)
ny=size(M,2);
minval=inf;
for c=1:ny
    M(:,c)=M(:,c)-(x+y(c));
    minval = min(minval,min(M(:,c)));
end
[rIdx,cIdx]=find(M==minval);

author: Sable

repo: mcbench-benchmarks


text:
function M = moveUDRVC(F, d)
%------------------------------------------------------------------------------
% Moves gridfunction F in vertical direction.
% Excess area is filled by Vertex-Centered (VC) Reflection across boundaries.
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: June 23, 2000.
% (c) 1999-2002 Stichting CWI, Amsterdam
%------------------------------------------------------------------------------
[n, m] = size(F);
if d > 0
  if (d+1) > n
    error(' moveUDRVC    d >= n ')
  else
    M = [flipud(F(2:(d+1),:)); F(1:(n-d),:)];
  end
elseif d < 0
  if (-d+1) > n
    error(' moveUDRVC   -d >= n ')
  else
    M = [F((-d+1):n,:); flipud(F((n+d):(n-1),:))];
  end
else
  M = F;
end
%------------------------------------------------------------------------------


function [ave, nvals] = calc_average(structure,field)
%CALC_AVERAGE Calculate the average of values in a field.
% Function CALC_AVERAGE calculates the average value
% of the elements in a particular field of a structure
% array.  It returns the average value and (optionally)
% the number of items averaged.

% Define variables:
%   arr       -- Array of values to average
%   ave       -- Average of arr
%   ii        -- Index variable
%
%  Record of revisions:
%      Date       Programmer          Description of change
%      ====       ==========          =====================
%    03/04/07    S. J. Chapman        Original code
%
% Check for a legal number of input arguments.
msg = nargchk(2,2,nargin);
error(msg);

% Create an array of values from the field
arr = [];
for ii = 1:length(structure)
   arr = [arr structure(ii).(field)];
end

% Calculate average
ave = mean(arr);

% Return number of values averaged
if nargout == 2
   nvals = length(arr);
end


text:
function p = paris(grains)
% Percentile Average Relative Indented Surface
%
% the paris (Percentile Average Relative Indented Surface) is shape
% parameter that measures the convexity of a grain
%
% Syntax
%   p = paris(grains)
%
% Input
%  grains - @grain2d
%
% Output
%  p - double
%
%

% store this in local variables for speed reasons
p = zeros(size(grains));

X = grains.V(:,1);
Y = grains.V(:,2);

poly = grains.poly;
% remove inclusions
incl = grains.inclusionId;
for i = find(incl>0).'
  poly{i} = poly{i}(1:end-incl(i));
end

% compute convex hull perimeters
for id = 1:length(grains)

  % for small grains there is no difference
  if length(poly{id}) <= 7, continue; end

  % extract coordinates
  xGrain = X(poly{id});
  yGrain = Y(poly{id});

  % compute perimeter
  perimeterGrain = sum(sqrt(...
    (xGrain(1:end-1) - xGrain(2:end)).^2 + ...
    (yGrain(1:end-1) - yGrain(2:end)).^2));

  % compute convex hull
  ixy = convhull(xGrain,yGrain);

  % compute perimenter
  perimeterHull = sum(sqrt(...
    (xGrain(ixy(1:end-1)) - xGrain(ixy(2:end))).^2 + ...
    (yGrain(ixy(1:end-1)) - yGrain(ixy(2:end))).^2));

  % paris is the relative difference between convex hull perimenter and true
  % perimeter
  p(id) = 200*(perimeterGrain - perimeterHull)./perimeterHull;

end

author: mtex-toolbox

repo: mtex

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


text:
function [score alignment] = ScanMatch(seq1, seq2, ScanMatchInfo, varargin)
%SCANMATCH compares two sequences using the Needleman-Wunsch
% algorithm providing a GUI is needed to display the results.
%   SCORE = SCANMATCH(SEQ1, SEQ2, SCANMATCHINFO) where SEQ1 and SEQ2 are
%   two double strings (where a single fixation is represented by two
%   characters). SCANMATCHINFO is a structure containing all the parameters
%   needed to perform the comparison and must contain the following fields:
%
%           .Xres
%           .Yres
%           .Xbin
%           .Ybin
%     .RoiModulus
%      .Threshold
%       .GapValue
%        .TempBin
%      .SubMatrix
%           .mask
%
%   As an extra input parameter 'ShowViewer'm can be set to 1 to display the
%   GUI. If not set at all or set to zero the Gui will be not shown.
%
%   SCORE is the normalised result of the alignment. The closest SCORE is
%   to one the more related the two sequences are. Alignment output the
%   alignment of the two sequences. This is quite computationally expensive
%   hence if not needed call the function as follow:
%           score = ScanMatch(seq1, seq2, ScanMatchInfo)
%
%   Part of the ScanMatch toolbox
%   Written by Filipe Cristino
%   $Version: 1.00 $  $Date: 10/09/2009

% Check the input parameters
okargs = {'ShowViewer'};
if nargin == 3
    ShowViewer = 0; % if not specified then don't show
else
    if strcmp(varargin(1), okargs)
        ShowViewer = varargin{2};
    else
        error('ScanMatch:InvalidInput','Invalid parameter: %s', char(varargin(1)));
    end
end

% Check if the ScanMatchInfo structure is of the right form
Ok = ScanMatch_CheckStructure(ScanMatchInfo);

% use numerical arrays for easy indexing
if ischar(seq1)
    intseq1 = uint8(ScanMatch_DoubleStrToNum(seq1,ScanMatchInfo.RoiModulus));
else
    intseq1 = double(seq1);
    seq1 = ScanMatch_NumToDoubleStr(intseq1,ScanMatchInfo.RoiModulus);
end
if ischar(seq2)
    intseq2 = uint8(ScanMatch_DoubleStrToNum(seq2,ScanMatchInfo.RoiModulus));
else
    intseq2 = double(seq2);
    seq2 = ScanMatch_NumToDoubleStr(intseq2,ScanMatchInfo.RoiModulus);
end

% Check for the sequences lenghts >0 and that they are even

m = length(seq1);
n = length(seq2);

if  ~n||~m || mod(m,2)||mod(n,2)
    error('ScanMatch:InvalidLengthSequences','Length of input sequences must be greater than 0 and even');
end

m = length(intseq1);
n = length(intseq2);

% Perform the Needleman-Wunsch alogorithm
[score path step F] = ScanMatch_nwAlgo(intseq1,m,intseq2,n,ScanMatchInfo.SubMatrix,ScanMatchInfo.GapValue);

% Normalise output score
max_sub = max(ScanMatchInfo.SubMatrix(:));
scale = max_sub * max(m,n);
score = score / scale;

% if alignment output is not needed then exit now (faster)
if nargout<=1 && ~ShowViewer
    return
end

path(step:end,:) = []; % step-1 is the length of the alignment
path = flipud(path);

%double the path to take into account the double string
p1 = interp1(1:step-1,double(path(:,1)),fix(1:0.5:step-1+0.5),'nearest');
p2 = interp1(1:step-1,double(path(:,2)),fix(1:0.5:step-1+0.5),'nearest');
path = [p1' p2'];

% setting the size of the alignment
n_l = (step-1)*2;
alignment = repmat(('- -')',1,n_l);

% adding sequence to alignment
alignment(1,path(:,1)>0) = seq1;
alignment(3,path(:,2)>0) = seq2;

% find locations where there are no gaps
h=find(all(path>0,2));
noGaps1=ScanMatch_DoubleStrToNum(alignment(1,h),ScanMatchInfo.RoiModulus);
noGaps2=ScanMatch_DoubleStrToNum(alignment(3,h),ScanMatchInfo.RoiModulus);

% score pairs with no gap
value = ScanMatchInfo.SubMatrix(sub2ind(size(ScanMatchInfo.SubMatrix),double(noGaps1),double(noGaps2)));
value = interp1(1:length(value),value,fix(1:0.5:length(value)+0.5),'nearest');

% insert symbols of the match string into the alignment
alignment(2,h(value>=0)) = ':';
alignment(2,h(value==max(ScanMatchInfo.SubMatrix(:)))) = '|';

if ShowViewer == 1 % display the Gui viewer

    % Set some variable
    res_size = [800 380];

    % Create and setup figure
    fh = figure('Position',[200 200 res_size],'Toolbar','none', ...
                'name', 'String Alignment Viewer','NumberTitle', 'off',...
                'MenuBar', 'none', 'Color', [0.95,0.95,0.95],'wvisual','27');

    uicontrol('Style','text', ...
                    'String','Alignment Results',...
                    'Position', [res_size(1)/2 - 100 res_size(2)-25 200 25], 'Fontsize', 16);

    ah2 = axes('Parent',fh,'units','pixels',...
               'Position',[40 40 280 180 ],...
               'Visible', 'off',...
               'DataAspectRatioMode','auto'...
           );

    % Display the Score matrix and winning path
    F=scale.*max(F(2:end,2:end,:),[],3);
    clim=max(max(max(abs(F(~isinf(F))))),eps);
    imagesc(F,[-clim clim]);
    colormap(privateColorMap(1));
    set(colorbar,'YLim',[min([F(:);-eps]) max([F(:);eps])])
    title('Scoring Space and Winning Path')
    xlabel('Sequence 1')
    ylabel('Sequence 2')
    hold on
    plot(path(all(path>0,2),1),path(all(path>0,2),2),'k.')

    % Create display for the strings
    S_r = 20;
    y_p = 290;

    ah3 = axes('Parent',fh,'units','pixels',...
               'Position',[0 y_p res_size(1) 61 ],...
               'Visible', 'on',...
               'DataAspectRatioMode','auto'...
           );
    axis off
    set(ah3,'xlim',[0 800]);
    set(ah3,'ylim',[0 60]);

    [th,ws] = create_display_string(alignment,ah3,0, S_r, res_size);

    % Do slider for roling window
    SliderW = uicontrol('Parent', fh, 'Style', 'slider',...
                   'Position', [20 y_p-30 res_size(1)-40  20],...
                   'Min', 0, ...
                   'Max', ws - 800, ...
                   'Tag', 'SliderWindow', ...
                   'SliderStep', [0.05 0.2] ,...
                   'Callback',{@popup_menu_Callback} );

    % show scores
    uicontrol('Style','text', ...
                    'String','Score:',...
                    'Position', [330 120 160 40],'FontSize', 18);
    Scorea = uicontrol('Style','text', ...
                    'String',score,...
                    'Position', [330 40 160 60], 'ForegroundColor', [1 0 0], 'FontSize', 24);

    % Display ScanMatchInfo parameters
    ah1 = uipanel('units','pixels',...
               'Position',[520 0 279 239]...
           );

    uicontrol(ah1,'Style','text', ...
                    'String','Parameters',...
                    'Position', [0 210 270 25], 'FontSize', 16);

    uicontrol(ah1,'Style','text', ...
                    'String','Alignment algorithm:',...
                    'Position', [5 170 130 20], 'FontSize', 10, ...
                    'HorizontalAlignment', 'right');

    uicontrol(ah1,'Style','text', ...
                    'String','Needleman-Wunsch',...
                    'Position', [145 170 130 20], 'FontSize', 10, ...
                    'FontWeight', 'bold', 'HorizontalAlignment', 'left');

    uicontrol(ah1,'Style','text', ...
                    'String','Gap value:',...
                    'Position', [5 150 130 20], 'FontSize', 10, ...
                    'HorizontalAlignment', 'right');

    uicontrol(ah1,'Style','text', ...
                    'String',num2str(ScanMatchInfo.GapValue),...
                    'Position', [145 150 130 20], 'FontSize', 10, ...
                    'FontWeight', 'bold','HorizontalAlignment', 'left');

    uicontrol(ah1,'Style','text', ...
                    'String','Substitution threshold:',...
                    'Position', [5 130 130 20], 'FontSize', 10, ...
                    'HorizontalAlignment', 'right');

    uicontrol(ah1,'Style','text', ...
                    'String',num2str(ScanMatchInfo.Threshold),...
                    'Position', [145 130 130 20], 'FontSize', 10, ...
                    'FontWeight', 'bold','HorizontalAlignment', 'left');

    uicontrol(ah1,'Style','text', ...
                    'String','RoI modulus:',...
                    'Position', [5 110 130 20], 'FontSize', 10, ...
                    'HorizontalAlignment', 'right');

    uicontrol(ah1,'Style','text', ...
                    'String',num2str(ScanMatchInfo.RoiModulus),...
                    'Position', [145 110 130 20], 'FontSize', 10, ...
                    'FontWeight', 'bold','HorizontalAlignment', 'left');

    uicontrol(ah1,'Style','text', ...
                    'String','Temporal Bin size:',...
                    'Position', [5 90 130 20], 'FontSize', 10, ...
                    'HorizontalAlignment', 'right');

    uicontrol(ah1,'Style','text', ...
                    'String',[num2str(ScanMatchInfo.TempBin) ' ms'],...
                    'Position', [145 90 130 20], 'FontSize', 10, ...
                    'FontWeight', 'bold','HorizontalAlignment', 'left');

    uicontrol(ah1,'Style','text', ...
                    'String','Stimuli Resolution:',...
                    'Position', [5 70 130 20], 'FontSize', 10, ...
                    'HorizontalAlignment', 'right');

    uicontrol(ah1,'Style','text', ...
                    'String',[num2str(ScanMatchInfo.Xres) ' x ' num2str(ScanMatchInfo.Yres) ' pixels'],...
                    'Position', [145 70 130 20], 'FontSize', 10, ...
                    'FontWeight', 'bold','HorizontalAlignment', 'left');

    uicontrol(ah1,'Style','text', ...
                    'String','Number of bins:',...
                    'Position', [5 50 130 20], 'FontSize', 10, ...
                    'HorizontalAlignment', 'right');

    uicontrol(ah1,'Style','text', ...
                    'String',[num2str(ScanMatchInfo.Xbin) ' x ' num2str(ScanMatchInfo.Ybin) ],...
                    'Position', [145 50 130 20], 'FontSize', 10, ...
                    'FontWeight', 'bold','HorizontalAlignment', 'left');
end

function popup_menu_Callback(source,eventdata)
% This function is Nested hence all variables are globals.
      pop_menu = get(source,'Tag');
      switch pop_menu
          case 'SliderGap'
              set(fh,'CurrentAxes',ah2)
              GapValue = get(source,'Value');
              [score, data]=nwalignds(StringA, StringB,'ScoringMatrix',Transition_dist, 'GapOpen', GapValue,'SHOWSCORE', true);
              set(Scorea, 'string', num2str(score));
              set(Gapa, 'string', ['Gap value: ' num2str(GapValue)]);
              set(fh,'CurrentAxes',ah3)
              delete(th); % delete old display
              [th ws] = create_display_string(data,ah3,0, S_r,info); % Create new display for the strings
              set(SliderW, 'Value', 0);
              set(SliderW, 'Max', ws-800);
          case 'Transition_mat'

          case 'SliderWindow'
              dx = get(source,'Value');
              set(ah3,'xlim',[dx 800+dx]);

      end
end

end

function [th, ws] = create_display_string(data,ah3,y_p, S_r, res_size)

n_string = size(data,2)/2; % IF MORE THAN 30 YOU SHOULD SHOW A SCROOLING BAR!
if n_string > 25
   S_c = res_size(1)/25;
else
    S_c = res_size(1)/n_string-1;
end
ws = n_string * S_c;
th = zeros(3,n_string);
tic
for i=1:3
    x_p = S_c/2;
    ind=1;
    for j=1:n_string
    switch data(2,ind)
        case '|'
            col = [0 0 0];
        case ':'
            col = [1 0 0];
        case ' '
            col = [0 0 1];
    end

    th(i,j)=text('parent', ah3', 'Position', [x_p y_p], ...
            'string', data(i,ind:ind+1), 'color', col,'HorizontalAlignment', 'center' );
    x_p = x_p + S_c;
    ind = ind +2;
    end
    y_p = y_p + S_r;
end

toc
end

function pcmap = privateColorMap(selection)
%PRIVATECOLORMAP returns a custom color map
switch selection
    case 1, pts = [0 0 .3 20;
            0 .1 .8 25;
            0 .9 .5 15;
            .9 1 .9 8;
            1 1 0 26;
            1 0 0 26;
            .4 0 0 0];
    otherwise, pts = [0 0 0 128; 1 1 1 0];
end
xcl=1;
for i=1:size(pts,1)-1
    xcl=[xcl,i+1/pts(i,4):1/pts(i,4):i+1];
end
pcmap = interp1(pts(:,1:3),xcl);
end

author: Allopart

repo: rbpf-gmapping


text:
function call_price=american_call_baw(S, X, r, b, sigma, time, accuracy)

%--------------------------------------------------------------------------
%
% DESCRIPTION:
%
% Barone-Adesi and Whaley (1987) quadratic approximation to the price
% of a call option.
%
%
% Reference:
%
% Giovanni Barone-Adesi and Robert E. Whaley,
% "Efficient analytic approximation of American option values",
% Journal of Finance, 42(2):301-320, June 1987.
%
%
%--------------------------------------------------------------------------
%
% INPUTS:
%
%   S:          spot price
%   X:          exercise price
%   r:          interest rate
%   b:          dividend yield
%   sigma:      volatility
%   time:       time to maturity
%   accuracy:   approximation accuracy
%
%--------------------------------------------------------------------------
%
% OUTPUT:
%
% call_price: price of a call option
%
%--------------------------------------------------------------------------
%
% Author:  Paolo Z., February 2012
%
%--------------------------------------------------------------------------

sigma_sqr = sigma*sigma;
time_sqrt = sqrt(time);

nn = 2.0*b/sigma_sqr;

m = 2.0*r/sigma_sqr;
K = 1.0-exp(-r*time);

q2 = ( -(nn-1) + sqrt( (nn-1)^2 + (4*m/K) ) )*0.5;
q2_inf = 0.5 * ( -(nn-1) + sqrt( (nn-1)^2 + 4.0*m ) );

S_star_inf = X / (1.0 - 1.0/q2_inf);
h2 = -(b*time+2.0*sigma*time_sqrt)*(X/(S_star_inf-X));

S_seed = X + (S_star_inf-X)*(1.0-exp(h2));

no_iterations=0;

Si     = S_seed;
g      = 1;
gprime = 1;

while ((abs(g) > accuracy) && (abs(gprime)>accuracy) && ( no_iterations<500) && (Si>0.0))
    c = european_call_contpay(Si,X,r,b,sigma,time);

    d1     = (log(Si/X)+(b+0.5*sigma_sqr)*time)/(sigma*time_sqrt);
    g      =(1.0-1.0/q2)*Si-X-c+(1.0/q2)*Si*exp((b-r)*time)*normcdf(d1);

    gprime =( 1.0-1.0/q2)*(1.0-exp((b-r)*time)*normcdf(d1))+(1.0/q2)*exp((b-r)*time)*normcdf(d1)*(1.0/(sigma*time_sqrt));
    Si     = Si-(g/gprime);
end

S_star = 0;

if (abs(g)>accuracy)
    S_star = S_seed; %// did not converge
else
    S_star = Si;
end

C=0;
c = european_call_contpay(S,X,r,b,sigma,time);

if (S>=S_star)
    C  = S-X;
else
    d1 = (log(S_star/X)+(b+0.5*sigma_sqr)*time)/(sigma*time_sqrt);
    A2 = (1.0-exp((b-r)*time)*normcdf(d1))* (S_star/q2);
    C  = c+A2*((S/S_star)^q2);
end

call_price = max(C,c);