text:
function value = r8_gmit ( a, x, algap1, sgngam, alx )

%*****************************************************************************80
%
%% R8_GMIT: Tricomi's incomplete gamma function for small X.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real A, the parameter.
%
%    Input, real X, the argument.
%
%    Input, real ALGAP1, the logarithm of Gamma ( A + 1 ).
%
%    Input, real SGNGAM, the sign of Gamma ( A + 1 ).
%
%    Input, real ALX, the logarithm of X.
%
%    Output, real VALUE, the Tricomi incomplete gamma function.
%
  persistent bot
  persistent eps

  if ( isempty ( eps ) )
    eps = 0.5 * r8_mach ( 3 );
    bot = log ( r8_mach ( 1 ) );
  end

  if ( x <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_GMIT - Fatal error!\n' );
    fprintf ( 1, '  X <= 0.\n' );
    error ( 'R8_GMIT - Fatal error!' )
  end

  if ( a < 0.0 )
    ma = r8_aint ( a - 0.5 );
  else
    ma = r8_aint ( a + 0.5 );
  end

  aeps = a - ma;

  if ( a < - 0.5 )
    ae = aeps;
  else
    ae = a;
  end

  t = 1.0;
  te = ae;
  s = t;
  converged = 0;
  for k = 1 : 200
    fk = k;
    te = - x * te / fk;
    t = te / ( ae + fk );
    s = s + t;
    if ( abs ( t ) < eps * abs ( s ) )
      converged = 1;
      break
    end
  end

  if ( ~ converged )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_GMIT - Fatal error!\n' );
    fprintf ( 1, '  No convergence in 200 iterations.\n' );
    error ( 'R8_GMIT - Fatal error!' )
  end

  if ( - 0.5 <= a )
    algs = - algap1 + log ( s );
    value = exp ( algs );
    return
  end

  algs = - r8_lngam ( 1.0 + aeps ) + log ( s );
  s = 1.0;
  m = - ma - 1;
  t = 1.0;
  for k = 1 : m
    t = x * t / ( aeps - ( m + 1 - k ) );
    s = s + t;
    if ( abs ( t ) < eps * abs ( s ) )
      break
    end
  end

  value = 0.0;
  algs = - ma * log ( x ) + algs;

  if ( s == 0.0 || aeps == 0.0 )
    value = exp ( algs );
    return
  end

  sgng2 = sgngam * r8_sign ( s );
  alg2 = - x - algap1 + log ( abs ( s ) );

  if ( bot < alg2 )
    value = sgng2 * exp ( alg2 );
  end

  if ( bot < algs )
    value = value + exp ( algs );
  end

  return
end


text:
function value = f_1_2d ( x, y )

%*****************************************************************************80
%
%% F_1_2D evaluates the function 1 in 2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the arguments.
%
%    Output, real F_1_2D, the value of the function.
%
  value = 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [X,rho,eta] = mr2(A,b,k,reorth)
%MR2 Solution of symmetric indefinite problems by the MR-II algorithm
%
% [X,rho,eta] = mr2(A,b,k,reorth)
%
% MR-II is a variant of the MINRES algorithm for symmetric indefinite linear
% systems A x = b, with starting vector A*b (instead of b as in MINRES).
% This function returns all k iterates, stored as the columns of the
% matrix X.  The solution norm and residual norm are returned in eta and
% rho, respectively.
%
% Reorthogonalization is controlled by means of reorth:
%    reorth = 0 : no reorthogonalization (default),
%    reorth = 1 : reorthogonalization by means of MGS.

% Reference: M. Hanke, "Conjugate Gradient Methods for Ill-Posed Problems",
% Longman Scientific and Technical, Essex, 1995.

% Per Christian Hansen, IMM, September 1, 2007.
% Based on the function mr2 from Restore Tools by James G. Nagy.

% Initialization.
if (k < 1), error('Number of steps k must be positive'), end
if (nargin==3), reorth = 0; end
[m,n] = size(A);
if (m ~= n || norm(A-A','fro')), error('The matrix must be symmetric'), end

% Allocate space.
X = zeros(n,k);
if reorth
  W = zeros(n,k);
  if (k>=n), error('No. of iterations must satisfy k < n'), end
end
if (nargout > 1)
  eta = zeros(k,1); rho = eta;
end

% Prepare for interation.
x = zeros(n,1); r = b;
vold = 0; v = A*r;
wold = 0; w = A*v;
beta = norm(w);
v = v./beta; w = w./beta;
if reorth, W(:,1) = w; end

% Perform k iterations.
for i=1:k

  rrho = r'*w;
  x = x + rrho*v;
  r = r - rrho*w;
  Aw = A*w;
  alpha = w'*Aw;
  vnew  =  w - alpha*v - beta*vold;
  wnew  = Aw - alpha*w - beta*wold;
  vold = v; wold = w; v = vnew; w = wnew;
  if reorth
    for j=1:i, w = w - (W(:,j)'*w)*W(:,j); end
  end;
  beta = norm(w);
  v = v./beta;	w = w./beta;
  if reorth, W(:,i+1) = w; end;

  X(:,i) = x;
  if (nargout>1), rho(i) = norm(r); end
  if (nargout>2), eta(i) = norm(x); end

end

repo: REST


text:
function h = vl_hikmeanshist(tree,path)
% VL_HIKMEANSHIST  Compute histogram of quantized data
%  H = VL_HIKMEANSHIST(TREE,PATH) computes the histogram of the HIKM tree
%  nodes activated by the root-to-leaf paths PATH. PATH is usually
%  obtained by quantizing data by means of VL_HIKMEANSPUSH().
%
%  The histogram H has one bin for each node of the HIKM tree TREE.
%  The tree has K = TREE.K nodes and depth D = TREE.DEPTH.  Therefore
%  there are M = (K^(D+1) - 1) / (K - 1) nodes in the tree (not
%  counting the root which carries no information). Nodes are stacked
%  into a vector of bins in breadth first order.
%
%  Example::
%    The following relations illustrate the structure of PATH:
%      H(1)   = # of paths such that PATH(1,:) = 1
%      H(K)   = # of paths such that PATH(1,:) = K
%      H(K+1) = # of paths such that PATH(1:2,:) = [1 ; 1]
%      H(K+K) = # of paths such that PATH(1:2,:) = [1 ; K]
%
%  See also: VL_HIKMEANS(), VL_HIKMEANSPUSH(), VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

% PATH(:,k) is colum of subscripts i1 i2 ... id identifying a path in
% the tree. In a bread first enumeration of the tree nodes (starting
% from one and not counting the root), the node of subscripts
% i1,i2,...id has index
%
% idx = i1 K^{d-1} + i2 K^{d-2} + ... + id
%
% where we assumed the indeces i1,i2,... start from 1. This formula
% can be easily computed recursively. Since we also have a root
% node, we need to add one.

K = tree.K ;
D = tree.depth ;
M = (K^(D+1) - 1) / (K - 1) ;

h = zeros(M, 1) ;
p = zeros(1,size(path,2)) ;

h(1) = size(path,2) ;

for d=1:D
  p = p * K + double(path(d,:))  ;
  h = vl_binsum(h, 1, p + 1) ;
end

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
% 1D Simple Bar
% clear memory
clear all

% E: modulus of elasticity
% A: area of cross section
% L: length of bar
E=30e6; A=1; L=[30 30 30];

% numberElements: number of elements
numberElements=3;
% numberNodes: number of nodes
numberNodes=4;
% generation of coordinates and connectivities
elementNodes=[1 2;2 3;3 4];
nodeCoordinates=[0 30 60 90];
% for structure:
   % displacements: displacement vector
   % force : force vector
   % stiffness: stiffness matrix
force=zeros(numberNodes,1);
stiffness=zeros(numberNodes,numberNodes);
% applied load at node 2
force(2)=3000.0;

% computation of the system stiffness matrix
for e=1:numberElements;
  % elementDof: element degrees of freedom (Dof)
  elementDof=elementNodes(e,:) ;
  k(e)=E*A/L(e) ;
  stiffness(elementDof,elementDof)=...
       stiffness(elementDof,elementDof)+k(e)*[1 -1;-1 1];
end
% boundary conditions and solution
% prescribed dofs
prescribedDof=[1;4];
% solution
GDof=numberNodes;
displacements=solution(GDof,prescribedDof,stiffness,force);
% output displacements/reactions
outputDisplacementsReactions(displacements,stiffness, ...
    numberNodes,prescribedDof,force)

repo: Aero-matlab


text:
function [fidout] = improvfid(ecg,fidin,fs,fidtype,fidhelp,N2)
%
% OVERVIEW, This function improves the fiducial point detection for the
% Q,R,S,Toff points in the beat variable. The improved fiducial points are
% stored in the beato variable.
%
% INPUTS    MANDATORY       DESCRIPTION
%
%           ecg             input ecg data in num of samples by 1
%                           dimensions
%           fidin           fiducial point input in samples number.
%           fs              sampling frequency of the input ecg data.
%           fidtype         type of fiducial point input. i.e. Q, R, S or T.
%           fidhelp         R peak annotations used for improving S and Q point detections.
%           N2              window size of 2*N2 is used to find improved fiducial point detection for fiducial point specified as fidtype.
%
% OUTPUTS
%           fidout          Adjusted fiducial point location after
%                           improving fiducial point detection.
%
%
% REPO:
%       https://github.com/cliffordlab/PhysioNet-Cardiovascular-Signal-Toolbox
%   ORIGINAL SOURCE AND AUTHORS:
%       Written by Shamim Nemati
%       editted by Ismail Sadiq on 10/26/2019.
%	COPYRIGHT (C) 2019
%   LICENSE:
%       This software is offered freely and without warranty under
%       the GNU (v3 or later) public license. See license file for
%       more information. The license may be found in
%       the Documents folder of the Physionet-Cardiovascular-Signal-Toolbox.

fidout=zeros(1,size(fidin,2));

if (nargin<5 || isempty(fidhelp))
    Rp=[]; Sp=[];
elseif (strcmp(fidtype,'S'))
    Rp=fidhelp; Sp=[];
elseif (strcmp(fidtype,'T'))
    Sp = fidhelp; Rp=[];
else
    Rp=[]; Sp=[];
end

if (nargin<4)
    direc=mean(ecg(fidin(:,1)))-mean(ecg);
else
    if (strcmp(fidtype,'R'))
        direc = 1;
    elseif (strcmp(fidtype,'T') & ~isempty(fidhelp))
        direc = median(ecg(fidin))-median(ecg(Sp));
    else
        direc=-1;
    end
end

if (nargin<6)
    N2 = round(0.04*fs);
end

for i=1:size(fidout,2)
    if (~isempty(Rp) & fidin(i)>Rp(i))%S
        upper = min(fidin(i)+N2,Rp(i)+2*N2+10);
        %upper=fidin(i)+N2;
        v=ecg(Rp(i):upper)';% if there's a R right before use it
    else
        v=ecg(fidin(i)-N2:fidin(i)+N2)';
    end
    vi=v;
    if direc>0
        [val,ind]=max(vi);
        fidout(i)=(fidin(i)-N2-1)+ind;
    else
        [val,ind]=min(vi);
        if (~isempty(Rp) & fidin(i)>Rp(i))%S
            fidout(i)=Rp(i)-1+ind;
        else
            fidout(i)=(fidin(i)-N2-1)+ind;
        end
    end
end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function cartPoints=pol2Cart(z,systemType,useHalfRange,zTx,zRx,M)
%%POL2CART Convert points from bistatic polar coordinates to 2D Cartesian
%          coordinates. The angle can measured either counterclockwise
%          from the x-axis, which is standard in mathematics, or clockwise
%          from the y axis, which is more common in navigation.
%
%INPUTS: z A 2XN matrix of points in polar coordinates. Each column of the
%          matrix has the format [r;azimuth], with azimuth given in
%          radians. Alternatively, if one just wishes to have unit vectors
%          returned, this can be a 1XN matrix of just azimuthal angles, in
%          which case the useHalfRange, zTx,and zRx inputs are ignored.
% systemType An optional parameter specifying the axis from which the
%          angles are measured. Possible values are:
%          0 (The default if omitted or an empty matrix is passed) The
%            azimuth angle is counterclockwise from the x axis.
%          1 The azimuth angle is measured clockwise from the y axis.
% useHalfRange A boolean value specifying whether the bistatic range value
%           should be divided by two. This normally comes up when operating
%           in monostatic mode, so that the range reported is a one-way
%           range. The default if this parameter is not provided, or an
%           empty matrix is passed, is true.
%       zTx The 2XN [x;y] location vectors of the transmitters in global
%           Cartesian coordinates. If this parameter is omitted or an
%           empty matrix is passed, then the transmitters are assumed to
%           be at the origin. If only a single vector is passed, then the
%           transmitter location is assumed the same for all of the target
%           states being converted. zTx can have more than 2 rows;
%           additional rows are ignored.
%       zRx The 2XN [x;y] location vectors of the receivers in Cartesian
%           coordinates.  If this parameter is omitted or an empty matrix
%           is passed, then the receivers are assumed to be at the origin.
%           If only a single vector is passed, then the receiver location
%           is assumed the same for all of the target states being
%           converted. zRx can have more than 2 rows; additional rows are
%           ignored.
%         M A 2X2XN hypermatrix of the rotation matrices to go from the
%           alignment of the global coordinate system to that at the
%           receiver. If omitted or an empty matrix is passed, then it is
%           assumed that the local coordinate system is aligned with the
%           global and M=eye(2) --the identity matrix is used. If only a
%           single 2X2 matrix is passed, then it is assumed to be the same
%           for all of the N conversions.
%
%OUTPUTS: cartPoints A 2XN or matrix of the points transformed into
%                    Cartesian coordinates. Each columns of cartPoints is
%                    of the format [x;y].
%
%The conversion utilizing bistatic polar measurements in 2D is similar to
%that using bistatic r-u-v measurements in 3D, which is discussed in [1].
%In both instances, one turns the direction components into a unit vector
%and the multiplied it by a one-way monostatic range that has to be
%computed.
%
%REFERENCES:
%[1] David F. Crouse , "Basic tracking using nonlinear 3D monostatic and
%    bistatic measurements," IEEE Aerospace and Electronic Systems
%    Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.
%
%February 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(isempty(z))
    cartPoints=zeros(2,0);
    return
end

N=size(z,2);

if(nargin<6||isempty(M))
    M=repmat(eye(2),[1,1,N]);
elseif(size(M,3)==1)
    M=repmat(M,[1,1,N]);
end

if(nargin<5||isempty(zRx))
    zRx=zeros(2,N);
elseif(size(zRx,2)==1)
    zRx=repmat(zRx,[1,N]);
end

if(nargin<4||isempty(zTx))
    zTx=zeros(2,N);
elseif(size(zTx,2)==1)
    zTx=repmat(zTx,[1,N]);
end

if(nargin<3||isempty(useHalfRange))
    useHalfRange=true;
end

if(nargin<2||isempty(systemType))
    systemType=0;
end

%Extract the components.
if(size(z,1)==1)
    %No ranges, so just get unit vectors.
    cartPoints=zeros(2,N);
    switch(systemType)
        case 0
            cartPoints(1,:)=cos(z);
            cartPoints(2,:)=sin(z);
        case 1
            cartPoints(1,:)=sin(z);
            cartPoints(2,:)=cos(z);
        otherwise
            error('Invalid system type specified.')
    end

    %Rotate the unit vectors into the global coordinate system.
    for k=1:N
        %Convert to global Cartesian coordinates.
        cartPoints(:,k)=M(:,:,k)\cartPoints(:,k);
    end
    return;
else
    rB=z(1,:);
    azimuth=z(2,:);
end

%The bistatic range is used in the conversions below.
if(useHalfRange)
   rB=2*rB;
end

%Get unit vectors from the azimuth.
u=zeros(2,N);
switch(systemType)
    case 0
        u(1,:)=cos(azimuth);
        u(2,:)=sin(azimuth);
    case 1
        u(1,:)=sin(azimuth);
        u(2,:)=cos(azimuth);
    otherwise
        error('Invalid system type specified.')
end

cartPoints=zeros(2,N);
for curPoint=1:N
    if(all(zRx(:,curPoint)==zTx(:,curPoint)))
        %If it is monostatic, set r1 directly. This avoids NaNs with
        %zero-ranges in the monostatic case.
        r1=rB(curPoint)/2;
        %This is the Cartesian location in the local coordinate system of
        %the receiver.
        zL=r1*u(:,curPoint);
        %Convert to global Cartesian coordinates.
        cartPoints(:,curPoint)=M(:,:,curPoint)\zL+zRx(1:2,curPoint);
    else
        %The transmitter location in the receiver's local coordinate
        %system.
        zTxL=M(:,:,curPoint)*(zTx(1:2,curPoint)-zRx(1:2,curPoint));

        r1=(rB(curPoint)^2-norm(zTxL)^2)/(2*(rB(curPoint)-dot(u(:,curPoint),zTxL)));
        %This is the Cartesian location in the local coordinate system of
        %the receiver.
        zL=r1*u(:,curPoint);
        %Convert to global Cartesian coordinates.
        cartPoints(:,curPoint)=M(:,:,curPoint)\zL+zRx(1:2,curPoint);
    end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [data,grd]=nj_tslice(ncRef,var,iTime,level)
% NJ_TSLICE - Get data and coordinates from CF-compliant dataset at specific time and level
% Usage:
%   [data,grd_temp]=nj_tslice(ncRef,var,[iTime],[level]);
% Inputs:
%   ncRef = OpenDAP Data URL, Local file, or existing 'ncgeodataset' object
%   var - variable to slice
%   iTime = time step to extract data. Use inf or -1 for last time step.
%   level = vertical level (if not supplied, volume data is retrieved). Use
%           inf or -1 for last level.
% Outputs:
%   data = data  - matlab array
%   grd_temp  = structure containing lon,lat,z,time (Matlab datenum)
% Example:
%   uri ='http://geoport.whoi.edu/thredds/dodsC/examples/bora_feb.nc';% NetCDF/OpenDAP/NcML file
%   [data,grd_temp]=nj_tslice(uri,'temp',2, 14); % get level 14 at time step 2
%   [data,grd_temp]=nj_tslice(uri,'temp',2); % get 3D data at time step 2
%   [data,grd_temp]=nj_tslice(uri,'h'); % Retrieve 2D non time dependent array
%
% NCTOOLBOX (https://github.com/nctoolbox/nctoolbox)

if nargin < 2, help(mfilename), return, end
try
  if (isa(ncRef, 'ncgeodataset')) %check for ncgeodataset object
    nc = ncRef;
  else
    % open CF-compliant dataset as ncgeodataset object
    nc = ncgeodataset(ncRef);
  end
  gvar = nc.geovariable(var); %get geovariable object
  switch nargin
    case 2  % return 2D or 3D data (z,lat,lon) for all times (could be huge)
      data = squeeze(gvar.data(:,:,:,:));
      if nargout==2,
        grd_temp = gvar.grid_interop(:,:,:,:);
        grd_temp.z = squeeze(grd_temp.z);
      end
    case 3  % return 2D or 3D data at a specific time
      a=size(gvar);
      if (isinf(iTime) || iTime == -1)
        iTime=a(1);  % assumes time is 1st dimension
      end
      data = squeeze(gvar.data(iTime,:,:,:));
      if nargout==2,
        grd_temp = gvar.grid_interop(iTime,:,:,:);
        grd_temp.z = squeeze(grd_temp.z);
      end
    case 4 % return a single 2D at a specific time and layer
      a=size(gvar);
      if (isinf(iTime) || iTime == -1)
        iTime=a(1);  % assumes time is 1st dimension
      end
      if (isinf(level) || level == -1)
        level=a(2);  % assumes level is 2nd dimension
      end
      data = squeeze(gvar.data(iTime,level,:,:));
      if nargout==2,
        grd_temp = gvar.grid_interop(iTime,level,:,:);
        grd_temp.z = squeeze(grd_temp.z);
      end
  end
  try
      grd.lat = grd_temp.lat;
  end
  try
      grd.lon = grd_temp.lon;
  end
  try
      grd.z = grd_temp.z;
  end
  try
      grd.time = grd_temp.time;
  end
catch ME
  ME.message
end

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
%EIMM_SMOOTH  EKF based fixed-interval IMM smoother using two IMM-EKF filters.
%
% Syntax:
%   [X_S,P_S,X_IS,P_IS,MU_S] = EIMM_SMOOTH(MM,PP,MM_i,PP_i,MU,p_ij,mu_0j,ind,dims,A,a,a_param,Q,R,H,h,h_param,Y)
%
% In:
%   MM    - Means of forward-time IMM-filter on each time step
%   PP    - Covariances of forward-time IMM-filter on each time step
%   MM_i  - Model-conditional means of forward-time IMM-filter on each time step
%   PP_i  - Model-conditional covariances of forward-time IMM-filter on each time step
%   MU    - Model probabilities of forward-time IMM-filter on each time step
%   p_ij  - Model transition probability matrix
%   ind   - Indices of state components for each model as a cell array
%   dims  - Total number of different state components in the combined system
%   A     - Dynamic model matrices for each linear model and Jacobians of each
%           non-linear model's measurement model function as a cell array
%   a     - Cell array containing function handles for dynamic functions
%           for each model having non-linear dynamics
%   a_param - Parameters of a as a cell array.
%   Q     - Process noise matrices for each model as a cell array.
%   R     - Measurement noise matrices for each model as a cell array.
%   H     - Measurement matrices for each linear model and Jacobians of each
%           non-linear model's measurement model function as a cell array
%   h     - Cell array containing function handles for measurement functions
%           for each model having non-linear measurements
%   h_param - Parameters of h as a cell array.
%   Y     - Measurement sequence
%
% Out:
%   X_S  - Smoothed state means for each time step
%   P_S  - Smoothed state covariances for each time step
%   X_IS - Model-conditioned smoothed state means for each time step
%   P_IS - Model-conditioned smoothed state covariances for each time step
%   MU_S - Smoothed model probabilities for each time step
%
% Description:
%   EKF based two-filter fixed-interval IMM smoother.
%
% See also:
%   EIMM_UPDATE, EIMM_PREDICT

% History:
%   09.01.2008 JH The first official version.
%
% Copyright (C) 2007,2008 Jouni Hartikainen
%
% $Id: imm_update.m 111 2007-11-01 12:09:23Z jmjharti $
%
% 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.

function [x_sk,P_sk,x_sik,P_sik,mu_sk] = eimm_smooth(MM,PP,MM_i,PP_i,MU,p_ij,mu_0j,ind,dims,A,a,a_param,Q,R,H,h,h_param,Y)
    % Default values for mean and covariance
    MM_def = zeros(dims,1);
    PP_def = diag(ones(dims,1));

    % Number of models
    m = length(A);

    % Number of measurements
    n = size(Y,2);

    % The prior model probabilities for each step
    p_jk = zeros(m,n);
    p_jk(:,1) = mu_0j;
    for i1 = 2:n
        for i2 = 1:m
            p_jk(i2,i1) = sum(p_ij(:,i2).*p_jk(:,i1-1));
        end
    end

    % Backward-time transition probabilities
    p_ijb = cell(1,n);
    for k = 1:n
        for i1 = 1:m
            % Normalizing constant
            b_i = sum(p_ij(:,i1).*p_jk(:,k));
            for j = 1:m
                p_ijb{k}(i1,j) = 1/b_i.*p_ij(j,i1).*p_jk(j,k);
            end
        end
    end

    % Space for overall smoothed estimates
    x_sk = zeros(dims,n);
    P_sk = zeros(dims,dims,n);
    mu_sk = zeros(m,n);

    % Values of smoothed estimates at the last time step.
    x_sk(:,end)   = MM(:,end);
    P_sk(:,:,end) = PP(:,:,end);
    mu_sk(:,end)  = MU(:,end);

    % Space for model-conditioned smoothed estimates
    x_sik = cell(m,n);
    P_sik = cell(m,n);

    % Values for last time step
    x_sik(:,end) = MM_i(:,end);
    P_sik(:,end) = PP_i(:,end);

    % Backward-time estimated model probabilities
    mu_bp = MU(:,end);

    % Space for model-conditioned backward-time updated means and covariances
    x_bki = cell(1,m);
    P_bki = cell(1,m);

    % Initialize with default values
    for i1 = 1:m
       x_bki{i1} = MM_def;
       x_bki{i1}(ind{i1}) = MM_i{i1,end};
       P_bki{i1} = PP_def;
       P_bki{i1}(ind{i1},ind{i1}) = PP_i{i1,end};
    end

    % Space for model-conditioned backward-time predicted means and covariances
    x_kp = cell(1,m);
    P_kp = cell(1,m);

    % Initialize with default values
    for i1 = 1:m
       x_kp{i1} = MM_def;
       P_kp{i1} = PP_def;
    end

    for k = n-1:-1:1
        % Space for normalizing constants and conditional model probabilities
        a_j = zeros(1,m);
        mu_bijp = zeros(m,m);

        for i2 = 1:m
            % Normalizing constant
            a_j(i2) = sum(p_ijb{k}(:,i2).*mu_bp(:));
             % Conditional model probability
            mu_bijp(:,i2) = 1/a_j(i2).*p_ijb{k}(:,i2).*mu_bp(:);

            % Retrieve the transition matrix or the Jacobian of the dynamic model
            if isnumeric(A{i2})
                A2 = A{i2};
            elseif isstr(A{i2}) | strcmp(class(A{i2}),'function_handle')
                A2 = feval(A{i2},x_bki{i2}(ind{i2}),a_param{i2});
            else
                A2 = A{i2}(x_bki{i2}(ind{i2}),a_param{i2});
            end

            % Backward-time EKF prediction step
            [x_kp{i2}(ind{i2}), P_kp{i2}(ind{i2},ind{i2})] = ekf_predict1(x_bki{i2}(ind{i2}),...
                                                              P_bki{i2}(ind{i2},ind{i2}),...
                                                              inv(A2),Q{i2},a{i2},[],a_param{i2});

        end

        % Space for mixed predicted mean and covariance
        x_kp0 = cell(1,m);
        P_kp0 = cell(1,m);

        % Space for measurement likelihoods
        lhood_j = zeros(1,m);

        for i2 = 1:m
            % Initialize with default values
            x_kp0{i2} = MM_def;
            P_kp0{i2} = PP_def;
            P_kp0{i2}(ind{i2},ind{i2}) = zeros(length(ind{i2}),length(ind{i2}));

            % Mix the mean
            for i1 = 1:m
                x_kp0{i2}(ind{i2}) = x_kp0{i2}(ind{i2}) + mu_bijp(i1,i2)*x_kp{i1}(ind{i2});
            end

            % Mix the covariance
            for i1 = 1:m
                P_kp0{i2}(ind{i2},ind{i2}) = P_kp0{i2}(ind{i2},ind{i2}) + mu_bijp(i1,i2)*(P_kp{i1}(ind{i2},ind{i2})+(x_kp{i1}(ind{i2})-x_kp0{i2}(ind{i2}))*(x_kp{i1}(ind{i2})-x_kp0{i2}(ind{i2}))');
            end

            % Backward-time EKF update
            %
            % If the measurement model is linear don't pass h and h_param to ekf_update1
            if isempty(h) | isempty(h{i2})
                [x_bki{i2}(ind{i2}), P_bki{i2}(ind{i2},ind{i2}),K,MUP,S,lhood_j(i2)] = ekf_update1(x_kp0{i2}(ind{i2}),P_kp0{i2}(ind{i2},ind{i2}),Y(:,k),H{i2},R{i2},[],[],[]);
            else
                [x_bki{i2}(ind{i2}), P_bki{i2}(ind{i2},ind{i2}),K,MUP,S,lhood_j(i2)] = ekf_update1(x_kp0{i2}(ind{i2}),P_kp0{i2}(ind{i2},ind{i2}),Y(:,k),H{i2},R{i2},h{i2},[],h_param{i2});
            end
        end

        % Normalizing constant
        a_s = sum(lhood_j.*a_j);
        % Updated model probabilities
        mu_bp = 1/a_s.*a_j.*lhood_j;

        % Space for conditional measurement likelihoods
        lhood_ji = zeros(m,m);
        for i1 = 1:m
            for i2 = 1:m
                d_ijk = MM_def;
                D_ijk = PP_def;
                d_ijk = d_ijk + x_kp{i1};
                d_ijk(ind{i2}) = d_ijk(ind{i2}) - MM_i{i2,k};
                PP2 = zeros(dims,dims);
                PP2(ind{i2},ind{i2}) = PP_i{i2,k};
                D_ijk = P_kp{i1} + PP2;

                % Calculate the (approximate) conditional measurement likelihoods
                lhood_ji(i2,i1) = gauss_pdf(d_ijk,0,D_ijk);
            end
        end

        d_j = zeros(m,1);
        for i2 = 1:m
           d_j(i2) = sum(p_ij(i2,:).*lhood_ji(i2,:));
        end
        d = sum(d_j.*MU(:,k));

        mu_ijsp = zeros(m,m);
        for i1 = 1:m
            for i2 = 1:m
                mu_ijsp(i1,i2) = 1./d_j(i2)*p_ij(i2,i1)*lhood_ji(i2,i1);
            end
        end

        mu_sk(:,k) = 1/d.*d_j.*MU(:,k);

        % Space for two-step conditional smoothing distributions p(x_k^j|m_{k+1}^i,y_{1:N}),
        % which are a products of two Gaussians
        x_jis = cell(m,m);
        P_jis = cell(m,m);
        for i2 = 1:m
            for i1 = 1:m
                MM1 = MM_def;
                MM1(ind{i2}) = MM_i{i2,k};

                PP1 = PP_def;
                PP1(ind{i2},ind{i2}) = PP_i{i2,k};

                iPP1 = inv(PP1);
                iPP2 = inv(P_kp{i1});

                % Covariance of the Gaussian product
                P_jis{i2,i1} = inv(iPP1+iPP2);
                % Mean of the Gaussian product
                x_jis{i2,i1} = P_jis{i2,i1}*(iPP1*MM1 + iPP2*x_kp{i1});
            end
        end

        % Mix the two-step conditional distributions to yield model-conditioned
        % smoothing distributions.
        for i2 = 1:m
            % Initialize with default values
            x_sik{i2,k} = MM_def;
            P_sik{i2,k} = PP_def;
            P_sik{i2,k}(ind{i2},ind{i2}) = zeros(length(ind{i2}),length(ind{i2}));

            % Mixed mean
            for i1 = 1:m
                x_sik{i2,k} = x_sik{i2,k} + mu_ijsp(i1,i2)*x_jis{i2,i1};
            end

            % Mixed covariance
            for i1 = 1:m
                P_sik{i2,k} = P_sik{i2,k} + mu_ijsp(i1,i2)*(P_jis{i2,i1} + (x_jis{i2,i1}-x_sik{i2,k})*(x_jis{i2,i1}-x_sik{i2,k})');
            end
        end

        % Mix the overall smoothed mean
        for i1 = 1:m
            x_sk(:,k) = x_sk(:,k) + mu_sk(i1,k)*x_sik{i1,k};
        end

        % Mix the overall smoothed covariance
        for i1 = 1:m
            P_sk(:,:,k) = P_sk(:,:,k) + mu_sk(i1,k)*(P_sik{i1,k} + (x_sik{i1,k}-x_sk(:,k))*(x_sik{i1,k}-x_sk(:,k))');
        end

    end

author: EEA-sensors

repo: ekfukf

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


text:
function [gmn, gm, dgm] = spm_eeg_inv_vbecd_getLF(s, sens, vol, step)
% Estimation of the leadfield matrix and its spatial derivative if required
% for a set of dipoles used in the VB-ECD solution
%% scales non-eeg data up by a fixed factor (1e8) for compatibility of
%% units
%
% FORMAT [gmn, gm, dgm] = spm_eeg_inv_vbecd_getLF(s, sens, vol, step)
%
% s      - location vector
% sens   - sensor locations (MNI [mm])
% vol    - volume structure needed by fieldtrip
% step   - stepsize to compute numerical derivatives
%
% gmn    - leadfields (three vectors for each dipole)
% gm     - vectorized leadfields
% dgm    - vectorized partials wrt locations
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Christophe Phillips & Stefan Kiebel
% $Id: spm_eeg_inv_vbecd_getLF.m 3833 2010-04-22 14:49:48Z vladimir $

%%% now does rank reduction (to 2) for non eeg data

MEGRANK=2;
MEGSCALE=1e10;  %% should be the same as in spm_eeg_inv_vbecd.m

if nargin<4
     step = 0;
end

gm = [];
for i = 1:length(s)/3

    % mean correction of LF, only for EEG data.
    if ft_senstype(sens, 'eeg')
       [tmp] = ft_compute_leadfield(s(1+(i-1)*3:i*3)', sens, vol);
        tmp = tmp - repmat(mean(tmp), size(tmp,1), 1);
    else %% reduce rank of leadfield for MEG- assume one direction (radial) is silent
       [tmp] = ft_compute_leadfield(s(1+(i-1)*3:i*3)', sens, vol,'reducerank',MEGRANK);
        tmp=tmp.*MEGSCALE;
    end
    gm = [gm tmp];
end

gmn = gm; % leadfield

[Nc, Np] = size(gmn);
if nargout >= 2
    gm = gmn(:); % vectorized leadfield
end

if all(step > 0) && nargout == 3,
%     if isempty(step),
%         step=randn(size(s));
%     end; % if isempty
    dgm = [];
    for j = 1:length(s)
        ds = s;
        ds(j) = s(j) + step(j);
        dtmp = [];
        for i = 1:length(s)/3
            if ceil(j/3) == i

                % mean correction of LF, only for EEG data.
                if ft_senstype(sens, 'eeg')
                    [tmp] = ft_compute_leadfield(ds(1+(i-1)*3:i*3)', sens, vol);
                    tmp = tmp - repmat(mean(tmp), size(tmp,1), 1);
                else  % MEG
                   [tmp] = ft_compute_leadfield(ds(1+(i-1)*3:i*3)', sens, vol,'reducerank',MEGRANK);
                    tmp=tmp.*MEGSCALE;
                end
                dtmp = [dtmp tmp];
            else
                dtmp = [dtmp gmn(:, 1+(i-1)*3:i*3)];
            end
        end
        dtmp = dtmp(:);
        dgm = [dgm (-gm + dtmp)./step(j)];
    end

    % correct order
    ind = reshape(1:Np^2 , Np, Np)';

    dgm = reshape(dgm, size(gmn, 1), Np^2);
    dgm = dgm(:, ind(:));
    dgm = reshape(dgm, Np*Nc, Np);
end

author: spm


text:
% extracts the center (cc,cr) and radius of the largest blob
function [cc,cr,radius,flag]=extractball(Imwork,Imback,index)%,fig1,fig2,fig3,fig15,index)

  cc = 0;
  cr = 0;
  radius = 0;
  flag = 0;
  [MR,MC,Dim] = size(Imback);

  % subtract background & select pixels with a big difference
  fore = zeros(MR,MC);          %image subtracktion
  fore = (abs(Imwork(:,:,1)-Imback(:,:,1)) > 10) ...
     | (abs(Imwork(:,:,2) - Imback(:,:,2)) > 10) ...
     | (abs(Imwork(:,:,3) - Imback(:,:,3)) > 10);

  % Morphology Operation  erode to remove small noise
  foremm = bwmorph(fore,'erode',2); %2 time

  % select largest object
  labeled = bwlabel(foremm,4);
  stats = regionprops(labeled,['basic']);%basic mohem nist
  [N,W] = size(stats);
  if N < 1
    return
  end

  % do bubble sort (large to small) on regions in case there are more than 1
  id = zeros(N);
  for i = 1 : N
    id(i) = i;
  end
  for i = 1 : N-1
    for j = i+1 : N
      if stats(i).Area < stats(j).Area
        tmp = stats(i);
        stats(i) = stats(j);
        stats(j) = tmp;
        tmp = id(i);
        id(i) = id(j);
        id(j) = tmp;
      end
    end
  end

  % make sure that there is at least 1 big region
  if stats(1).Area < 100
    return
  end
  selected = (labeled==id(1));

  % get center of mass and radius of largest
  centroid = stats(1).Centroid;
  radius = sqrt(stats(1).Area/pi);
  cc = centroid(1);
  cr = centroid(2);
  flag = 1;
  return

author: Sable

repo: mcbench-benchmarks


text:
function family=getclusterstat(c)

%GETCLUSTERSTAT gets detais of individual correlation clusters.
% family = GETCLUSTERSTAT(C) outputs a structure FAMILY which contains
% pertinent details of clusters that exist in C. This is a simple routine
% that assumes all cluster calculation has already been done. Note that in
% many cases it may be more appropriate to get the direct CLUST field from
% a correlation object. In some cases however, GETCLUSTERSTAT may be more
% expedient. The FAMILY structure includes the fields:
%   rank:    size order of the cluster (scalar)
%   numel:   number of events in the cluster (scalar)
%   index:   the index of the cluster events in the correlation object(Nx1)
%   trig:    trigger times of the cluster events (Nx1)
%   begin:   the eariest trigger time in the cluster (scalar)
%   finish:  the last trigger time in the cluster (scalar)
% The fields index and trig are cell arrays.
%
% See also correlation/cluster

% AUTHOR: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

if nargin <= 0
    error('Wrong number of inputs');
end

if ~isa(c,'correlation')
    error('First input must be a correlation object');
end

if isempty(get(c,'CLUST'))
    error('CLUSTER field must be filled in input argument. See HELP CLUSTER');
end

for n = 1:max(get(c,'CLUST'))
    f = find(c,'CLUST',n);
    c1 = subset(c,f);
    family.rank(n) = n;
    family.numel(n) = numel(f);
    family.index(n) = {f};
    family.begin(n) = min(get(c1,'TRIG'));
    family.finish(n) = max(get(c1,'TRIG'));
    family.trig(n) = {get(c1,'TRIG')};
end

author: geoscience-community-codes

repo: GISMO


text:
function [rpm] = THz2rpm(THz)
% Convert frequency from terahertz to revolutions per minute.
% Chad A. Greene 2012
rpm = THz*60*1e+12;


text:
function [errorsA, errorsLUB, stuckPoints] = verifyPoints(sampleStruct)
% Verify that a set of points are in the solutoin space of `sampleStruct`.
% Typically, this method would be called to check a set of warmup points or
% points generated via `gpSampler`. Also verifies if points moved from
% warmup points.
%
% USAGE:
%
%    [errorsA, errorsLUB, stuckPoints] = verifyPoints(sampleStruct)
%
% INPUT:
%    sampleStruct:    LPProblem containing points and warmup points
%
% OUTPUTS:
%    errorsA:         Row index of the constraint in `sampleStruct` that
%                     is not consistent with the given points
%    errorsLUB:       Upper and lower bounds of the constraint + tolerance
%    stuckPoints:     Index of points which did not move.
%
% .. Authors:
%       - Ellen Tsai 2007
%       - Richard Que 12/1/09 Combined with checkWP.m

[warmupPts, points] = deal(sampleStruct.warmupPts, sampleStruct.points);
%above a check to see if warmup points moved
n=size(points,2);
len=size(n,1);
stuckPoints=zeros(0,1);

for i=1:n
    len(i,1)=norm(warmupPts(:,i)-points(:,i));
    if len(i,1)<10
        stuckPoints=[stuckPoints; i];
    end
end
minL = min(len);

%check to see if points are within solution space
[A,b,csense]=deal(sampleStruct.A,sampleStruct.b,sampleStruct.csense);
t=.01; %tolerance
npoints = size(points, 2);

LHS = A*points;
EIndex = (csense == 'E');
LIndex = (csense == 'L');
GIndex = (csense == 'G');

[errorsL(:,1), errorsL(:,2)] = find(LHS(LIndex,:) - b(LIndex)* ones(1, npoints) > t);
[errorsG(:,1), errorsG(:,2)] = find(b(GIndex)* ones(1, npoints) - LHS(GIndex,:)   > t);
[errorsE(:,1), errorsE(:,2)] = find( abs(LHS(EIndex,:) - b(EIndex)* ones(1, npoints))   > t);

errorsUB = find(points > sampleStruct.ub*ones(1,npoints) + t);
errorsLB = find(points < sampleStruct.lb*ones(1,npoints) - t);

errorsA = [errorsL; errorsG; errorsE];
errorsLUB = [errorsLB; errorsUB];

end

author: opencobra

repo: cobratoolbox


text:
function [interpImg, interpImg3D] = rxInterpSlice(rx, rxSlice, ori);
%
% [interpImg, interpImg3D] = rxInterpSlice(rx, rxSlice, [ori=3]);
%
% Given the transform matrix and slice specified by
% the mrRx GUI, compute an interpolated slice.
%
% The optional third argument, ori, specifies the orientation of the slice
% to take relative to the prescription:
%	3 [default] - prescribed slice (corresponds to the slices shown in the
%					Rx window)
%	2 - slice corresponds to a column of the prescription. (e.g., if the
%					prescription is sagittal, this is coronal.)
%	1 - slice corresponds to a row of the prescription. (e.g., if the
%					prescription is sagittal, this is axial.)
%
%
% Returns two images: interpImg is a 2D matrix
% with the raw interpolated values; interpImg3D
% is a 3D TrueColor image, with the brightness
% and contrast settings of the GUI taken into
% account.
%
%
% ras 03/05.
% ras, 02/08/08: added orientation flag.
if ~exist('rx', 'var') | isempty(rx),
	rx = get(findobj('Tag','rxControlFig'), 'UserData');
end
if ~exist('ori', 'var') | isempty(ori), ori = 3;				 end
if ~exist('rxSlice', 'var') | isempty(rxSlice)
    rxSlice = get(rx.ui.rxSlice.sliderHandle, 'Value');
end

ysz = rx.rxDims(1);
xsz = rx.rxDims(2);
zsz = rx.rxDims(3);

%% get the sampling range for each dimension
switch ori
	case 1,	 % row
		yRange = rxSlice;
		xRange = 1:xsz;
		zRange = 1:zsz;
		imgSize = [xsz zsz];

	case 2,  % column
		yRange = 1:ysz;
		xRange = rxSlice;
		zRange = 1:zsz;
		imgSize = [ysz zsz];

	case 3,  % slice
		yRange = 1:ysz;
		xRange = 1:xsz;
		zRange = rxSlice;
		imgSize = [ysz xsz];

	otherwise, error('Invalid orientation flag.')
end

% build up coords for interp vol slice
% here we flip: x->cols, y->rows
[X Y Z] = meshgrid(xRange, yRange, zRange);
coords = [X(:) Y(:) Z(:) ones(size(Z(:)))]';
coords = double(coords);

% get new coords for the interp slice
newCoords = rx.xform * coords;

% interpolate
interpVals = myCinterp3(rx.vol, rx.volDims(1:2), rx.volDims(3), ...
						newCoords(1:3,:)', 0);
interpImg = reshape(interpVals, imgSize);

% convert to truecolor if requested
if nargout >= 2
	if ishandle(rx.ui.interpBright.sliderHandle)
		brightness = rx.ui.interpBright;
	else
		brightness = 0.5;
	end

	if ishandle(rx.ui.interpContrast.sliderHandle)
		contrast = rx.ui.interpContrast;
	else
		contrast = .7;
	end

	if ishandle(rx.ui.interpHistoThresh)
		histoThresh = rx.ui.interpHistoThresh;
	else
		histoThresh = 1;
	end

    interpImg3D = rxClip(interpImg, [], brightness, contrast, histoThresh);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function calSave7 (~, ~, xt, cumu2, as)
    % Save cum #  and z value to an interactively chosen file

    %msg.infodisp('  ','Save Data');

    [file1,path1] = uiputfile(fullfile(ZmapGlobal.Data.Directories.output, '*.dat'), 'Earthquake Datafile');
    if file1 && path1
        data = [xt', cumu2', as']';
        fid = fopen([path1 file1],'w') ;
        fprintf(fid, '%6.2f  %6.2f %6.2f\n' , data);
        fclose(fid) ;
    else
        %msg.dbdisp('  ','cancelled save');
    end

end

author: CelsoReyes


text:
function [e1, e2, E1_l, E2_l] = ahmLinEndpoints(l,t1,t2)

% AHMLINENDPOINTS  AHM line endpoints.
%   AHMLINENDPOINTS(L,T) returns the 3D Euclidean endpoint of the
%   AHM line L at abscissa T.
%
%   [E1,E2] = AHMLINENDPOINTS(L,T1,T2) returns two 3D Euclidean endpoints
%   at abscissas T1 and T2.
%
%   [e1,e2,E1_l,E2_l] = AHMLINENDPINTS(...) returns the Jacobians of the
%   endopints wrt L.
%
%   See also AHMLINSEGMENT, AHMLIN2SEG, AHMLIN2IDPPNTS.

%   Copyright 2009 Teresa Vidal.

if nargout <= 2

    % support segment
    s  = ahmLin2seg(l);

    % support points
    p1 = s(1:3);
    p2 = s(4:6);

    % endpoints
    e1 = (1-t1)*p1 + t1*p2;
    e2 = (1-t2)*p1 + t2*p2;

else

    [s, S_l] = ahmLin2seg(l);       % support segment

    p1 = s(1:3);                    % support points
    p2 = s(4:6);

    P1_l = S_l(1:3,:);              % jacobians of sup points
    P2_l = S_l(4:6,:);

    e1 = (1-t1)*p1 + t1*p2;         % endpoints
    e2 = (1-t2)*p1 + t2*p2;

    E1_p1 = 1-t1;                   % jacobians of endpoints
    E1_p2 = t1;
    E2_p1 = 1-t2;
    E2_p2 = t2;

    E1_l = E1_p1*P1_l + E1_p2*P2_l; % chain rule
    E2_l = E2_p1*P1_l + E2_p2*P2_l;

end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function [y, s2, fmu, fs2] = simulGPnaive(hyp, inf, mean, cov, lik, input, target, test, lag)
% simulGPnaive - 'Naive' (i.e. without propagation of variance) simulation of the GP
% ARX and AR model.
%
%% Syntax
%   [y, s2] = simulGPnaive(hyp, inf, mean, cov, lik, input, target, test, lag)
%
%% Description
% GP-NARX model is used to predict a further step ahead by replacing the data at
% present time instant with the data at one time instant before and using the mean
% value of prediction from the previous prediction step instead of the measured
% output value. This is then repeated as many times as there are test samples.
%
% Input:
% * hyp    ... the column vector of hyperparameters
% * inf    ... the function specifying the inference method
% * cov    ... the prior covariance function (see below)
% * mean   ... the prior mean function
% * lik    ... the likelihood function
% * input  ... the input part of the training data,  NxD matrix
% * target ... the output part of the training data (ie. target), Nx1 vector
% * test   ... the input matrix, kxD matrix, see construct.m for more info
% * lag    ... the order of the model (number of used lagged outputs)
%
% Output:
% * y    ... the mean predicted output
% * s2   ... associated variances
% * fmu  ... the mean predicted output without noise
% * fs2  ... associated variances without noise
%
% See also:
% gpx.m, simulGPmc, construct.m
%
% Examples:
% demo_example_gp_simulation.m
%
%%
% Written by K. Azman in J. Kocijan, 2007

[NN,D] = size(test); % D - input space dimension, NN - number or test samples

% allocate memory
y = zeros(1,D);
s2 = y;
fmu = y;
fs2 = y;

% first step: k=1
xt = test(1,:);
[y(1), s2(1) ,fmu(1), fs2(1), post] = gpx(hyp, inf, mean, cov, lik, input, target, xt);

% future steps ...
for k=2:NN

    if (mod(k,100) == 0)  % remark on every 100th step
        disp(['simulGPnaive, step: ', int2str(k), '/', int2str(NN)]);
    end
        % For the next prediction prepare the input vector ...
        if D > lag
            if (k>lag)
                xt = [y(k-lag:k-1) test(k, lag+1:end)];
            elseif (k<=lag)
                xt = [test(k, 1:lag-k+1) y(1:k-1) test(k, lag+1:end)];
            end
        else                                            % if D <= lag
            if (k>lag)
                xt = y(k-lag:k-1);
            elseif (k<=lag)
                xt = [test(k, 1:lag-k+1) y(1:k-1)];
            end
        end
        % ... and calculate the one-step-ahead prediction
    [y(k), s2(k) ,fmu(k), fs2(k), post] = gpx(hyp, inf, mean, cov, lik, input, target, xt, post);

end

% transform into column vectors
y = y';
s2 = s2';
fmu = fmu';
fs2 = fs2';

author: Dynamic-Systems-and-GP

repo: GPdyn


text:
function out=lagx(X,p)

% function out=lagX(X,p)
% returns lagged values of a matrix, sequentially lagging by 0,1,...,p periods and returning all the lags stacked in a single matrix
% inputs:  - matrix 'X': data matrix
%          - integer 'p': number of lags
% outputs: -matrix 'out': matrix of current and lagged series

% Compute the number of rows and columns of input matrix X, and save the values as r and c
[r,c]=size(X);
out = zeros(numel(p+1:r), size(X,2));
% start with the reformating of the original vectors
out(:,1:c)=X(p+1:r,:);

% treat each lag consecutively
for ii=1:p
out(:,c*ii+1:c*(ii+1))=X(p+1-ii:r-ii,:);
end

author: european-central-bank

repo: BEAR-toolbox