text:
function x=v_choosenk(n,k)
%V_CHOOSENK All choices of K elements taken from 1:N [X]=(N,K)
% The output X is a matrix of size (N!/(K!*(N-K)!),K) where each row
% contains a choice of K elements taken from 1:N without duplications.
% The rows of X are in lexically sorted order.
%
% To choose from the elements of an arbitrary vector V use
% V(CHOOSENK(LENGTH(V),K)).

% CHOOSENK(N,K) is the same as the MATLAB5 function NCHOOSEK(1:N,K) but is
% much faster for large N and most values of K.

%   Copyright (c) 1998 Mike Brookes,  mike.brookes@ic.ac.uk
%      Version: $Id: v_choosenk.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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

kk=min(k,n-k);
if kk<2
   if kk<1
      if k==n
         x=1:n;
      else
         x=[];
      end
   else
      if k==1
         x=(1:n)';
      else
         x=1:n;
         x=reshape(x(ones(n-1,1),:),n,n-1);
      end
   end
else
   n1=n+1;
   m=prod(n1-kk:n)/prod(1:kk);
   x=zeros(m,k);
   f=n1-k;
   x(1:f,k)=(k:n)';
   for a=k-1:-1:1
      d=f;
      h=f;
      x(1:f,a)=a;
      for b=a+1:a+n-k
         d=d*(n1+a-b-k)/(n1-b);
         e=f+1;
         f=e+d-1;
         x(e:f,a)=b;
         x(e:f,a+1:k)=x(h-d+1:h,a+1:k);
      end
   end
end

author: ImperialCollegeLondon

repo: sap-voicebox


x = [ 1 2 3 4 5 6];
y = [ 2 6 8 7 8 5];
stairs(x,y);
title('\bfExample of a Stair Plot');
xlabel('\bf\itx');
ylabel('\bf\ity');
axis([0 7 0 10]);


text:
function sp=mtimes(sp1,sp2)
% sympoly/mtimes: Matrix multiplication of sympoly objects or scalars
% usage: sp=sp1*sp2;
%
% arguments:
%  sp,sp1,sp2   - sympoly objects or numeric scalars or a numeric array

% are they of compatible sizes?
s1 = size(sp1);
s2 = size(sp2);
if (length(s1)>2) || (length(s2)>2)
  error 'Matrix multiplication is only defined for vectors and simple arrays.'
end

if (s1(2) == s2(1)) || (numel(sp1)==1) || (numel(sp2)==1)
  % they are compatible

  if (numel(sp1) == 1) && (numel(sp2) == 1)
    % both are scalars. Just use .*
    sp = sp1.*sp2;

  elseif (numel(sp1) == 1) && (numel(sp2) > 1)
    % sp1 is a scalar, but not sp2
    sp = sympoly(sp2);
    for i = 1:numel(sp2)
      sp(i) = sp1.*sp2(i);
    end

  elseif (numel(sp1) > 1) && (numel(sp2) == 1)
    % sp2 is a scalar, but not sp1
    sp = sympoly(sp1);
    for i = 1:numel(sp1)
      sp(i) = sp1(i).*sp2;
    end

  elseif (numel(sp1) > 1) && (numel(sp2) > 1)
    % both are arrays.
    sp = sympoly(zeros(s1(1),s2(2)));

    for i = 1:s1(1)
      for j = 1:s2(2)
        for k = 1:s1(2)
          sp(i,j) = sp(i,j) + sp1(i,k).*sp2(k,j);
        end
      end
    end

  end

else
  error 'sp1 and sp2 are of incompatible sizes for .* operation.'

end

author: Sable

repo: mcbench-benchmarks


text:
function dc = diffusivity_2d_elman ( a, cl, dc0, m_1d, omega, n1, n2, x, y )

%*****************************************************************************80
%
%% DIFFUSIVITY_2D_ELMAN evaluates a 2D stochastic diffusivity function.
%
%  Discussion:
%
%    The 2D diffusion equation has the form
%
%      - Del ( DC(X,Y) Del U(X,Y) ) = F(X,Y)
%
%    where DC(X,Y) is a function called the diffusivity.
%
%    In the stochastic version of the problem, the diffusivity function
%    includes the influence of stochastic parameters:
%
%      - Del ( DC(X,Y;OMEGA) Del U(X,Y;OMEGA) ) = F(X,Y).
%
%    In this function, the domain is assumed to be the square [-A,+A]x[-A,+A].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 December 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Howard Elman, Darran Furnaval,
%    Solving the stochastic steady-state diffusion problem using multigrid,
%    IMA Journal on Numerical Analysis,
%    Volume 27, Number 4, 2007, pages 675-688.
%
%    Roger Ghanem, Pol Spanos,
%    Stochastic Finite Elements: A Spectral Approach,
%    Revised Edition,
%    Dover, 2003,
%    ISBN: 0486428184,
%    LC: TA347.F5.G56.
%
%  Parameters:
%
%    Input, real A, the "radius" of the square region.  The region
%    is assumed to be [-A,+A]x[-A,+A].
%    0 < A.
%
%    Input, real CL, the correlation length.
%    0 < CL.
%
%    Input, real DC0, the constant term in the expansion of the
%    diffusion coefficient.  Take DC0 = 10.
%
%    Input, integr M_1D, the first and second dimensions of the stochastic
%    parameter array.
%
%    Input, real OMEGA(M_1D*M_1D), the stochastic parameters.
%
%    Input, integer N1, N2, the dimensions of the X and Y arrays.
%
%    Input, real X(N1,N2), Y(N1,N2), the points where the diffusion coefficient is to
%    be evaluated.
%
%    Output, real DC(N1,N2), the value of the diffusion coefficient at X.
%

%
%  Compute THETA.
%
  theta_1d = theta_solve ( a, cl, m_1d );
%
% m = m_1d * m_1d;
% theta_1d = theta_solve ( a, cl, m );
%
%  Compute LAMBDA_1D.
%
  lambda_1d(1:m_1d) = 2.0 * cl ./ ( 1.0 + cl^2 * theta_1d(1:m_1d).^2 );
%
%  Compute C_1DX(1:M1D)  and C_1DY(1:M1D) at (X,Y).
%
  c_1dx = zeros ( m_1d, n1, n2 );
  c_1dy = zeros ( m_1d, n1, n2 );

  k = 0;

  while ( 1 )

    if ( m_1d <= k )
      break
    end

    k = k + 1;
    c_1dx(k,1:n1,1:n2) = cos ( theta_1d(k) * a * x(1:n1,1:n2) ) ...
      / sqrt ( a + sin ( 2.0 * theta_1d(k) * a ) / ( 2.0 * theta_1d(k) ) );
    c_1dy(k,1:n1,1:n2) = cos ( theta_1d(k) * a * y(1:n1,1:n2) ) ...
      / sqrt ( a + sin ( 2.0 * theta_1d(k) * a ) / ( 2.0 * theta_1d(k) ) );

    if ( m_1d <= k )
      break
    end

    k = k + 1;
    c_1dx(k,1:n1,1:n2) = sin ( theta_1d(k) * a * x(1:n1,1:n2) ) ...
      / sqrt ( a - sin ( 2.0 * theta_1d(k) * a ) / ( 2.0 * theta_1d(k) ) );
    c_1dy(k,1:n1,1:n2) = sin ( theta_1d(k) * a * y(1:n1,1:n2) ) ...
      / sqrt ( a - sin ( 2.0 * theta_1d(k) * a ) / ( 2.0 * theta_1d(k) ) );

  end
%
%  Evaluate the diffusion coefficient DC at (X,Y).
%  This nonsense of fussy array shapes really frustrates me!
%
  k = 0;
  dc = dc0 * ones ( 1, n1, n2 );
  for j = 1 : m_1d
    for i = 1 : m_1d
      k = k + 1;
      dc(1,1:n1,1:n2) = dc(1,1:n1,1:n2) + sqrt ( lambda_1d(i) * lambda_1d(j) ) ...
        * c_1dx(i,1:n1,1:n2) .* c_1dy(j,1:n1,1:n2) * omega(k);
    end
  end

  dc = reshape ( dc, n1, n2 );

  return
end
function theta = theta_solve ( a, cl, m )

%*****************************************************************************80
%
%% THETA_SOLVE solves a pair of transcendental equations.
%
%  Discussion:
%
%    The vector THETA returned by this function is needed in order to define
%    the terms in a Karhunen-Loeve expansion of a diffusion coefficient.
%
%    The two equations are:
%
%      1/CL - THETA * TAN ( A * THETA ) = 0
%      THETA - 1/CL * TAN ( A * THETA ) = 0
%
%    A and CL are taken to be positive.  Over each open interval
%
%      ( n - 1/2 pi, n + 1/2 pi ) / A, for N = 0, 1, ...
%
%    the function TAN ( A * THETA ) monotonically rises from -oo to +00;
%    therefore, it can be shown that there is one root of each equation
%    in every interval of this form.  Moreover, because of the positivity
%    of A and CL, we can restrict our search to the interval
%
%      [ n pi, n + 1/2 pi ) / A, for N = 0, 1, ...
%
%    This function computes K such roots, starting in the first interval,
%    finding those two roots, moving to the next interval, and so on, until
%    the requested number of roots have been found.  Odd index roots will
%    correspond to the first equation, and even index roots to the second.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 July 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Howard Elman, Darran Furnival,
%    Solving the Stochastic Steady-State Diffusion Problem Using Multigrid,
%    University of Maryland Department of Computer Science,
%    Technical Report TR-4786.
%
%  Parameters:
%
%    Input, real A, the "radius" of the domain, D = (-A,A)x(-A,A).
%    0 < A.
%
%    Input, real CL, the correlation length.
%    0 < CL.
%
%    Input, integer M, the number of values to compute.
%
%    Output, real THETA(M), the values of Theta.
%
  k = 0;
  theta = zeros(m,1);
%
%  [ XA_INIT, XB_INIT] = [ n * pi, n+1/2 pi ] / a, n = 0, 1, 2, ...
%
  xa_init = 0.0;
  xb_init = ( pi / 2 ) / a;

  while ( 1 )
%
%  Seek root of equation 1 in interval.
%
    if ( m <= k )
      break
    end

    k = k + 1;
    xa = xa_init;
    fa = 1.0 / cl - xa * tan ( a * xa );
    ftol = eps * ( abs ( fa ) + 1.0 );
    xb = xb_init;
    fb = - fa;
    fc = fa;
    bmatol = 100.0 * eps * ( abs ( xa ) + abs ( xb ) );

    while ( bmatol < xb - xa )

      xc = ( xa + xb ) / 2;
      fc = 1.0 / cl - xc * tan ( a * xc );

      if ( abs ( fc ) <= ftol )
        break
      elseif ( 0 < fc )
        xa = xc;
      else
        xb = xc;
      end

    end

    theta(k) = xc;
%
%  Seek root of equation 2 in interval.
%
    if ( m <= k )
      break
    end

    k = k + 1;
%
%  In the first interval, we need to skip the zero root of equation 2.
%
    if ( k == 2 )

      k = k - 1;

    else

      xa = xa_init;
      fa = xa - tan ( a * xa ) / cl;
      ftol = eps * ( abs ( fa ) + 1.0 );
      xb = xb_init;
      fb = - fa;

      while ( bmatol < xb - xa )

        xc = ( xa + xb ) / 2;
        fc = xc - tan ( a * xc ) / cl;

        if ( abs ( fc ) <= ftol )
          break
        elseif ( 0 < fc )
          xa = xc;
        else
          xb = xc;
        end

      end

      theta(k) = xc;

    end
%
%  Advance the interval.
%
    xa_init = xa_init + pi / a;
    xb_init = xb_init + pi / a;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function varargout = model(x,parameters,nvp)
% model   A BERT model forward pass.
%
%   Z = model(X,parameters) performs inference with a BERT model on the
%   input X. X is a 1-by-numInputTokens-by-numObs array of encoded tokens.
%   The return is an array Z of size
%   (NumHeads*HeadSize)-by-numInputTokens-by-numObs. Each
%   Z(:,i,j) corresponds to the BERT embedding of input token X(1,i,j).
%
%   Z = model(X,parameters,'PARAM1', VAL1, 'PARAM2', VAL2, ...) specifies
%   the optional parameter name/value pairs:
%
%      'InputMask'             - A logical mask with the same size as X.
%                                The mask should be false at indices
%                                (i,j,k) for which X(i,j,k) corresponds to
%                                padding, and true elsewhere. The default
%                                is empty, for which the padding is
%                                inferred by the entries of X that match
%                                the PaddingCode name-value pair.
%
%      'DropoutProb'           - The probability of dropout for the output
%                                activation. It is standard to set this to a
%                                non-zero value during training, for example
%                                0.1. The default is 0.
%
%      'AttentionDropoutProb'  - The probability of dropout used in the
%                                attention layer. It is standard to set
%                                this to a non-zero value during training,
%                                for example 0.1. The default is 0.
%
%      'Outputs'               - Specify the indices of the layers to
%                                return outputs from as a vector of
%                                positive integers, or 'last' to specify
%                                the final encoder layer only. The default
%                                is 'last'.
%
%      'SeparatorCode'         - The positive integer corresponding to the
%                                separator token. The default is 103, as
%                                specified in the default BERT vocab.txt.
%
%      'PaddingCode'           - The positive integer corresponding to the
%                                padding token. The default is 1, as
%                                specified in the default BERT vocab.txt.
%
% References:
% [1] https://arxiv.org/abs/1810.04805
% [2] https://github.com/google-research/bert/

