text:
function [x,stat] = libqp_gsmo(H,f,a,b,LB,UB,x0,opt)
% LIBQP_GSMO Generalized SMO algorithm.
%
% Synopsis:
%  [x,stat] = libqp_gsmo(H,f,a,b,LB,UB)
%  [x,stat] = libqp_gsmo(H,f,a,b,LB,UB,x0)
%  [x,stat] = libqp_gsmo(H,f,a,b,LB,UB,[],opt)
%  [x,stat] = libqp_gsmo(H,f,a,b,LB,UB,x0,opt)
%
% Description:
%  This function implements the Generalized SMO algorithm which solves
%  the following convex quadratic programming task:
%
%   min 0.5*x'*H*x + f'*x
%    x
%
%   subject to   a'*x = b
%                LB(i) <= x(i) <= UB(i)   for all i=1:n
%
% Reference:
%  S.-S. Keerthi, E.G.Gilbert. Convergence of a Generalized SMO Algorithm for SVM
%   Classifier Design. Technical Report CD-00-01, Control Division, Dept. of Mechanical
%   and Production Engineering, National University of Singapore, 2000.
%   http://citeseer.ist.psu.edu/keerthi00convergence.html
%
% Input:
%  H [n x n] Symmetric positive semi-definite matrix.
%  f [n x 1] Vector.
%  a [n x 1] Vector which must not contain zero entries.
%  b [1 x 1] Scalar.
%  LB [n x 1] Lower bound; -inf is allowed.
%  UB [n x 1] Upper bound; inf is allowed.
%
% Optional inputs:
%  x0 [n x 1] Initial solution.
%
%  options [struct]
%    .TolKKT [1 x 1] Determines relaxed KKT conditions (default TolKKT=0.001);
%       it correspondes to $\tau$ in Keerthi's paper.
%    .verb [1 x 1] if > 0 then prints info every verb-th iterations (default 0)
%    .MaxIter [1 x 1] Maximal number of iterations (default inf).
%
% Output:
%  x [n x 1] Solution vector.
%  stat [struct]
%   .QP [1x1] Primal objective value.
%   .exitflag [1 x 1] Indicates which stopping condition was used:
%      nIter >= MaxIter                  ->  exitflag = 0
%      relaxed KKT conditions satisfied  ->  exitflag = 4
%   .nIter [1x1] Number of iterations.
%
% Example:
%  n=50; X = rand(n,n); H = X'*X; f = -10*rand(n,1);
%  a=rand(1,n)+0.1; b=rand; tmp = rand(n,1); LB = tmp-1; UB = tmp+1;
%
%  tic; [x1,stat1] = libqp_gsmo(H,f,a,b,LB,UB); fval1=stat1.QP, toc
%  tic; [x2,fval2] = quadprog(H,f,[],[],a,b,LB,UB); fval2, toc
%

% Copyright (C) 2006-2008 Vojtech Franc, xfrancv@cmp.felk.cvut.cz
% Center for Machine Perception, CTU FEL Prague

% options
if nargin < 6, error('At least six input arguments are required.'); end
if nargin < 7 || isempty(x0),
  % create feasible solution x0 if not given
  x0 = zeros(length(f),1);
  xa = 0; i = 0;
  while x0'*a(:) ~= b,
    i = i + 1;
    if i > length(a),
      error('No feasible solution exists.');
    end
    x0(i) = min(UB(i),max(LB(i),(b-xa)/a(i)));
    xa = xa + x0(i)*a(i);
  end
end

if nargin < 8, opt = []; end
if ~isfield(opt,'TolKKT'), opt.TolKKT = 0.001; end
if ~isfield(opt,'MaxIter'), opt.MaxIter = inf; end
if ~isfield(opt,'verb'), opt.verb = 0; end

[x,stat.QP,stat.exitflag,stat.nIter] = ...
    libqp_gsmo_mex(H,f,a,b,LB,UB,x0,opt.MaxIter,opt.TolKKT,opt.verb);

return;
% EOF

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function cs_demo3 (do_pause, matrixpath)
%CS_DEMO3 MATLAB version of the CSparse/Demo/cs_demo3.c program.
%   Cholesky update/downdate.
%
% Example:
%   cs_demo3
% See also: cs_demo

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

if (nargin < 2)
    matrixpath = [] ;
end

if (isempty (matrixpath))
    try
        % older versions of MATLAB do not have an input argument to mfilename
        p = mfilename ('fullpath') ;
        t = strfind (p, filesep) ;
        matrixpath = [ p(1:t(end)) '../../Matrix' ] ;
    catch
        % assume we are in the C*Sparse/MATLAB/CSparse/Demo directory
        matrixpath = '../../Matrix' ;
    end
end

matrices = { 'HB/bcsstk01', 'HB/bcsstk16' } ;

if (nargin < 1)
    do_pause = 1 ;
end

for i = 1:length(matrices)
    name = matrices {i} ;
    [C sym] = get_problem (matrixpath, name, 1e-14) ;
    demo3 (C, sym, name) ;
    if (do_pause)
        input ('Hit enter to continue: ') ;
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function ins=ins_mech(ins,imu)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INS mechanization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2020-2025, by Kai Chen, All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent old_dw old_dv

ins.time=imu.time;
if isempty(old_dw)
    old_dw=zeros(3,1);
end
if isempty(old_dv)
    old_dv=zeros(3,1);
end

% correct bias and scaling factor errors
dw0 = imu.dw';
dv0 = imu.dv';
dw = ins.Kg*dw0-ins.bg*ins.nt;
dv = ins.Ka*dv0-ins.ba*ins.nt;

% extrapolate velocity and position
vel_mid = ins.vel+ins.acc*(ins.nt/2);
pos_mid = ins.pos+ins.Mpv*(ins.vel+vel_mid)/2*ins.nt;

% update the related parameters
ins.eth = earth_update(pos_mid,vel_mid);
ins.wib = dw/ins.nt;
ins.fb  = dv/ins.nt;
ins.fn  = ins.Cnb*ins.fb;
ins.web = ins.wib-ins.Cnb'*ins.eth.wnie;

% update velocity
dv_rot  = 0.5*cross(dw0,dv0);
dv_scul = 1/12*(cross(old_dw,dv0)+cross(old_dv,dw0));
dv_sf   = (eye(3)-0.5*ins.nt*askew(ins.eth.wnin))*ins.Cnb*dv + ins.Cnb*(dv_rot+dv_scul);
dv_cor  = ins.eth.gcc*ins.nt;
vel_new = ins.vel+dv_sf+dv_cor;

% update position
ins.Mpv(2) = ins.eth.Mpv2;
ins.Mpv(4) = ins.eth.Mpv4;
pos_new    = ins.pos + ins.Mpv*(ins.vel+vel_new)/2*ins.nt;

