text:
% COST function.
% It shows how accurate our model is based on current model parameters.
function [cost] = cost_function(X, y, theta, lambda)
    % Input:
    % X - input features - (m x n) matrix.
    % theta - our model parameters - (n x 1) vector.
    % y - a vector of correct output - (m x 1) vector.
    % lambda - regularization parameter.
    %
    % Output:
    % cost - number that represents the cost (error) of our model with specified parameters theta.
    %
    % Where:
    % m - number of training examples,
    % n - number of features.

    % Get the size of the trainging set.
    m = size(X, 1);

    % Get the difference between predictions and correct output values.
    differences = hypothesis(X, theta) - y;

    % Calculate regularization parameter.
    % Remember that we should not regularize the parameter theta_zero.
    theta_cut = theta(2:end, 1);
    regularization_param = lambda * (theta_cut' * theta_cut);

    % Calculate current predictions cost.
    cost = (1 / 2 * m) * (differences' * differences + regularization_param);
end

author: trekhleb

repo: machine-learning-octave


text:
function value = r4_atanh ( x )

%*****************************************************************************80
%
%% R4_ATANH evaluates the arc-hyperbolic tangent of an R4 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the arc-hyperbolic tangent of X.
%
  persistent atnhcs;
  persistent dxrel;
  persistent nterms;
  persistent sqeps;

  if ( isempty ( nterms ) )
    atnhcs = [ ...
    0.094395102393195492E+00;
    0.049198437055786159E+00;
    0.002102593522455432E+00;
    0.000107355444977611E+00;
    0.000005978267249293E+00;
    0.000000350506203088E+00;
    0.000000021263743437E+00;
    0.000000001321694535E+00;
    0.000000000083658755E+00;
    0.000000000005370503E+00;
    0.000000000000348665E+00;
    0.000000000000022845E+00;
    0.000000000000001508E+00;
    0.000000000000000100E+00;
    0.000000000000000006E+00 ];
    nterms = r4_inits ( atnhcs, 15, 0.1 * r4_mach ( 3 ) );
    dxrel = sqrt ( r4_mach ( 4 ) );
    sqeps = sqrt ( 3.0 * r4_mach ( 3 ) );
  end

  y = abs ( x );

  if ( y <= sqeps )
    value = x;
  elseif ( y <= 0.5 )
    value = x * ( 1.0 ...
      + r4_csevl ( 8.0 * x * x - 1.0, atnhcs, nterms ) );
  elseif ( y < 1.0 )
    value = 0.5 * log ( ( 1.0 + x ) / ( 1.0 - x ) );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_ATANH - Fatal error!\n' );
    fprintf ( 1, '  1 <= |X|.\n' );
    error ( 'R4_ATANH - Fatal error!' )
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function f = bohach2_xy ( x, y )

%*****************************************************************************80
%
%% BOHACH2_XY evaluates the Bohachevsky function #2.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Zbigniew Michalewicz,
%    Genetic Algorithms + Data Structures = Evolution Programs,
%    Third Edition,
%    Springer Verlag, 1996,
%    ISBN: 3-540-60676-9,
%    LC: QA76.618.M53.
%
%  Parameters:
%
%    Input, real X, Y, the argument of the function.
%
%    Output, real F, the value of the function at X.
%
  f =       x * x ...
    + 2.0 * y * y ...
    - 0.3 * cos ( 3.0 * pi * x ) * cos ( 4.0 * pi * y ) ...
    + 0.3;

  return
end


text:
function [img,nrg] = rayImage(ray,mic,rad,rMax)
%+========================================================================+
%|                                                                        |
%|           OPENRAY - LIBRARY FOR TRI-DIMENSIONAL RAY TRACING            |
%|           openRay is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : rayImage.m                                    |
%|    #    |   VERSION    : 0.41                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 01.04.2018                                    |
%| ( === ) |   SYNOPSIS   : Image-sources from ray-tracing                |
%|  `---'  |                                                              |
%+========================================================================+

% Material properties (http://www.odeon.dk/material-manufactures)
load('odeon.mat')
mat  = odeon.mat(ray.msh.col,:);
Nfrq = size(mat,2);

% Air absorbing formulae for T = 20°C (bouquin Jouhaneau p 68-69)
rhm = 30;
air = 5.5 * (50/rhm) .* (odeon.frq/1000).^1.7 * 1e-4;

% Measures
[Isrc,src] = ray.measure(mic,rad,rMax);

% Initialization
img = cell(size(src));
nrg = cell(size(src));
rfl = ones(length(ray),Nfrq);

% Loop on images sources
for i = 1:length(src)
    % Unicity for images
    [~,Ia,Ic] = unique(round(src{i}*1e12),'rows','stable');

    % Image coordinates
    img{i} = src{i}(Ia,:) - ones(length(Ia),1)*mic;

    % Image Energy
    nrg{i} = zeros(length(Ia),Nfrq);
    for j = 1:Nfrq
        nrg{i}(:,j) = accumarray(Ic,rfl(Isrc{i},j),[length(Ia),1]);
    end

    % Update reflecting coeff
    ind         = (ray.iel(:,i+1) > 0);
    rfl(ind,:)  = rfl(ind,:) .* (1 - mat(ray.iel(ind,i+1),:));
    rfl(~ind,:) = 0;
end

% Vectoriel format
img = cell2mat(img);
nrg = cell2mat(nrg);

% Atmosphere dissipation
dst = sqrt(sum(img.^2,2));
nrg = nrg .* exp(-dst*air);

% Sort in phase
[~,ind] = sort(dst);
img     = img(ind,:);
nrg     = nrg(ind,:);

% Normalisation
nrg = nrg ./ max(max(nrg));
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [y,p_avg,p_std]=multinomrnd(p,m,n)
%Performs random sampling from a binomial distribution
%
% [y]=multinomrnd(p,m,n)
% where p=1-by-k vector of probabilities of occurrence
%       n=sample size
% and   m= number of trials
%       y=samples-matrix of size k-by-m
%
% for picking out one of k mixture components, set n=1;
%
if nargin<3
  n=1;
end

k=length(p);
x=rand(n,m);

if (sum(p)-1>100*eps)
  p(k+1)=1-sum(p);
  k=k+1;
end
p=cumsum(p);

y(1,:)=sum(x<=p(1),1);
for i=2:k
  y(i,:)=sum(x>p(i-1) & x<=p(i),1);
end

p_avg=mean(y'./n);
p_std=std(y'./n);

author: OHBA-analysis

repo: HMM-MAR


text:
% erpregout() - regress out the ERP from the data
%
% Usage:
%   newdata = erpregout(data);
%   [newdata erp factors] = erpregout(data, tlim, reglim);
%
% Inputs:
%   data    - [float] 2-D data (times x trials) or 3-D data
%             (channels x times x trials).
%
% Optional inputs:
%   tlim    - [min max] time limits in ms.
%   reglim  - [min max] regression time window in ms (by default
%             the whole time period is used
% Outputs:
%   newdata - data with ERP regressed out
%   erp     - data ERP
%   factors - factors used for regressing out the ERP (size is the same
%             as the number of trials or (channels x trials)
%
% Note: it is better to regress out the ERP about 4 times (launch the
%       function 4 times in a row) to really be able to regress out the
%       ERP and have a residual ERP close to 0.
%
% Author: Arnaud Delorme, Salk, SCCN, UCSD, CA, April 29, 2004

%123456789012345678901234567890123456789012345678901234567890123456789012

% Copyright (C) 2004 Arnaud Delorme
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

function [data, erp, factors] = erpregout(data, tlim, reglim);

    if nargin < 1
        help erpregout;
        return;
    end;
    if nargin < 2
        tlim = [0 1];
    end;
    if nargin < 3
        reglim = tlim;
    end;
    if ndims(data) == 2
        data = reshape(data, 1, size(data,1), size(data,2));
        redim = 1;
    else
        redim = 0;
    end;

    % find closest points
    % -------------------
    timevect = linspace(tlim(1), tlim(2), size(data,2));
    [tmp begpoint] = min( abs(timevect-reglim(1)) );
    [tmp endpoint] = min( abs(timevect-reglim(2)) );
    erp = mean(data, 3);

    % regressing out erp in channels and trials
    % -----------------------------------------
    for chan = 1:size(data,1)
        fprintf('Channel %d (trials out of %d):', chan, size(data,3));
        for trial = 1:size(data,3)
            if ~mod(trial, 10) , fprintf('%d ', trial); end;
            if ~mod(trial, 200), fprintf('\n', trial); end;
            [factors(chan, trial) tmpf exitflag] = fminbnd('erpregoutfunc', 0, 10, [], ...
                                   data(chan, begpoint:endpoint, trial), erp(chan, begpoint:endpoint));
            data(chan,:,trial) = data(chan,:,trial) - factors(chan, trial)*erp(chan, :);
        end;
        fprintf('\n');
    end;

    if redim
        data = squeeze(data);
    end;

author: PatternRecognition

repo: OpenBMI


text:
function[varargout]=windtrans(varargin)
%WINDTRANS  Ekman-like transfer-functions for the wind-driven response.
%
%   G=WINDTRANS(OMEGA,Z,FC,DELTA,MU,H) returns the no-slip transfer
%   function for the wind-driven currents evaluated at frequencies OMEGA
%   and at depths Z.  G will have LENGTH(OMEGA) rows and LENGTH(Z) columns.
%
%   Here FC is the local Coriolis frequency, DELTA is the Ekman depth.
%   MU is the Madsen depth, and H is the boundary layer depth.
%
%   The units of these quantitites are important.  For consistency with
%   other routines, OMEGA and FC are in radians per day, while Z, DELTA,
%   MU, and H are all in meters.  The units of G are then m^2 s / kg.
%
%   WINDTRANS(...,'free') instead returns the free slip transfer function.
%
%   For details on the expressions implemented by this function, see
%
%        Lilly, J. M. and S. Elipot (2021). A unifying perspective on
%            transfer function solutions to the unsteady Ekman problem.
%            Fluids, 6 (2): 85, 1--36.
%   __________________________________________________________________
%
%   Special forms
%
%   By default WINDTRANS uses the general transfer function.  WINDTRANS
%   will also employ limiting expressions for special cases, as follows.
%
%       H = Inf                --  Mixed Ekman / Madsen solution
%       H = Inf, MU = 0        --  Ekman solution
%       H = Inf, DELTA = 0     --  Madsen solution
%       MU = 0                 --  Finite-layer Ekman solution
%       DELTA = 0              --  Finite-layer Madsen solution
%
%   These have to be coded separately because the full solution is singular
%   in these cases.
%   __________________________________________________________________
%
%   Computational options
%
%   When computed over a wide range of parameter space, the transfer
%   function tends to be encounter numerical overflow when the arguments to
%   Bessel functions become large, causing its computation to fail.
%
%   To avoid this problems, by default WINDTRANS switches to using a highly
%   accurate thirty-term expansion about the large-argument exponential
%   behavior of the Bessel functions when their arguments exceed 10^2.9.
%
%   Two other options are available, primarily for testing purposes. Both
%   of these other algorithms lead to artifacts and are not recommended.
%
%   WINDTRANS(...'far',...) switches instead to use the (inferior) one-term
%   expansion, also known as the far-inertial limit.
%
%   WINDTRANS(...,'general',...) uses the general formula with no switch.
%
%   Note that these options only apply to no-slip solution.  The free-slip
%   solution, which is not deemed to be physically relevant, is only
%   computed with the general formula.
%
%   For details on these algorithms, see Lilly and Elipot (2021).
%   __________________________________________________________________
%
%   'windtrans --t' runs a some tests.
%
%   Usage: G=windtrans(omega,z,fc,delta,mu,h);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2019--2021 J.M. Lilly --- type 'help jlab_license' for details

%   WINDTRANS(OMEGA,Z,FC,DELTA,MU,A,'AMP') alternately parameterizes the
%   transfer function in terms of the magntide of the near-inertial peak A
%   in 1/m, together with the Ekman and Madsen depths DELTA and MU.

%   Gradient output
%
%   [G,DDELTA,DH]=WINDTRANS(OMEGA,Z,FC,DELTA,0,H) for the Ekman-like model
%   also returns the partial derivatives with respect to DELTA and H.

if strcmp(varargin{1}, '--t')
    windtrans_test,return
end

omega=varargin{1}(:)/24/3600;%convert to rad/second
z=varargin{2};
fc=varargin{3}/24/3600;%convert to rad/second
delta=varargin{4};
mu=varargin{5};
h=varargin{6};

tol=100;
rho=1027;
str='expansion';%use the 30-term tilde expansion by default
%str='two';%use the two-term tilde expansion by default
slipstr='noslip';

bool=false(2,1);
for i=7:nargin
    istr=varargin{i}(1:3);
    if strcmpi(istr,'gen')||strcmpi(istr,'eli')||strcmpi(istr,'far')||strcmpi(istr,'ada')||strcmpi(istr,'exp')||strcmpi(istr,'two')
        str=varargin{i};
        if strcmpi(str(1:3),'exp')
            bool(1)=true;
        end
    elseif strcmpi(istr,'fre')||strcmpi(istr,'nos')
        slipstr=varargin{i};
        if strcmpi(slipstr(1:3),'fre')
            bool(2)=true;
        end
    end
end

%this is just to keep track of an unsuitable combination of input arguments
if all(bool)
     disp('Sorry, the tilde expansion algorithm is not implemented for the free-slip transfer function.')
end
%these may have been either input, or reflecting the default behavior
if strcmpi(str(1:3),'exp')&&strcmpi(slipstr(1:3),'fre')
     str='gen';
end

omega1=omega;
if length(z)~=1
    [z,omega]=meshgrid(z,omega);
end

%if isinf(h)||delta==0||mu==0
%    str='lilly';%use the lilly version of expressions for special cases
%end

%str,slipstr
G=zeros(size(omega));
if strcmpi(slipstr(1:3),'nos')
    %----------------------------------------------------------------------
    %The noslip forms
        %[G,xiz,xih,xi0]=windtrans_general_noslip(delta,fc,rho,z,mu,omega,h,str);
   %     G=windtrans_general_noslip(delta,fc,rho,z,mu,omega,h,str);
    if strcmpi(str(1:3),'gen')|| strcmpi(str(1:3),'far')||strcmpi(str(1:3),'exp')||strcmpi(str(1:3),'two')
        [G,varargout{2},varargout{3}]=windtrans_lilly_noslip(delta,fc,rho,z,mu,omega,h,str);
    elseif strcmpi(str(1:3),'eli')%asymptotic forms, elipot edition
        G=windtrans_elipot_noslip(delta,fc,rho,z,mu,omega,h);
    end
elseif strcmpi(slipstr(1:3),'fre')
    %----------------------------------------------------------------------
    %The freeslip forms

    K0=frac(1,2)*delta.^2.*abs(fc);
    K1=frac(1,2)*mu*abs(fc);

    s=sign(fc).*sign(1+omega./fc);
    zo=delta.^2./mu;
    [xiz,xih,xi0]=xis(s,zo,delta,z,omega,fc,h);

    coeff=frac(sqrt(2).*rot(-s.*pi/4),delta.*rho.*abs(fc).*sqrt(abs(1+omega./fc)));
    [k0z,i0z,k1h,i1h,k10,i10]=bessels_freeslip(xiz,xih,xi0);

    if strcmpi(str(1:3),'gen')||(strcmpi(str(1:3),'lil')&&((K0~=0)&&(K1~=0)))
        numer=i0z.*k1h+i1h.*k0z;
        denom=i1h.*k10-i10.*k1h;
        G=coeff.*frac(numer,denom);
    elseif strcmpi(str(1:3),'lil')
        %    [h,K0,K1]
        if mu==0
            coeff=frac(sqrt(2).*rot(-s*pi/4),delta.*abs(fc).*rho.*sqrt(abs((1+omega./fc))));
            cosharg=sqrt(2).*rot(s*pi/4).*frac(h-z,delta).*sqrt(abs((1+omega./fc)));
            sinharg=sqrt(2).*rot(s*pi/4).*frac(h,delta).*sqrt(abs((1+omega./fc)));
            G=coeff.*frac(cosh(cosharg),sinh(sinharg));
        elseif delta==0
            coeff=frac(2,rho*K1);
            k0z=besselk(0,2*rot(s*pi/4).*sqrt(frac(z,K1./abs(fc)).*abs(1+omega./fc)));
            k1h=besselk(1,2*rot(s*pi/4).*sqrt(frac(h,K1./abs(fc)).*abs(1+omega./fc)));
            i0z=besseli(0,2*rot(s*pi/4).*sqrt(frac(z,K1./abs(fc)).*abs(1+omega./fc)));
            i1h=besseli(1,2*rot(s*pi/4).*sqrt(frac(h,K1./abs(fc)).*abs(1+omega./fc)));
            G=coeff.*(k0z+frac(k1h.*i0z,i1h));
        end
    elseif strcmpi(str(1:3),'eli')

        delta1=sqrt(2.*K0./(omega+fc));
        delta2=K1./(omega+fc);
        xize=2*sqrt(1i*(zo+z)./delta2);
        xi0e=2*sqrt(1i*zo./delta2);
        xihe=2*sqrt(1i*(zo+h)./delta2);

        if mu==0
            coeff=frac(1,rho.*sqrt(1i*(omega+fc).*K0));
            numer=cosh((1+1i).*(h-z)./delta1);
            denom=sinh((1+1i).*h./delta1);
            G=coeff.*numer./denom;
        elseif delta==0
            coeff=(2./rho./K1);
            k0z=besselk(0,2.*sqrt(1i.*z./delta2));
            k1h=besselk(1,2*sqrt(1i.*h./delta2));
            i0z=besseli(0,2*sqrt(1i.*z./delta2));
            i1h=besseli(1,2*sqrt(1i.*h./delta2));
            G=coeff.*(k0z+frac(k1h.*i0z,i1h));
        else
            coeff=frac(1,rho.*sqrt(1i.*(omega+fc).*K0));
            k0z=besselk(0,xize);
            k10=besselk(1,xi0e);
            i10=besseli(1,xi0e);
            i0z=besseli(0,xize);
            k1h=besselk(1,xihe);
            i1h=besseli(1,xihe);
            numer=i0z.*k1h+i1h.*k0z;
            denom=i1h.*k10-i10.*k1h;
            G=coeff.*frac(numer,denom);
        end
    end
end
G(z>h)=nan;

varargout{1}=G;
%varargout{2}=xih;
%varargout{3}=xiz;
%varargout{4}=xi0;

function[xiz,xih,xi0]=xis(s,zo,delta,z,omega,fc,h)

xiz=2*sqrt(2).*rot(s.*pi/4).*(zo./delta).*sqrt((1+z./zo).*abs((1+omega./fc)));
xih=2*sqrt(2).*rot(s.*pi/4).*(zo./delta).*sqrt((1+h./zo).*abs((1+omega./fc)));
xi0=2*sqrt(2).*rot(s.*pi/4).*(zo./delta).*sqrt(abs((1+omega./fc)));

function[k0z,i0z,k0h,i0h,k10,i10]=bessels_noslip(varargin)

argz=varargin{1};
argh=varargin{2};
if length(varargin)==3
    arg0=varargin{3};
end

k0z=besselk(0,argz);
i0z=besseli(0,argz);
k0h=besselk(0,argh);
i0h=besseli(0,argh);

if nargout>4
    k10=besselk(1,arg0);
    i10=besseli(1,arg0);
end

function[k0z,i0z,k0h,i0h,k10,i10]=besseltildes_noslip(argz,argh,arg0,nterms)

k0z=besselktilde(0,argz,nterms);
i0z=besselitilde(0,argz,nterms);
k0h=besselktilde(0,argh,nterms);
i0h=besselitilde(0,argh,nterms);

k10=besselktilde(1,arg0,nterms);
i10=besselitilde(1,arg0,nterms);

function[k0z,i0z,k1h,i1h,k10,i10]=bessels_freeslip(argz,argh,arg0)

k0z=besselk(0,argz);
i0z=besseli(0,argz);
k1h=besselk(1,argh);
i1h=besseli(1,argh);

if nargout>4
    k10=besselk(1,arg0);
    i10=besseli(1,arg0);
end

function[G]=windtrans_expansion_noslip(delta,fc,rho,z,mu,omega,h,str)

zo=delta.^2./mu;
s=sign(fc).*sign(1+omega./fc);
[xiz,xih,xi0]=xis(s,zo,delta,z,omega,fc,h);
coeff=frac(sqrt(2).*rot(-s.*pi/4),delta.*rho.*abs(fc).*sqrt(abs(1+omega./fc)));

if strcmpi(str(1:3),'two')
    [k0z,i0z,k0h,i0h,k10,i10]=besseltildes_noslip(xiz,xih,xi0,2);
elseif strcmpi(str(1:3),'exp')
    [k0z,i0z,k0h,i0h,k10,i10]=besseltildes_noslip(xiz,xih,xi0,30);
end

numer=exp(xi0-xiz).*i0h.*k0z-exp(xi0+xiz-2*xih).*k0h.*i0z;
denom=i0h.*k10+exp(2*xi0-2*xih).*k0h.*i10;
G=coeff.*frac(numer,denom);

bool=(omega==-fc);
if ~isinf(h)&&~isinf(zo)
    if length(z)==1
        G(bool)=frac(4*zo,rho*abs(fc)*delta.^2).*frac(sqrt(1+h./zo)-sqrt(1+z./zo),(1+z./zo).^(1/4));
    else
        G(bool)=frac(4*zo,rho*abs(fc)*delta.^2).*frac(sqrt(1+h./zo)-sqrt(1+z(bool)./zo),(1+z(bool)./zo).^(1/4));
    end
elseif isinf(h)
    G(bool)=inf;
end

function[G,xiz,xih,xi0]=windtrans_general_noslip(delta,fc,rho,z,mu,omega,h,str)

zo=delta.^2./mu;
s=sign(fc).*sign(1+omega./fc);
[xiz,xih,xi0]=xis(s,zo,delta,z,omega,fc,h);
coeff=frac(sqrt(2).*rot(-s.*pi/4),delta.*rho.*abs(fc).*sqrt(abs(1+omega./fc)));

[k0z,i0z,k0h,i0h,k10,i10]=bessels_noslip(xiz,xih,xi0);

%numerically, better to do it like this to avoid having huge numbers dominate
%tests shows this removes the large-argument overflow when only h (not z) is large

numer=k0z./k10-(k0h./k10).*(i0z./i0h);
denom=1+(k0h./k10).*(i10./i0h);
G=coeff.*frac(numer,denom);

if strcmpi(str(1:3),'exp')||strcmpi(str(1:3),'two')
    bool=log10(abs(xiz))>2.9;
    if ~isempty(find(bool,1))
        if length(z)==1
            G(bool)=windtrans_expansion_noslip(delta,fc,rho,z,mu,omega(bool),h,str);
        else
            G(bool)=windtrans_expansion_noslip(delta,fc,rho,z(bool),mu,omega(bool),h,str);
        end
    end
elseif strcmpi(str(1:3),'far')
    bool=log10(abs(xiz))>2.9;
    if ~isempty(find(bool,1))
        if length(z)==1
            G(bool)=windtrans_farinertial_noslip(delta,fc,rho,z,mu,omega(bool),h);
        else
            G(bool)=windtrans_farinertial_noslip(delta,fc,rho,z(bool),mu,omega(bool),h);
        end
    end
end
G=windtrans_inertiallimit(G,delta,fc,rho,z,mu,omega,h);

function[G]=windtrans_farinertial_noslip(delta,fc,rho,z,mu,omega,h)

zo=delta.^2./mu;
s=sign(fc).*sign(1+omega./fc);
[xiz,xih,xi0]=xis(s,zo,delta,z,omega,fc,h);
coeff=frac(sqrt(2).*rot(-s.*pi/4),delta.*abs(fc).*rho);

numer=exp(-xiz+xi0)-exp(xiz+xi0-2*xih);
denom=1+exp(2*xi0-2*xih);
G=coeff.*frac(1,sqrt(abs(1+omega./fc)).*(1+z./zo).^(1/4)).*frac(numer,denom);
%Don't apply the inertial limit!

function[G]=windtrans_inertiallimit(G,delta,fc,rho,z,mu,omega,h)

zo=delta.^2./mu;
bool=(omega==-fc);

if ~isempty(find(bool,1))
    if ~isinf(h)&&~isinf(zo)
        if length(z)==1
            %G(bool)=frac(2*zo,rho*abs(fc)*delta.^2).*log(frac(1+h./zo,1+z./zo));
            G(bool)=frac(2,rho*abs(fc)*mu).*log(frac(1+h./zo,1+z./zo));
        else
            %G(bool)=frac(2*zo,rho*abs(fc)*delta.^2).*log(frac(1+h./zo,1+z./zo));
            G(bool)=frac(2,rho*abs(fc)*mu).*log(frac(1+h./zo,1+z(bool)./zo));
        end
    elseif ~isinf(h)&&isinf(zo)
        if length(z)==1
            G(bool)=frac(2,rho*abs(fc)*delta.^2).*(h-z);
        else
            G(bool)=frac(2,rho*abs(fc)*delta.^2).*(h-z(bool));
        end
    else
        G(bool)=inf;
    end
end

function[G,ddelta,dh]=windtrans_lilly_noslip(delta,fc,rho,z,mu,omega,h,str)

zo=delta.^2./mu;
s=sign(fc).*sign(1+omega./fc);
[xiz,xih,xi0]=xis(s,zo,delta,z,omega,fc,h);
%[delta,zo,K0,K1]

if h==inf
    if mu==0
        %Ekman solution
        coeff=frac(sqrt(2).*rot(-s.*pi/4),delta.*abs(fc).*rho);
        G=coeff.*frac(exp(-(1+s.*1i).*(z./delta).*sqrt(abs((1+omega./fc)))),sqrt(abs((1+omega./fc))));
    elseif delta==0
        %Madsen solution
        coeff=frac(4,rho*abs(fc)*mu);
        G=coeff.*besselk(0,2*sqrt(2)*rot(s.*pi/4).*sqrt(frac(z,mu).*abs(1+omega./fc)));
    else
        %Mixed solution
        k0z=besselk(0,xiz);
        k10=besselk(1,xi0);
        coeff=frac(sqrt(2).*rot(-s.*pi/4),delta.*abs(fc).*rho.*sqrt(abs((1+omega./fc))));
        G=coeff.*frac(k0z,k10);
    end
else
    if mu==0
        %finite-layer Ekman
        coeff=frac(sqrt(2).*rot(-s*pi/4),delta.*abs(fc).*rho.*sqrt(abs((1+omega./fc))));
        %sinharg=sqrt(2).*rot(s*pi/4).*frac(h-z,delta).*sqrt(abs((1+omega./fc)));
        %cosharg=sqrt(2).*rot(s*pi/4).*frac(h,delta).*sqrt(abs((1+omega./fc)));
        %G=coeff.*frac(sinh(sinharg),cosh(cosharg));
        %can't do it this way because you divide two huge numbers

        argh=sqrt(2).*rot(s*pi/4).*frac(h,delta).*sqrt(abs((1+omega./fc)));
        argz=sqrt(2).*rot(s*pi/4).*frac(z,delta).*sqrt(abs((1+omega./fc)));

        numer=exp(-argz)-exp(argz).*exp(-2*argh);
        denom=1+exp(-2*argh);
        G=coeff.*frac(numer,denom);

        bool=(omega==-fc);
        if length(z)==1
            G(bool)=frac(2,rho*abs(fc)*delta.^2).*(h-z);
        else
            G(bool)=frac(2,rho*abs(fc)*delta.^2).*(h-z(bool));
        end
    elseif delta==0
        %         if z==0
        %             coeff=frac(4,rho*abs(fc)*mu);
        %             argz=2*sqrt(2)*rot(s*pi/4).*sqrt(frac(z,mu).*abs(1+omega./fc));
        %             argh=2*sqrt(2)*rot(s*pi/4).*sqrt(frac(h,mu).*abs(1+omega./fc));
        %
        %             [k0z,i0z,k0h,i0h]=bessels_noslip(argz,argh);
        %             length(find(isnan(k0z)))
        %             length(find(isnan(i0z)))
        %             length(find(isnan(k0h)))
        %             length(find(isnan(i0h)))
        %             length(find(isnan(argh)))
        %
        %             %figure,plot(abs(argh))
        %
        %             G=-coeff.*(log(xi0)+frac(1,2)*(z./zo)+frac(k0h,i0h));
        %             G=-coeff.*(log(xi0./xih)+frac(1,2)*(z./zo));
        %             length(find(isnan(xi0./xih)))
        %
        %         else
        %finite-layer Madsen
        coeff=frac(4,rho*abs(fc)*mu);
        argz=2*sqrt(2)*rot(s*pi/4).*sqrt(frac(z,mu).*abs(1+omega./fc));
        argh=2*sqrt(2)*rot(s*pi/4).*sqrt(frac(h,mu).*abs(1+omega./fc));
        [k0z,i0z,k0h,i0h]=bessels_noslip(argz,argh);
        G=coeff.*(k0z-i0z.*frac(k0h,i0h));

        bool=(omega==-fc);
        if length(z)==1
            G(bool)=frac(1,2)*coeff.*log(h./z);
        else
            G(bool)=frac(1,2)*coeff.*log(h./z(bool));
        end
        %end
    else
        %General solution
        G=windtrans_general_noslip(delta,fc,rho,z,mu,omega,h,str);
    end
end

if mu==0
    s=sign(fc).*sign(1+omega./fc);
    Gamma=sqrt(2).*rot(s*pi/4).*sqrt(abs((1+omega./fc)));
    ddelta1=(Gamma.*frac(h,delta).*tanh(Gamma.*frac(h,delta))-1).*G./delta;
    numer=exp(Gamma.*frac(-z,delta))+exp(-Gamma.*frac(2*h-z,delta));
    denom=1+exp(-Gamma.*frac(2*h,delta));
    ddelta2=-frac(2,delta.^2.*abs(fc).*rho).*frac(h-z,delta).*frac(numer,denom);
    ddelta=ddelta1+ddelta2;
    dh1=-Gamma.*frac(1,delta).*tanh(Gamma.*frac(h,delta)).*G;
    dh2=frac(2,delta.^2.*abs(fc).*rho).*frac(numer,denom);
    dh=dh1+dh2;
else
    ddelta=[];
    dh=[];
end

function[G]=windtrans_elipot_noslip(delta,fc,rho,z,mu,omega,h)

zo=delta.^2./mu;
K0=frac(1,2)*delta.^2.*abs(fc);
K1=frac(1,2)*mu*abs(fc);

delta1=sqrt(2.*K0./(omega+fc));
delta2=K1./(omega+fc);
xiz=2*sqrt(1i*(zo+z)./delta2);
xi0=2*sqrt(1i*zo./delta2);
xih=2*sqrt(1i*(zo+h)./delta2);

if h==inf
    if K1==0
       %Ekman solution
        coeff=frac(1,rho.*sqrt(1i*(omega+fc).*K0));
        G=coeff.*exp(-z.*(1+1i)./delta1);
    elseif K0==0
        %Madsen solution
        coeff=(2./rho./K1);
        k0z=besselk(0,2.*sqrt(1i.*z./delta2));
        G=coeff.*k0z;
    else
        %Mixed solution
        coeff=frac(1,rho.*sqrt(1i.*(omega+fc).*K0));
        G=coeff.*besselk(0,xiz)./besselk(1,xi0);
    end
else
    if K1==0
        %finite-layer Ekman
        coeff=frac(1,rho.*sqrt(1i*(omega+fc).*K0));
        numer=sinh((1+1i).*(h-z)./delta1);
        denom=cosh((1+1i).*h./delta1);
        G=coeff.*numer./denom;
    elseif K0==0
        %finite-layer Madsen
        coeff=(2./rho./K1);
        k0z=besselk(0,2.*sqrt(1i.*z./delta2));
        k0h=besselk(0,2*sqrt(1i.*h./delta2));
        i0z=besseli(0,2*sqrt(1i.*z./delta2));
        i0h=besseli(0,2*sqrt(1i.*h./delta2));
        G=coeff.*(k0z-frac(k0h.*i0z,i0h));
    else
        %General solution
        coeff=frac(1,rho.*sqrt(1i.*(omega+fc).*K0));
        [k0z,i0z,k0h,i0h,k10,i10]=bessels_noslip(xiz,xih,xi0);
        numer=i0h.*k0z-k0h.*i0z;
        denom=i10.*k0h+k10.*i0h;
        G=coeff.*frac(numer,denom);
    end
end

function[]=windtrans_test
windtrans_test_gradient;
windtrans_test_limits;

function[]=windtrans_test_gradient

delta=10.^(-1:0.05:3);
h=10.^(log10(15.15):.05:5);
[dg,hg]=meshgrid(delta,h);

[G,G1,G2,G3,G4,dg1,dg2]=vzeros(size(dg));
ddelta=1e-6;dh=1e-6;
omega=-1e-4;
%omega=-1.5e-4;
%omega=-5e-4;
for j=1:length(delta)
    for i=1:length(h)
        %G=windtrans(omega,z,fc,delta,mu,h);
       [G(i,j),dg1(i,j),dg2(i,j)]=windtrans(omega,15,1e-4,dg(i,j),0,hg(i,j));
        G1(i,j)=windtrans(omega,15,1e-4,dg(i,j)-ddelta/2,0,hg(i,j));
        G2(i,j)=windtrans(omega,15,1e-4,dg(i,j)+ddelta/2,0,hg(i,j));
        G3(i,j)=windtrans(omega,15,1e-4,dg(i,j),0,hg(i,j)-dh/2);
        G4(i,j)=windtrans(omega,15,1e-4,dg(i,j),0,hg(i,j)+dh/2);
    end
end

dg1hat=frac(G2-G1,ddelta);
dg2hat=frac(G4-G3,dh);

bool=dg1hat~=0;
eps1=maxmax(log10(frac(abs(dg1hat(bool)-dg1(bool)),sqrt(squared(dg1hat(bool))+squared(dg1(bool))))));
bool=dg2hat~=0;
eps2=maxmax(log10(frac(abs(dg2hat(bool)-dg2(bool)),sqrt(squared(dg2hat(bool))+squared(dg2(bool))))));

%jpcolor(log10(delta),log10(h),real(dg1))
%figure,jpcolor(log10(delta),log10(h),log10(frac(abs(dg1hat-dg1),sqrt(squared(dg1hat)+squared(dg1)))))
%jpcolor(log10(delta),log10(h),log10(frac(abs(dg2hat-dg2),sqrt(squared(dg2hat)+squared(dg2)))))

reporttest(['WINDTRANS analytic and numerical gradients match for Ekman case'],eps1<-4&&eps2<-4)

function[]=windtrans_test_limits

N=1000;
z=[0.1:1:100];
h=[inf 200];
K0=[0 1/10 1/10];
K1=[1 0 1];
fc=1e-4;
delta=sqrt(frac(2*K0,fc));
mu=frac(2*K1,fc);

[bool1,bool2,bool3,bool4]=vzeros(length(K0),2);
[Ge,Go,Gl,Ga]=vzeros(N,length(z),length(K0),2);
slipstr='noslip';
omega=fourier(1,N,'two');

for s=[1 -1]
    for i=1:length(K0)
        for j=1:2
            Ge(:,:,i,j)=windtrans(omega,z,s/2,delta(i),mu(i),h(j),'elipot',slipstr);
            Gl(:,:,i,j)=windtrans(omega,z,s/2,delta(i),mu(i),h(j),slipstr);
            Go(:,:,i,j)=windtrans(omega,z,s/2,delta(i)+1e-14,mu(i)+0.5,h(j),slipstr,'general');
            Ga(:,:,i,j)=windtrans(omega,z,s/2,delta(i)+1e-14,mu(i)+0.5,h(j),slipstr,'two');
            Ge(:,:,i,j)=Ge(:,:,i,j)./maxmax(abs(Gl(:,:,i,j)));
            Go(:,:,i,j)=Go(:,:,i,j)./maxmax(abs(Gl(:,:,i,j)));
            Ga(:,:,i,j)=Ga(:,:,i,j)./maxmax(abs(Gl(:,:,i,j)));
            Gl(:,:,i,j)=Gl(:,:,i,j)./maxmax(abs(Gl(:,:,i,j)));

            bool1(i,j)=aresame(Ge(:,:,i,j),Gl(:,:,i,j),1e-8);
            bool2(i,j)=aresame(Ge(:,:,i,j),Go(:,:,i,j),1e-1);
            bool3(i,j)=aresame(Gl(:,:,i,j),Go(:,:,i,j),1e-1);
            bool4(i,j)=aresame(Ga(:,:,i,j),Go(:,:,i,j),1e-20);
        end
    end

    if s==1
        reporttest(['WINDTRANS Elipot and Lilly forms match for f>0 and ' slipstr ' condition'],allall(bool1))
        reporttest(['WINDTRANS general and special forms match for f>0 and ' slipstr ' condition'],allall(bool2).*allall(bool3))
        reporttest(['WINDTRANS general and expansion forms match for f>0 and ' slipstr ' condition'],allall(bool4).*allall(bool3))
    else
        reporttest(['WINDTRANS Elipot and Lilly forms match for f<0 and ' slipstr ' condition'],allall(bool1))
        reporttest(['WINDTRANS general and special forms match for f<0 and ' slipstr ' condition'],allall(bool2).*allall(bool3))
        reporttest(['WINDTRANS general and expansion forms match for f>0 and ' slipstr ' condition'],allall(bool4).*allall(bool3))
    end
end

% older test for freeslip forms; not currently  particularly relevant
% h=20;
% z=[0.1:1:20];
% [bool1,bool2,bool3]=vzeros(length(K0),1);
% [Ge,Go,Gl]=vzeros(N,length(z),length(K0));
% slipstr='freeslip';
% for s=[1 -1]
%     for i=1:length(K0)
%         for j=1:length(h)
%             i,j,s
%             Ge(:,:,i,j)=windtrans(omega,z,s/2,delta(i),mu(i),h(j),'elipot',slipstr);
%             Gl(:,:,i,j)=windtrans(omega,z,s/2,delta(i),mu(i),h(j),'lilly',slipstr);
%             Go(:,:,i,j)=windtrans(omega,z,s/2,delta(i)+1e-10,mu(i)+0.5,h(j),slipstr,'general');
%             Ge(:,:,i,j)=Ge(:,:,i,j)./maxmax(abs(Gl(:,:,i,j)));
%             Go(:,:,i,j)=Go(:,:,i,j)./maxmax(abs(Gl(:,:,i,j)));
%             Gl(:,:,i,j)=Gl(:,:,i,j)./maxmax(abs(Gl(:,:,i,j)));
%             bool1(i,j)=aresame(Ge(:,:,i,j),Gl(:,:,i,j),1e-8);
%             bool2(i,j)=aresame(Ge(:,:,i,j),Go(:,:,i,j),1e-1);
%             bool3(i,j)=aresame(Gl(:,:,i,j),Go(:,:,i,j),1e-1);
%         end
%     end
%
%     if s==1
%         reporttest(['WINDTRANS Elipot and Lilly forms match for f>0 and ' slipstr ' condition'],allall(bool1))
%         reporttest(['WINDTRANS general and special forms match for f>0 and ' slipstr ' condition'],allall(bool2).*allall(bool3))
%     else
%         reporttest(['WINDTRANS Elipot and Lilly forms match for f<0 and ' slipstr ' condition'],allall(bool1))
%         reporttest(['WINDTRANS general and special forms match for f<0 and ' slipstr ' condition'],allall(bool2).*allall(bool3))
%     end
% end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function tips = correct_rv_crescent_tips(tips, valve, dthr)
% CORRECT_RV_CRESCENT_TIPS  Correct the segmentation of the Right
% Ventricle's crescent tips using the segmentation of the tricuspid and
% pulmonary annula
%
% TIPS = CORRECT_RV_CRESCENT_TIPS(TIPS, VALVE)
%
%   TIPS is a 3-column matrix with the coordinates of all the points in the
%   Right Ventricle crest, that corresponds to the crescent shape tips when
%   you cut the ventricle by an axial plane.
%
%   VALVE is a 3-colum matrix with the coordinates of either the tricuspid
%   valve or pulmonary valve segmetations.
%
%   The points in the segmentations are assumed to be consecutive.
%
% TIPS = CORRECT_RV_CRESCENT_TIPS(TIPS, VALVE, DTHR)
%
%   DTHR is the threshold distance value between the crest line and the
%   valve.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2010 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
error( nargchk( 2, 3, nargin, 'struct' ) );
error( nargoutchk( 0, 1, nargout, 'struct' ) );

% number of crest points
P = size(tips, 1);

% compute distance from each crest point to the valve
d = inf(P, 1);
for I = find(~isnan(sum(tips, 2)));
    d(I) = min(dmatrix(valve', tips(I,:)'));
end

% default value of dthr
if (nargin < 3 || isempty(dthr))
    dthr = min(d) * 5;
end

% find the first and last times the crest gets within a small distance of
% the valve
idx = find(d <= dthr);
if isempty(idx)
    error('The crest does not get close enough to the valve')
end
t1 = idx(1);
t2 = idx(end);

% % find one of the intersection points of the valve with the crest (crest
% % point index)
% [aux, t1] = min(d);
%
% % we are going to ignore a bunch of points either side of the intersection
% % point
% L = round(size(valve, 1)/4);
% d2 = d;
% d2(t1-L:t1+L) = Inf;
%
% % find second intersection point (crest point index)
% [aux, t2] = min(d2);
% aux = sort([t1 t2]);
% t1 = aux(1);
% t2 = aux(2);

% find intersection points (valve point index)
[aux, v1] = min(dmatrix(valve', tips(t1,:)'));
[aux, v2] = min(dmatrix(valve', tips(t2,:)'));

% find which half of the valve is closer to the centroid
aux = sort([v1 v2]);
v1 = aux(1);
v2 = aux(2);
if (mean(dmatrix(valve(v1:v2, :)', [0;0;0])) < ...
        mean(dmatrix(valve([v2:end 1:v1], :)', [0;0;0])))
    valveseg = valve(v1:v2, :);
else
    valveseg = valve([v2:end 1:v1], :);
end

% if the first point we find in the valve segment is closer to the _last_
% point we find in the crest gap, then we have to invert the segment's
% order
if (dmatrix(valveseg(1, :)', tips(t1,:)') > ...
        dmatrix(valveseg(end, :)', tips(t1,:)'))
    valveseg = valveseg(end:-1:1, :);
end

% replace the incorrect crest segment by the valve segment
tips = [tips(1:t1, :); valveseg; tips(t2:end, :)];

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [solversToUse] = prepareTest(varargin)
% Checks the prerequisites of the test, and returns solvers depending on
% the input parameters. If the requirements are NOT met, a
% COBRA:RequirementsNotMet error is thrown.
%
% USAGE:
%    [solversToUse] = prepareTest(varargin)
%
% OPTIONAL INPUTS:
%    varagin:      `ParameterName` value pairs with the following options:
%
%                   - `toolboxes` or `requiredToolboxes`: Names of required toolboxes (the license feature name).(default: {})
%                   - `minimalMatlabSolverVersion`: Minimal version of the optimization toolbox required to use the matlab solver.(default: 0)
%                   - `requiredSolvers`: Names of all solvers that MUST be available. If not empty, the resulting solvers struct will contain cell arrays (default: {})
%                   - `useSolversIfAvailable`: Names of solvers that should be used if available. If not empty, the resulting solvers struct will contain cell arrays (will not throw an error if not). (default: {})
%                   - `requireOneSolverOf`: Names of solvers, at least one of which has to be available
%                   - `excludeSolvers`: Names of solvers which should never be used for the test (because they fail)
%                   - `useMinimalNumberOfSolvers`: Always use only one solver. This option allows tests which only use FBA to generate input, but where no expicit linear programming code is present, to only validate on one solver (default: false).
%                   - `needsLP`: Whether a LP solver is required (default: false)
%                   - `needsMILP`: Whether a MILP solver is required (default: false)
%                   - `needsQP`: Whether a QP solver is required (default: false)
%                   - `needsMIQP`: Whether a MIQP solver is required (default: false)
%                   - `needsNLP`: Whether a NLP solver is required (default: false)
%                   - `needsEP`: Whether a EP solver is required (default: false)
%                   - `needsUnix`: Whether the test only works on a Unix system (macOS or Linux) (default: false)
%                   - `needsWindows`: Whether the test only works on a Windows system (default: false)
%                   - `needsMac`: Whether the test only works on a Mac system (default: false)
%                   - `needsLinux`: Whether the test only works on a Linux system (default: false)
%                   - `needsWebAddress`: Tests, whether the supplied url exists (default: '')
%                   - `needsWebRead`: Tests, whether webread can be used with the given url
%
% OUTPUTS:
%
%    solversToUse:  A struct with one field per solver type listing
%                   the solvers to use for that type of problem in a cell array.
%                   If neither of the 'useIfAvailable' nor the 'reqSolvers'
%                   parameter is provided, only at most one solver
%                   per type will be returned (i.e. the default
%                   solver for that type). See the examples below
%
% EXAMPLE:
%
%      % request a check for the parallel processing toolbox
%      >> solvers = prepareTest('requiredToolboxes', {'distrib_computing_toolbox'})
%      solvers =
%
%                struct with fields:
%
%                      LP: {'gurobi'}
%                    MILP: {'gurobi'}
%                      QP: {'gurobi'}
%                    MIQP: {'gurobi'}
%                     NLP: {'matlab'}
%                      EP: {'pdco'}
%
%      % request gurobi, ibm_cplex and tomlab if available
%      >> solvers = prepareTest('useIfAvailable', {'tomlab', 'ibm_cplex', 'gurobi'})
%      solvers =
%
%                struct with fields:
%
%                      LP: {2×1 cell}
%                    MILP: {2×1 cell}
%                      QP: {2×1 cell}
%                    MIQP: {'gurobi'}
%                     NLP: {'matlab'}
%                      EP: {'pdco'}
%
%

global CBT_MISSING_REQUIREMENTS_ERROR_ID
global OPT_PROB_TYPES

% some Matlab Toolboxes currently in use in the COBRA Toolbox.
% This might have to be extended in the future.
toolboxInfo = struct('statistics_toolbox', {{'Statistics and Machine Learning Toolbox', 'Statistics Toolbox'}}, ...
                     'bioinformatics_toolbox', {{'Bioinformatics Toolbox'}}, ...
                     'distrib_computing_toolbox', {{'Parallel Computing Toolbox'}}, ...
                     'optimization_toolbox', {{'Optimization Toolbox'}}, ...
                     'global_optimization_toolbox', {{'Global Optimization Toolbox'}}, ...
                     'image_toolbox', {{'Image Processing Toolbox'}}, ...
                     'gads_toolbox', {{'Global Optimization Toolbox'}});

%generate this on the fly rather than using a persistient variable - Ronan
availableSolvers = getAvailableSolversByType();

parser = inputParser();
parser.addParamValue('toolboxes', {}, @iscell);
parser.addParamValue('requiredToolboxes', {}, @iscell);
parser.addParamValue('minimalMatlabSolverVersion',0,@isnumeric);

parser.addParamValue('requiredSolvers', {}, @iscell);
parser.addParamValue('useSolversIfAvailable', {}, @iscell);
parser.addParamValue('requireOneSolverOf', {}, @iscell);
parser.addParamValue('excludeSolvers', {}, @(x) iscell(x) || ischar(x) );
parser.addParamValue('useMinimalNumberOfSolvers', false, @(x) islogical(x) || x == 1 || x == 0);
parser.addParamValue('needsLP', false, @(x) islogical(x) || x == 1 || x == 0);
parser.addParamValue('needsMILP', false, @(x) islogical(x) || x == 1 || x == 0);
parser.addParamValue('needsNLP', false, @(x) islogical(x) || x == 1 || x == 0);
parser.addParamValue('needsQP', false, @(x) islogical(x) || x == 1 || x == 0);
parser.addParamValue('needsMIQP', false, @(x) islogical(x) || x == 1 || x == 0);
parser.addParamValue('needsEP', false, @(x) islogical(x) || x == 1 || x == 0);
parser.addParamValue('needsUnix', false, @(x) islogical(x) || x == 1 || x == 0);
parser.addParamValue('needsLinux', false, @(x) islogical(x) || x == 1 || x == 0);
parser.addParamValue('needsWindows', false, @(x) islogical(x) || x == 1 || x == 0);
parser.addParamValue('needsMac', false, @(x) islogical(x) || x == 1 || x == 0);
parser.addParamValue('needsWebAddress', '', @ischar);
parser.addParamValue('needsWebRead', false, @(x) islogical(x) || x == 1 || x == 0);

parser.parse(varargin{:});

useQP = parser.Results.needsQP;
useLP = parser.Results.needsLP;
useMIQP = parser.Results.needsMIQP;
useNLP = parser.Results.needsNLP;
useMILP = parser.Results.needsMILP;
useEP = parser.Results.needsEP;

macOnly = parser.Results.needsMac;
windowsOnly = parser.Results.needsWindows;
unixOnly = parser.Results.needsUnix;
linuxOnly = parser.Results.needsLinux;

toolboxes = union(parser.Results.toolboxes, parser.Results.requiredToolboxes);
requiredSolvers = parser.Results.requiredSolvers;
possibleSolvers = parser.Results.requireOneSolverOf;
excludedSolvers = parser.Results.excludeSolvers;
if ischar(excludedSolvers)
    excludedSolvers = {excludedSolvers};
end
preferredSolvers = parser.Results.useSolversIfAvailable;

needsWebAddress = parser.Results.needsWebAddress;
needsWebRead = parser.Results.needsWebRead;
useMinimalNumberOfSolvers = parser.Results.useMinimalNumberOfSolvers;
runtype = getenv('CI_RUNTYPE');

minimalMatlabSolverVersion = parser.Results.minimalMatlabSolverVersion;

errorMessage = {};
infoMessage = {};

% first, check whether the OS is applicable
if macOnly
    if ~ismac
        errorMessage{end + 1} = 'This test only works on macOS';
    end
end

if windowsOnly
    if ~ispc
        errorMessage{end + 1} = 'This test only works on Windows';
    end
end

if linuxOnly
    if ~isunix
        errorMessage{end + 1} = 'This test only works on Linux Systems';
    else
        if ~strcmp(computer('arch'), 'glnx64')
            errorMessage{end + 1} = 'This test only works on Linux Systems';
        end
    end
end

if ~isempty(needsWebAddress)
    [status_curl, result_curl] = system(['curl -s -k ' needsWebAddress]);
    if status_curl ~= 0 || isempty(result_curl)
        errorMessage{end + 1} = sprintf('This function needs to connect to %s and was unable to do so.',needsWebAddress);
    end
    if needsWebRead
        if verLessThan('MATLAB','9.3') && isunix && strncmp(needsWebAddress,'https',5)
            errorString = sprintf(['This function needs to connect to a ''https'' address using webread.\n', ...
                                   'Your MATLAB version is shipped with an invalid libssl.so.1.0.0 \n',...
                                   'which will cause MATLAB to crash if webread is called with an \n',...
                                   '''https'' website.\n',...
                                   'To fix this, you can replace your MATLAB library with the system library \n',...
                                   'by running the following commands in the console:\n',...
                                   '$ sudo mv %s/bin/glnxa64/libssl.so.1.0.0 %s/bin/glnxa64/libssl.so.1.0.0.old\n',...
                                   '$ sudo cp /lib/x86_64-linux-gnu/libssl.so.1.0.0 %s/bin/glnxa64/libssl.so.1.0.0\n',...
                                   'Please note that this test will not be able to run on your system,\n',...
                                   'regardless on whether you fixed the library or not. If you want to run it,',...
                                   'you will have to remove the ''needsWebRead'' flag from the ''prepareTest''',...
                                   'statement in the test and run it again.'],...
                                   matlabroot,matlabroot,matlabroot);
            errorMessage{end + 1} = errorString;
        end
    end
end

if unixOnly
    if ~isunix
        errorMessage{end + 1} = 'This test only works on Unix Systems (Mac and Linux)';
    end
end

% restrict the solvers available for this test
solversForTest = availableSolvers;

if any(ismember(solversForTest.LP,'matlab'))
    boxes = ver();
    optBox = find(ismember({boxes.Name},'Optimization Toolbox'));
    if ~isempty(optBox)
        optVer = boxes(optBox).Version;
        if str2double(optVer) < minimalMatlabSolverVersion
            excludedSolvers = [columnVector(excludedSolvers);'matlab'];
        end
    end
end

if ~isempty(excludedSolvers)
    solverTypes = fieldnames(availableSolvers);
    for i = 1:numel(solverTypes)
        excludedPos = ismember(solversForTest.(solverTypes{i}),excludedSolvers);
        solversForTest.(solverTypes{i}) = solversForTest.(solverTypes{i})(~excludedPos);
    end
    if numel(excludedSolvers) == 1
        infoMessage{end+1} = sprintf('%s is not compatible with the tested function and thus excluded.',excludedSolvers{1});
    else
        infoMessage{end+1} = sprintf('The following solvers are not compatible with this test and therefore excluded:\n%s and %s', strjoin(excludedSolvers(1:end-1),' and '),excludedSolvers{end});
    end
end

% then, check the required Solvers
if ~isempty(requiredSolvers) && ~all(ismember(requiredSolvers, solversForTest.ALL))
    % we have required solvers and some are missing
    missing = ~ismember(requiredSolvers, solversForTest.ALL);
    if sum(missing) == 1
        misssolver = requiredSolvers{missing};
            errorMessage{end + 1} = sprintf('%s is a required solver for the test and not available on your system.', misssolver);
    else
        errorMessage{end + 1} = sprintf('%s are missing solvers required for the test but not available on your system.', strjoin(requiredSolvers(missing), ' and '));
    end
else
    % otherwise add the required Solvers to the preferred solvers.
    preferredSolvers = union(preferredSolvers, requiredSolvers);
end

if ~isempty(possibleSolvers) && ~any(ismember(possibleSolvers,solversForTest.ALL))
    if numel(possibleSolvers) == 1
        errorMessage{end + 1} = sprintf('The test requires that the following solver is installed:\n%s', strjoin(possibleSolvers, ' or '));
    else
        errorMessage{end + 1} = sprintf('The test requires that at least one of the following solvers is installed:\n%s or %s', strjoin(possibleSolvers(1:end-1), ', '),possibleSolvers{end});
    end
else
    if ~isempty(possibleSolvers)
        %We have a set of possible solvers.
        %So we restrict the preferredSolvers to those
        % if there are preferred solvers.
        if ~isempty(preferredSolvers)
            preferredSolvers = intersect(preferredSolvers,possibleSolvers);
        else
            solverOptions = intersect(possibleSolvers,solversForTest.ALL);
            preferredSolvers = solversForTest.ALL(find(ismember(solversForTest.ALL,possibleSolvers),1));
        end

    end
end

% check the Toolboxes
res = ver;
missingTBs = struct('License', {{}}, 'Installation', {{}});
for i = 1:numel(toolboxes)
    tbstring = lower(toolboxes{i});
    licpres = license('test', tbstring);
    if any(ismember(tbstring, fieldnames(toolboxInfo)))
        tbpres = any(ismember(toolboxInfo.(lower(toolboxes{i})), {res.Name}));
    else
        % rely on the license
        tbpres = licpres;
    end
    if ~tbpres
        missingTBs.Installation{end + 1} = tbstring;
    end
    if ~licpres
        missingTBs.License{end + 1} = tbstring;
    end
end

% append the error message
if ~isempty(missingTBs.License)
    errorMessage{end + 1} = sprintf('The test requires licenses for the following Toolboxes: %s', strjoin(missingTBs.License, ' and '));
end
if ~isempty(missingTBs.Installation)
    errorMessage{end + 1} = sprintf('The test requires the following Toolboxes to be installed: %s', strjoin(missingTBs.Installation, ' and '));
end

% set up default solvers and test whether the test is useable
if isempty(solversForTest.LP)
    if useLP
        errorMessage{end + 1} = 'The test requires at least one LP solver but no solver is installed';
    end
else
    if ~isempty(solversForTest.LP)
        defaultLPSolver = solversForTest.LP{1};
    else
        defaultLPSolver = '';
    end
end

if isempty(solversForTest.QP)
    if useQP
        errorMessage{end + 1} = 'The test requires at least one QP solver but no solver is installed';
    end
else
    if ~isempty(solversForTest.QP)
        defaultQPSolver = solversForTest.QP{1};
    else
        defaultQPSolver = '';
    end
end

if isempty(solversForTest.MILP)
    if useMILP
        errorMessage{end + 1} = 'The test requires at least one MILP solver but no solver is installed';
    end
else
    if ~isempty(solversForTest.MILP)
        defaultMILPSolver = solversForTest.MILP{1};
    else
        defaultMILPSolver = '';
    end

end

if isempty(solversForTest.MIQP)
    if useMIQP
        errorMessage{end + 1} = 'The test requires at least one MIQP solver but no solver is installed';
    end
else
    if ~isempty(solversForTest.MIQP)
        defaultMIQPSolver = solversForTest.MIQP{1};
    else
        defaultMIQPSolver = '';
    end

end

if isempty(solversForTest.NLP)
    if useNLP
        errorMessage{end + 1} = 'The test requires at least one NLP solver but no solver is installed';
    end
else
    if ~isempty(solversForTest.NLP)
        defaultNLPSolver = solversForTest.NLP{1};
    else
        defaultNLPSolver = '';
    end
end

if isempty(solversForTest.EP)
    if useEP
        errorMessage{end + 1} = 'The test requires at least one EP solver but no solver is installed';
    end
else
    if ~isempty(solversForTest.EP)
        defaultEPSolver = solversForTest.EP{1};
    else
        defaultEPSolver = '';
    end
end

if ~isempty(errorMessage)
    errorString = strjoin(errorMessage, '\n');
    infoString = strjoin(infoMessage, '\n');
    error(CBT_MISSING_REQUIREMENTS_ERROR_ID, strjoin({errorString,infoString}, '\n'));
end

% collect the Used Solvers.
solversToUse = struct();
problemTypes = OPT_PROB_TYPES;
% if this is the extensive test suite, and the solver use not just about
% testing whether the actual work succeeded.
if strcmpi(runtype, 'extensive') && ~useMinimalNumberOfSolvers
    solversToUse = solversForTest;
    if 0
        % exclude pdco if not explicitly requested and available, as it has issues at the moment.
        if ~any(ismember('pdco',preferredSolvers)) && any(ismember('pdco',solversToUse.LP))
            solversToUse.LP(ismember(solversToUse.LP,'pdco')) = [];
            solversToUse.QP(ismember(solversToUse.LP,'pdco')) = [];
        end
    end
    if ~isempty(possibleSolvers)
        %Restrict to the possibleSolvers
        for i = 1:numel(problemTypes)
            solversToUse.(problemTypes{i}) = intersect(solversToUse.(problemTypes{i}),possibleSolvers);
        end
    end
else
    for i = 1:numel(problemTypes)
        solversToUse.(problemTypes{i}) = intersect(preferredSolvers, solversForTest.(problemTypes{i}));
        if isempty(solversToUse.(problemTypes{i})) && ~isempty(solversForTest.(problemTypes{i}))
            if isempty(preferredSolvers)
                eval(['solversToUse.' problemTypes{i} ' = {default' problemTypes{i} 'Solver};']);
            else
                if isempty(solversForTest.(problemTypes{i}))
                    % no solver exists, the cell array is empty.
                    solversToUse.(problemTypes{i}) = {};
                else
                    % a solver exists, provide it.
                    eval(['solversToUse.' problemTypes{i} ' = {default' problemTypes{i} 'Solver};']);
                end
            end
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function [M0,q0,X,x,f,M1,L] = spm_fp(M,x,u)
% Fokker-Planck operators and equilibrium density for dynamic systems
% FORMAT [M0,q0,X,x,f,M1,L] = spm_fp(M,x,u)
%--------------------------------------------------------------------------
% M   - model specification structure
% Required fields:
%    M.f   - dx/dt    = f(x,u,P) or f(x,u,a,P)  {function string or m-file}
%    M.g   - y(t)     = g(x,u,P)                {function string or m-file}
%    M.m   - m inputs
%    M.n   - n states
%    M.l   - l outputs
%    M.x   - (n x 1) = x(0) = expansion point
%    M.W   - (n x n) - precision matrix of state noise
% x    - cell array of vectors specifying evaluation grid
% u    - expansion point for inputs or causes;
%
% M0   - 1st order FP operator dq/dt = M0*q + u*M1*q,  q = p(X);
% q0   - stable or equilibrium mode: M0*q0 = 0
% X    - evaluation points of state space
% x    - cell array of vectors specifying evaluation grid
% f    - flow
% M1   - 2nd order FP operator
% L    - output matrix                 <y> = L*q;
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

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

% default: first level of hierarchical model
%--------------------------------------------------------------------------
M   = M(1);

% default expansion point for inputs
%--------------------------------------------------------------------------
if nargin < 3
    u0 = zeros(M.m,1);
else
    u0 = u;
end

% event space: get or create X - coordinates of evaluation grid
%--------------------------------------------------------------------------
n  = length(M.x);
if nargin < 2

    % use M.X
    %----------------------------------------------------------------------
    X     = M.X;
    for i = 1:size(X,2)
        x{i} = unique(X(:,i));
    end
else

    % use x
    %----------------------------------------------------------------------
    [X,x] = spm_ndgrid(x);
end

% f(x,0)
%--------------------------------------------------------------------------
N     = length(X);
f     = sparse(n,N);
try
    for i = 1:N
        f(:,i) = feval(M.f,X(i,:)',u0,M.pE);
    end
catch
    for i = 1:N
        f(:,i) = feval(M.f,X(i,:)',u0,0,M.pE);
    end
end

% Fokker-Planck operator; i.e., Jacobian J - M0
%==========================================================================
for i = 1:n
    d(i)  = length(x{i});
    dx(i) = x{i}(2) - x{i}(1);
    I{i}  = speye(d(i),d(i));
end
J     = sparse(N,N);

% cycle over dimensions of state space
%--------------------------------------------------------------------------
for i = 1:length(x)

    % differential operators (for positive and negative flows)
    %----------------------------------------------------------------------
    j          = find(f(i,:) < 0);
    u          = sparse(j,1,f(i,j),N,1);
    j          = find(f(i,:) > 0);
    v          = sparse(j,1,f(i,j),N,1);

    dp         = spdiags(ones(d(i),1)*[1 -1],[ 0;1],d(i),d(i))/dx(i);
    dp(:,1)    = 0;
    dq         = spdiags(ones(d(i),1)*[1 -1],[-1;0],d(i),d(i))/dx(i);
    dq(:,end)  = 0;

    % Kronecker tensor products
    %----------------------------------------------------------------------
    Dp         = 1;
    Dq         = 1;
    for j = 1:(i - 1)
        Dp = spm_kron(I{j},Dp);
        Dq = spm_kron(I{j},Dq);
    end
    Dp    = spm_kron(dp,Dp);
    Dq    = spm_kron(dq,Dq);
    for j = (i + 1):n
        Dp = spm_kron(I{j},Dp);
        Dq = spm_kron(I{j},Dq);
    end
    DP{i} = Dp;

    % augment Jacobian
    %----------------------------------------------------------------------
    J     = J + Dp*spdiags(u,0,N,N) + Dq*spdiags(v,0,N,N);

end

% dispersion
%--------------------------------------------------------------------------
C     = inv(M.W);
if length(C) ~= n
    C = C(1)*speye(n,n);
end
for i = 1:n
    for j = 1:n
        if C(i,j)
            J  = J - C(i,j)*DP{i}*DP{j}'/2;
        end
    end
end

% return if only M0 is required
%--------------------------------------------------------------------------
M0    = J;
if nargout == 1, return, end

% stable mode - stochastic iteration
%==========================================================================
q     = sparse(N,1) + 1/N;
for i = 1:1024
    dq = M0*q;
    if norm(dq,1) < exp(-4), break, end
    while min(q + dq) < 0
        dq = dq - dq/2;
    end
    q  = q + dq;
end
for i = 1:n
    nx{i} = length(x{i});
end
q0    = reshape(q,nx{:});

if nargout < 6,  return,  end

% input induced changes to Jacobian dJ/du - M1
%==========================================================================
M1    = spm_diff('spm_fp',M,x,u0,3);
if nargout < 6,  return,  end

% Output matrix L - M1
%==========================================================================

% l(x,0)
%--------------------------------------------------------------------------
L     = sparse(M.l,N);
for i = 1:N
    try
        L(:,i) = feval(M.g,X(i,:)',u0,M.pE);
    catch
        L(:,i) = feval(M.g,X(i,:)',u0,[],M.pE);
    end
end

return

% NOTES: analytic solution for q0
%--------------------------------------------------------------------------
q   = length(M0);
J   = [M0; ones(1,q)];
q0  = sum(inv(J'*J),2);

author: spm


text:
function [J, grad] = cofiCostFunc(params, Y, R, num_users, num_movies, ...
                                  num_features, lambda)
%COFICOSTFUNC Collaborative filtering cost function
%   [J, grad] = COFICOSTFUNC(params, Y, R, num_users, num_movies, ...
%   num_features, lambda) returns the cost and gradient for the
%   collaborative filtering problem.
%

% Unfold the U and W matrices from params
X = reshape(params(1:num_movies*num_features), num_movies, num_features);
Theta = reshape(params(num_movies*num_features+1:end), ...
                num_users, num_features);

% You need to return the following values correctly
J = 0;
X_grad = zeros(size(X));
Theta_grad = zeros(size(Theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost function and gradient for collaborative
%               filtering. Concretely, you should first implement the cost
%               function (without regularization) and make sure it is
%               matches our costs. After that, you should implement the
%               gradient and use the checkCostFunction routine to check
%               that the gradient is correct. Finally, you should implement
%               regularization.
%
% Notes: X - num_movies  x num_features matrix of movie features
%        Theta - num_users  x num_features matrix of user features
%        Y - num_movies x num_users matrix of user ratings of movies
%        R - num_movies x num_users matrix, where R(i, j) = 1 if the
%            i-th movie was rated by the j-th user
%
% You should set the following variables correctly:
%
%        X_grad - num_movies x num_features matrix, containing the
%                 partial derivatives w.r.t. to each element of X
%        Theta_grad - num_users x num_features matrix, containing the
%                     partial derivatives w.r.t. to each element of Theta
%

% FIXME(SaveTheRbtz@): Not optical: preforms calculations on cells with R(i,j) == 0
J = sum(sum((R==1) .* ((X * Theta' - Y) .^ 2))) / 2;

X_grad = (R==1) .* (X * Theta' - Y) * Theta + lambda * X;
Theta_grad = (R==1)' .* (X * Theta' - Y)' * X + lambda * Theta;

Regularization = lambda * (sum(sum(Theta .^ 2)) + sum(sum(X .^ 2))) / 2;
J += Regularization;

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

grad = [X_grad(:); Theta_grad(:)];

end

author: SaveTheRbtz

repo: ml-class


text:
function [muX,SigmaX,suffStat] = VBA_EKF(y,u,posterior,dim,options,flag)
% standard EKF & computation of the predictive density
% function [muX,SigmaX,suffStat] = VBA_EKF(y,u,posterior,dim,options,flag)
%
% This function inverts any nonlinear state-space model of the form:
%   y_t = g( x_t,u_t,phi ) + e_t
%   x_t+1 = f( x_t,u_t,theta ) + f_t
% using a standard extended Kalman filter (EKF).
%
% IN :          [ see VBA_NLStateSpaceModel.m ]
%   - y: pxn_t mesurements matrix
%   - u: mxn_t known input matrix (which is required as an argument in
%   the obvservation/evolution functions)
%
%   - posterior: posterior pdf structure (see VBA_NLStateSpaceModel.m).
%   This is used to extract first and second order parameters of the
%   generative model (theta, phi, alpha, sigma), as well as the initial
%   conditions of the hidden states (X0).
%   - dim: a structure variable containing the dimensions of the 3 sets of
%   model unknown variables (see VBA_NLStateSpaceModel.m)
%   - options: user-defined structure containing specific informations
%   regarding the model, ie (see VBA_check.m):
%       .f_fname (resp. g_fname): name/handle of the function that outputs
%       the evolution (resp. observation) of the hidden states.
%       .u0: the mx1 initial value of the input {0}
%       .inF: a (possibly structure) variable containing the additional
%       (internal) fixed parameters which may have to be sent to the
%       evolution function (eg pointing to different variants) {[]}
%       .inG: idem for the observation function {[]}
%   -  flag: a switch for computing just the mode of the deterministic
%   predictive density (flag=0), the 1st and 2d order statistics of the EKF
%   ({flag=1}), or the 1st and 2d order statistics of the predictive
%   density (flag = 2).
%
% OUT:
%   - muX: posterior mean of the hidden states X (nxn_t matrix)
%   - SigmaX: covariance matrices of the variational posterior pdf of
%       the dynamic hidden-states.

if numel(options.sources) > 1
    error('*** EKF is not yet compatible with multisource observations');
end

if options.sources.type > 1
    error('*** EKF is not yet compatible with multinomial observations');
end

% By default, this function implements an EKF:
if ~exist('flag','var') || isempty(flag)
    flag = 1;
end

% This checks and fills in required dummy variables
if isempty(u)
    u = zeros(1,dim.n_t);
end
[dim.p,dim.n_t] = size(y);
try
    dim.u;
catch
    dim.u = size(u,1);
end
if isfield(options,'microU') && options.microU
    u = VBA_getU(u,options,dim,'2macro');
end
if ~isfield(options,'nout_f')
    options.nout_f = nargout(options.f_fname);
end
if ~isfield(options,'nout_g')
    options.nout_g = nargout(options.g_fname);
end
if ~isfield(options,'OnLine')
    options.OnLine = 0;
end
try
    X0 = posterior.muX0;
    SigmaX0 = posterior.SigmaX0;
catch
    X0 = zeros(dim.n,1);
    SigmaX0 = zeros(dim.n,dim.n);
end
try
    theta = posterior.muTheta;
catch
    theta = [];
end
try
    phi = posterior.muPhi;
catch
    phi = [];
end
try
    iQx = options.priors.iQx;
    iQy = options.priors.iQy;
catch
    iQx = cell(dim.n_t,1);
    iQy = cell(dim.n_t,1);
    for t= 1:dim.n_t
        iQx{t} = eye(dim.n);
        iQy{t} = eye(dim.p);
    end
end

switch flag
    case 0
        suffStat = [];
        str = 'deterministic time series';
    case {1,2}
        try
            alpha = posterior.a_alpha(end)./posterior.b_alpha(end);
            if options.sources.type == 0
                sigma = posterior.a_sigma(end)./posterior.b_sigma(end);
            end
        catch
            error('Not enough info in posterior structure!')
        end
        mStar = zeros(dim.n,dim.n_t);
        SigmaX = cell(dim.n_t,1);
        %--- Initialize sufficient statistics time-series ---%
        suffStat = VBA_getSuffStat(options);
        if isequal(flag,1)
            str = 'standard EKF';
        else
            str = 'predictive density';
        end
end
muX = zeros(dim.n,dim.n_t);
gx = zeros(dim.p,dim.n_t);

% First time iteration (from initial conditions)
if ~options.OnLine && options.verbose
    fprintf(1,['Deriving ',str,' ...'])
end
if flag>=1
    %--- Prediction
    [fx0,dF_dX0] = VBA_evalFun('f',X0,theta,u(:,1),options,dim,1);
    mStar(:,1) = fx0;
    Rp = dF_dX0'*SigmaX0*dF_dX0 + 1./alpha.*VBA_inv(iQx{1});
    if flag == 1 % EKF update
        [gx(:,1),dG_dX] = VBA_evalFun('g',mStar(:,1),phi,u(:,1),options,dim,1);
        iRp = pinv(Rp);
        C =  dG_dX*iQy{1}*dG_dX';
        iSX = iRp + sigma*C;
        SigmaX{1} = pinv( iSX );
        muX(:,1) = mStar(:,1) + sigma.*SigmaX{1}*dG_dX*iQy{1}* (y(:,1)-gx(:,1));
    else % Predictive density
        muX(:,1) = mStar(:,1);
        SigmaX{1} = Rp;
    end
    % get predicted observation at the mode
    [gx(:,1),dG_dX] = VBA_evalFun('g',muX(:,1),phi,u(:,1),options,dim,1);
    suffStat.dy(:,1) = y(:,1) - gx(:,1);
    if options.sources.type == 0
        suffStat.vy(:,1) = diag( sigma.^-1.*pinv(iQy{1}) + dG_dX'*SigmaX{1}*dG_dX );
        suffStat.dy2 = suffStat.dy2 + suffStat.dy(:,1)'*iQy{1}*suffStat.dy(:,1);
    else
        suffStat.vy(:,1) = gx(:,1).*(1-gx(:,1));
        suffStat.logL = y(:,1)'*log(gx(:,1)) + (1-y(:,1))'*log(1-gx(:,1));
    end
    suffStat.dx(:,1) = muX(:,1) - fx0;
    suffStat.dx2 = suffStat.dx2 + suffStat.dx(:,1)'*iQx{1}*suffStat.dx(:,1);
else
    muX(:,1) = VBA_evalFun('f',X0,theta,u(:,1),options,dim);
end

% Loop over time samples
if ~options.OnLine && options.verbose
    fprintf(1,'%6.2f %%',0)
end
for t = 1:dim.n_t-1
    if flag >= 1
        %-- Prediction
        [fx,dF_dX] = VBA_evalFun('f',muX(:,t),theta,u(:,t+1),options,dim,t+1);
        mStar(:,t+1) = fx;
        Rp = dF_dX'*SigmaX{t}*dF_dX + 1./alpha.*VBA_inv(iQx{t+1});
        if flag == 1    % EKF update
            [gx(:,t+1),dG_dX] = VBA_evalFun('g',mStar(:,t+1),phi,u(:,t+1),options,dim,t+1);
            C =  dG_dX*iQy{t+1}*dG_dX';
            iRp = pinv(Rp);
            iSX = iRp + sigma*C;
            SigmaX{t+1} = pinv( iSX );
            muX(:,t+1) = mStar(:,t+1) + sigma.*SigmaX{t+1}*dG_dX*iQy{t+1}* (y(:,t+1)-gx(:,t+1));
        else
            muX(:,t+1) = mStar(:,t+1);
            SigmaX{t+1} = Rp;
        end
        % get predicted observation at the mode
        [gx(:,t+1),dG_dX] = VBA_evalFun('g',muX(:,t+1),phi,u(:,t+1),options,dim,t+1);
        suffStat.dy(:,t+1) = y(:,t+1) - gx(:,t+1);
        if options.sources.type == 0
            suffStat.vy(:,t+1) = diag( sigma.^-1.*pinv(iQy{t+1}) + dG_dX'*SigmaX{t+1}*dG_dX );
            suffStat.dy2 = suffStat.dy2 + suffStat.dy(:,t+1)'*iQy{t+1}*suffStat.dy(:,t+1);
        else
            suffStat.vy(:,t+1) = gx(:,t+1).*(1-gx(:,t+1));
            suffStat.logL = suffStat.logL + y(:,t+1)'*log(gx(:,t+1)) + (1-y(:,t+1))'*log(1-gx(:,t+1));
        end
        suffStat.dx(:,t+1) = muX(:,t+1) - fx;
        suffStat.dx2 = suffStat.dx2 + suffStat.dx(:,t+1)'*iQx{t+1}*suffStat.dx(:,t+1);
    else
        muX(:,t+1) = VBA_evalFun('f',muX(:,t),theta,u(:,t+1),options,dim);
    end
    if ~options.OnLine && isequal(mod(t,32),0) && options.verbose
        fprintf(1,repmat('\b',1,8))
        fprintf(1,'%6.2f %%',100*t/dim.n_t)
    end
end
if ~options.OnLine  && options.verbose
    fprintf(1,repmat('\b',1,8))
    fprintf(' OK.')
    fprintf('\n')
end
if flag >= 1
    suffStat.gx = gx;
end

author: MBB-team

repo: VBA-toolbox


text:
function [ n_data, x, fx ] = bessel_i1_spherical_values ( n_data )

%*****************************************************************************80
%
%% BESSEL_I1_SPHERICAL_VALUES returns some values of the Spherical Bessel function i1.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      Sqrt[Pi/(2*x)] * BesselI[3/2,x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    LC: QA47.A34,
%    ISBN: 0-486-61272-4.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 21;

  fx_vec = [ ...
    0.03336667857363341E+00, ...
    0.06693371456802954E+00, ...
    0.1354788933285401E+00, ...
    0.2072931911031093E+00, ...
    0.2841280857128948E+00, ...
    0.3678794411714423E+00, ...
    0.4606425870674146E+00, ...
    0.5647736480096238E+00, ...
    0.6829590627779635E+00, ...
    0.8182955028627777E+00, ...
    0.9743827435800610E+00, ...
    1.155432469636406E+00, ...
    1.366396525527973E+00, ...
    1.613118767572064E+00, ...
    1.902515460838681E+00, ...
    2.242790117769266E+00, ...
    2.643689828630357E+00, ...
    3.116811526884873E+00, ...
    3.675968313148932E+00, ...
    4.337627987747642E+00, ...
    5.121438384183637E+00 ];

  x_vec = [ ...
     0.1E+00, ...
     0.2E+00, ...
     0.4E+00, ...
     0.6E+00, ...
     0.8E+00, ...
     1.0E+00, ...
     1.2E+00, ...
     1.4E+00, ...
     1.6E+00, ...
     1.8E+00, ...
     2.0E+00, ...
     2.2E+00, ...
     2.4E+00, ...
     2.6E+00, ...
     2.8E+00, ...
     3.0E+00, ...
     3.2E+00, ...
     3.4E+00, ...
     3.6E+00, ...
     3.8E+00, ...
     4.0E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    fx = 0.0;
  else
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function test_bug2770

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY eeglab2fieldtrip

% the *.set file is actually a MATLAB file with the EEG structure in it
% but it require EEGLAB to read it together with the ftd (which has the binary data)

% I imported the data in EEGLAB with the GUI, and then saved it to a *.mat file
load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug2770/164_MIST_prac.mat'));

% the eeglab2fieldtrip function is maintained
ft_hastoolbox('eeglab', 1);

data = eeglab2fieldtrip(EEG, 'preprocessing', 'none');

cfg = [];
cfg.method = 'mtmfft';
cfg.taper = 'hanning';
cfg.foilim = [4 7];
cfg.output = 'pow';
freq = ft_freqanalysis(cfg,data); % this failed

cfg = [];
cfg.trials = 'all';
cfg.channel = 'all';
data1 = ft_selectdata(cfg, data);
assert(isfield(data1, 'trial')); % this failed

% this is actually at the core of the problem
assert(numel(data.label)==size(data.trial{1},1));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function cmap = revjetCmap(numGrays,numColors)
%
% cmap = revjetCmap(numGrays,numColors)
%
% Makes colormap array with:
%   gray scale - 1:numGrays
%   jet colors - numGrays+1:numGrays+numColors
%
% djh 1/98

if ~exist('numGrays','var')
    numGrays = 128;
end
if ~exist('numColors','var')
    numColors = 128;
end

cmap = zeros(numGrays+numColors,3);
cmap(1:numGrays+numColors,:) = [gray(numGrays);1-jet(numColors)];

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [y,deriv] = cllr_obj(w,T,weights,logit_prior)
% This is an MV2DF. See MV2DF_API_DEFINITION.readme.
%
% Weighted binary classifier cross-entropy objective, based on logarithmic
% cost function.
%
%  Differentiable inputs:
%   w: is vector of N detection scores (in log-likelihood-ratio format)
%
%  Fixed parameters:
%   T: is vector of N labels: 1 for target and -1 for non-target.
%   weights: is N-vector of objective function weights, one per trial.
%   logit_prior: is logit(prior), this controls the region of interest

if nargin==0
    test_this();
    return;
end

if isempty(w)
    y = @(w)cllr_obj(w,T,weights,logit_prior);
    return;
end

if isa(w,'function_handle')
    outer = cllr_obj([],T,weights,logit_prior);
    y = compose_mv(outer,w,[]);
    return;
end

w = w(:);
scores = w.';
arg = bsxfun(@plus,scores,logit_prior).*T;
neglogp1 = neglogsigmoid(arg); % 1*N           p1 = p(tar)
y = neglogp1*weights(:);

if nargout>1
    neglogp2 = neglogsigmoid(-arg); % 1*N      p2 = 1-p1 = p(non)
    deriv = @(dy) deriv_this(dy,-neglogp1(:),-neglogp2(:),T(:),weights(:));
end

function [g,hess,linear] = deriv_this(dy,logp1,logp2,T,weights)
g0 = -exp(logp2).*weights.*T;
g = dy*g0;
linear = false;
hess = @(d) hessianprod(d,dy,g0,logp1,logp2,weights);

function [h,Jv] = hessianprod(d,dy,g0,logp1,logp2,weights)

h = dy*(exp(logp1+logp2).*weights(:).*d(:));

if nargout>1
    Jv = d.'*g0;
end

function test_this()
N = 30;
T = [ones(1,N/3),-ones(1,N/3),zeros(1,N/3)];
W = randn(1,N);
weights = [rand(1,2*N/3),zeros(1,N/3)];
f = @(w) cllr_obj(w,T,weights,-2.23);
test_MV2DF(f,W(:));

author: nesl


text:
%% Predict Masked Tokens Using BERT
% This example shows how to predict masked tokens using a pretrained BERT
% model.
%
% BERT models are trained to perform various tasks. One of the tasks is
% known as masked language modeling which is the task of predicting tokens
% in text that have been replaced by a mask value.
%
% This example shows how to predict masked tokens for text data and
% calculate the token probabilities using a pretrained BERT model.

%% Load Pretrained BERT Model
% Load a pretrained BERT model using the |bert| function. The model
% consists of a tokenizer that encodes text as sequences of integers, and a
% structure of parameters.
mdl = bert

%%
% View the BERT model tokenizer. The tokenizer encodes text as sequences of
% integers and holds the details of padding, start, separator and mask
% tokens.
tokenizer = mdl.Tokenizer

%% Predict Masked Token
% Create a string containing a piece of text and replace a single word with
% the tokenizer mask token.
str = "Text Analytics Toolbox includes tools for preprocessing raw text from sources such as equipment logs, news feeds, surveys, operator reports, and social media.";
strMasked = replace(str,"sources",tokenizer.MaskToken)

%%
% Predict the masked token using the |predictMaskedToken| function. The
% function returns the original string with the mask tokens replaced.
sentencePred = predictMaskedToken(mdl,strMasked)

%% Calculate Prediction Scores
% To get the prediction scores for each word in the model vocabulary, you
% can evaluate the language model directly using the |bert.languageModel|
% function.

%%
% First, tokenize the input sentence with the BERT model tokenizer using
% the |tokenize| function. Note that the tokenizer may split single words
% and also prepends a [CLS] token and appends a [SEP] token to the input.
tokens = tokenize(tokenizer,str);
tokens{1}

%%
% Replace one of the tokens with the mask token.
idx = 16;
tokens{1}(idx) = tokenizer.MaskToken

%%
% Encode the tokens using the BERT model tokenizer using the |encodeTokens|
% function.
X = encodeTokens(tokenizer,tokens);
X{1}

%%
% To get the predictions scores from for the encoded tokens, evaluate the
% BERT language model directly using the |bert.languageModel| function. The
% language model output is a VocabularySize-by-SequenceLength array.
scores = bert.languageModel(X{1},mdl.Parameters);

%%
% View the tokens of the BERT model vocabulary corresponding to the top 10
% prediction scores for the mask token.
[~,idxTop] = maxk(extractdata(scores(:,idx)),10);
tbl = table;
tbl.Token = arrayfun(@(x) decode(tokenizer,x), idxTop);
tbl.Score = scores(idxTop,idx)

author: matlab-deep-learning

repo: transformer-models


text:
%SETDATA Dummy routine for setdata for doubles
%
%   OUT = SETDATA(A,B,PARS)
%
%DESCRIPTION
%This routine catches calls to SETDATA in case A is not a dataset but a
%double. In that case just B will be returned instead of a dataset with all
%annotations as defined by A but with data B.

function b = setdata(a,b,varargin)

return

author: marianux

repo: ecg-kit

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


text:
function y = sparbwslv(L,b)
% y = sparbwslv(L,b)
%
% SPARBWSLV Solves block sparse upper-triangular system.
%    y = sparbwslv(L,b) yields the same result as
%              y(L.perm,:) = L.L'\b
%    However, SPARBWSLV is faster than the built-in operator "\",
%    because it uses dense linear algebra and loop-unrolling on
%    supernodes.
%
%    Typical use, with X sparse m x m positive definite and b is m x n:
%            L = sparchol(symbchol(X),X);
%            L.d(L.dep) = inf;
%            y = sparbwslv(L,sparfwslv(L,b) ./ L.d);
%    Then y solves X*y=b.
%
% See also  symbchol, sparchol, sparfwslv, mldivide, mrdivide.

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

y = bwblkslv(L,b);
%y(L.perm,:) = L.L'\b;

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function erfx = r8_erf ( x )

%*****************************************************************************80
%
%% R8_ERF evaluates the error function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 May 2007
%
%  Author:
%
%    W J Cody,
%    Mathematics and Computer Science Division,
%    Argonne National Laboratory,
%    Argonne, Illinois, 60439.
%
%  Reference:
%
%    W J Cody,
%    "Rational Chebyshev approximations for the error function",
%    Mathematics of Computation,
%    1969, pages 631-638.
%
%  Parameters:
%
%    Input, real X, the argument of the error function.
%
%    Output, real ERFX, the value of the error function.
%
  a = [ ...
    3.16112374387056560E+00, ...
    1.13864154151050156E+02, ...
    3.77485237685302021E+02, ...
    3.20937758913846947E+03, ...
    1.85777706184603153E-01 ];
  b = [ ...
    2.36012909523441209E+01, ...
    2.44024637934444173E+02, ...
    1.28261652607737228E+03, ...
    2.84423683343917062E+03 ];
  c = [ ...
    5.64188496988670089E-01, ...
    8.88314979438837594E+00, ...
    6.61191906371416295E+01, ...
    2.98635138197400131E+02, ...
    8.81952221241769090E+02, ...
    1.71204761263407058E+03, ...
    2.05107837782607147E+03, ...
    1.23033935479799725E+03, ...
    2.15311535474403846E-08 ];
  d = [ ...
    1.57449261107098347E+01, ...
    1.17693950891312499E+02, ...
    5.37181101862009858E+02, ...
    1.62138957456669019E+03, ...
    3.29079923573345963E+03, ...
    4.36261909014324716E+03, ...
    3.43936767414372164E+03, ...
    1.23033935480374942E+03 ];
  p = [ ...
    3.05326634961232344E-01, ...
    3.60344899949804439E-01, ...
    1.25781726111229246E-01, ...
    1.60837851487422766E-02, ...
    6.58749161529837803E-04, ...
    1.63153871373020978E-02 ];
  q = [ ...
    2.56852019228982242E+00, ...
    1.87295284992346047E+00, ...
    5.27905102951428412E-01, ...
    6.05183413124413191E-02, ...
    2.33520497626869185E-03 ];
  sqrpi = 0.56418958354775628695E+00;
  thresh = 0.46875E+00;
  xbig = 26.543E+00;
  xsmall = 1.11E-16;
%
  xabs = abs ( x );
%
%  Evaluate ERF(X) for |X| <= 0.46875.
%
  if ( xabs <= thresh )

    if ( xsmall < xabs )
      xsq = xabs * xabs;
    else
      xsq = 0.0;
    end

    xnum = a(5) * xsq;
    xden = xsq;
    for i = 1 : 3
      xnum = ( xnum + a(i) ) * xsq;
      xden = ( xden + b(i) ) * xsq;
    end

    erfx = x * ( xnum + a(4) ) / ( xden + b(4) );
%
%  Evaluate ERFC(X) for 0.46875 <= |X| <= 4.0.
%
  elseif ( xabs <= 4.0 )

    xnum = c(9) * xabs;
    xden = xabs;
    for i = 1 : 7
      xnum = ( xnum + c(i) ) * xabs;
      xden = ( xden + d(i) ) * xabs;
    end

    erfx = ( xnum + c(8) ) / ( xden + d(8) );
    xsq = floor ( xabs * 16.0 ) / 16.0;
    del = ( xabs - xsq ) * ( xabs + xsq );
    erfx = exp ( - xsq * xsq ) * exp ( - del ) * erfx;

    erfx = ( 0.5 - erfx ) + 0.5;

    if ( x < 0.0 )
      erfx = -erfx;
    end
%
%  Evaluate ERFC(X) for 4.0 < |X|.
%
  else

    if ( xbig <= xabs )

      if ( 0.0 < x )
        erfx = 1.0;
      else
        erfx = -1.0;
      end

    else

      xsq = 1.0 / ( xabs * xabs );

      xnum = p(6) * xsq;
      xden = xsq;
      for i = 1 : 4
        xnum = ( xnum + p(i) ) * xsq;
        xden = ( xden + q(i) ) * xsq;
      end

      erfx = xsq * ( xnum + p(5) ) / ( xden + q(5) );
      erfx = ( sqrpi -  erfx ) / xabs;
      xsq = floor ( xabs * 16.0 ) / 16.0;
      del = ( xabs - xsq ) * ( xabs + xsq );
      erfx = exp ( - xsq * xsq ) * exp ( - del ) * erfx;

      erfx = ( 0.5 - erfx ) + 0.5;
      if ( x < 0.0 )
        erfx = -erfx;
      end

    end

  end

  return
end


text:
%IMESHGRID Domain matrices for image
%
% [U,V] = IMESHGRID(IM) are matrices that describe the domain of image IM (HxW)
% and are each HxW.  These matrices are used for the evaluation of functions
% over the image. The element U(R,C) = C and V(R,C) = R.
%
% [U,V] = IMESHGRID(W, H) as above but the domain is WxH.
%
% [U,V] = IMESHGRID(S) as above but the domain is described by S which can
% be a scalar SxS or a 2-vector S=[W, H].
%
% See also MESHGRID, INTERP2.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB 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.
%
% MVTB 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 MVTB.  If not, see <http://www.gnu.org/licenses/>.

function [U,V] = imeshgrid(a1, a2)

    if nargin == 1
        if length(a1) == 1
            % we specified a size for a square output image
            [U,V] = meshgrid(1:a1, 1:a1);
        elseif numel(a1) == 2
            % we specified a size for a rectangular output image (w,h)
            [U,V] = meshgrid(1:a1(1), 1:a1(2));
        elseif ndims(a1) >= 2
            [U,V] = meshgrid(1:numcols(a1), 1:numrows(a1));
        else
            error('incorrect argument');
        end
    elseif nargin == 2
        [U,V] = meshgrid(1:a1, 1:a2);
    end

author: petercorke

repo: machinevision-toolbox-matlab