% Copyright 2021 The MathWorks, Inc.

arguments
    x dlarray {mustBeNumericDlarray,mustBeNonempty}
    parameters {mustBeA(parameters,'struct')}
    nvp.InputMask {mustBeNumericOrLogical} = logical.empty()
    nvp.DropoutProb (1,1) {mustBeNonnegative,mustBeLessThanOrEqual(nvp.DropoutProb,1),mustBeNumeric} = 0
    nvp.AttentionDropoutProb (1,1) {mustBeNonnegative,mustBeLessThanOrEqual(nvp.AttentionDropoutProb,1),mustBeNumeric} = 0
    nvp.Outputs {mustBePositive,mustBeLessThanOrEqualNumLayers(nvp.Outputs,parameters),mustBeInteger,mustBeNumeric} = parameters.Hyperparameters.NumLayers
    nvp.PaddingCode (1,1) {mustBePositive,mustBeInteger,mustBeNumeric} = 1
    nvp.SeparatorCode (1,1) {mustBePositive,mustBeInteger,mustBeNumeric} = 103
end

nvp.Outputs = nvp.Outputs(:);
nargoutchk(0,numel(nvp.Outputs));

maxLayer = max(nvp.Outputs);

w = parameters.Weights;
hyperparameters = parameters.Hyperparameters;

% Identify padding tokens - if InputMask is set just use it.
if isempty(nvp.InputMask)
    inputMask = x~=nvp.PaddingCode;
else
    assert(isequal(size(nvp.InputMask),size(x)),"bert:model:InvalidMaskSize","Expected InputMask to have same size as input X.");
    inputMask = logical(nvp.InputMask);
end

% Assuming CTB format of x.
xsz = size(x);
seq_len = xsz(2);

% Apply embeddings
types = dlarray(bert.internal.inferTypeID(x,nvp.SeparatorCode));
positions = dlarray(1:seq_len);
z = bert.layer.embedding(x,types,positions,w.embeddings,nvp.DropoutProb);

% Transformer layers
num_layers = min(hyperparameters.NumLayers,maxLayer);
varargout = cell(numel(nvp.Outputs),1);
for i = 1:num_layers
    z = bert.layer.block(z,w.encoder_layers.("layer_"+i),hyperparameters,'InputMask',inputMask);
    toAssign = nvp.Outputs==i;
    varargout(toAssign) = repelem({z},sum(toAssign));
end
end

function mustBeLessThanOrEqualNumLayers(x,params)
mustBeLessThanOrEqual(x,params.Hyperparameters.NumLayers);
end

function mustBeALogicalOrDlarrayLogical(val)
if isa(val,'dlarray')
    val = extractdata(val);
end
mustBeA(val,'logical');
end

function mustBeNumericDlarray(val)
mustBeA(val,'dlarray');
mustBeNumeric(extractdata(val));
end

author: matlab-deep-learning

repo: transformer-models


text:
function [ n_data, n, c ] = moebius_values ( n_data )