% update attitude
dw_cone  = 1/12*cross(old_dw,dw0);
phi_b_ib = dw+dw_cone;
phi_n_in = ins.eth.wnin*ins.nt;
Cbb = rvec2mat(phi_b_ib);
Cnn = rvec2mat(phi_n_in)';
Cnb_new = Cnn*ins.Cnb*Cbb;
att_new = Cnb2att(Cnb_new);

% update INS result
ins.Cnb = Cnb_new;
ins.att = att_new;
ins.vel = vel_new;
ins.pos = pos_new;
ins.x = [ins.att;ins.vel;ins.pos;ins.bg;ins.ba];

old_dw=dw0;
old_dv=dv0;

return

repo: GINav


text:
%% housekeeping
clear
close all
clc()

%% Create dataset
clc

do_plot=false;

scale=100;

[db,varlist0]=create_dataset(scale,do_plot);

varlist=fieldnames(varlist0);
%% Choose a model type: see cell "create the structural VAR model" below

for imod=1:7

model_type=imod-1;

% set up the restrictions
%---------------------------
close()

% create restrictions on parameters as well as markov chains
%------------------------------------------------------------
switch model_type
    case 0
        % constant-parameter model
        [restrictions,markov_chains,switch_prior]=create_restrictions_and_markov_chains0();
    case 1
        % Coefficients are switching regimes across all equations
        % (synchronized case)
        [restrictions,markov_chains,switch_prior]=create_restrictions_and_markov_chains1();
    case 2
        % Coefficients and variances have different chains, different
        % regimes, and different durations
        [restrictions,markov_chains,switch_prior]=create_restrictions_and_markov_chains2();
    case 3
        % Only coefficients in monetary policy equation are changing
        [restrictions,markov_chains,switch_prior]=create_restrictions_and_markov_chains3();
    case 4
        % Only variance in monetary policy equation is changing
        [restrictions,markov_chains,switch_prior]=create_restrictions_and_markov_chains4();
    case 5
        % Both coefficients and variances in monetary policy equation
        % change with two independent Markov processes
        [restrictions,markov_chains,switch_prior]=create_restrictions_and_markov_chains5();
    case 6 % ok
        % Only variances in ALL three equations switch
        [restrictions,markov_chains,switch_prior]=create_restrictions_and_markov_chains6();
    otherwise
        error('the coded model types are 0, 1, 2, 3, 4 and 6')
end

% Create the VAR
%----------------

clc

nlags=2;

exog={};

panel=[];

constant=true;

% first we create a template structure
% ------------------------------------
sv0=svar(varlist,exog,nlags,constant,panel,markov_chains);

% set priors
%-----------

var_prior=svar.prior_template();

var_prior.type='sz';

prior=struct('var',var_prior,'nonvar',switch_prior);

is_prior=true;

if ~is_prior

    prior=rmfield(prior,'var');

end

% Find posterior mode
%-------------------------

sv=sv0;

sv=estimate(sv,db,{'1960Q1','2015Q2'},prior,restrictions);

end

author: jmaih


text:
function table_io_test ( )

%*****************************************************************************80
%
%% TABLE_IO_TEST tests the TABLE_IO library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2006
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TABLE_IO_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the TABLE_IO library.\n' );

  table_io_test01 ( );
  table_io_test02 ( );
  table_io_test03 ( );
  table_io_test04 ( );
  table_io_test05 ( );
  table_io_test06 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TABLE_IO_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function E = ShortTimeEnergy(signal, windowLength,step)
signal = signal / max(max(signal));
curPos = 1;
L = length(signal);
numOfFrames = floor((L-windowLength)/step) + 1;
%H = hamming(windowLength);
E = zeros(numOfFrames,1);
for (i=1:numOfFrames)
    window = (signal(curPos:curPos+windowLength-1));
    E(i) = (1/(windowLength)) * sum(abs(window.^2));
    curPos = curPos + step;
end
%Max = max(E);
%med = median(E);
%m = mean(E);
%a = length(find(E>2*med))/numOfFrames;
%b = length(find(E<(med/2)))/numOfFrames;

%S_EnergyM = length(find(E>2*med))/numOfFrames;
%S_EnergyV = length(find(E<(med/2)))/numOfFrames;
%S_EnergyM = std(E);
%S_EnergyV = a;

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function [P3, S_bar, V, Z, RO, Tr] = random_nr_motion(T, J, K, state)
% [P3, S_bar, V, Z, RO, Tr] = random_nr_motion(T, J, K, state)
%
%  INPUT:
%
%  T           - number of frames
%  J           - number of points
%  K           - number of deformation basis
%
%  OUTPUT:
%
%  P3          - (3*T) x J 3D-motion matrix:    P3([t t+T t+2*T],:) contains the 3D coordinates of the J points at time t
%  S_bar       - shape average:                 3 x J matrix
%  V           - deformation shapes:            (3*K) x J matrix     ( V((n-1)*3+[1:3],:) contains the n-th deformation basis )
%  Z           - deformation weights:           T x K matrix
%  RO          - rotation:                      cell array           ( RO{t} gives the rotation matrix at time t )
%  Tr          - translation:                   T x 2 matrix

rand('state', state);
randn('state', state);

S_bar = rand(3, J);

[q,r] = qr(rand(3*J));

V = zeros(3*K, J);
for kk=1:K,
   V(1+(kk-1)*3:3*kk, :) = reshape(q(:,kk), 3, J);
end

Z = randn(T,K);
Tr = randn(T,2);

a = (rand(T,1)-0.5)*2*pi;
b = (rand(T,1)-0.5)*2*pi;
c = (rand(T,1)-0.5)*2*pi;
P3 = zeros(3*T, J);
for t=1:T,
   R1 = [1 0 0; 0 cos(a(t)) -sin(a(t)); 0 sin(a(t)) cos(a(t))];
   R2 = [cos(b(t)) 0 sin(b(t)); 0 1 0; -sin(b(t)) 0 cos(b(t))];
   R3 = [cos(c(t)) -sin(c(t)) 0; sin(c(t)) cos(c(t)) 0; 0 0 1];

   RO{t} = R1*R2*R3;

   Sdef = S_bar;
   for kk = 1:K,
      Sdef = Sdef+Z(t,kk)*V((kk-1)*3+[1:3],:);
   end;
   Sdef = RO{t}*Sdef +[Tr(t,:)'; 0]*ones(1, J);

   P3([t t+T t+2*T], :) = Sdef;
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function [output] = F_power_spectrum(input_layer)
% assume the input is a DxTxN matrix of complex spectrum, where D is the dimension of the feature vector,
% T is the number of frames in the minibatch or utterance, and N is the
% number of sentences
ComplexSpectrum = input_layer.a;
output = abs( ComplexSpectrum .* conj(ComplexSpectrum) );

end

author: singaxiong

repo: SignalGraph


text:
function [X,Zf] = filter_fast(B,A,X,Zi,dim)
% Like filter(), but faster when both the filter and the signal are long.
% [Y,Zf] = filter_fast(B,A,X,Zi,Dim)
%
% Uses FFT convolution. The function is faster than filter when approx. length(B)>256 and
% size(X,Dim)>1024, otherwise slower (due size-testing overhead).
%
% See also:
%   filter, fftfilt
%
%                           Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                           2010-07-09
%
%                           contains fftfilt.m from Octave:
%                           Copyright (C) 1996-1997 John W. Eaton

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

if nargin <= 4
    dim = find(size(X)~=1,1); end
if nargin <= 3
    Zi = []; end

lenx = size(X,dim);
lenb = length(B);
if lenx == 0
    % empty X
    Zf = Zi;
elseif lenb < 256 || lenx<1024 || lenx <= lenb || lenx*lenb < 4000000 || ~isequal(A,1)
    % use the regular filter
    if nargout > 1
        [X,Zf] = filter(B,A,X,Zi,dim);
    else
        X = filter(B,A,X,Zi,dim);
    end
else
    % fftfilt can be used, determine which algorithm to use
    try
        fftfilt(1,1);
        do_fftfilt = @fftfilt;
    catch
        do_fftfilt = @oct_fftfilt;
    end
    was_single = isa(X,'single');

    if isempty(Zi)
        % no initial conditions to take care of
        if nargout < 2
            % and no final ones
            X = unflip(do_fftfilt(B,flip(double(X),dim)),dim);
        else
            % final conditions needed
            X = flip(X,dim);
            [dummy,Zf] = filter(B,1,X(end-length(B)+1:end,:),Zi,1); %#ok<ASGLU>
            X = do_fftfilt(B,double(X));
            X = unflip(X,dim);
        end
    else
        % initial conditions available
        X = flip(X,dim);
        % get a Zi-informed piece
        tmp = filter(B,1,X(1:length(B),:),Zi,1);
        if nargout > 1
            % also need final conditions
            [dummy,Zf] = filter(B,1,X(end-length(B)+1:end,:),Zi,1); %#ok<ASGLU>
        end
        X = do_fftfilt(B,double(X));
        % incorporate the piece
        X(1:length(B),:) = tmp;
        X = unflip(X,dim);
    end
    if was_single
        X = single(X); end
end

function X = flip(X,dim)
if dim ~= 1
    order = 1:ndims(X);
    order = order([dim 1]);
    X = permute(X,order);
end

function X = unflip(X,dim)
if dim ~= 1
    order = 1:ndims(X);
    order = order([dim 1]);
    X = ipermute(X,order);
end

function y = oct_fftfilt(b, x, N)
% Copyright (C) 1996, 1997 John W. Eaton
%
% This file is part of Octave.
%
% Octave 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, or (at your option)
% any later version.
%
% Octave 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 Octave; see the file COPYING.  If not, write to the Free
% Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA.
%
% -*- texinfo -*-
% @deftypefn {Function File} {} fftfilt (@var{b}, @var{x}, @var{n})
%
% With two arguments, @code{fftfilt} filters @var{x} with the FIR filter
% @var{b} using the FFT.
%
% Given the optional third argument, @var{n}, @code{fftfilt} uses the
% overlap-add method to filter @var{x} with @var{b} using an N-point FFT.
%
% If @var{x} is a matrix, filter each column of the matrix.
% @end deftypefn
%
% Author: Kurt Hornik <Kurt.Hornik@wu-wien.ac.at>
% Created: 3 September 1994
% Adapted-By: jwe

% If N is not specified explicitly, we do not use the overlap-add
% method at all because loops are really slow.  Otherwise, we only
% ensure that the number of points in the FFT is the smallest power
% of two larger than N and length(b).  This could result in length
% one blocks, but if the user knows better ...
transpose = (size(x,1) == 1);

if transpose
    x = x.'; end

[r_x,c_x] = size(x);
[r_b,c_b] = size(b);
if min([r_b, c_b]) ~= 1
    error('octave:fftfilt','fftfilt: b should be a vector'); end

l_b = r_b*c_b;
b = reshape(b,l_b,1);

if nargin == 2
    % Use FFT with the smallest power of 2 which is >= length (x) +
    % length (b) - 1 as number of points ...
    N = 2^(ceil(log(r_x+l_b-1)/log(2)));
    B = fft(b,N);
    y = ifft(fft(x,N).*B(:,ones(1,c_x)));
else
    % Use overlap-add method ...
    if ~isscalar(N)
        error ('octave:fftfilt','fftfilt: N has to be a scalar'); end
    N = 2^(ceil(log(max([N,l_b]))/log(2)));
    L = N - l_b + 1;
    B = fft(b, N);
    B = B(:,ones(c_x,1));
    R = ceil(r_x / L);
    y = zeros(r_x, c_x);
    for r = 1:R
        lo = (r - 1) * L + 1;
        hi = min(r * L, r_x);
        tmp = zeros(N, c_x);
        tmp(1:(hi-lo+1),:) = x(lo:hi,:);
        tmp = ifft(fft(tmp).*B);
        hi = min(lo+N-1, r_x);
        y(lo:hi,:) = y(lo:hi,:) + tmp(1:(hi-lo+1),:);
    end
end

y = y(1:r_x,:);
if transpose
    y = y.'; end

% Final cleanups: if both x and b are real respectively integer, y
% should also be
if isreal(b) && isreal(x)
    y = real(y); end
if ~any(b - round(b))
    idx = ~any(x - round(x));
    y(:,idx) = round(y(:,idx));
end

repo: REST


text:
function y = construct_harmonic_signal( freq,amp,phase,dc,t )
% y = construct_harmonic_signal( freq,amp,phase,dc,t )
% reconstructs a periodical signal from given harmonics
% Written by Dr. Yoash Levron, April 2013.
%
% The input is the individual harmonics of the signal.
% the function reconstructs the signal y(t).
%
% for example, given the input :
% freq        = [50   100  150   200 ...]
% amp       = [3     0      0        0     ...]
% phase    = [0     0      0        0     ...]
% dc  = 2
% The function will output the signal:
% y = 2 + 3*cos(2*pi*50*t)
%
% Generally, the output signal y is:
% y = dc + amp(1)*cos(2*pi*freq(1)*t + phase(1))
%             +amp(2)*cos(2*pi*freq(2)*t + phase(2)) + ...
%
% This function is designed compatible to
% function 'fourier_series_real'. The two functions are inversible:
% The following code will generate a signal y(t) identical to x(t) :
% [ freq,amp,phase,dc ] = fourier_series_real( t,x );
% y = construct_harmonic_signal( freq,amp,phase,dc,t);
%
% inputs:
% freq - [Hz] frequencies of the fourier series, not including zero.
% amp - amplitudes vector, not including the DC component.
% phase - [rad/sec] . phases, not including the DC component.
% dc - the DC value (average of the signal).
% t - [sec] time vector.
%
% outputs:
% y - output signal. same length as t.

y = dc * t.^0;
for ii = 1:length(freq)
    y = y + amp(ii)*cos(2*pi*freq(ii)*t + phase(ii));
end

end

author: Sable

repo: mcbench-benchmarks


text:
function rxnFBS = diffexprs2rxnFBS(model, diff_exprs, Vref, varargin)
% Returns Forward - Backward - Unchanged (+1;0;-1) for each reaction.
% (+1)R_f    (-1)R_b     (0)unchanged
%
% USAGE:
%
%    rxnFBS = diffexprs2rxnFBS(model, diff_exprs, Vref, 'SeparateTranscript', '', 'logFC', 0, 'pval',0.05)
%
% INPUT:
%    model:             The COBRA Model structure
%    diff_exprs:        MATLAB Table including the information of the
%                       differentially expressed genes.
%                       Required columns (with theses names):
%                           - gene ( ID of gene, same as in the meabolic model)
%                           - logFC (SOURCE VS TARGET)
%                           - pval (p-value or adjusted-p-value)
%    Vref:              Reference flux of the model
%
% OPTIONAL INPUTS:
%    varargin:          `ParameterName` value pairs with the following options:
%
%                           - `SeparateTranscript`: Character used to separate different transcripts of a gene. (default: '')
%                           - `logFC`: minimum log2 (fold change) requiered (default = 0)
%                           - `pval`: maximum p-value admited (default = 0.05)
%
% OUTPUT:
%    rxnFBS:             array containting the information of altered
%                       reactions: Forward - Backward - Unchanged
%                       (+1;0;-1);
%
% .. Note:
%    It is highly recommended to use as diff_exprs the TopTable result from
%    limma in R and change names if neccesary.
%    In the tutorial there is further information to load data from GEO or
%    using R to conduct absolute and differential gene expression analysis.
%
% .. Authors:
%       - Luis V. Valcarcel, 25/06/2015, University of Navarra, CIMA & TECNUN School of Engineering.
%       - Luis V. Valcarcel, 26/10/2018, University of Navarra, CIMA & TECNUN School of Engineering.
%       - Francisco J. Planes, 26/10/2018, University of Navarra, TECNUN School of Engineering.

p = inputParser; % check input information
p.CaseSensitive = false;
addParameter(p, 'SeparateTranscript', '', @(x)ischar(x));
addParameter(p, 'logFC', 0, @(x)isnumeric(x)&&isscalar(x));
addParameter(p, 'pval', 0.05, @(x)isnumeric(x)&&isscalar(x));
parse(p, varargin{:});
SeparateTranscript = p.Results.SeparateTranscript;
logFC = p.Results.logFC;
pval = p.Results.pval;

% Traslate omics data from gene level to reaction level
if ~isempty(SeparateTranscript)
    aux_table = table(strtok(model.genes,SeparateTranscript),model.genes,'VariableNames',{'gene' 'transcript'});
    diff_exprs = innerjoin(aux_table, diff_exprs, 'Keys','gene');
	diff_exprs.gene = diff_exprs.transcript;
end

% Check that genes are in the model
assert(sum(ismember(diff_exprs.gene,model.genes))>0,...
    'Gene ID are not in the model. Revise gene ID in input table or SeparateTranscript option parameter')

% Minimun Fold Change to admit the changes
idx = abs(diff_exprs.logFC) > logFC;
diff_exprs = diff_exprs(idx,:);

% Maximum p-value to admit the changes
idx = abs(diff_exprs.pval) <= pval;
diff_exprs = diff_exprs(idx,:);

pos_up = diff_exprs.logFC > 0;     % source state bigger than target state
pos_down = diff_exprs.logFC < 0;   % source state smaller than target state

% As we cannot map the reaction to genes, we can discretize thess values as
% (+1, 0 -1), to create geneFBS
geneFBS_aux = zeros(size(diff_exprs,1),1);
geneFBS_aux(pos_up) = +1;
geneFBS_aux(pos_down) = -1;

% Generate geneFBS for all genes in the model
geneFBS = zeros(length(model.genes),1);
[~,idx] = ismember(diff_exprs.gene,model.genes);
geneFBS_aux = geneFBS_aux(idx~=0);
idx= idx(idx~=0);
geneFBS(idx) = geneFBS_aux;

% If a gene is down-regulated in the source state,the flux activity should
% be increased. Similarly, if a gene is up-regulated in the source state,
% the flux activity should be decreased.
% For this reason we change the sign of geneFBS (target vs source here).
geneFBS = -geneFBS;

% Change to Cobra format
geneFBS = struct('gene',{model.genes},'value',geneFBS);

if isempty(SeparateTranscript)
    fprintf('\tGene expression changes calculated\n');
    fprintf('\tThere are %u genes that are differentially expressed\n',sum(geneFBS.value~=0));
else
    fprintf('\tGene expression changes calculated\n');
    fprintf('\tThere are %u trainscripts that are differentially expressed\n',sum(geneFBS.value~=0));
    [~,idx] = unique(strtok(geneFBS.gene,SeparateTranscript));
    fprintf('\tThere are %u genes that are differentially expressed\n',sum(geneFBS.value(idx)~=0));
end

% Transform geneFBS into rxnFBS
% in order to produce changes, all changes must be in the same sense
minSum = false;
parsedGPR = GPRparser(model,minSum);

% Calculate changes to produce more flux
geneFBS_plus = geneFBS;
geneFBS_plus.value(geneFBS_plus.value<=0) = 0;
geneFBS_plus.value(geneFBS_plus.value>0) = +1;
[gene_id_plus, gene_expr_plus] = findUsedGenesLevels(model,geneFBS_plus);
rxnFBS_plus = selectGeneFromGPR(model, gene_id_plus, gene_expr_plus, parsedGPR, minSum);

% Calculate changes to produce less flux
geneFBS_minus = geneFBS;
geneFBS_minus.value(geneFBS_minus.value>=0) = 0;
geneFBS_minus.value(geneFBS_minus.value<0) = +1;
[gene_id_minus, gene_expr_minus] = findUsedGenesLevels(model,geneFBS_minus);
rxnFBS_minus = selectGeneFromGPR(model, gene_id_minus, gene_expr_minus, parsedGPR, minSum);

% all changes in same sense or they neglect each other
rxnFBS_plus(rxnFBS_plus<=0) = 0;
rxnFBS_minus(rxnFBS_minus<=0) = 0;
rxnFBS = rxnFBS_plus - rxnFBS_minus;
rxnFBS(strcmp(model.grRules,'')) = 0; % non-mapped reactions have no change in activity

% We have considered all reactions as Forward reactions. GeneFBS gives the
% information whatever a reaction should increase activity in absolute
% values.
% For Forward reactions, (+1) increase flux, (-1) decrease flux
% For Backward reactions, (+1) decrease flux (increase in absolute values),
% (-1) increase flux, (decrease in absolute values)
rxnFBS(Vref < 0) = - rxnFBS(Vref < 0);
fprintf('\tReaction expression changes calculated\n');
fprintf('\tThere are %u reactions that are differentially expressed\n',sum(rxnFBS~=0));

end

author: opencobra

repo: cobratoolbox


text:
function Dvect=entropyGradDist(npd)
%
% Compute entropy estimate using nearest neighbor estimate
%

% Copyright (C) 2003 Alexander Ihler; distributable under GPL -- see README.txt

pts = getPoints(npd);
Ce = .57721566490153286;
[N1,N2] = size(pts);
[I,D] = knn(npd,pts,2);
I = I(2,:);
Dvect = pts - pts(:,I);

%Sr = N1* pi^(N1/2) / gamma((N1/2) + 1);
%h = N1/N2 * sum( log(D) ) + log(Sr * (N2-1)/N1 ) + Ce;

Dvect = N1/N2 * Dvect ./ repmat(D.^2,[N1,1]);   % find gradient direction
%Dvect = .1 * Dvect / max(max(Dvect));        % scale for epsilon steps

author: ShapeNet

repo: RenderForCNN


text:
function install_eq_toolbox(arg)
%INSTALL_EQ_TOOLBOX Install using Toolbox Installer, with sensible defaults
%
%Syntax
% install_eq_toolbox(arg);
%
%Description
% INSTALL_EQ_TOOLBOX uses Toolbox Installer to install this toolbox. It must
% be run from the top level directory of this toolbox, which must contain the
% file info.ins.
%
% If INSTALL_EQ_TOOLBOX is called with no argument, it asks the user for
% confirmation before taking any action.
%
% INSTALL_EQ_TOOLBOX first checks if the file pathname.m exists on the current
% Matlab path and can be appended to.
% If so, it asks the user for confirmation and then simply calls the INSTALL
% function of Toolbox Installer to update the file pathname.m.
% If not, it tries to choose a sensible directory in which to create a new
% pathdef.m file, and suggests this to the user before proceeding.
%
% INSTALL_EQ_TOOLBOX DIRNAME, that is INSTALL_EQ_TOOLBOX('DIRNAME'),
% calls INSTALL('DIRNAME'), using DIRNAME as the directory for pathdef.m.
%
% INSTALL_EQ_TOOLBOX -, that is INSTALL_EQ_TOOLBOX('-'), calls INSTALL('-'),
% which displays a dialog box so the user can select the directory to be used
% for pathdef.m.
%
%Notes
% If it is necessary to create a new pathdef.m file, the directory which
% INSTALL_EQ_TOOLBOX suggests to the user for this file defaults to
% the directory which contains startup.m, if this exists.
%
% Otherwise INSTALL_EQ_TOOLBOX suggests the following directory to the user:
% For Unix systems including Linux and Mac: $HOME/matlab
% For Windows systems: [matlabroot '\toolbox\local'].
% If this directory does not exist, INSTALL_EQ_TOOLBOX tries to create it.
%
% For Toolbox Installer 2.2 by B. Rasmus Anthin, see
% Matlab Central File Exchange
% http://www.mathworks.com/matlabcentral/fileexchange/
% Home page of B. Rasmus Anthin
% http://www.etek.chalmers.se/~e8rasmus/
%
%Examples
% > install_eq_toolbox
% Info: Installer will create the file /c/leopardi/matlab/pathdef.m
% Proceed (y/n)? y
% Info: Installer has created the directory /c/leopardi/matlab
% Info: Installed.
%
% > install_eq_toolbox
% Info: Installer will update the file /c/leopardi/matlab/pathdef.m
% Proceed (y/n)? y
% Info: Installed.
%
% > install_eq_toolbox
% Info: Installer will update the file /c/leopardi/matlab/pathdef.m
% Info: Proceed (y/n)? n
%
% > cd ~
% > install_eq_toolbox
% ??? Error using ==> install_eq_toolbox
% Please run install_eq_toolbox from toolbox directory, eg. /c/leopardi/[...]
%
%See also
% UNINSTALL_EQ_TOOLBOX, PATHDEF, PATH, Toolbox Installer 2.2

% Copyright 2004-2005 Paul Leopardi for the University of New South Wales.
% $Revision 1.10 $ $Date 2005-06-01 $
% Partial overhaul
% o Restructure tests for existing pathdef.m file and guesses for pathdef_dir
% o Simplify tests after mkdir
% o On Windows, do not use [getenv('USERPROFILE') '\matlab'] directory
% o Clearly distinguish between info and error messages
% Documentation files renamed
% $Revision 1.00 $ $Date 2005-02-12 $
%
% For licensing, see COPYING.
% For references, see AUTHORS.
% For revision history, see CHANGELOG.

if (nargin == 1) && strcmp(arg,'-')
    install(arg);
    return;
end

command = mfilename;
[command_dir name ext ver] = fileparts(mfilename('fullpath'));
wd = pwd;
if ~strcmp(command_dir,wd)
    error(['Please run ' command ' from toolbox directory, eg. ' command_dir]);
end

if nargin
    %
    % If user has used an argument, just call install with this argument.
    %
    install(arg);
    return;
end

estr = 'Error: ';
wstr = 'Warning: ';
istr = 'Info: ';
was_not_installed_msg = 'Installer cannot continue.';
try_to_create_msg =     'Try to create the directory outside of Matlab';

current_pathdef_name = which('pathdef');
pathdef_exists = (exist('pathdef') == 2);
if pathdef_exists
    pathdef_name = current_pathdef_name;
    %
    % We have been told that pathdef.m exists.
    % Try to open it for updating
    %
    pathdefid = fopen(pathdef_name,'a');
    if pathdefid > 0
        %
        % File open for pathdef.m succeeded.
        % Ask whether to proceed with update.
        %
        fclose(pathdefid);
        reply = input([istr 'Installer will update the file ' strrep(pathdef_name,'\','\\') ...
            '\nProceed (y/n)? '],'s');
        if strcmpi(reply(1),'y')
            install;
        end
        return;
    elseif (exist('startup') == 2)
        %
        % File open for pathdef.m failed.
        % We have been told that startup.m exists.
        % Guess that we should use pathdef.m in the same directory.
        %
        startup_name = which('startup');
        [pathdef_dir name ext ver] = fileparts(startup_name);
        pathdef_name = fullfile(pathdef_dir,'pathdef.m');
    else
        %
        % File open for pathdef.m failed.
        % We have been told that no startup.m exists.
        %
        if isunix
            %
            % We are using Unix.
            % Guess that we should use $HOME/matlab/pathdef.m
            %
            [status,dir_str] = system('echo $HOME');
            home_dir = deblank(dir_str);
            pathdef_name = fullfile(home_dir,'matlab','pathdef.m');
        else
            %
            % We are using Windows.
            % Guess that we should use $matlabroot\toolbox\local\pathdef.m
            %
            pathdef_name = fullfile(matlabroot,'toolbox','local','pathdef.m');
        end
    end
    [pathdef_dir name ext ver] = fileparts(pathdef_name);
end

if strcmp(pathdef_name, current_pathdef_name)
    %
    % If we have not changed pathdef_name,
    % then this is the same file that we have been told we cannot write to,
    % and we have not been able to guess another possiblity for pathdef_name.
    %
    disp([wstr 'Cannot write to the pathdef file ' ...
        pathdef_name ' which was found on the Matlab path.']);
    reply = input([istr 'Do you want to select an existing pathdef.m file to update (y/n)? '], 's');
    if strcmpi(reply(1),'y')
        %
        % The user has chosen to update an existing pathdef.m file.
        % Use install -
        %
        install('-');
        return;
    end
    reply = input([istr 'Do you want to select a directory and create a new pathdef.m file (y/n)? '], 's');
    if ~strcmpi(reply(1),'y')
        %
        % The user has chosen not to create a new pathdef.m file.
        % There are no other possibilities, so return.
        %
        return;
    end
    %
    % The user has chosen to create a new pathdef.m file.
    % Ask where it should go. The directory must exist.
    %
    pname=uigetdir('','Please choose an existing directory for the new ''pathdef.m'' file.');
    if isequal(pname,0)
        %
        % The user has cancelled the directory selection.
        % There are no other possibilities, so return.
        %
        return;
    end
    pathdef_name = fullfile(pname,'pathdef.m');
    [pathdef_dir name ext ver] = fileparts(pathdef_name);
end
%
% Determine if Matlab thinks that pathdef_name already exists as a file.
%
if exist(pathdef_name) == 2
    willstr = 'update';
else
    willstr = 'create';
end
%
% Ask before attempting to perform the appropriate action.
%
reply = input([istr 'Installer will ' willstr ' the file ' strrep(pathdef_name,'\','\\') ...
    '\nProceed (y/n)? '],'s');
if ~strcmpi(reply(1),'y')
    return;
end
%
% Determine if Matlab thinks that pathdef_dir already exists as a directory.
%
pathdef_dir_exists = (exist(pathdef_dir) == 7);
%
% Even if Matlab says that pathdef_dir exists, it might not exist as a directory.
% Try to create it anyway. This should not hurt if it really does exist.
%
[parent dirname ext ver] = fileparts(pathdef_dir);
[status message] = mkdir(parent,dirname);
if status < 0
    created_pathdef_dir = false;
    if ~pathdef_dir_exists
        %
        % Matlab said pathdef_dir does not exist and that we cannot not create it.
        %
        disp([estr was_not_installed_msg]);
        disp([estr 'Cannot create ' pathdef_dir ': ' message]);
        error([try_to_create_msg '\n and then try "' command pathdef_dir '"']);
    end
else
    created_pathdef_dir = true;
    if ~pathdef_dir_exists
        %
        % Matlab has said that pathdef_dir did not exist and that we have now created it.
        %
        disp([istr command ' has created the directory ' pathdef_dir]);
    end
end
%
% At this point, pathdef_dir should exist as a directory.
% Check again.
%
pathdef_dir_exists = (exist(pathdef_dir) == 7);
if ~pathdef_dir_exists
    %
    % Matlab said pathdef_dir does not exist.
    %
    disp([estr was_not_installed_msg]);
    disp([estr 'Cannot create ' pathdef_dir ': ' message]);
    error([try_to_create_msg '\n and then try "' command pathdef_dir '"']);
end
%
% Finally, we can ask install to use pathdef_dir.
%
install(pathdef_dir);
%
% If we created pathdef_dir, we should add it to the path,
% so that the current Matlab session can see the new pathdef.m
%
if created_pathdef_dir
    addpath(pathdef_dir);
end
rehash pathreset;
%
%end function

repo: gspbox


text:
% datasetCandi = {'siftsmall'};
datasetCandi = {'sift','gist'};
% datasetCandi = {'sift'};
% datasetCandi = {'gist'};
% datasetCandi = {'imagenet'};

% methodCandi = {'LSH','ITQ'};

% methodCandi = {'USPLH'};
% methodCandi = {'CPH'};
% methodCandi = {'SH'};
% methodCandi = {'ITQ'};
% methodCandi = {'IsoH' 'SpH'};
% methodCandi = {'LSH'};

% methodCandi = {'DSH'};
% methodCandi = {'BRE'};
% methodCandi = {'SpH'};

% methodCandi = {'IsoH1'};

methodCandi = {'LSH','SpH','BRE','USPLH','ITQ','IsoH','SH'};
% methodCandi = {'KLSH','AGH1','AGH2','CH'};

codelengthCandi = [32, 64, 128, 256, 512, 1024];

for d=1:length(datasetCandi)
    dataset = datasetCandi{d};

    for m=1:length(methodCandi)
        method = methodCandi{m};

        for c=1:length(codelengthCandi)

            codelength = codelengthCandi(c);

            ResultFile = ['../',dataset,'/hashingCodeLong.',num2str(codelength),'/',method,'table',upper(dataset),'32b_1'];
            bFound = checkFILEmkDIR(ResultFile);
            if(bFound)
                error('Result file exists!');
            end

            disp('==============================');
            disp([method,' ',num2str(codelength),'bit ',dataset]);
            disp('==============================');

            trainset = double(fvecs_read (['../',dataset,'/',dataset,'_base.fvecs']));
            testset = fvecs_read (['../',dataset,'/',dataset,'_query.fvecs']);
            trainset = trainset';
            testset = testset';

            trainStr = ['[model, trainB ,train_elapse] = ',method,'_learn(trainset, codelength);'];
            testStr = ['[testB,test_elapse] = ',method,'_compress(testset, model);'];
            eval(trainStr);
            eval(testStr);

            ntrain=size(trainB,1);
            ntest=size(testB,1);

            realcodelength = size(trainB,2);
            disp([method,' ',num2str(realcodelength),'bit learned. ',dataset]);

            nFile = fix(realcodelength/32);
            nRemain = mod(realcodelength,32);

            disp('==============================');
            disp(['Total training time: ',num2str(train_elapse)]);
            disp(['Total testing time: ',num2str(test_elapse)]);

            %continue;

            for j =1:nFile
                ResultFile = ['../',dataset,'/hashingCodeLong.',num2str(codelength),'/',method,'table',upper(dataset),'32b_',num2str(j)];

                fid = fopen(ResultFile,'w');
                fwrite(fid, 1, 'uint32');
                fwrite(fid, 32, 'uint32');
                fwrite(fid, ntrain, 'uint32');
                for i = 1 : ntrain
                    tmpV = uint32(0);
                    for k = 32*(j-1)+1:32*j
                        tmpV = bitshift(tmpV,1);
                        tmpV = tmpV + uint32(trainB(i,k));
                    end
                    fwrite(fid, tmpV, 'uint32');
                end
                fclose(fid);

                ResultFile = ['../',dataset,'/hashingCodeLong.',num2str(codelength),'/',method,'query',upper(dataset),'32b_',num2str(j)];
                fid = fopen(ResultFile,'w');
                fwrite(fid, 1, 'uint32');
                fwrite(fid, 32, 'uint32');
                fwrite(fid, ntest, 'uint32');
                for i = 1 : ntest
                    tmpV = uint32(0);
                    for k = 32*(j-1)+1:32*j
                        tmpV = bitshift(tmpV,1);
                        tmpV = tmpV + uint32(testB(i,k));
                    end
                    fwrite(fid, tmpV, 'uint32');
                end
                fclose(fid);
            end

            if nRemain > 0
                ResultFile = ['../',dataset,'/hashingCodeLong.',num2str(codelength),'/',method,'table',upper(dataset),'32b_',num2str(nFile+1)];

                fid = fopen(ResultFile,'w');
                fwrite(fid, 1, 'uint32');
                fwrite(fid, 32, 'uint32');
                fwrite(fid, ntrain, 'uint32');
                for i = 1 : ntrain
                    tmpV = uint32(0);
                    for k = 32*nFile+1:size(trainB,2)
                        tmpV = bitshift(tmpV,1);
                        tmpV = tmpV + uint32(trainB(i,k));
                    end
                    fwrite(fid, tmpV, 'uint32');
                end
                fclose(fid);

                ResultFile = ['../',dataset,'/hashingCodeLong.',num2str(codelength),'/',method,'query',upper(dataset),'32b_',num2str(nFile+1)];
                fid = fopen(ResultFile,'w');
                fwrite(fid, 1, 'uint32');
                fwrite(fid, 32, 'uint32');
                fwrite(fid, ntest, 'uint32');
                for i = 1 : ntest
                    tmpV = uint32(0);
                    for k = 32*nFile+1:size(testB,2)
                        tmpV = bitshift(tmpV,1);
                        tmpV = tmpV + uint32(testB(i,k));
                    end
                    fwrite(fid, tmpV, 'uint32');
                end
                fclose(fid);

            end

            disp('binary file saved!');
        end

    end
end

author: ZJULearning

repo: MatlabFunc


text:
function [A,xHyp,PHyp,GateMat]=makeStandardCartOnlyLRMatHyps(xPred,SPred,zCart,SRCart,SRCartMax,PD,lambda,rPred,gammaVal,measJacobDet)
%%MAKESTANDARDCARTONLYLRMATHYPS Create the likelihood ratio matrix used for
%         2D assignment under standard Gaussian approximations in the
%         coordinate system of the target states using Gaussian-
%         approximated Cartesian position-only measurements and provide
%         updated state estimates when tracks gate with the measurements.
%         Brute-force gating is used. This function is appropriate for use
%         in Cartesian-converted measurement tracking without range rate
%         information. This assumes that the false alarms occur due to a
%         Poisson point process with a constant density lambda that is
%         either in the coordinate system of the states (Cartesian) or in
%         the coordinate system of a converted measurement when appropriate
%         parameters are passed. This implements the unitless ratio of [1],
%         except the missed detection probability with gating is 1-PD*PG
%         and not just 1-PD due to Section 3.5.3 of [2] (for multiple
%         gammaVal values, the larged PG is chosen). Additionally, the
%         covariance inflation due to gating applied to a missed detection
%         hypothesis (due to [3], Equation 28) can be used, if desired.
%
%INPUTS: xPred A numDimXnumTar set of predicted target states. The
%              first posDim Components are position.
%        SPred A numDimXnumDimXnumTar set of lower-triangular square-root
%              state covariance matrices.
%        zCart The zDimXnumMeas set of measurements in Cartesian
%              coordinates.
%       SRCart A zDimXzDimXnumMeas set of lower-triangular square-root
%              measurement covariance matrices. If all of the covariance
%              matrices are the same, then SRCart can be a single zDimXzDim
%              lower triangular matrix.
%    SRCartMax A zDimXzDimXnumTar set of maximum measurement covariance
%              matrices that can be taken into consideration for gating and
%              will be used to compute the covariance inflation of [3]
%              (Equation 28) for missed detection hypotheses. If only a
%              single zDimXzDim matrix is passed, it is assumed the same
%              for all targets. If an empty matrix is passed, then the
%              covariance inflation is not performed. The inflation is
%              insignificant if the gate probability PG is almost 1 and
%              PD is less than 1 by a decent amount. On the other hand, if
%              PD is high but PG is low, then this can be the difference
%              between track divergence or not.
%           PD The target detection probabilities. This can either be a
%              numTarX1 vector if the targets have different probabilities,
%              or this can be a scalar value if all of the target detection
%              probabilities are the same. This value does not include the
%              effects of gating (set by gammaVal) eliminating measurements
%              from consideration.
%       lambda The false alarm density. When given in Cartesian
%              coordinates, this has units of # false alarm/volume, for
%              example, false alarms/m^3. When given in the coordinate
%              system of the measurement (anticipating a conversion at a
%              single point using a Jacobian), this has units of inverse of
%              the "volume" in the measurement coordinate system. For
%              example, for polar measurements, this might be inverse
%              meters*radians. If this is not in Cartesian coordinates,
%              then the inputs measJacobDet must be provided.
%        rPred When using a single-scan tracking algorithm with integrated
%              target existence probabilities, this is the numTarX1 set of
%              target existence probabilities. If omitted or an empty
%              matrix is passed, then rPred=1 is used (everything exists).
%     gammaVal This is the threshold for gating. Gating is performed
%              by a brute-force evaluation of Mahalanobis distances between
%              predicted target locations and the measurements. gammaVal is
%              related to the probability region about the target as
%              described in Chatper 2.3.2 of [2]. For a zDim-dimensional
%              measurement, one can obtain the value of gammaVal for a
%              99.97% probability region using
%              ChiSquareD.invCDF(0.9997,zDim). The default if this
%              parameter is omitted or an empty matrix is passed is Inf,
%              meaning that everything gates.
% measJacobDet If this input is omitted, then it is assumed that lambda is
%              given in the same coordinate system as the target state
%              (typically Cartesian). Otherwise, this is a length numMeas
%              array such that measJacobDet(i) is
%              det(measJacob(zNative(:,i))) where measJacob is a function
%              that computes the Jacobian matrix of the measurement and
%              zNative is the measurement in the original coordinate system
%              of the measurement, not the coordinate system of the state.
%              For example, if the original measurements were in spherical
%              coordinates, measJacob could be calcSpherConvJacob(z,0), and
%              thus measJacobDet would be determinates of the outputs of
%              calcSpherConvJacob evaluated at all of the measurements.
%
%OUTPUTS: A A numTarX(numMeas+numTar) matrix of target-measurement and
%           missed detection likelihood ratios computed as in [1] for
%           Gaussian states/(converted) measurements. Columns > numMeas
%           hold missed-detection likelihoods. Thus, off-diagonal terms
%           for columns > numMeas are set to 0 and the diagonal terms
%           set to the costs of a missed detection for each given target.
%           When given rPred, PD is just multiplied by these values before
%           computing the assignment matrix. If a measurement does not
%           gate, then the corresponding entry in A is explicitly set to
%           0.
%      xHyp A numDimXnumTarX(numMeas+1) set of conditionally updated target
%           states (using the prior, the measurement, and the
%           sqrtKalmanMeasPred and sqrtKalmanUpdateWithPred functions,
%           which are equivalent to the sqrtKalmanUpdate function), one for
%           each association hypothesis in A. For hypotheses that do not
%           gate, then the corresponding entry in xHyp is all zeros. The
%           last hypothesis (missed detection) is the same as the
%           prediction value, but the covariance might be inflated as in
%           [3] if the SRCartMax input is provided.
%      PHyp A numDimXnumDimXnumTarX(numMeas+1) set of conditionally updated
%           target covariance matrices corresponding to the state estimate
%           in xHyp. Hypotheses are not provided in square root form so as
%           to simplify their usage in JPDAF-style filters.
%   GateMat A numTarXnumMeas boolean matrix indicating whether each target
%           gates with each measurement.
%
%We are using the dimensionless score function from [1]. This formulation
%in terms of likelihood ratios eliminates the need to compute the volumes
%of detection gates. Of course, if these likelihood ratios are computed in
%Cartesian coordinates, but lambda is in a measurement coordinate system,
%this is just an approximation.
%
%The assignment matrix returned by this function is commonly used in
%functions such as assign2D, singleScanUpdate, and
%singleScanUpdateWithExistence. This function is useful when performing
%tracking using Cartesian-converted measurements.
%
%Developments of the score function are given in [1] and relevant parts of
%Section 3.5.3 and other sections of [2].
%
%REFERENCES:
%[1] Y. Bar-Shalom, S. S. Blackman, and R. J. Fitzgerald, "Dimensionless
%    score function for multiple hypothesis tracking," IEEE Transactions on
%    Aerospace and Electronic Systems, vol. 43, no. 1, pp. 392-400, Jan.
%    2007.
%[2] Y. Bar-Shalom, P. K. Willett, and X. Tian, Tracking and Data Fusion.
%    Storrs, CT: YBS Publishing, 2011.
%[3] X. R. Li, "Tracking in clutter with strongest neighbor measurements -
%    Part i: Theoretical analysis," IEEE Transactions on Automatic Control,
%    vol. 43, no. 11, pp. 1560-1578, Nov. 1998.
%
%June 2018 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

xDim=size(xPred,1);

if(nargin<10)
    measJacobDet=[];
else
    %They must all be positive.
    measJacobDet=abs(measJacobDet);
end

if(nargin<9||isempty(gammaVal))
    gammaVal=Inf;
end

numTar=size(xPred,2);
zDim=size(zCart,1);
numMeas=size(zCart,2);

if(nargin<8||isempty(rPred))
    %Assume all targets exist with probability 1.
    rPred=ones(numTar,1);
end

%If only one is given, assume they are all the same.
if(size(SRCart,3)==1)
    SRCart=repmat(SRCart,[1,1,numMeas]);
end

%If only one is given, assume they are all the same.
if(~isempty(SRCartMax)&&size(SRCartMax,3)==1)
    SRCartMax=repmat(SRCartMax,[1,1,numTar]);
end

if(~isempty(SRCartMax)&&isempty(zCart))
    %Get the correct measurement dimensionality to use to update the missed
    %detection covariance matrix with PG for the possibility that the
    %measurement is actually outside of the gate.
    zDim=size(SRCartMax,1);
end

if(isscalar(PD))
    PD=PD*ones(numTar,1);
end

%The gate probability (from an inverse chi-squared PDF). Choose PG to be
%the maximum gate when multiple gammaVal values are given.
PG=ChiSquareD.CDF(max(gammaVal(:)),zDim);

%Whether or not the target is detected is affected by whether it exists
%(the rPred).
PD=PD(:).*rPred(:);

%The measurement matrix.
H=[eye(zDim,zDim),zeros(zDim,xDim-zDim)];

numHyp=numMeas+1;%The extra one is the missed detection hypothesis.

%The likelihood ratio matrix.
A=zeros(numTar,numMeas+numTar);
%The gating matrix
GateMat=false(numTar,numMeas);

%Conditional hypothesis updates.
xHyp=zeros(xDim,numTar,numHyp);
PHyp=zeros(xDim,xDim,numTar,numHyp);
for curTar=1:numTar
    xPredCur=xPred(:,curTar);
    SPredCur=SPred(:,:,curTar);

    %The missed detection likelihood (includes the gating probability).
    A(curTar,numMeas+curTar)=(1-PD(curTar)*PG);
    xHyp(:,curTar,numHyp)=xPredCur;
    PHyp(:,:,curTar,numHyp)=SPredCur*SPredCur';

    if(numMeas>0||~isempty(SRCartMax))
        %Get the filter measurement prediction once (not for all
        %measurements). This is independent of the measurement covariance
        %matrix.
        [zPred,PzPred,otherInfo]=sqrtKalmanMeasPred(xPredCur,SPredCur,H);

        %If upper bounds on the measurement covariances to consider for
        %each target are provided.
        if(~isempty(SRCartMax))
            %The gain corresponding to a gate with the maximum considered
            %Cartesian covariance matrix for that target.
            W=calcSqrtKalmanGain(SRCartMax(:,:,curTar),otherInfo);
            Pzz=PzPred+(SRCartMax(:,:,curTar)*SRCartMax(:,:,curTar)');
            PMissed=calcMissedGateCov(PHyp(:,:,curTar,numHyp),Pzz,W,PD(curTar),gammaVal,PG);
            PHyp(:,:,curTar,numHyp)=PMissed;
        end
    end

    for curMeas=1:numMeas
        %First, we evaluate the Mahalanobis distance necessary for gating,
        %then we evaluate the measurement update and the likelihood ratio
        %if the target gates.
        innov=zCart(:,curMeas)-zPred;

        Pzz=PzPred+(SRCart(:,:,curMeas)*SRCart(:,:,curMeas)');
        mahabDist=innov'*inv(Pzz)*innov;
        if(mahabDist>gammaVal)
            %If it does not gate, then assign a zero likelihood ratio.
            A(curTar,curMeas)=0;
            continue;
        end

        %It gates.
        GateMat(curTar,curMeas)=true;

        %Perform the measurement update.
        [xUpdate, SUpdate,innov,Szz]=sqrtKalmanUpdateWithPred(zCart(:,curMeas),SRCart(:,:,curMeas),zPred,otherInfo);

        xHyp(:,curTar,curMeas)=xUpdate;
        PHyp(:,:,curTar,curMeas)=SUpdate*SUpdate';

        if(~isempty(measJacobDet))
            %If the clutter density is given in some local measurement
            %coordinate system, then use the Jacobian to do a single-
            %point transformation.
            JDet=measJacobDet(curMeas);
        else
            JDet=1;
        end

        %Evaluate the likelihood ratio.
        A(curTar,curMeas)=PD(curTar)/(JDet*lambda)*GaussianD.PDFS(innov,zeros(zDim,1),Szz);
    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