%*****************************************************************************80
%
%% MOEBIUS_VALUES returns some values of the Moebius function.
%
%  Discussion:
%
%    MU(N) is defined as follows:
%
%      MU(N) = 1 if N = 1;
%              0 if N is divisible by the square of a prime;
%              (-1)^K, if N is the product of K distinct primes.
%
%    In Mathematica, the function can be evaluated by:
%
%      MoebiusMu[n]
%
%  First values:
%
%     N  MU(N)
%
%     1    1
%     2   -1
%     3   -1
%     4    0
%     5   -1
%     6    1
%     7   -1
%     8    0
%     9    0
%    10    1
%    11   -1
%    12    0
%    13   -1
%    14    1
%    15    1
%    16    0
%    17   -1
%    18    0
%    19   -1
%    20    0
%
%    As special cases, MU(N) is -1 if N is a prime, and MU(N) is 0
%    if N is a square, cube, etc.
%
%  Formula:
%
%    The Moebius function is related to Euler's totient function:
%
%      PHI(N) = Sum ( D divides N ) MU(D) * ( N / D ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    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, integer N, the argument of the Moebius function.
%
%    Output, integer C, the value of the Moebius function.
%
  n_max = 20;

  c_vec = [ ...
      1,  -1,  -1,   0,  -1,   1,  -1,   0,   0,   1, ...
     -1,   0,  -1,   1,   1,   0,  -1,   0,  -1,   0 ];

  n_vec = [ ...
      1,   2,   3,   4,   5,   6,   7,   8,   9,  10, ...
     11,  12,  13,  14,  15,  16,  17,  18,  19,  20 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    c = 0;
  else
    n = n_vec(n_data);
    c = c_vec(n_data);
  end

  return
end


text:
function b = r8lt_vxm ( m, n, a, x )

%*****************************************************************************80
%
%% R8LT_VXM multiplies a vector by a R8LT matrix.
%
%  Discussion:
%
%    The R8LT storage format is used for an M by N lower triangular matrix,
%    and sets aside storage even for the entries that must be zero.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows of the matrix.
%    M must be positive.
%
%    Input, integer N, the number of columns of the matrix.
%    N must be positive.
%
%    Input, real A(M,N), the R8LT matrix.
%
%    Input, real X(M), the vector to be multiplied by A.
%
%    Output, real B(N), the product A * x.
%
  for i = 1 : n
    b(i) = x(i:m) * a(i:m,i);
  end

  return
end


text:
function PSNRdb = PSNR(x, y)

err = x - y;
err = err(:);
PSNRdb = 20 * log10(256/sqrt(mean(err .^2)));

author: thomaskuestner


text:
function []=TVEmavardisp(beta_median,beta_std,beta_lbound,beta_ubound,theta_median,psi_std,psi_lbound,psi_ubound,sigma_median,X,Z,Y,n,m,p,k1,k3,q1,q2,T,ar,lambda1,lambda2,lambda3,lambda4,lambda5,IRFt,beta_gibbs,endo,exo,psi_gibbs,startdate,enddate,stringdates1,decimaldates1,datapath,data_endo,TVEH,indH,pref)

% function []=mavardisp(beta_median,beta_std,beta_lbound,beta_ubound,psi_median,psi_std,psi_lbound,psi_ubound,sigma_median,X,Z,Y,n,m,p,k1,k3,q1,q2,T,ar,lambda1,lambda2,lambda3,lambda4,lambda5,IRFt,beta_gibbs,endo,exo,psi_gibbs,startdate,enddate,stringdates1,decimaldates1,datapath)
% displays estimation results for the MABVAR model on Matlab prompt, creates a copy of these results on the text file results.txt
% inputs:  - vector 'beta_median': median value of the posterior distribution of beta
%          - vector 'beta_std': standard deviation of the posterior distribution of beta
%          - vector 'beta_lbound': lower bound of the credibility interval of beta
%          - vector 'beta_ubound': upper bound of the credibility interval of beta
%          - vector 'psi_median': median value of the posterior distribution of psi
%          - vector 'psi_std': standard deviation of the posterior distribution of psi
%          - vector 'psi_lbound': lower bound of the credibility interval of psi
%          - vector 'psi_ubound': upper bound of the credibility interval of psi
%          - vector 'sigma_median': median value of the posterior distribution of sigma (vectorised)
%          - matrix 'X': matrix of regressors for the VAR model (defined in 3.5.10)
%          - matrix 'Z': matrix of exogenous regressors for the MABVAR model (defined in 3.5.10)
%          - matrix 'Y': matrix of regressands for the VAR model (defined in 3.5.10)
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'm': number of exogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'p': number of lags included in the model (defined p 7 of technical guide)
%          - integer 'k1': number of endogenous coefficients to estimate for each equation in the MABVAR model (defined p 77 of technical guide)
%          - integer 'k3': number of exogenous coefficients to estimate for each equation in the reformulated MABVAR model (defined p 77 of technical guide)
%          - integer 'q1': total number of endogenous coefficients to estimate in the MABVAR model (defined p 77 of technical guide)
%          - integer 'q2': total number of exogenous coefficients to estimate in the MABVAR model (defined p 77 of technical guide)
%          - integer 'T': number of sample time periods (defined p 7 of technical guide)
%          - scalar 'lambda1': overall tightness hyperparameter (defined p 16 of technical guide)
%          - scalar 'lambda2': cross-variable weighting hyperparameter(defined p 16 of technical guide)
%          - scalar 'lambda3': lag decay hyperparameter (defined p 16 of technical guide)
%          - scalar 'lambda4': exogenous variable tightness hyperparameter (defined p 17 of technical guide)
%          - scalar 'lambda5': block exogeneity shrinkage hyperparameter (defined p 32 of technical guide)
%          - integer 'IRFt': determines which type of structural decomposition to apply (none, Choleski, triangular factorization)
%          - matrix 'beta_gibbs': record of the gibbs sampler draws for the beta vector
%          - cell 'endo': list of endogenous variables of the model
%          - matrix 'psi_gibbs': record of the gibbs sampler draws for the psi vector
%          - cell 'exo': list of exogenous variables of the model
%          - string 'startdate': start date of the sample
%          - string 'enddate': end date of the sample
%          - cell 'stringdates1': date strings for the sample period
%          - vector 'decimaldates1': dates converted into decimal values, for the sample period
%          - string 'datapath': user-supplied path to excel data spreadsheet
% outputs: none

% before displaying and saving the results, start estimating the evaluation measures for the model

% obtain first a point estimate betatilde of the VAR coefficients related to endogenous variables
% this is simply beta_median, which is both the mean and the median of the normal distribution
betatilde=beta_median;
Btilde=reshape(betatilde,k1,n);
% generate U, using Btilde, from (3.5.15)
% U=eye(n*m);
%    for jj=1:p
%    U=[U;kron(eye(m),Btilde((jj-1)*n+1:jj*n,:)')];
%    end

% obtain then a point estimate psitilde of the VAR coefficients related to exogenous variables
% this is simply psi_median, which is both the mean and the median of the normal distribution
theta=theta_median;

% % combine the two to obtain vec(DELTA'), using (3.5.14)
% vecdeltap=U*psitilde;
% % recover delta
% deltap=reshape(vecdeltap,n,k3);
% DELTA=deltap';

for it=1:T+p
    eq(it,:)=(squeeze(TVEH(:,:,indH(it),it))*theta)'; % compute the equilibrium values given theta
end

temp2=data_endo-eq;
temp3=bear.lagx(temp2,p);
Yhat=temp3(:,1:n);
Xhat=temp3(:,n+1:end);

% then produce the corresponding residuals, using (3.5.9)
EPStilde=Yhat-Xhat*Btilde;
Ytilde=Y-EPStilde;

% check first whether the model is stationary, using (1.9.1)
[stationary eigmodulus]=mabear.checkstable(Btilde,n,p);

% Compute then the sum of squared residuals
% compute first the RSS matrix, defined in (1.9.5)
RSS=EPStilde'*EPStilde;
% retain only the diagonal elements to get the vector of RSSi values
rss=diag(RSS);

% Go on calculating R2
% generate Mbar
Mbar=eye(T)-ones(T,T)/T;
% then compute the TSS matrix, defined in (1.9.8)
TSS=Y'*Mbar*Y;
% generate the R2 matrix in (1.9.9)
R2=eye(n)-RSS./TSS;
% retain only the diagonal elements to get the vector of R2 values
r2=diag(R2);

% then calculate the adjusted R2, using (1.9.11)
R2bar=eye(n)-((T-1)/(T-(k1+m)))*(eye(n)-R2);
% retain only the diagonal elements to get the vector of R2bar values
r2bar=diag(R2bar);

% now start displaying and saving the results

% preliminary task: create and open the txt file used to save the results

% filelocation=datapath;
% filelocation=[filelocation '\results.txt'];
% fid=fopen(filelocation,'wt');

filelocation=fullfile(pref.results_path, [pref.results_sub '.txt']);
fid=fopen(filelocation,'wt');

% print toolbox header

fprintf('%s\n','');
fprintf(fid,'%s\n','');

fprintf('%s\n','%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
fprintf(fid,'%s\n','%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
fprintf('%s\n','%                                                                                %%');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    BAYESIAN ESTIMATION, ANALYSIS AND REGRESSION (BEAR) TOOLBOX                                         %');
fprintf(fid,'%s\n','%    BAYESIAN ESTIMATION, ANALYSIS AND REGRESSION (BEAR) TOOLBOX                                     %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    This statistical package has been developed by the external developments division of the ECB.       %');
fprintf(fid,'%s\n','%    This statistical package has been developed by the external developments division of the ECB.   %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    Authors:                                                                                %');
fprintf(fid,'%s\n','%    Authors:                                                                                %');
fprintf('%s\n','%    Romain Legrand  (Romain Legrand <b00148883@essec.edu>)                                              %');
fprintf(fid,'%s\n','%    Romain Legrand  (Romain Legrand <b00148883@essec.edu>)                                          %');
fprintf('%s\n','%    Alistair Dieppe (adieppe@worldbank.org)                                                             %');
fprintf(fid,'%s\n','%    Alistair Dieppe (adieppe@worldbank.org)                                                         %');
fprintf('%s\n','%    Björn van Roye  (Bjorn.van_Roye@ecb.europa.eu)                                                      %');
fprintf(fid,'%s\n','%    Björn van Roye  (Bjorn.van_Roye@ecb.europa.eu)                                                  %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    Version 4.4                                                                                %');
fprintf(fid,'%s\n','%    Version 4.4                                                                                %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    The authors are grateful to Paolo Bonomolo, Marta Banbura, Martin Bruns, Fabio Canova,              %');
fprintf(fid,'%s\n','%    The authors are grateful to Paolo Bonomolo, Marta Banbura, Martin Bruns, Fabio Canova,          %');
fprintf('%s\n','%    Matteo Ciccarelli, Marek Jarocinski, Niccolo Battistini, Gabriel Bobeica                            %');
fprintf(fid,'%s\n','%    Matteo Ciccarelli, Marek Jarocinski, Niccolo Battistini, Gabriel Bobeica                        %');
fprintf('%s\n','%    Michele Lenza, Chiara Osbat, Mirela Miescu, Gary Koop, Giorgio Primiceri                            %');
fprintf(fid,'%s\n','%    Michele Lenza, Chiara Osbat, Mirela Miescu, Gary Koop, Giorgio Primiceri,                       %');
fprintf('%s\n','%    Michal Rubaszek, Barbara Rossi, Ben Schumann, Peter Welz, Hugo Vega de la Cruz and Francesca Loria. %');
fprintf(fid,'%s\n','%  Michal Rubaszek, Barbara Rossi, Ben Schumann, Peter Welz, Hugo Vega de la Cruz and Francesca Loria%');
fprintf('%s\n','%    valuable input and advice which contributed to improve the quality of this work.                    %');
fprintf(fid,'%s\n','%  valuable input and advice which contributed to improve the quality of this work.                  %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%   These programmes are the responsibilities of the authors and not of the ECB and the Worldbank.       %');
fprintf(fid,'%s\n','%   These programmes are the responsibilities of the authors and not of the ECB and the Worldbank.   %');
fprintf('%s\n','%   Errors and ommissions remain those of the authors.                                                   %');
fprintf(fid,'%s\n','%   Errors and ommissions remain those of the authors.                                               %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%    Please do not use or quote this work without permission.                                            %');
fprintf(fid,'%s\n','%    Please do not use or quote this work without permission.                                        %');
fprintf('%s\n','%                                                                                %');
fprintf(fid,'%s\n','%                                                                                %');
fprintf('%s\n','%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');
fprintf(fid,'%s\n','%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%');

% print then estimation results

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

toolboxinfo='BEAR toolbox estimates';
fprintf('%s\n',toolboxinfo);
fprintf(fid,'%s\n',toolboxinfo);

time=clock;
datestring=datestr(time);
dateinfo=['Date: ' datestring(1,1:11) '   Time: ' datestring(1,13:17)];
fprintf('%s\n',dateinfo);
fprintf(fid,'%s\n',dateinfo);

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

VARtypeinfo='Mean-adjusted BVAR';
fprintf('%s\n',VARtypeinfo);
fprintf(fid,'%s\n',VARtypeinfo);

if IRFt==1
SVARinfo='structural decomposition: none';
elseif IRFt==2
SVARinfo='structural decomposition: choleski factorisation';
elseif IRFt==3
SVARinfo='structural decomposition: triangular factorisation';
elseif IRFt==4
SVARinfo='structural decomposition: sign restrictions';
end
fprintf('%s\n',SVARinfo);
fprintf(fid,'%s\n',SVARinfo);

temp='endogenous variables: ';
for ii=1:n
temp=[temp ' ' endo{ii,1} ' '];
end
endoinfo=temp;
fprintf('%s\n',endoinfo);
fprintf(fid,'%s\n',endoinfo);

temp='exogenous variables: constant';
% if m>1
%    for ii=1:m-1
%    temp=[temp ' ' exo{ii,1} ' '];
%    end
% end
exoinfo=temp;
fprintf('%s\n',exoinfo);
fprintf(fid,'%s\n',exoinfo);

sampledateinfo=['estimation sample: ' startdate '-' enddate];
fprintf('%s\n',sampledateinfo);
fprintf(fid,'%s\n',sampledateinfo);

samplelengthinfo=['sample size (omitting initial conditions): ' num2str(T)];
fprintf('%s\n',samplelengthinfo);
fprintf(fid,'%s\n',samplelengthinfo);

laginfo=['number of lags included in regression: ' num2str(p)];
fprintf('%s\n',laginfo);
fprintf(fid,'%s\n',laginfo);

hyperparam1='hyperparameters:';
fprintf('%s\n',hyperparam1);
fprintf(fid,'%s\n',hyperparam1);

%hyperparam2=['autoregressive coefficient (ar):          ' num2str(ar)];
%fprintf('%s\n',hyperparam2);
%fprintf(fid,'%s\n',hyperparam2);

hyperparam3=['overall tightness (lambda1):              ' num2str(lambda1)];
fprintf('%s\n',hyperparam3);
fprintf(fid,'%s\n',hyperparam3);

hyperparam4=['cross-variable weighting (lambda2):       ' num2str(lambda2)];
fprintf('%s\n',hyperparam4);
fprintf(fid,'%s\n',hyperparam4);

hyperparam5=['lag decay (lambda3):                      ' num2str(lambda3)];
fprintf('%s\n',hyperparam5);
fprintf(fid,'%s\n',hyperparam5);

%hyperparam6=['exogenous variable tightness (lambda4):   ' num2str(lambda4)];
%fprintf('%s\n',hyperparam6);
%fprintf(fid,'%s\n',hyperparam6);

hyperparam7=['block exogeneity shrinkage (lambda5):     ' num2str(lambda5)];
fprintf('%s\n',hyperparam7);
fprintf(fid,'%s\n',hyperparam7);

% display coefficient estimates

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

coeffinfo=['VAR coefficients (beta): posterior estimates'];
fprintf('%s\n',coeffinfo);
fprintf(fid,'%s\n',coeffinfo);

for ii=1:n

fprintf('%s\n','');
fprintf(fid,'%s\n','');
if ii~=1
fprintf('%s\n','');
fprintf(fid,'%s\n','');
end

endoinfo=['Endogenous: ' endo{ii,1}];
fprintf('%s\n',endoinfo);
fprintf(fid,'%s\n',endoinfo);

coeffheader=fprintf('%25s %15s %15s %15s %15s\n','','Median','St.dev','Low.bound','Upp.bound');
coeffheader=fprintf(fid,'%25s %15s %15s %15s %15s\n','','Median','St.dev','Low.bound','Upp.bound');

% handle the endogenous
   for jj=1:n
      for kk=1:p
      values=[beta_median((ii-1)*k1+n*(kk-1)+jj,1) beta_std((ii-1)*k1+n*(kk-1)+jj,1) beta_lbound((ii-1)*k1+n*(kk-1)+jj,1) beta_ubound((ii-1)*k1+n*(kk-1)+jj,1)];
      fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n',strcat(endo{jj,1},'(-',int2str(kk),')'),values);
      fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n',strcat(endo{jj,1},'(-',int2str(kk),')'),values);
      end
   end

% handle the exogenous
% display the results related to the constant
values=[theta_median(ii,1) psi_std(ii,1) psi_lbound(ii,1) psi_ubound(ii,1)];
fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n','Constant',values);
fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n','Constant',values);
   % if there is no other exogenous, stop here
   if m==1
   % if there are other exogenous, display their results
   elseif m>1
      for jj=1:m-1
      %values=[theta_median(n*jj+ii,1) psi_std(n*jj+ii,1)
      %psi_lbound(n*jj+ii,1) psi_ubound(n*jj+ii,1)]; BVR, commeted out as
      %it was unused!
      %fprintf('%25s %15.3f %15.3f %15.3f %15.3f\n',exo{jj,1},values);
     %fprintf(fid,'%25s %15.3f %15.3f %15.3f %15.3f\n',exo{jj,1},values);
     end
   end

fprintf('%s\n','');
fprintf(fid,'%s\n','');

% display evaluation measures
rssinfo=['Sum of squared residuals: ' num2str(rss(ii,1),'%.2f')];
fprintf('%s\n',rssinfo);
fprintf(fid,'%s\n',rssinfo);

r2info=['R-squared: ' num2str(r2(ii,1),'%.3f')];
fprintf('%s\n',r2info);
fprintf(fid,'%s\n',r2info);

adjr2info=['adj. R-squared: ' num2str(r2bar(ii,1),'%.3f')];
fprintf('%s\n',adjr2info);
fprintf(fid,'%s\n',adjr2info);

end

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% display VAR stability results
eigmodulus=reshape(eigmodulus,p,n);
stabilityinfo1=['Roots of the characteristic polynomial (modulus):'];
fprintf('%s\n',stabilityinfo1);
fprintf(fid,'%s\n',stabilityinfo1);
for ii=1:p
temp=num2str(eigmodulus(ii,1),'%.3f');
   for jj=2:n
   temp=[temp,'  ',num2str(eigmodulus(ii,jj),'%.3f')];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end
if stationary==1;
stabilityinfo2=['No root lies outside the unit circle.'];
stabilityinfo3=['The estimated VAR model satisfies the stability condition'];
fprintf('%s\n',stabilityinfo2);
fprintf(fid,'%s\n',stabilityinfo2);
fprintf('%s\n',stabilityinfo3);
fprintf(fid,'%s\n',stabilityinfo3);
else
stabilityinfo2=['Warning: at leat one root lies on or outside the unit circle.'];
stabilityinfo3=['The estimated VAR model will not be stable'];
fprintf('%s\n',stabilityinfo2);
fprintf(fid,'%s\n',stabilityinfo2);
fprintf('%s\n',stabilityinfo3);
fprintf(fid,'%s\n',stabilityinfo3);
end

fprintf('%s\n','');
fprintf(fid,'%s\n','');
fprintf('%s\n','');
fprintf(fid,'%s\n','');

% display posterior for sigma
sigmainfo=['sigma (residual covariance matrix): posterior estimates'];
fprintf('%s\n',sigmainfo);
fprintf(fid,'%s\n',sigmainfo);
% calculate the (integer) length of the largest number in sigma, for formatting purpose
width=length(sprintf('%d',floor(max(abs(bear.vec(sigma_median))))));
% add a separator, a potential minus sign, and three digits (total=5) to obtain the total space for each entry in the matrix
width=width+5;
for ii=1:n
temp=[];
   for jj=1:n
   % convert matrix entry into string
   number=num2str(sigma_median(ii,jj),'% .3f');
      % pad potential missing blanks
      while numel(number)<width
      number=[' ' number];
      end
   number=[number '  '];
   temp=[temp number];
   end
fprintf('%s\n',temp);
fprintf(fid,'%s\n',temp);
end

fclose(fid);

% Finally, display the results in terms of graph
if pref.plot

% graph the empirical posterior distribution of the VAR parameters
% first reshape the gibbs sampler output to be compatible with the way Matlab uses subplot
% first concatenate psi_gibbs on beta_gibbs to obtain a vector similar to beta in the classical BVAR model
temp=[];
for ii=1:n
temp=[temp;beta_gibbs((ii-1)*k1+1:ii*k1,:)];
   for jj=1:m
   temp=[temp;psi_gibbs((jj-1)*m+ii,:)];
   end
end
beta_reshape=temp;
[beta_swap]=bear.betaswap(beta_reshape,n,m,p,k1+m);
plotvar=1;
plotlag=1;
plotexo=1;
plotconst=1;
% then plot the figure
%{
postdis=figure;
set(postdis,'Color',[0.9 0.9 0.9]);
set(postdis,'name','posterior ditribution of VAR coefficients');
for ii=1:(q1+q2)
subplot(n,k1+m,ii)
hist(beta_swap(ii,:),20);
histproperties=findobj(gca,'Type','patch');
set(histproperties,'FaceColor',[0.7 0.78 1],'EdgeColor',[0 0 0],'LineWidth',0.1);
set(gca,'FontName','Times New Roman');
   % top labels for endogenous
   if ii<=n*p
   temp=[endo{plotvar,1} '(-' num2str(plotlag) ')'];
   title(temp,'FontWeight','normal');
      if plotlag<p
      plotlag=plotlag+1;
      elseif plotlag==p
      plotlag=1;
      plotvar=plotvar+1;
      end
   end
   % top labels for exogenous
   %if ii==k1+1
   %title('Constant','FontWeight','normal');
   %end
   %if ii>k1+1 && ii<=k1+m
   %title(exo{plotexo,1},'FontWeight','normal');
   %plotexo=plotexo+1;
   %end
   % side labels
   %if rem((ii-1)/(k1+m),1)==0
   ylabel(endo{(ii-1)/(k1+m)+1,1},'FontWeight','normal');
   %end
end
%}

% plot actual vs. fitted
actualfitted=figure;
set(actualfitted,'Color',[0.9 0.9 0.9]);
set(actualfitted,'name','model estimation: actual vs fitted')
ncolumns=ceil(n^0.5);
nrows=ceil(n/ncolumns);
for ii=1:n
subplot(nrows,ncolumns,ii)
hold on
plot(decimaldates1,Y(:,ii),'Color',[0 0 0],'LineWidth',2);
plot(decimaldates1,Ytilde(:,ii),'Color',[1 0 0],'LineWidth',2);
hold off
set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
title(endo{ii,1},'FontName','Times New Roman','FontSize',10);
   if ii==1
   plotlegend=legend('actual','fitted');
   set(plotlegend,'FontName','Times New Roman');
   end
end

% plot the residuals
residuals=figure;
set(residuals,'Color',[0.9 0.9 0.9]);
set(residuals,'name','model estimation: residuals')
for ii=1:n
subplot(nrows,ncolumns,ii)
plot(decimaldates1,EPStilde(:,ii),'Color',[0 0 0],'LineWidth',2)
set(gca,'XLim',[decimaldates1(1,1) decimaldates1(end,1)],'FontName','Times New Roman');
title(endo{ii,1},'FontName','Times New Roman','FontSize',10);
end

end

% finally, save the results on excel
bear.data.excelrecord2fcn(T, n, endo, stringdates1, Y, Ytilde, pref, EPStilde)

author: european-central-bank

repo: BEAR-toolbox


text:
% DEMROBOTWIRELESSFGPLVM3 Wireless Robot data from University of Washington with dynamics and no back constraints.

% FGPLVM

% Fix seeds
randn('seed', 1e5);
rand('seed', 1e5);

dataSetName = 'robotWireless';
experimentNo = 3;

% load data
[Y, lbls] = lvmLoadData(dataSetName);

% Set up model
options = fgplvmOptions('ftc');

latentDim = 2;
d = size(Y, 2);

model = fgplvmCreate(latentDim, d, Y, options);

% Add dynamics model.
options = gpOptions('ftc');
options.kern = kernCreate(model.X, {'rbf', 'white'});
options.kern.comp{1}.inverseWidth = 0.2;
% This gives signal to noise of 0.1:1e-3 or 100:1.
options.kern.comp{1}.variance = 0.1^2;
options.kern.comp{2}.variance = 1e-3^2;
model = fgplvmAddDynamics(model, 'gp', options);

% Optimise the model.
iters = 1000;
display = 1;

model = fgplvmOptimise(model, display, iters);

% Save the results.
modelWriteResult(model, dataSetName, experimentNo);

if exist('printDiagram') & printDiagram
  lvmPrintPlot(model, lbls, dataSetName, experimentNo);
end

% Load the results and display dynamically.
lvmResultsDynamic(model.type, dataSetName, experimentNo, 'vector')

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% Copyright 2009 - 2010 The MathWorks, Inc.
function plot_trajectory(z,y)
%#eml
eml.extrinsic('title','xlabel','ylabel','plot','axis','pause');
title('Trajectory of object [blue] its Kalman estimate[green]');
xlabel('horizontal position');
ylabel('vertical position');
plot(z(1), z(2), 'bx-');
plot(y(1), y(2), 'go-');
axis([-1.1, 1.1, -1.1, 1.1]);
pause(0.02);
end


text:
function [I,x,W,F] = sphericalradial(f,m,P,param)
% SPHERICALRADIAL - ND scaled spherical-radial cubature rule
%
% Syntax:
%   [I,x,W,F] = sphericalradial(f,m,P[,param])
%
% In:
%   f - Function f(x,param) as inline, name or reference
%   m - Mean of the d-dimensional Gaussian distribution
%   P - Covariance of the Gaussian distribution
%   param - Parameters for the function (optional)
%
% Out:
%   I - The integral
%   x - Evaluation points
%   W - Weights
%   F - Function values
%
% Description:
%   Apply the spherical-radial cubature rule to integrals of form:
%     int f(x) N(x | m,P) dx

% History:
%   Aug 5, 2010 - Renamed from 'cubature' to 'sphericalradial' (asolin)

% Copyright (c) 2010 Arno Solin
%
% This software is distributed under the GNU General Public
% Licence (version 2 or later); please refer to the file
% Licence.txt, included with the software, for details.

%% Spherical-radial cubature rule

  % The dimension of m is
  n = size(m,1);

  % Evaluation points (nx2n)
  x = [eye(n) -eye(n)];

  % Scaling
  x = sqrt(n)*x;
  x = chol(P)'*x + repmat(m,1,2*n);

  % Evaluate the function at the points
    if ischar(f) || strcmp(class(f),'function_handle')
      if nargin < 4
        F = feval(f,x);
      else
        F = feval(f,x,param);
      end
    elseif isnumeric(f)
         F = f*x;
    else
      if nargin < 4
         F = f(x);
      else
         F = f(x,param);
      end
    end

  % The weights are
  W = 1/(2*n);

  % Return integral value
  I = W*sum(F,2);

  % Return weights
  W = W*ones(1,2*n);

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% SABR European option pricer comparison
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Descritpion: Script to Compare European Option pricers in SABR Model
% Author:      Justin Kirkby
% References:  (1) General Valuation Framework for SABR and Stochastic Local Volatility
%                   Models. SIAM J. Financial Mathematics, 2018.
%              (2) Full Fledge SABR with Markov Chains. Wimott, 2019.
% Disclaimer: this is research code, not production code. The parameter settings etc
%            should not be expected to work as is in all cases. Determining the right parameter
%            settings will depend on your usecase. This is left up to the user.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[folder, name, ext] = fileparts(which( mfilename('fullpath')));
cd(folder);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Step 1) CHOOSE CONTRACT/GENERAL PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T    = 1;         % time to matruity
r    = 0.05;      % interest rate
F_0  = 1.1;       % Futures price
K    = F_0*0.95;  % strike
call = 1;         % 1 = call option, else put

% -----------------
% Model Params
% -----------------
ModParams.beta   = .6;
ModParams.alpha  = 0.08;
ModParams.v0     = 0.2;
ModParams.rho    = 0.0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Reference Price (Using Obloj Asymptotic Formula with correction)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath('../../Approx/SABR')
addpath('../../Analytical/BlackScholes')

tic
IV = SABR_Hagan_Obloj_ImpliedVol( K, F_0, ModParams.v0, T, ModParams.alpha, ModParams.beta, ModParams.rho );
price_obloj = BSM_Greeks( 0, F_0/exp(r*T), IV, r, 0, T, K, call)
time_oblog=toc;

price_ref = price_obloj;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Antonov Approximation (Using Asymptotic Formula)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
price_ant = SABR_European_AntonovApprox(F_0,K,T,call,r,ModParams)
time_ant = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Monte Carlo (Euler)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath('../../Monte_Carlo')
addpath('../../Monte_Carlo/European')
M = 400;
N_sim = 5*10^4;

tic;
Spath = Simulate_SLV_func( N_sim, M, T, F_0, 0, 0, 2, ModParams);
disc = exp(-r*T);
[price_MC, stdErr] = Price_MC_European_Strikes_func(Spath, disc, call, K );

price_MC_L = price_MC - 2*stdErr; price_MC_U = price_MC + 2*stdErr;
time_MC = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  CTMC-Double Layer Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath('../../PROJ/SABR/European_American_Barrier/')
addpath('../../PROJ/SABR/Helper_Functions/')

% -----------------
% Numerical Params
% -----------------
M    = 100;   % number of time steps
contract_type = 1;    % 1 = European, 2 = American, 3 = Down and Out Barrier
L    = 0.6*F_0;     % For barrier contract, this is the barrier

CTMCParams.m_0 = 26;
CTMCParams.N = 90;
CTMCParams.gridMult_v = 0.5;
CTMCParams.gridMult_s = 0.05;  %Grid mult param for S
CTMCParams.gamma = 6;  %Grid width param for variance grid

tic
price_ctmc = SABR_EurBarAmer_func(call, M, T, F_0, K, r, CTMCParams, ModParams, contract_type, L)
time_ctmc = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% COMPARE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n---------------------------------------------\n')
fprintf('Method      |    Price    |    Err   |  CPU \n')
fprintf('---------------------------------------------\n')
fprintf('Ref (Obloj) | %.8f  |          |       \n', price_ref)
fprintf('---------------------------------------------\n')
fprintf('Antonov     | %.8f  | %.2e | %.4f \n', price_ant, abs(price_ref-price_ant), time_ant)
fprintf('CTMC        | %.8f  | %.2e | %.4f \n', price_ctmc, abs(price_ref-price_ctmc), time_ctmc)
fprintf('MC-Euler    |[%.3f,%.3f]| %.2e | %.4f \n', price_MC_L, price_MC_U, abs(price_ref-price_MC), time_MC)
fprintf('---------------------------------------------\n')


text:
%  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 St, Fifth Floor, Boston, MA  02110-1301  USA
%function tfrgabot
%TFRGABT Unit test for the function tfrgabor.

%	O. Lemoine - February 1996.

N1=256;
M=16;
N=16;

% Test of the biorthogonality between h and gam (dual frame window)
Nh=65;					% length of window h
h=amgauss(Nh); h=h/norm(h);		% h must be of odd length
n0=1;
m0=1;
sig=zeros(N1,1);
sig(((n0-1)*N+1):((n0-1)*N+Nh))=h.*fmconst(Nh,(m0-1)/M);
tfr=tfrgabor(sig,N,1,h);		% Critical sampling case
if abs(tfr(1,1)-1)>sqrt(eps),
  error('tfrgabor test 1 failed');
end
errors=find(any(tfr>eps));
if length(errors)~=1,
  error('tfrgabor test 2 failed');
end

% Localization
Nh=37;					% length of window h
h=amgauss(Nh); h=h/norm(h);		% h must be of odd length
n0=8;
m0=8;
sig=zeros(N1,1);
sig(((n0-1)*N+1):((n0-1)*N+Nh))=h.*fmconst(Nh,(m0-1)/M);
tfr=tfrgabor(sig,N,1,h);		% Critical sampling case
if abs(tfr(m0,n0)-1)>sqrt(eps),		% C(n0,m0)=1
  error('tfrgabor test 3 failed');
end
errors=find(any(tfr>eps));		% Only one non-zero coeff
if length(errors)~=1,
  error('tfrgabor test 4 failed');
end

% For another window and position of sig
Nh=19;					% length of window h
h=amexpo2s(Nh);	h=h/norm(h);		% h must be of odd length
n0=13;
m0=5;
sig=zeros(N1,1);
sig(((n0-1)*N+1):((n0-1)*N+Nh))=h.*fmconst(Nh,(m0-1)/M);
tfr=tfrgabor(sig,N,1,h);		% Critical sampling case
if abs(tfr(m0,n0)-1)>sqrt(eps),		% C(n0,m0)=1
  error('tfrgabor test 5 failed');
end
errors=find(any(tfr>eps));		% Only one non-zero coeff
if length(errors)~=1,
  error('tfrgabor test 6 failed');
end

% For another window and position of sig
Nh=17;					% length of window h
h=amexpo2s(Nh);	h=h/norm(h);		% h must be of odd length
n0=15;
m0=8;
sig=zeros(N1,1);
sig(((n0-1)*N+1):((n0-1)*N+Nh))=h.*fmconst(Nh,(m0-1)/M);
tfr=tfrgabor(sig,N,1,h);		% Critical sampling case
if abs(tfr(m0,n0)-1)>sqrt(eps),		% C(n0,m0)=1
  error('tfrgabor test 7 failed');
end
errors=find(any(tfr>eps));		% Only one non-zero coeff
if length(errors)~=1,
  error('tfrgabor test 8 failed');
end

% Synthesis for q=1
Nh=33;q=1;
h=amgauss(Nh); h=h/norm(h);
sig=zeros(N1,1);
sig=amgauss(N1);
sigr=zeros(N1,1);
[tfr,dgr]=tfrgabor(sig,N,q,h);
alpha=round((2*N1/N-1-Nh)/(2*q));
hN1=zeros(N1,1);
hN1((N1-(Nh-1))/2-alpha:(N1+Nh-1)/2-alpha)=h;
n=1:N;
som=M*ifft(dgr);
for k=1:N1,
  indice=modulo(k-round(n*M/q),N1);
  sigr(k)=som(modulo(k,M),:)*fftshift(hN1(indice));
end
if any(abs(sig-sigr)>sqrt(eps))~=0,
  error('tfrgabor test 9 failed');
end

% Synthesis for q=4
Nh=33;q=4;N=32;M=32;
h=amgauss(Nh); h=h/norm(h);
sig=zeros(N1,1);
sig=amgauss(N1);
sigr=zeros(N1,1);
[tfr,dgr]=tfrgabor(sig,N,q,h);
alpha=round((2*N1/N-1-Nh)/(2*q));
hN1=zeros(N1,1);
hN1((N1-(Nh-1))/2-alpha:(N1+Nh-1)/2-alpha)=h;
n=1:N;
som=M*ifft(dgr);
for k=1:N1,
  indice=modulo(k-round(n*M/q),N1);
  sigr(k)=som(modulo(k,M),:)*fftshift(hN1(indice));
end
if any(abs(sig-sigr)>sqrt(eps))~=0,
  error('tfrgabor test 10 failed');
end

author: HeLiangHIT


text:
% Fig. 6.46   Feedback Control of Dynamic Systems, 5e
%             Franklin, Powell, Emami
%

clear all;
close all;

num=1;
den=[1 0 0];
w=logspace(-2,2,100);
[m,p]=bode(num,den,w);
loglog(w,m);
grid;
xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Fig. 6.46 Spacecraft frequency-response magnitude');


text:
function [intersects, inds] = intersectEdgePolygon(edge, poly, varargin)
%INTERSECTEDGEPOLYGON  Intersection point of an edge with a polygon
%
%   INTER = intersectEdgePolygon(EDGE, POLY)
%   Computes intersection(s) point(s) between the edge EDGE and the polygon
%   POLY. EDGE is given by [x1 y1 x2 y2]. POLY is a N-by-2 array of vertex
%   coordinates.
%   INTER is a M-by-2 array containing coordinates of intersection(s). It
%   can be empty if no intersection is found.
%
%   [INTER, INDS] = intersectEdgePolygon(EDGE, POLY)
%   Also returns index/indices of edge(s) involved in intersections.
%
%   Example
%   % Intersection of an edge with a square
%     poly = [0 0;10 0;10 10;0 10];
%     edge = [9 2 9+3*1 2+3*2];
%     exp = [10 4];
%     inter = intersectEdgePolygon(edge, poly)
%     ans =
%         10   4
%
%   See also
%   edges2d, polygons2d, intersectLinePolygon, intersectRayPolygon
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2012-02-24,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2012 INRA - Cepia Software Platform.

% get computation tolerance
tol = 1e-14;
if ~isempty(varargin)
    tol = varargin{1};
end

% get supporting line of edge
line = edgeToLine(edge);

% compute all intersections of supporting line with polygon
[intersects, inds] = intersectLinePolygon(line, poly, tol);

% keep only intersection points located on the edge
if ~isempty(intersects)
    pos = linePosition(intersects, line);
    keep = pos >= -tol & pos <= (1+tol);
    intersects = intersects(keep, :);
    inds = inds(keep);
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function plot_circbootfitF()
    % plot_circbootfitF Selects earthquakes in the radius ra around a grid node and calculates the forecast
    % by using calc_bootfitF.m
    %
    % Jochen Woessner
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals
    report_this_filefun();

    ShapeGeneral.clearplot(); % was axes h1, tag plos1

    % interactively get the circle of interest
    shape=ShapeCircle();
    [ZG.newt2, max_km] = selectCircle(newa, shape.toStruct());

    % Select radius in time
    newt3=ZG.newt2;
    vSel = (ZG.newt2.Date <= ZG.maepi.Date+days(time));
    ZG.newt2 = ZG.newt2.subset(vSel);
    R2 = ZG.newt2.Count;
    messtext = ['Number of selected events: ' num2str(ZG.newt2.Count)  ];
    disp(messtext)

    % Sort the catalog
    ZG.newt2.sort('Date');
    R2 = ra;

    % Plot selected earthquakes
    shape.plot([],ZG.newt2); % linespec was xk, tag was plos1

    % Compute and Plot the forecast
    calc_bootfitF(newt3.Date,time,timef,bootloops,ZG.maepi.Date)

    ZG.newcat = ZG.newt2;
    ctp=CumTimePlot(ZG.newt2);
    ctp.plot();

end

author: CelsoReyes


text:
function [vert, fac] = voxel_image( pts, vox_sz, color, alpha, edgec )
%VOXEL_IMAGE Creates a 3D voxel image
%   Parameters:
%   pts    - n x 3 matrix with 3D points
%   vox_sz - 1 x 3 vector with voxel size - if vox_sz is a scalar, all
%            edges will have the same length
%   color  - face color
%   alpha  - face alpha (opacity)
%   edgec  - edge color
%
%   Return values:
%   vert   - 8n x 3 matrix containing the vertices of the voxels
%   fac    - 6n x 4 matrixes containing the indexes of vert that form a
%            face

% Example:
%
% pts = [1,1,1; 2,2,2];
% vs = [0.5, 0.5, 0.5];
% voxel_image(pts, vs);
% view([-37.5, 30]);
%
% This example creates two voxels (at (1,1,1) and (2,2,2)) with edges of
% length 0.5
%
% Author: Stefan Schalk, 11 Feb 2011

if (nargin < 1)
    error('No input arguments given');
end
if (nargin < 2)
    vox_sz = [1,1,1];
end
if (nargin < 3)
    color = 'b';
end
if (nargin < 4)
    alpha = 1;
end
if (nargin < 5)
    edgec = 'k';
end

if (size(pts,2) ~= 3)
    error('pts should be an n x 3 matrix');
end
if (isscalar(vox_sz))
    vox_sz = vox_sz*ones(1,3);
end
if (size(vox_sz,1) ~= 1 || size(vox_sz,2) ~= 3)
    error('vox_sz should be an 1 x 3 vector');
end

np = size(pts,1);
vert = zeros(8*np,3);
fac = zeros(6*np,4,'uint32');
vert_bas = [...
        -0.5,-0.5,-0.5;
        0.5,-0.5,-0.5;
        0.5,0.5,-0.5;
        -0.5,0.5,-0.5;
        -0.5,-0.5,0.5;
        0.5,-0.5,0.5;
        0.5,0.5,0.5;
        -0.5,0.5,0.5];
vert_bas = vert_bas.*([vox_sz(1).*ones(8,1), vox_sz(2).*ones(8,1), vox_sz(3).*ones(8,1)]);
fac_bas = [...
        1,2,3,4;
        1,2,6,5;
        2,3,7,6;
        3,4,8,7;
        4,1,5,8;
        5,6,7,8];
for vx = 1:np
    a = ((vx-1)*8+1):vx*8;
    for dim = 1:3
        vert( a,dim ) = vert_bas(:,dim) + pts(vx,dim);
    end
    fac ( ((vx-1)*6+1):vx*6,: ) = (vx - 1)*8*ones(6,4) + fac_bas;
end
patch('Vertices',vert,'Faces',fac,'FaceColor',color,'FaceAlpha',alpha,'Edgecolor',edgec);

end


text:
function roiColorRgba = dtiRoiGetColor(roi,alpha)
%
%  roiColorRgba = dtiRoiGetColor(roi,[defaultAlpha=0.75])
%
% If the roi does not specify an alpha channel (opacity), then defaultAlpha
% is appended to the rgb triplet.
%
% HISTORY:
% 2006.08.18 RFD: extracted from dtiSaveImageSlicesOverlay.
% 2006.11.06 RFD: added csome comments about alpha.

% Default transparency (alpha) channel
if ~exist('alpha','var'), alpha = 0.75; end

if(isnumeric(roi.color))
  roiColorRgba = roi.color(:)';
  if(length(roiColorRgba)==3) roiColorRgba(4) = alpha;
  elseif(length(roiColorRgba)~=4) error('ROI color must be 1x3 or 1x4!'); end
else
  switch(roi.color)
   case 'r'
    roiColorRgba = [1 0 0 alpha];
   case 'g'
    roiColorRgba = [0 1 0 alpha];
   case 'b'
    roiColorRgba = [0 0 1 alpha];
   case 'm'
    roiColorRgba = [1 0 1 alpha];
   case 'w'
    roiColorRgba = [1 1 1 alpha];
   case 'y'
    roiColorRgba = [1 1 0 alpha];
   case 'c'
    roiColorRgba = [0 1 1 alpha];
   case 'k'
    roiColorRgba = [0 0 0 alpha];
   otherwise
    error('Unknown ROI color');
  end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function this = subsasgn(this, subs, A)
% Subscript assignment for GIfTI objects
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: subsasgn.m 6513 2015-08-05 17:52:13Z guillaume $

switch subs(1).type
    case '.'
        if ~ismember(subs(1).subs, ...
                {'vertices' 'faces' 'normals' 'cdata','mat','indices','private'})
            error('Reference to non-existent field ''%s''.',subs(1).subs);
        end
        % TODO % handle cases when length(subs) > 1
        [i,n] = isintent(this,subs(1).subs);
        if isempty(i) && ~strcmp(subs(1).subs,'private')
            n = length(this.data) + 1;
            if n==1, this.data = {}; end
            % TODO % Initialise data field appropriately
            this.data{n}.metadata = struct([]);
            this.data{n}.space    = [];
            this.data{n}.attributes.Dim = size(A);
            % TODO % set DataType according to intent type
            this.data{n}.data = [];
            switch subs(1).subs
                case {'vertices','mat'}
                    in = 'NIFTI_INTENT_POINTSET';
                    dt = 'NIFTI_TYPE_FLOAT32';
                    this.data{n}.space.DataSpace = 'NIFTI_XFORM_UNKNOWN';
                    this.data{n}.space.TransformedSpace = 'NIFTI_XFORM_UNKNOWN';
                    this.data{n}.space.MatrixData = eye(4);
                case 'faces'
                    in = 'NIFTI_INTENT_TRIANGLE';
                    dt = 'NIFTI_TYPE_INT32';
                case 'indices'
                    in = 'NIFTI_INTENT_NODE_INDEX';
                    dt = 'NIFTI_TYPE_INT32';
                case 'normals'
                    in = 'NIFTI_INTENT_VECTOR';
                    dt = 'NIFTI_TYPE_FLOAT32';
                case 'cdata'
                    in = 'NIFTI_INTENT_NONE';
                    dt = 'NIFTI_TYPE_FLOAT32';
                otherwise
                    error('This should not happen.');
            end
            this.data{n}.attributes.Intent = in;
            this.data{n}.attributes.DataType = dt;
        end

        switch subs(1).subs
            %- .private
            %--------------------------------------------------------------
            case 'private'
                this = builtin('subsasgn',this,subs(2:end),A);

            % .mat
            %--------------------------------------------------------------
            case 'mat'
                if length(subs) > 1
                    this.data{n}.space(1).MatrixData = builtin('subsasgn',...
                        this.data{n}.space(1).MatrixData,subs(2:end),A);
                else
                    if ~isequal(size(A),[4 4])
                        error('Invalid Coordinate System Transform Matrix.');
                    end
                    this.data{n}.space(1).MatrixData = A;
                end

            %- .faces
            %--------------------------------------------------------------
            case 'faces'
                if length(subs) > 1
                    this.data{n}.data = int32(builtin('subsasgn',this.data{n}.data,subs(2:end),A-1));
                else
                    this.data{n}.data = int32(A - 1);
                    this.data{n}.attributes.Dim = size(A);
                end

            %- .indices
            %--------------------------------------------------------------
            case 'indices'
                if n ~= 1
                    this.data = this.data([n setdiff(1:numel(this.data),n)]);
                    n = 1;
                end
                if length(subs) > 1
                    this.data{n}.data = int32(builtin('subsasgn',this.data{n}.data,subs(2:end),A-1));
                else
                    A = A(:);
                    this.data{n}.data = int32(A - 1);
                    this.data{n}.attributes.Dim = size(A);
                end

            %- .vertices, .normals, .cdata
            %--------------------------------------------------------------
            otherwise
                if length(subs) > 1
                    if numel(n) == 1
                        this.data{n}.data = single(builtin('subsasgn',this.data{n}.data,subs(2:end),A));
                        this.data{n}.attributes.Dim = size(this.data{n}.data);
                    else
                        if numel(subs(2).subs) == 1
                            error('Linear indexing not supported: use multiple subscripts.');
                        end
                        idx = subs(2).subs{2};
                        if isequal(idx,':'), idx = 1:numel(this.data); end
                        for k=1:numel(idx)
                            s = subs(2);
                            s.subs{2} = 1;
                            if numel(A) == 1
                                this.data{idx(k)}.data = single(builtin('subsasgn',this.data{idx(k)}.data,s,A));
                            else
                                this.data{idx(k)}.data = single(builtin('subsasgn',this.data{idx(k)}.data,s,A(:,k)));
                            end
                            this.data{idx(k)}.attributes.Dim = size(this.data{idx(k)}.data);
                        end
                    end
                else
                    if numel(n) == 1
                        if isa(A,'file_array')
                            this.data{n}.data = A;
                            this.data{n}.attributes.Dim = A.dim;
                        else
                            this.data{n}.data = single(A);
                            this.data{n}.attributes.Dim = size(A);
                        end
                    else
                        error('Syntax not implemented.');
                    end
                end
        end

    case '()'
    case '{}'
    otherwise
        error('This should not happen.');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [sel1, sel2] = match_str(a, b);

% MATCH_STR looks for matching labels in two listst of strings
% and returns the indices into both the 1st and 2nd list of the matches.
% They will be ordered according to the first input argument.
%
% [sel1, sel2] = match_str(strlist1, strlist2)
%
% The strings can be stored as a char matrix or as an vertical array of
% cells, the matching is done for each row.

% Copyright (C) 2000, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: match_str.m 2865 2011-02-12 19:24:57Z roboos $

% ensure that both are cell-arrays
if isempty(a)
  a = {};
elseif ~iscell(a)
  a = cellstr(a);
end
if isempty(b)
  b = {};
elseif ~iscell(b)
  b = cellstr(b);
end

% regardless of what optimizations are implemented, the code should remain
% functionally compatible to the original, which is
% for i=1:length(a)
%   for j=1:length(b)
%     if strcmp(a(i),b(j))
%       sel1 = [sel1; i];
%       sel2 = [sel2; j];
%     end
%   end
% end

% ensure that both are column vectors
a = a(:);
b = b(:);
Na = numel(a);
Nb = numel(b);

% replace all unique strings by a unique number and use the fact that
% numeric comparisons are much faster than string comparisons
[dum1, dum2, c] = unique([a; b]);
a = c(1:Na);
b = c((Na+1):end);

sel1 = [];
sel2 = [];
for i=1:length(a)
  % s = find(strcmp(a(i), b));  % for string comparison
  s = find(a(i)==b);            % for numeric comparison
  sel2 = [sel2; s];
  s(:) = i;
  sel1 = [sel1; s];
end

author: PatternRecognition

repo: OpenBMI


text:
function baseName=baseNameFromFileName(fileName)

[~,baseName]=fileparts(fileName);

end

author: kristinbranson

repo: JAABA


text:
function callBoxDeleteCallback(rectangle,evt)

% Only do anything if the face was selected, not the markers
if ~strcmp(evt.SelectedPart, 'face')
    return
end

hObject = get(rectangle,'Parent');
handles = guidata(hObject);
clicked_tag = str2double(get(rectangle, 'Tag'));

switch  evt.SelectionType
    case 'right' % delete if right click
        handles.data.calls(clicked_tag,:) = [];
        SortCalls(hObject, [], handles, 'time', 0, clicked_tag - 1);
    case 'left' % Make it the current call if left clicked
        if clicked_tag ~= handles.data.currentcall
            handles.data.currentcall = clicked_tag;
            update_fig(hObject, [], handles)
        end
    case {'ctrl', 'double'}
        choosedialog(rectangle,handles);
end
end

function choosedialog(rectangle,handles)

hObject = handles.figure1;
axes_position = get(handles.focusWindow, 'Position');

axes_x_lim = xlim(handles.focusWindow);
axes_y_lim = ylim(handles.focusWindow);
position = get(rectangle, 'Position');
tag = get(rectangle, 'Tag');

dialog_width = 354;
dialog_heigth = 154;

rect_x_relative_to_axes = (position(1) + position(3) - axes_x_lim(1)) / (axes_x_lim(2) - axes_x_lim(1));
rect_y_relative_to_axes = (position(2) + position(4) - axes_y_lim(1)) / (axes_y_lim(2) - axes_y_lim(1));

set(hObject,'Units','Pixels');
fig_pos = get(hObject, 'Position');

x_position = fig_pos(1) + fig_pos(3) * (axes_position(1) + axes_position(3)*rect_x_relative_to_axes);
y_position = fig_pos(2) + fig_pos(4) * (axes_position(2) + axes_position(4)*rect_y_relative_to_axes) - dialog_heigth;

%Make sure the dialog stays within the figure
x_position = min(x_position,fig_pos(3)-dialog_width);
y_position = min(y_position,fig_pos(4)-dialog_heigth);

d = dialog('Position',[x_position y_position dialog_width dialog_heigth],'Name',['Change label for call ' tag],'Units','normal');

% Get the number of rows to columns
n_columns = floor(sqrt(length(handles.data.settings.labels)));
n_rows = ceil(length(handles.data.settings.labels) / n_columns);

button_width = 1 ./ n_columns;
button_height = 1 ./ n_rows;
label_index = 0;
padding = .1;

for c = 1:n_columns
    for r = n_rows:-1:1
        label_index = label_index + 1;
        if label_index > length(handles.data.settings.labels)
            continue
        end
        x_position = (c-1) / n_columns + button_width*padding/2;
        y_position = (r-1) / n_rows + button_height*padding/2;
        uicontrol('Parent',d,...
            'Units', 'normalized',...
            'Position',[x_position, y_position, button_width*(1-padding), button_height*(1-padding)],...
            'String',handles.data.settings.labels{label_index},...
            'Callback',{@label_callback,handles.data.settings.labels{label_index}},...
            'BackgroundColor',[0.302,0.502,0.302],...
            'ForegroundColor',[1.0 1.0 1.0]);
    end
end

uiwait(d);

    function label_callback(~,~,new_type)
        clicked_tag = str2double(tag);
        if ~isempty(i)
            handles.data.calls{clicked_tag,'Type'} = {new_type};
            delete(gcf);
            update_fig(hObject, [], handles);
        end
    end

end

author: DrCoffey

repo: DeepSqueak

save_path: github-repos/MATLAB/DrCoffey-DeepSqueak


text:
% StackExchange Signal Processing Q76344
% https://dsp.stackexchange.com/questions/76344
% Generate the Matrix Form of 1D Convolution Kernel
% References:
%   1.
% Remarks:
%   1.  B
% TODO:
% 	1.  C
% Release Notes
% - 1.0.000     19/07/2021
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

CONVOLUTION_SHAPE_FULL         = 1;
CONVOLUTION_SHAPE_SAME         = 2;
CONVOLUTION_SHAPE_VALID        = 3;

%% Simulation Parameters

vK              = [1; 4; 6; 4; 1] / 16; %<! Regularization
numSamples      = 50;
vConvShape      = [1; 2; 3];
cConvShapeStr   = {['Full'], ['Same'], ['Valid']};

%% Generate Data

cConvMtx = cell(3, 1);

for ii = 1:length(vConvShape)
    cConvMtx{ii} = CreateConvMtx1D(vK, numSamples, vConvShape(ii));
end

%% Display Results

figureIdx = figureIdx + 1;

hFigure     = figure('Position', [100, 100, 1200, 600]);
for ii = 1:length(vConvShape)
    hAxes       = subplot(1, 3, ii);
    hScatterObj = ScatterSparse(cConvMtx{ii}, 'fill');
    set(hAxes, 'YDir', 'reverse');
    set(hAxes, 'DataAspectRatio', [1, 1, 1]);
    set(hAxes, 'YLim', [0, numSamples + size(vK, 1) - 1]);
    % set(hLineObj, 'LineWidth', lineWidthNormal);
    set(get(hAxes, 'Title'), 'String', {['Convolution Type: ', cConvShapeStr{ii}]}, ...
        'FontSize', fontSizeTitle);
    % set(get(hAxes, 'XLabel'), 'String', {['Sample Index']}, ...
    %     'FontSize', fontSizeAxis);
    % set(get(hAxes, 'YLabel'), 'String', {['Sample Value']}, ...
    %     'FontSize', fontSizeAxis);
    % hLegend = ClickableLegend({['Ground Truth Signal'], ['Measured Signal']});
end

if(generateFigures == ON)
    % saveas(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
    print(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png'], '-dpng', '-r0'); %<! Saves as Screen Resolution
end

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function [val table]= muller(strEqn,strVar,p0,p1,p2,tol,n)
%MULLER uses the Muller's method to solve for a root of f(x)
%
% Find a solution to f(x) = 0 given three approximations p0, p1
% and p2 (Note that Muller's method may return complex roots)
%
%
%  INPUT:
%         strEqn: is the string representation of the equation f(x)
%         strVar: is the varialbe in strEqn
%     p0, p1, p2: initial guesses at the solution
%            tol: is the solution tolerence
%              n: is the maximum number of iterations of the algorithum.
%
%  OUTPUT:   val: is the approximate solution (can be complex)
%          table: is the iteration table
%
%
%  USAGE:   val = muller(strEqn,strVar,x0,x1,x2,tol,n)
%          [val table] = muller(...) to display a table of the iterations
%           val = muller('x^3+2*x^2-5','x',0,1,2,10^-4,100)
%
%     BY: Jonathan Lister (jlister1@utk.edu)

%% Additonal Examples
%
%           strEqn = '16*x^4-40*x^3+5*x^2+20*x+6';
%           strVar = 'x';
%           tol = 10^-5;
%
%           disp('p0 = 0.5, p1 = -0.5, p2 = 0.0')
%           [root1 table1]= muller(strEqn,strVar,0.5,-0.5,0,tol,100)
%
%           disp('p0 = 0.5, p1 = 1.0, p2 = 1.5')
%           [root2 table2]= muller(strEqn,strVar,0.5,1.0,1.5,tol,100)
%
%           disp('p0 = 2.5, p1 = 2.0, p2 = 2.25')
%           [root3 table3]= muller(strEqn,strVar,2.5,2.0,2.25,tol,100)
%

%% Input Error Checking
if nargin < 7
    error('Incorrect number of input arguments type "help muller"')
end

if tol < 0
    error('tolerance must be a positive number')
end

if ~ischar(strEqn)
    error('Input 1 Error, a string of the equation of interest is required')
end

if ~ischar(strVar)
    error('Input 3 Error, a string representation of the variable is required')
end

%parse equation
F = vectorize(inline(strEqn,strVar));

%intitialize output variables
val = [];
table ={};

%% Algorithm 2.8 from text.
%Step 1
h1 = p1 - p0;
h2 = p2 - p1;
del1 = (F(p1)-F(p0))./h1;
del2 = (F(p2)-F(p1))./h2;
d = (del2-del1)./(h2+h1);
I = 3;

%Step 2
while I <= n
    %Step 3
    b = del2+h2.*d;
    D = sqrt(b.^2-4.*F(p2).*d); % could be complex
    %Step 4
    if abs(b - D) < abs(b + D)
        E = b + D;
    else
        E = b - D;
    end
    %Step 5
    h = -2.*F(p2)./E;
    p = p2 + h;

    if I == 3
        table{1} = 'Muller''s Method Iterations';
        table{2}='  I             P                     f(P)           ';
        table{3}='-----------------------------------------------------';
    end
    str = sprintf('%3u:  % 6.6f + %6.6fi % 6.6f + %6.6fi',I,real(p),imag(p),real(F(p)),imag(F(p)));
    table{I + 1} = str; %#ok<*AGROW>

    %Step 6
    if abs(h) < tol
        val = p;
        table = char(table);
        break
    end
    p0 = p1;
    p1 = p2;
    p2 = p;
    h1 = p1 - p0;
    h2 = p2 - p1;
    del1 = (F(p1)-F(p0))./h1;
    del2 = (F(p2)-F(p1))./h2;
    d = (del2-del1)./(h2+h1);
    I = I + 1;
end

if isempty(val)
    disp('The procedure was unsuccessful.')
    table = char(table);
end


text:
function [uout, tout] = spin2(varargin)
%SPIN2  Solve stiff PDEs in 2D periodic domains, Fourier spectral method and
%exponential integrators.
%
%   UOUT = SPIN2(PDECHAR) solves the PDE specified by the string PDECHAR, and
%   plays a movie of the solution. Possible strings include 'GL' and 'GS' for
%   the Ginzburg-Landau and Gray-Scott equations. Other PDEs are available, see
%   Remark 1 and Examples 1-4. The output UOUT is a CHEBFUN2 corresponding to
%   the solution at the final time (a CHEBMATRIX for systems of equations, each
%   row representing one variable).
%
%   UOUT = SPIN2(S, N, DT) solves the PDE specified by the SPINOP2 S with N grid
%   points in each direction and time-step DT, and plays a movie of the solution.
%   See HELP/SPINOP2 and Example 5.
%
%   [UOUT, TOUT] = SPIN2(...) also returns the times chunks TOUT at which UOUT
%   was computed.
%
%   Users of SPIN2 will quickly find they want to vary aspects of the plotting.
%   The fully general syntax for this involves using preferences specified by
%   a SPINPREF2 object PREF. See HELP/SPINPREF2 and Example 6. However for many
%   purposes it is most convenient to use the syntax
%
%   UOUT = SPIN2(..., 'PREF1', VALUE1, 'PREF2', VALUE2, ...)
%
%   For example:
%
%   UOUT = SPIN2(..., 'Clim', [a b]) changes colorbar limits to [a b]
%   UOUT = SPIN2(..., 'colormap', 'jet') changes the colormap to 'jet'
%   UOUT = SPIN2(..., 'dataplot', 'abs') plots absolute value
%   UOUT = SPIN2(..., 'iterplot', 4) plots only every 4th time step
%   UOUT = SPIN2(..., 'Nplot', 256) plays a movie at 256x256 resolution
%   UOUT = SPIN2(..., 'plot', 'off') for no movie
%   UOUT = SPIN2(..., 'view', [a b]) changes the view angle to [a b]
%
% Remark 1: List of PDEs (case-insensitive)
%
%    - 'GL' for the Ginzburg-Landau equation,
%    - 'GS' for the Gray-Scott equations (stripes),
%    - 'GSspots' for the Gray-Scott equations (spots),
%    - 'SCHNAK' for the Schnakenberg equations,
%    - 'SH' for the Swift-Hohenberg equation.
%
% Example 1: Ginzburg-Landau equation (spiral waves)
%
%       u = spin2('GL');
%
%    solves the Ginzburg-Landau equation
%
%       u_t = laplacian(u) + u - (1+1.5i)*u*|u|^2,
%
%    on [0 100]^2 from t=0 to t=100, with a RANDNFUN2 initial condition.
%    The movie shows the real part of u.  For a movie of the absolute
%    value of u rather than the real part, execute
%
%       S = spinop2('GL');
%       u = spin2(S, 128, 1e-1, 'dataplot', 'abs');
%
% Example 2: Gray-Scott equations (pattern formation - stripes)
%
%       u = spin2('GS');
%
%    solves the Gray-Scott equations
%
%       u_t = 2e-5*laplacian(u) + F*(1-u) - u*v^2,
%       v_t = 1e-5*laplacian(v) - (F+K)*v + u*v^2,
%
%    on [0 1]^2 from t=0 to t=5000, with initial condition
%
%       u0(x,y) = 1 - exp(-100*((x-1/2.05)^2 + (y-1/2.05)^2)),
%       v0(x,y) = exp(-100*((x-1/2)^2 + 2*(y-1/2)^2)),
%
%    with F=0.030 and K=0.057.
%
% Example 2 (bis): Gray-Scott equations (pattern formation - spots)
%
%       u = spin2('GSspots');
%
%    solves the Gray-Scott equations
%
%       u_t = 2e-5*laplacian(u) + F*(1-u) - u*v^2,
%       v_t = 1e-5*laplacian(v) - (F+K)*v + u*v^2,
%
%    on [0 1]^2 from t=0 to t=5000, with initial condition
%
%       u0(x,y) = 1 - exp(-100*((x-1/2.05)^2 + (y-1/2.05)^2)),
%       v0(x,y) = exp(-100*((x-1/2)^2 + 2*(y-1/2)^2)),
%
%    with F=0.026 and K=0.059.
%
% Example 3: Schnakenberg equations (pattern formation - spots)
%
%       u = spin2('SCHNAK');
%
%    solves the Schnakenberg equations
%
%       u_t = laplacian(u) + 3*(a - u + u^2*v),
%       v_t = 10*laplacian(v) + 3*(b - u^2*v),
%
%    on [0 50]^2 from t=0 to t=500, with initial condition
%
%       u0(x,y) = (a+b) - exp(-2*((x-G/2.15)^2 + (y-G/2.15)^2)),
%       v0(x,y) = b/(a+b)^2 + exp(-2*((x-G/2)^2 + 2*(y-G/2)^2)),
%           with G=50, a=0.1 and b=0.9.
%
% Example 4: Swift-Hohenberg equation (Rayleigh-Benard convection rolls)
%
%       u = spin2('SH');
%
%    solves the Swift-Hohenberg equation
%
%       u_t = -2*laplacian(u) - biharmonic(u) - .9*u - u^3,
%
%    on [0 50]^2 from t=0 to t=800, with a RANDNFUN2 initial condition.
%
% Example 5: PDE specified by a SPINOP2
%
%       dom = [0 100 0 100]; tspan = [0 100];
%       S = spinop2(dom, tspan);
%       S.lin = @(u) lap(u);
%       S.nonlin = @(u) u - (1 + 1.5i)*u.*(abs(u).^2);
%       S.init = randnfun2(4, dom, 'trig');
%       S.init = S.init/norm(S.init, inf);
%       u = spin2(S, 128, 1e-1);
%
%   is equivalent to u = spin2('GL');
%
% Example 6: Using preferences
%
%       pref = spinpref2('plot', 'off', 'scheme', 'pecec433');
%       S = spinop2('SH');
%       u = spin2(S, 128, 5e-1, pref);
%   or simply,
%       u = spin2(S, 128, 5e-1, 'plot', 'off', 'scheme', 'pecec433');
%
%   solves the Swift-Hohenberg equation using N=128 grid points in each
%   direction, a time-step dt=5e-1, doesn't play any movie and uses the
%   time-stepping scheme PECEC433.
%
% See also SPINOP2, SPINPREF2, EXPINTEG.

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

% We are going to parse the inputs and call SOLVEPDE in the following ways,
%
%       SPINOPERATOR.SOLVEPDE(S, N, dt)
%  or
%       SPINOPERATOR.SOLVEPDE(S, N, dt, pref)
%
% where S is a SPINOP2 object, N is the number of grid points in each direction,
% DT is the time-step and PREF is a SPINPREF2 oject.

% CASE 1. U = SPIN2('GL'):
if ( nargin == 1 )

    try spinop2(varargin{1});
    catch
        error('Unrecognized PDE. See HELP/SPINSPHERE for the list of PDEs.')
    end
    [S, N, dt, pref] = parseInputs(varargin{1});
    varargin{1} = S;
    varargin{2} = N;
    varargin{3} = dt;
    varargin{4} = pref;

% CASE 2. U = SPIN2('GL', 'PREF1', VALUE1) or U = SPINSPHERE(S, N, DT):
elseif ( nargin == 3 )

    % CASE 2.1. U = SPIN2('GL', 'PREF1', VALUE1):
    if ( isa(varargin{1}, 'char') == 1 && isa(varargin{2}, 'char') == 1 )
        [S, N, dt, pref] = parseInputs(varargin{1});
        pref.(varargin{2}) = varargin{3};
        varargin{1} = S;
        varargin{2} = N;
        varargin{3} = dt;
        varargin{4} = pref;

    % CASE 2.2. U = SPIN2(S, N, DT):
    else
        % Nothing to do here.
    end

% CASE 3. U = SPIN2(S, N, DT, PREF)
elseif ( nargin == 4 )
    % Nothing to do here.

% CASE 4.
elseif ( nargin >= 5 )

    % CASE 4.1. U = SPIN2('GL', 'PREF1', VALUE1, 'PREF2', VALUE2, ...)
    if ( isa(varargin{1}, 'char') == 1 && isa(varargin{2}, 'char') == 1 )
        [S, N, dt, pref] = parseInputs(varargin{1});
        j = 2;
        while j < nargin
            pref.(varargin{j}) = varargin{j+1};
            varargin{j} = [];
            varargin{j+1} = [];
            j = j + 2;
        end
        varargin{1} = S;
        varargin{2} = N;
        varargin{3} = dt;
        varargin{4} = pref;
        varargin = varargin(~cellfun(@isempty, varargin));

    % CASE 4.2. U = SPIN2(S, N, DT, 'PREF1', VALUE1, 'PREF2', VALUE2, ...)
    else
        pref = spinpref2();
        j = 4;
        while j < nargin
            pref.(varargin{j}) = varargin{j+1};
            varargin{j} = [];
            varargin{j+1} = [];
            j = j + 2;
        end
        varargin{4} = pref;
        varargin = varargin(~cellfun(@isempty, varargin));
    end

end

% SPIN2 is a wrapper for SOLVPDE:
[uout, tout] = spinoperator.solvepde(varargin{:});

end

function [S, N, dt, pref] = parseInputs(pdechar)
%PARSEINPUTS   Parse the inputs.

pref = spinpref2();
S = spinop2(pdechar);
if ( strcmpi(pdechar, 'GL') == 1 )
    dt = 1e-1;
    N = 128;
    pref.Clim = [-1 1];
    pref.iterplot = 2;
    pref.Nplot = 256;
elseif ( strcmpi(pdechar, 'GS') == 1 )
    dt = 5;
    N = 64;
    pref.Clim = [.25 .8 0 .35];
    pref.iterplot = 10;
    pref.Nplot = 128;
elseif ( strcmpi(pdechar, 'GSspots') == 1 )
    dt = 2;
    N = 64;
    pref.Clim = [.25 .85 0 .35];
    pref.iterplot = 25;
    pref.Nplot = 128;
elseif ( strcmpi(pdechar, 'SCHNAK') == 1 )
    dt = 5e-1;
    N = 64;
    pref.Clim = [.7 1.7 .65 1.05];
    pref.iterplot = 10;
    pref.Nplot = 128;
elseif ( strcmpi(pdechar, 'SH') == 1 )
    dt = 1;
    N = 128;
    pref.Clim = [-.4 .5];
    pref.iterplot = 4;
    pref.Nplot = 256;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [E, D] = pcamat(vectors, firstEig, lastEig, s_interactive, ...
    s_verbose);
%PCAMAT - Calculates the pca for data
%
% [E, D] = pcamat(vectors, firstEig, lastEig, ...
%                 interactive, verbose);
%
% Calculates the PCA matrices for given data (row) vectors. Returns
% the eigenvector (E) and diagonal eigenvalue (D) matrices containing the
% selected subspaces. Dimensionality reduction is controlled with
% the parameters 'firstEig' and 'lastEig' - but it can also be done
% interactively by setting parameter 'interactive' to 'on' or 'gui'.
%
% ARGUMENTS
%
% vectors       Data in row vectors.
% firstEig      Index of the largest eigenvalue to keep.
%               Default is 1.
% lastEig       Index of the smallest eigenvalue to keep.
%               Default is equal to dimension of vectors.
% interactive   Specify eigenvalues to keep interactively. Note that if
%               you set 'interactive' to 'on' or 'gui' then the values
%               for 'firstEig' and 'lastEig' will be ignored, but they
%               still have to be entered. If the value is 'gui' then the
%               same graphical user interface as in FASTICAG will be
%               used. Default is 'off'.
% verbose       Default is 'on'.
%
%
% EXAMPLE
%       [E, D] = pcamat(vectors);
%
% Note
%       The eigenvalues and eigenvectors returned by PCAMAT are not sorted.
%
% This function is needed by FASTICA and FASTICAG

% For historical reasons this version does not sort the eigenvalues or
% the eigen vectors in any ways. Therefore neither does the FASTICA or
% FASTICAG. Generally it seams that the components returned from
% whitening is almost in reversed order. (That means, they usually are,
% but sometime they are not - depends on the EIG-command of matlab.)

% @(#)$Id: pcamat.m,v 1.5 2003/12/15 18:24:32 jarmo Exp $

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Default values:
if nargin < 5, s_verbose = 'on'; end
if nargin < 4, s_interactive = 'off'; end
if nargin < 3, lastEig = size(vectors, 1); end
if nargin < 2, firstEig = 1; end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check the optional parameters;
switch lower(s_verbose)
 case 'on'
  b_verbose = 1;
 case 'off'
  b_verbose = 0;
 otherwise
  error(sprintf('Illegal value [ %s ] for parameter: ''verbose''\n', s_verbose));
end

switch lower(s_interactive)
 case 'on'
  b_interactive = 1;
 case 'off'
  b_interactive = 0;
 case 'gui'
  b_interactive = 2;
 otherwise
  error(sprintf('Illegal value [ %s ] for parameter: ''interactive''\n', ...
		s_interactive));
end

oldDimension = size (vectors, 1);
if ~(b_interactive)
  if lastEig < 1 | lastEig > oldDimension
    error(sprintf('Illegal value [ %d ] for parameter: ''lastEig''\n', lastEig));
  end
  if firstEig < 1 | firstEig > lastEig
    error(sprintf('Illegal value [ %d ] for parameter: ''firstEig''\n', firstEig));
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate PCA

% Calculate the covariance matrix.
if b_verbose, fprintf ('Calculating covariance...\n'); end
covarianceMatrix = cov(vectors', 1);

% Calculate the eigenvalues and eigenvectors of covariance
% matrix.
[E, D] = eig (covarianceMatrix);

% The rank is determined from the eigenvalues - and not directly by
% using the function rank - because function rank uses svd, which
% in some cases gives a higher dimensionality than what can be used
% with eig later on (eig then gives negative eigenvalues).
rankTolerance = 1e-7;
maxLastEig = sum (diag (D) > rankTolerance);
if maxLastEig == 0,
  fprintf (['Eigenvalues of the covariance matrix are' ...
	    ' all smaller than tolerance [ %g ].\n' ...
	    'Please make sure that your data matrix contains' ...
	    ' nonzero values.\nIf the values are very small,' ...
	    ' try rescaling the data matrix.\n'], rankTolerance);
  error ('Unable to continue, aborting.');
end

% Sort the eigenvalues - decending.
eigenvalues = flipud(sort(diag(D)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Interactive part - command-line
if b_interactive == 1

  % Show the eigenvalues to the user
  hndl_win=figure;
  bar(eigenvalues);
  title('Eigenvalues');

  % ask the range from the user...
  % ... and keep on asking until the range is valid :-)
  areValuesOK=0;
  while areValuesOK == 0
    firstEig = input('The index of the largest eigenvalue to keep? (1) ');
    lastEig = input(['The index of the smallest eigenvalue to keep? (' ...
                    int2str(oldDimension) ') ']);
    % Check the new values...
    % if they are empty then use default values
    if isempty(firstEig), firstEig = 1;end
    if isempty(lastEig), lastEig = oldDimension;end
    % Check that the entered values are within the range
    areValuesOK = 1;
    if lastEig < 1 | lastEig > oldDimension
      fprintf('Illegal number for the last eigenvalue.\n');
      areValuesOK = 0;
    end
    if firstEig < 1 | firstEig > lastEig
      fprintf('Illegal number for the first eigenvalue.\n');
      areValuesOK = 0;
    end
  end
  % close the window
  close(hndl_win);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Interactive part - GUI
if b_interactive == 2

  % Show the eigenvalues to the user
  hndl_win = figure('Color',[0.8 0.8 0.8], ...
    'PaperType','a4letter', ...
    'Units', 'normalized', ...
    'Name', 'FastICA: Reduce dimension', ...
    'NumberTitle','off', ...
    'Tag', 'f_eig');
  h_frame = uicontrol('Parent', hndl_win, ...
    'BackgroundColor',[0.701961 0.701961 0.701961], ...
    'Units', 'normalized', ...
    'Position',[0.13 0.05 0.775 0.17], ...
    'Style','frame', ...
    'Tag','f_frame');

b = uicontrol('Parent',hndl_win, ...
	'Units','normalized', ...
	'BackgroundColor',[0.701961 0.701961 0.701961], ...
	'HorizontalAlignment','left', ...
	'Position',[0.142415 0.0949436 0.712077 0.108507], ...
	'String','Give the indices of the largest and smallest eigenvalues of the covariance matrix to be included in the reduced data.', ...
	'Style','text', ...
	'Tag','StaticText1');
e_first = uicontrol('Parent',hndl_win, ...
	'Units','normalized', ...
	'Callback',[ ...
          'f=round(str2num(get(gcbo, ''String'')));' ...
          'if (f < 1), f=1; end;' ...
          'l=str2num(get(findobj(''Tag'',''e_last''), ''String''));' ...
          'if (f > l), f=l; end;' ...
          'set(gcbo, ''String'', int2str(f));' ...
          ], ...
	'BackgroundColor',[1 1 1], ...
	'HorizontalAlignment','right', ...
	'Position',[0.284831 0.0678168 0.12207 0.0542535], ...
	'Style','edit', ...
        'String', '1', ...
	'Tag','e_first');
b = uicontrol('Parent',hndl_win, ...
	'Units','normalized', ...
	'BackgroundColor',[0.701961 0.701961 0.701961], ...
	'HorizontalAlignment','left', ...
	'Position',[0.142415 0.0678168 0.12207 0.0542535], ...
	'String','Range from', ...
	'Style','text', ...
	'Tag','StaticText2');
e_last = uicontrol('Parent',hndl_win, ...
	'Units','normalized', ...
	'Callback',[ ...
          'l=round(str2num(get(gcbo, ''String'')));' ...
          'lmax = get(gcbo, ''UserData'');' ...
          'if (l > lmax), l=lmax; fprintf([''The selected value was too large, or the selected eigenvalues were close to zero\n'']); end;' ...
          'f=str2num(get(findobj(''Tag'',''e_first''), ''String''));' ...
          'if (l < f), l=f; end;' ...
          'set(gcbo, ''String'', int2str(l));' ...
          ], ...
	'BackgroundColor',[1 1 1], ...
	'HorizontalAlignment','right', ...
	'Position',[0.467936 0.0678168 0.12207 0.0542535], ...
	'Style','edit', ...
        'String', int2str(maxLastEig), ...
        'UserData', maxLastEig, ...
	'Tag','e_last');
% in the first version oldDimension was used instead of
% maxLastEig, but since the program would automatically
% drop the eigenvalues afte maxLastEig...
b = uicontrol('Parent',hndl_win, ...
	'Units','normalized', ...
	'BackgroundColor',[0.701961 0.701961 0.701961], ...
	'HorizontalAlignment','left', ...
	'Position',[0.427246 0.0678168 0.0406901 0.0542535], ...
	'String','to', ...
	'Style','text', ...
	'Tag','StaticText3');
b = uicontrol('Parent',hndl_win, ...
	'Units','normalized', ...
	'Callback','uiresume(gcbf)', ...
	'Position',[0.630697 0.0678168 0.12207 0.0542535], ...
	'String','OK', ...
	'Tag','Pushbutton1');
b = uicontrol('Parent',hndl_win, ...
	'Units','normalized', ...
	'Callback',[ ...
          'gui_help(''pcamat'');' ...
          ], ...
	'Position',[0.767008 0.0678168 0.12207 0.0542535], ...
	'String','Help', ...
	'Tag','Pushbutton2');

  h_axes = axes('Position' ,[0.13 0.3 0.775 0.6]);
  set(hndl_win, 'currentaxes',h_axes);
  bar(eigenvalues);
  title('Eigenvalues');

  uiwait(hndl_win);
  firstEig = str2num(get(e_first, 'String'));
  lastEig = str2num(get(e_last, 'String'));

  % close the window
  close(hndl_win);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% See if the user has reduced the dimension enought

if lastEig > maxLastEig
  lastEig = maxLastEig;
  if b_verbose
    fprintf('Dimension reduced to %d due to the singularity of covariance matrix\n',...
           lastEig-firstEig+1);
  end
else
  % Reduce the dimensionality of the problem.
  if b_verbose
    if oldDimension == (lastEig - firstEig + 1)
      fprintf ('Dimension not reduced.\n');
    else
      fprintf ('Reducing dimension...\n');
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Drop the smaller eigenvalues
if lastEig < oldDimension
  lowerLimitValue = (eigenvalues(lastEig) + eigenvalues(lastEig + 1)) / 2;
else
  lowerLimitValue = eigenvalues(oldDimension) - 1;
end

lowerColumns = diag(D) > lowerLimitValue;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Drop the larger eigenvalues
if firstEig > 1
  higherLimitValue = (eigenvalues(firstEig - 1) + eigenvalues(firstEig)) / 2;
else
  higherLimitValue = eigenvalues(1) + 1;
end
higherColumns = diag(D) < higherLimitValue;

% Combine the results from above
selectedColumns = lowerColumns & higherColumns;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% print some info for the user
if b_verbose
  fprintf ('Selected [ %d ] dimensions.\n', sum (selectedColumns));
end
if sum (selectedColumns) ~= (lastEig - firstEig + 1),
  error ('Selected a wrong number of dimensions.');
end

if b_verbose
  fprintf ('Smallest remaining (non-zero) eigenvalue [ %g ]\n', eigenvalues(lastEig));
  fprintf ('Largest remaining (non-zero) eigenvalue [ %g ]\n', eigenvalues(firstEig));
  fprintf ('Sum of removed eigenvalues [ %g ]\n', sum(diag(D) .* ...
    (~selectedColumns)));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Select the colums which correspond to the desired range
% of eigenvalues.
E = selcol(E, selectedColumns);
D = selcol(selcol(D, selectedColumns)', selectedColumns);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Some more information
if b_verbose
  sumAll=sum(eigenvalues);
  sumUsed=sum(diag(D));
  retained = (sumUsed / sumAll) * 100;
  fprintf('[ %g ] %% of (non-zero) eigenvalues retained.\n', retained);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function newMatrix = selcol(oldMatrix, maskVector);

% newMatrix = selcol(oldMatrix, maskVector);
%
% Selects the columns of the matrix that marked by one in the given vector.
% The maskVector is a column vector.

% 15.3.1998

if size(maskVector, 1) ~= size(oldMatrix, 2),
  error ('The mask vector and matrix are of uncompatible size.');
end

numTaken = 0;

for i = 1 : size (maskVector, 1),
  if maskVector(i, 1) == 1,
    takingMask(1, numTaken + 1) = i;
    numTaken = numTaken + 1;
  end
end

newMatrix = oldMatrix(:, takingMask);

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
function line_fekete_rule_test02 ( m )

%*****************************************************************************80
%
%% LINE_FEKETE_RULE_TEST02 seeks Fekete points in [-1,+1].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 March 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    L Bos, N Levenberg,
%    On the calculation of approximate Fekete points: the univariate case,
%    Electronic Transactions on Numerical Analysis,
%    Volume 30, pages 377-397, 2008.
%
%  Parameters:
%
%    Input, integer M, the dimension of the polynomial space.
%
  if ( nargin < 1 )
    m = 5;
  end
  n = 1000;
  a = -1.0;
  b = +1.0;
  x = linspace ( a, b, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'LINE_FEKETE_RULE_TEST02:\n' );
  fprintf ( 1, '  Seek Fekete points in [%g,%g]\n', a, b );
  fprintf ( 1, '  using %d equally spaced sample points\n', n );
  fprintf ( 1, '  for polynomials of degree M = %d\n', m );
  fprintf ( 1, '  with the Chebyshev basis\n' );
  fprintf ( 1, '  and weight 1/sqrt(1-x^2).\n' );

  [ nf, xf, wf, vf ] = line_fekete_chebyshev ( m, a, b, n, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  NF = %d\n', nf );
  r8vec_print ( nf, xf, '  Estimated Fekete points XF:' );

  yf = ones ( size ( xf ) );
  plot ( xf, yf, '*' );
  title ( 'Estimated Fekete points' );
  grid on
  filename = 'line_fekete_rule_test02.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saved plotfile as "%s"\n', filename );

  wf_sum = sum ( wf(1:nf) );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sum(WF) = %g\n', wf_sum );

  return
end

author: johannesgerer


text:
%%%% calculation of field spatial distributions, starting from the eigenvectors (Fourier coefficients)
%%%% omega = normalized frequency
%%%% eta = inverse of matrix containing the Fourier coeff. of dielectric
%%%% function
%%%% Phi = matrix of column eigenvectors [hx,hy]'
%%%% kz = vector of real eigenvalues
%%%% u = index of selected eigenvalue/eigenvector
function [Ex,Ey,Ez,Hx,Hy,Hz] = rfields(omega,eta,kGx,kGy,kz,Phi,N1,N2,u)
N=N1*N2;
hx=Phi(1:N,u); hy=Phi((N+1):2*N,u);
hz=-(1/kz(u))*(kGx*hx+kGy*hy);
ex=-(1/omega)*eta*(kGy*hz-kz(u)*hy);
ey=-(1/omega)*eta*(kz(u)*hx-kGx*hz);
ez=-(1/omega)*eta*(kGx*hy-kGy*hx);

ex=reshape(ex,N1,N2); ey=reshape(ey,N1,N2); ez=reshape(ez,N1,N2);
hx=reshape(hx,N1,N2); hy=reshape(hy,N1,N2); hz=reshape(hz,N1,N2);
Ex=(N1*N2)*ifft2(ifftshift(ex'));
Ey=(N1*N2)*ifft2(ifftshift(ey'));
Ez=(N1*N2)*ifft2(ifftshift(ez'));

Hx=(N1*N2)*ifft2(ifftshift(hx'));
Hy=(N1*N2)*ifft2(ifftshift(hy'));
Hz=(N1*N2)*ifft2(ifftshift(hz'));

author: Sable


text:
function z = EulerToQuaternion(phi, theta, psi)
%
%
sinPhi   = sin(phi/2);    cosPhi   = cos(phi/2);
sinTheta = sin(theta/2);  cosTheta = cos(theta/2);
sinPsi   = sin(psi/2);    cosPsi   = cos(psi/2);

z = [ cosPhi*cosTheta*cosPsi + sinPhi*sinTheta*sinPsi;
      sinPhi*cosTheta*cosPsi - cosPhi*sinTheta*sinPsi;
      cosPhi*sinTheta*cosPsi + sinPhi*cosTheta*sinPsi;
      cosPhi*cosTheta*sinPsi - sinPhi*sinTheta*cosPsi;
    ];

author: philbooks

repo: Kalman-Filter-for-Beginners


text:
clear; close all; clc;

Data = load('data.txt');
x = Data(:, [1, 2]); Y = Data(:, 3);

fprintf(['Plotting data with + indicating (Y = 1) examples and o ' ...
         'indicating (Y = 0) examples.\n']);

plotdata(x, Y);

hold on;

xlabel('Exam 1 score');
ylabel('Exam 2 score');

legend('Admitted', 'Not admitted');
hold off;

fprintf('\nProgram paused, press enter to continue.\n');
pause;

[m, n] = size(x);

X = [ones(m, 1) x];

Initial_Theta = zeros(n + 1, 1);

[Cost, Grad] = costfunction(Initial_Theta, X, Y);

fprintf('Cost at initial theta (zeros): %f\n', Cost);
fprintf('Expected cost (approx): 0.693\n');
fprintf('Gradient at initial theta (zeros): \n');
fprintf(' %f \n', Grad);
fprintf('Expected gradients (approx):\n -0.1000\n -12.0092\n -11.2628\n');

Test_Theta = [-24; 0.2; 0.2];
[Cost, Grad] = costfunction(Test_Theta, X, Y);

fprintf('\nCost at test theta: %f\n', Cost);
fprintf('Expected cost (approx): 0.218\n');
fprintf('Gradient at test theta: \n');
fprintf(' %f \n', Grad);
fprintf('Expected gradients (approx):\n 0.043\n 2.566\n 2.647\n');

fprintf('\nProgram paused, press enter to continue.\n');
pause;

Options = optimset('GradObj', 'on', 'MaxIter', 400);

[Theta, Cost] = ...
	fminunc(@(t)(costfunction(t, X, Y)), Initial_Theta, Options);

fprintf('Cost at theta found by fminunc: %f\n', Cost);
fprintf('Expected cost (approx): 0.203\n');
fprintf('theta: \n');
fprintf(' %f \n', Theta);
fprintf('Expected theta (approx):\n');
fprintf(' -25.161\n 0.206\n 0.201\n');

plotdecisionboundary(Theta, X, Y);

hold on;

xlabel('Exam 1 score');
ylabel('Exam 2 score');

legend('Admitted', 'Not admitted');
hold off;

fprintf('\nProgram paused, press enter to continue.\n');
pause;

Prob = Sigmoid([1 45 85] * Theta);
fprintf(['For a student with scores 45 and 85, we predict an admission ' ...
         'probability of %f\n'], Prob);
fprintf('Expected value: 0.775 +/- 0.002\n\n');

P = Predict(Theta, X);

fprintf('Train accuracy: %f\n', mean(double(P == Y)) * 100);
fprintf('Expected accuracy (approx): 89.0\n');
fprintf('\n');

author: TheAlgorithms

repo: MATLAB-Octave