text:
function [ t, dtdr, dtds ] = shape_t6 ( r, s )

%*****************************************************************************80
%
%% SHAPE_T6 evaluates shape functions for a 6 node reference triangle.
%
%  Reference Element T6:
%
%    |
%    1  3
%    |  |\
%    |  | \
%    S  6  5
%    |  |   \
%    |  |    \
%    0  1--4--2
%    |
%    +--0--R--1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R, S, the reference coordinates of a point.
%
%    Output, real T(6), the basis functions at the point.
%
%    Output, real DTDR(6), the R basis derivatives at the point.
%
%    Output, real DTDS(6), the S basis derivatives at the point.
%
  t(1) = 2.0 *     ( 1.0 - r - s ) * ( 0.5 - r - s );
  t(2) = 2.0 * r * ( r - 0.5 );
  t(3) = 2.0 * s * ( s - 0.5 );
  t(4) = 4.0 * r * ( 1.0 - r - s );
  t(5) = 4.0 * r * s;
  t(6) = 4.0 * s * ( 1.0 - r - s );

  dtdr(1) = - 3.0 + 4.0 * r + 4.0 * s;
  dtdr(2) = - 1.0 + 4.0 * r;
  dtdr(3) =   0.0;
  dtdr(4) =   4.0 - 8.0 * r - 4.0 * s;
  dtdr(5) =                   4.0 * s;
  dtdr(6) =                 - 4.0 * s;

  dtds(1) = - 3.0 + 4.0 * r + 4.0 * s;
  dtds(2) =   0.0;
  dtds(3) = - 1.0           + 4.0 * s;
  dtds(4) =       - 4.0 * r;
  dtds(5) =         4.0 * r;
  dtds(6) =   4.0 - 4.0 * r - 8.0 * s;

  return
end


text:
% Copyright (C) 2010 Quan Wang <wangq10@rpi.edu>
% Signal Analysis and Machine Perception Laboratory
% Department of Electrical, Computer, and Systems Engineering
% Rensselaer Polytechnic Institute, Troy, NY 12180, USA

%% The testing of binary Fisher's linear discriminant
%    f: n*m feature matrix, each row being a data point
%    w: projection weights
%    t: threshold, obtained by
%    l0: n*1 ground truth of binary label vector, each element being 0 or 1
%    l: n*1 resulting binary label vector, each element being 0 or 1
%    precision, recall, accuracy, F1: measurement of performance
%    display: a flag indicating whether to display results, 0 for silent
function [l,precision,recall,accuracy,F1]=fisher_testing(f,w,t,l0,display)

fp=f*w;
l=(fp>t);

if nargin<5
    display=1;
end

if nargin<4
    precision=NaN;
    recall=NaN;
    accuracy=NaN;
    F1=NaN;
else
    tp= sum( (l0==1) & (l==1) );
    fn= sum( (l0==1) & (l==0) );
    fp= sum( (l0==0) & (l==1) );
    tn= sum( (l0==0) & (l==0) );
    precision=tp/(tp+fp);
    recall=tp/(tp+fn);
    accuracy=(tp+tn)/(tp+tn+fp+fn);
    F1=2*precision*recall/(precision+recall);
    if display>0
        fprintf(['Fisher --- precision: %f, recall: %f, accuracy: %f, ' ...
            'F1-measure: %f\n'],precision,recall,accuracy,F1);
    end
end

author: Sable

repo: mcbench-benchmarks


text:
% OP_DIV_V_Q_TP: assemble the matrix B = [b(i,j)], b(i,j) = (q_i, div v_j), exploiting the tensor product structure.
%
%   mat = op_div_v_q_tp (spv, spq, msh);
%   [rows, cols, values] = op_div_v_q_tp (spv, spq, msh);
%
% INPUT:
%
%   spv:     object representing the vector-valued space of trial functions for the velocity (see sp_vector)
%   spq:     object representing the scalar-valued space of test functions for the pressure (see sp_scalar)
%   msh:     object that defines the domain partition and the quadrature rule (see msh_cartesian)
%
% OUTPUT:
%
%   mat:    assembled matrix
%   rows:   row indices of the nonzero entries
%   cols:   column indices of the nonzero entries
%   values: values of the nonzero entries
%
% Copyright (C) 2009, 2010 Carlo de Falco
% Copyright (C) 2011 Andrea Bressan
% Copyright (C) 2011, 2017 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

function varargout = op_div_v_q_tp (spv, spq, msh)

  for idim = 1:msh.ndim
    size2 = size (spq.sp_univ(idim).connectivity);
    for icomp = 1:spv.ncomp_param
      size1 = size (spv.scalar_spaces{icomp}.sp_univ(idim).connectivity);
      if (size1(2) ~= size2(2) || size1(2) ~= msh.nel_dir(idim))
        error ('One of the discrete spaces is not associated to the mesh')
      end
    end
  end

  A = spalloc (spq.ndof, spv.ndof, 3*spv.ndof);

  for iel = 1:msh.nel_dir(1)
    msh_col = msh_evaluate_col (msh, iel);
    spv_col = sp_evaluate_col (spv, msh_col, 'divergence', true, ...
			       'value', false);
    spq_col = sp_evaluate_col (spq, msh_col);

    A = A + op_div_v_q (spv_col, spq_col, msh_col);
  end

  if (nargout == 1)
    varargout{1} = A;
  elseif (nargout == 3)
    [rows, cols, vals] = find (A);
    varargout{1} = rows;
    varargout{2} = cols;
    varargout{3} = vals;
  end

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
%> see function mfcc.m from Slaneys Auditory Toolbox
function [H] = ToolMfccFb (iMagSpecLength, f_s)

    % initialization
    f_start         = 133.3333;

    iNumLinFilters  = 13;
    iNumLogFilters  = 27;
    iNumFilters     = iNumLinFilters + iNumLogFilters;

    linearSpacing   = 66.66666666;
    logSpacing      = 1.0711703;

    % compute band frequencies
    f = f_start + (0:iNumLinFilters-1)*linearSpacing;
    f(iNumLinFilters+1:iNumFilters+2) = ...
		      f(iNumLinFilters) * logSpacing.^(1:iNumLogFilters+2);

    % sanity check
    f = f (f<f_s/2);
    iNumFilters = length(f)-2;

    f_l = f(1:iNumFilters);
    f_c = f(2:iNumFilters+1);
    f_u = f(3:iNumFilters+2);

    % allocate memory for filters and set max amplitude
    H = zeros(iNumFilters,iMagSpecLength);
    afFilterMax = 2./(f_u-f_l);
    f_k = (0:iMagSpecLength-1)/(iMagSpecLength-1)*f_s/2;

    % compute the transfer functions
    for c = 1:iNumFilters
        H(c,:)      = ...
            (f_k > f_l(c) & f_k <= f_c(c)).* ...
            afFilterMax(c).*(f_k-f_l(c))/(f_c(c)-f_l(c)) + ...
            (f_k > f_c(c) & f_k < f_u(c)).* ...
            afFilterMax(c).*(f_u(c)-f_k)/(f_u(c)-f_c(c));
    end
end

author: alexanderlerch

repo: ACA-Code


text:
function X = applyOBM(x, OBM, N_frame, NFFT, NUMBAND)

x_hat = stdft(x, N_frame, N_frame/2, NFFT);
x_hat = x_hat(:, 1:(NFFT/2+1)).';

X = zeros(NUMBAND, size(x_hat, 2));
for i = 1:size(x_hat, 2)
    X(:, i)	= sqrt(OBM*abs(x_hat(:, i)).^2);
end

author: mpariente

repo: pystoi

save_path: github-repos/MATLAB/mpariente-pystoi


text:
function m = size(t,idx)
%SIZE Sparse tensor dimensions.
%
%   D = SIZE(T) returns the size of the tensor.
%
%   I = size(T,DIM) returns the sizes of the dimensions specified by DIM.
%
%   See also SPTENSOR, SPTENSOR/NDIMS.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if exist('idx','var')
    m = t.size(idx);
else
    m = t.size;
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
% FileSeries
% Version 1.50 02-Jun-2010
%
%   cdw               - Change directory (like CD), allowing wildcard (*).
%   lsw               - List directory (like LS), allowing wildcard (*).
%   rdir              - Recursive list directory.
%   rdelete           - Delete files recursively.
%   rrmdir            - Delete directories recursively.
%   renamefile        - Rename a series of files.
%   renumberfile      - Re-number the indices of a series of files
%   getfilenum        - Get the index of a series of files.


text:
function [dStates, contactForces] = dynamics_flight(States, Actuators, Parameters)
% function [dStates, Actuators] = dynamics_flight(States, Actuators, Parameters)
%
% Computer Generated File -- DO NOT EDIT
%
% This function was created by the function Write_ContinuousDynamics()
% 10-Dec-2013 19:40:47
%
% Dymanics Model: retractable double pendulum biped
% Motion Phase: Flight
%
% Matthew Kelly
% Cornell University
%

x = States(1,:); % (m) Foot One horizontal position
y = States(2,:); % (m) Foot One vertical position
th1 = States(3,:); % (rad) Leg One absolute angle
th2 = States(4,:); % (rad) Leg Two absolute angle
L1 = States(5,:); % (m) Leg One length
L2 = States(6,:); % (m) Leg Two length
dx = States(7,:); % (m/s) Foot One horizontal velocity
dy = States(8,:); % (m/s) Foot One vertical velocity
dth1 = States(9,:); % (rad/s) Leg One absolute angular rate
dth2 = States(10,:); % (rad/s) Leg Two absolute angular rate
dL1 = States(11,:); % (m/s) Leg One extension rate
dL2 = States(12,:); % (m/s) Leg Two extensioin rate

F1 = Actuators(1,:); % (N) Compresive axial force in Leg One
F2 = Actuators(2,:); % (N) Compresive axial force in Leg Two
T1 = 0; %Foot One not in contact with ground!
T2 = 0; %Foot Two not in contact with ground!
Thip = Actuators(5,:); % (Nm) Torque acting on Leg Two from Leg One

m1 = Parameters.m1; % (kg) Foot One mass
m2 = Parameters.m2; % (kg) Foot Two mass
M = Parameters.M; % (kg) Hip mass
g = Parameters.g; % (m/s^2) Gravity

% Constraints for this phase:
H1 = 0; %(N) Foot One, horizontal contact force
V1 = 0; %(N) Foot One, vertical contact force
H2 = 0; %(N) Foot Two, horizontal contact force
V2 = 0; %(N) Foot Two, vertical contact force

dStates = zeros(size(States));
dStates(1:6,:) = States((1+6):(6+6),:);
dStates(7,:) = (T1.*sin(th1) - Thip.*sin(th1) + H1.*L1 - F1.*L1.*cos(th1))./(L1.*m1);
dStates(8,:) = -(T1.*cos(th1) - L1.*V1 - Thip.*cos(th1) + L1.*g.*m1 + F1.*L1.*sin(th1))./(L1.*m1);
dStates(9,:) = -(L2.*Thip.*m1 - L2.*T1.*m1 + L1.*T2.*m1.*cos(th1 - th2) + L1.*Thip.*m1.*cos(th1 - th2) - L2.*M.*T1.*cos(th1).^2 + L2.*M.*Thip.*cos(th1).^2 - L2.*M.*T1.*sin(th1).^2 + L2.*M.*Thip.*sin(th1).^2 + L1.*L2.*M.*V1.*cos(th1) - H1.*L1.*L2.*M.*sin(th1) - F2.*L1.*L2.*m1.*sin(th1 - th2) + 2.*L1.*L2.*M.*dL1.*dth1.*m1)./(L1.^2.*L2.*M.*m1);
dStates(10,:) = (L1.*M.*T2.*m1 + L1.*M.*Thip.*m1 + L1.*T2.*m1.*m2.*cos(th1 - th2).^2 + L1.*Thip.*m1.*m2.*cos(th1 - th2).^2 + L1.*T2.*m1.*m2.*sin(th1 - th2).^2 + L1.*Thip.*m1.*m2.*sin(th1 - th2).^2 - L2.*T1.*m1.*m2.*cos(th1 - th2) + L2.*Thip.*m1.*m2.*cos(th1 - th2) - L2.*M.*T1.*m2.*cos(th1 - th2).*sin(th1).^2 + L2.*M.*Thip.*m2.*cos(th1 - th2).*sin(th1).^2 - F1.*L1.*L2.*m1.*m2.*sin(th1 - th2) + L2.*M.*T1.*m2.*cos(th1).*cos(th2) - L2.*M.*Thip.*m2.*cos(th1).*cos(th2) + L2.*M.*T1.*m2.*sin(th1).*sin(th2) - L2.*M.*Thip.*m2.*sin(th1).*sin(th2) - L1.*L2.*M.*V1.*m2.*cos(th2) + L1.*L2.*M.*V2.*m1.*cos(th2) + H1.*L1.*L2.*M.*m2.*sin(th2) - H2.*L1.*L2.*M.*m1.*sin(th2) - L2.*M.*T1.*m2.*cos(th1 - th2).*cos(th1).^2 + L2.*M.*Thip.*m2.*cos(th1 - th2).*cos(th1).^2 + L1.*L2.*M.*V1.*m2.*sin(th1 - th2).*sin(th1) - F1.*L1.*L2.*M.*m2.*sin(th1 - th2).*cos(th1).^2 - F1.*L1.*L2.*M.*m2.*sin(th1 - th2).*sin(th1).^2 - F1.*L1.*L2.*M.*m2.*cos(th1).*sin(th2) + F1.*L1.*L2.*M.*m2.*cos(th2).*sin(th1) - 2.*L1.*L2.*M.*dL2.*dth2.*m1.*m2 + L1.*L2.*M.*V1.*m2.*cos(th1 - th2).*cos(th1) - H1.*L1.*L2.*M.*m2.*cos(th1 - th2).*sin(th1) + H1.*L1.*L2.*M.*m2.*sin(th1 - th2).*cos(th1))./(L1.*L2.^2.*M.*m1.*m2);
dStates(11,:) = -(T2.*m1.*sin(th1 - th2) + Thip.*m1.*sin(th1 - th2) - F1.*L2.*m1 + H1.*L2.*M.*cos(th1) + L2.*M.*V1.*sin(th1) + F2.*L2.*m1.*cos(th1 - th2) - F1.*L2.*M.*cos(th1).^2 - F1.*L2.*M.*sin(th1).^2 - L1.*L2.*M.*dth1.^2.*m1)./(L2.*M.*m1);
dStates(12,:) = (T1.*m1.*m2.*sin(th1 - th2) - Thip.*m1.*m2.*sin(th1 - th2) + F2.*L1.*M.*m1 + M.*T1.*m2.*cos(th1).*sin(th2) - M.*T1.*m2.*cos(th2).*sin(th1) - M.*Thip.*m2.*cos(th1).*sin(th2) + M.*Thip.*m2.*cos(th2).*sin(th1) + F2.*L1.*m1.*m2.*cos(th1 - th2).^2 + F2.*L1.*m1.*m2.*sin(th1 - th2).^2 - H1.*L1.*M.*m2.*cos(th2) + H2.*L1.*M.*m1.*cos(th2) - L1.*M.*V1.*m2.*sin(th2) + L1.*M.*V2.*m1.*sin(th2) + M.*T1.*m2.*sin(th1 - th2).*cos(th1).^2 - M.*Thip.*m2.*sin(th1 - th2).*cos(th1).^2 + M.*T1.*m2.*sin(th1 - th2).*sin(th1).^2 - M.*Thip.*m2.*sin(th1 - th2).*sin(th1).^2 - F1.*L1.*m1.*m2.*cos(th1 - th2) + F1.*L1.*M.*m2.*cos(th1).*cos(th2) + F1.*L1.*M.*m2.*sin(th1).*sin(th2) + H1.*L1.*M.*m2.*cos(th1 - th2).*cos(th1) + L1.*M.*V1.*m2.*cos(th1 - th2).*sin(th1) - L1.*M.*V1.*m2.*sin(th1 - th2).*cos(th1) + H1.*L1.*M.*m2.*sin(th1 - th2).*sin(th1) + L1.*L2.*M.*dth2.^2.*m1.*m2 - F1.*L1.*M.*m2.*cos(th1 - th2).*cos(th1).^2 - F1.*L1.*M.*m2.*cos(th1 - th2).*sin(th1).^2)./(L1.*M.*m1.*m2);

% contactForces(1,:) == H1 == (N) Foot One, horizontal contact force
% contactForces(2,:) == V1 == (N) Foot One, vertical contact force
% contactForces(3,:) == H2 == (N) Foot Two, horizontal contact force
% contactForces(4,:) == V2 == (N) Foot Two, vertical contact force
contactForces = zeros(4,size(States,2));
contactForces(1,:) = zeros(1,size(States,2));
contactForces(2,:) = zeros(1,size(States,2));
contactForces(3,:) = zeros(1,size(States,2));
contactForces(4,:) = zeros(1,size(States,2));

end

author: MatthewPeterKelly

repo: dscTutorials


text:
%% The Class spinTensor
%
%%
%

author: mtex-toolbox

repo: mtex

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


text:
%it is a state obsever of v1
function [sys,x0,str,ts]=ctrl5_2_3ESO(t,x,u,flag)
switch flag,
    case 0,
        [sys,x0,str,ts]=mdlInit();
    case 1,
        sys=mdlDer(t,x,u);
    case 3,
        sys=mdlOutput(t,x,u);
    otherwise,
        sys=[];
end;
function [sys,x0,str,ts]= mdlInit()
size=simsizes;
size.NumContStates=3;
size.NumDiscStates=0;
size.NumOutputs=3;
size.NumInputs=2;
size.DirFeedthrough=0;
size.NumSampleTimes=1;
sys=simsizes(size);
x0=[0,0,0];
str=[];
ts=[0,0];

function sys=mdlDer(t,x,u)
global b0;
b=[50,800,10000];
e=x(1)-u(1);
fe=fal(e,0.5,0.01);
fe1=fal(e,0.25,0.01);
sys(1)=x(2)-b(1)*e;
sys(2)=x(3)-b(2)*fe+b0*u(2);
sys(3)=-b(3)*fe1;
function sys=mdlOutput(t,x,u)
sys=x;

author: TianfaYao

repo: ADRC

save_path: github-repos/MATLAB/TianfaYao-ADRC


text:
%% -------------------- only variances switch model -------------------- %%
%
% Only variances in ALL three equations switch

%% housekeeping
clear
close all
clc()

%% Create dataset
clc

do_plot=true;

scale=100;

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

varlist=fieldnames(varlist0);

%% set up Markov chains

markov_chains=struct('name','syncvol',...
    'number_of_states',3,...
    'controlled_parameters',{{'s'}},...
    'endogenous_probabilities',[],...
    'probability_parameters',[]);

%% Create the VAR

clc

nlags=2;

exog={};

constant=true;

panel=[];

sv0=svar(varlist,exog,nlags,constant,panel,markov_chains);

%% set up restrictions

% syntax is alag(eqtn,vname)
%-------------------------------
lin_restr={
    % first equation or "FFR" equation
    %----------------------------------
    'a1(1,pi)=0'
    'a2(1,pi)=0'
    'a1(1,ygap)=0'
    'a2(1,ygap)=0'
    'a2(1,FFR)=0'
    % second equation or "pi" equation
    %----------------------------------
    'a0(2,FFR)=0'
    'a1(2,FFR)=0'
    'a2(2,FFR)=0'
    'a1(2,ygap)=0'
    'a2(2,ygap)=0'
    % third equation or "ygap" equation
    %-----------------------------------
    'a1(3,FFR)=0'
    'a2(3,FFR)=0'
    'a1(3,pi)=0'
    'a2(3,pi)=0'
    'a0(3,pi)+a0(3,FFR)=0'
    };
nonlin_restr={
    'a0(3,FFR)>=0'
    'a1(1,FFR)>=0'
    'a1(1,FFR)<=1'
    };

restrictions=[lin_restr;nonlin_restr];

%% set priors

% priors for the VAR coefficients
%--------------------------------
var_prior=svar.prior_template();
var_prior.type='sz';
% priors for the syncvol transition probabilities
%------------------------------------------------
switch_prior=struct();
switch_prior.dirichlet_1={0.1,'syncvol_tp_1_2',0.2,'syncvol_tp_1_3',0.2};
switch_prior.dirichlet_2={0.1,'syncvol_tp_2_1',0.2,'syncvol_tp_2_3',0.2};
switch_prior.dirichlet_3={0.1,'syncvol_tp_3_1',0.2,'syncvol_tp_3_2',0.2};

prior=struct();

prior.var=var_prior;

prior.nonvar=switch_prior;

%% Find posterior mode
clc

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

%% Find posterior mode : known regimes
clc

db2=db;

nobs=get(db.FFR,'NumberOfObservations');

h=sv0.nregs;

db2.hist_regimes=ts(get(db.FFR,'start'),randi(h,nobs,1));

is_fixed_regime=true;

sv2=estimate(sv0,db2,{'1960Q1','2015Q2'},prior,restrictions,'fmincon',is_fixed_regime);

%% estimates
clc

pmode=posterior_mode(sv)

%% Printing estimates
clc

print_structural_form(sv)

%% Printing solution
clc

print_solution(sv)

%% plot smoothed state and regime probabilities
clc

close all

plot_probabilities(sv)

%% plots probabilities against data
close all

plot_data_against_probabilities(sv,'regime')

%% Impulse responses

myirfs=irf(sv);

snames=fieldnames(myirfs);

figure('name','Simple impulse response functions')

iter=0;

jter=0;

for ishk=1:numel(snames)

    shk=snames{ishk};

    for iv=1:numel(varlist)

        v=varlist{iv};

        iter=iter+1;

        subplot(3,3,iter)

        plot(myirfs.(shk).(v),'linewidth',2)

        if iter<=3,title(v),end

        if any(iter==[1,4,7])

            jter=jter+1;

            ylabel([shk,'(',varlist{jter},' shock)'])

        end

    end

end

%% Generalized impulse responses

shock_names=[];

irf_periods=[];

params=[];

girf_setup=struct();% girf_setup=struct('nsims',300);

mygirfs=irf(sv,shock_names,irf_periods,params,girf_setup);

figure('name','Generalized impulse response functions')

iter=0;

jter=0;

for ishk=1:numel(snames)

    shk=snames{ishk};

    for iv=1:numel(varlist)

        v=varlist{iv};

        iter=iter+1;

        subplot(3,3,iter)

        plot(mygirfs.(shk).(v),'linewidth',2)

        if iter<=3,title(v),end

        if any(iter==[1,4,7])

            jter=jter+1;

            ylabel([shk,'(',varlist{jter},' shock)'])

        end

    end

end

%% Posterior sampling

[ff,lb,ub,x0,vcov,self]=pull_objective(sv);

options=struct();
options.alpha=0.234;
options.thin=10;
options.burnin=10^3;
options.N=2*10^4;
options.nchain=2;

lb(~isfinite(lb))=-500;
ub(~isfinite(ub))=500;

results=mh_sampler(ff,lb,ub,options,x0,vcov);

%% MCMC diagnostics
clc

pnames=fieldnames(pmode);

a2tilde_to_a=sv.estim_.linres.a2tilde_to_a;

mcmcobj=mcmc(results.pop,pnames,0,1,1,a2tilde_to_a);

myList={'syncvol_tp_1_2','syncvol_tp_1_3','syncvol_tp_2_1','syncvol_tp_2_3',...
    'syncvol_tp_3_1','syncvol_tp_3_2','s_1_1_syncvol_1','s_1_1_syncvol_2',...
    's_1_1_syncvol_3','s_2_2_syncvol_1','s_2_2_syncvol_2','s_2_2_syncvol_3',...
    's_3_3_syncvol_1','s_3_3_syncvol_2','s_3_3_syncvol_3'};

chain_id=[];

for ii=1:numel(myList)

    v=myList{ii};

    figure('name',['Univariate convergence diagnostics for ',v]);

    subplot(3,2,1),autocorrplot(mcmcobj,v,chain_id,40); title('autocorrelogram')

    subplot(3,2,2),densplot(mcmcobj,v,chain_id,250); title('density')

    subplot(3,2,3),meanplot(mcmcobj,v,chain_id); title('recursive mean')

    subplot(3,2,4),psrf_plot(mcmcobj,v); title('pot. scale red. factor')

    subplot(3,2,5),traceplot(mcmcobj,v,chain_id,20); title('trace')

    pause

    close

end

%% Marginal data density
clc

mddobj=mdd(results,ff,lb,ub);

fprintf('Importance sampling::%0.2f\n',is(mddobj,[],mdd.global_options))

fprintf('Reciprocal Importance sampling::%0.2f\n',ris(mddobj,[],mdd.global_options))

fprintf('Meng and Wong''s bridge::%0.2f\n',bridge(mddobj,true,mdd.global_options))

fprintf('Ulrich Mueller::%0.2f\n',mueller(mddobj,[],mdd.global_options))

fprintf('Laplace::%0.2f\n',laplace(mddobj))

fprintf('Sims, Waggoner and Zha::%0.2f\n',swz(mddobj,[],mdd.global_options))

fprintf('Laplace MCMC::%0.2f\n',laplace_mcmc(mddobj))

fprintf('Chib and Jeliazkov::%0.2f\n',cj(mddobj,[],mdd.global_options))

%% Out-of sample forecasts at the mode
shock_uncertainty=false;

nsteps=12;

mycast=forecast(sv,[],[],[],nsteps,shock_uncertainty);

do_plot_unconditional_forecasts(mycast,sv)

%% Conditional forecast on ygap

% conditional information
%-------------------------
ygap=scale*(-0.025:0.005:8*0.005).';
cond_db=struct('ygap',ts('2015Q3',ygap));

% options for the exercise
%--------------------------
myoptions=struct('cbands',[10,20,50,80,90],'do_plot',true,'nsteps',20,...
    'param_uncertainty',false,'shock_uncertainty',true,'ndraws',200);

% do it
%-------
tic
[fkst,bands,hdl]=do_conditional_forecasts(mest,db,cond_db,Results.pop,myoptions);
fprintf('\n\n Computing conditional forecasts took %0.4f minutes\n\n',toc/60);
%% Median conditional forecasts
% we take advantage of the fact that in the bands above we specified 50 in
% the bands above
%--------------------------------------------------------------------------
figure('name','Median Conditional Forecasts');
nvars=mest.endogenous.number;
for ivar=1:nvars
    thisname=mest.endogenous.name{ivar};
    subplot(nvars,1,ivar)
    plot(bands.(thisname)('ci_50'),'linewidth',2)
    title(mest.endogenous.tex_name{ivar})
end

%% irfs distributions

clc

shock_names=[];

params=[results.pop.x];

myirfs2=irf(sv,shock_names,[],params);

%% plotting of fan chargs

for ireg=1:sv.nregs

    reg=sprintf('regime_%d',ireg);

    figure('name',['Simple impulse response functions(',reg,')'])

    iter=0;

    jter=0;

    for ishk=1:numel(snames)

        shk=snames{ishk};

        for iv=1:numel(varlist)

            v=varlist{iv};

            iter=iter+1;

            subplot(3,3,iter)

            ffr=fanchart(myirfs2.(shk).(v),[30,50,70,90]);

            plot_fanchart(ffr.(reg),'r')

            if iter<=3,title(v),end

            if any(iter==[1,4,7])

                jter=jter+1;

                ylabel([shk,'(',varlist{jter},' shock)'])

            end

        end

    end

end

author: jmaih


text:
%USING_HG2 Determine if the HG2 graphics pipeline is used
%
%   tf = using_hg2(fig)
%
%IN:
%   fig - handle to the figure in question.
%
%OUT:
%   tf - boolean indicating whether the HG2 graphics pipeline is being used
%        (true) or not (false).

function tf = using_hg2(fig)
try
    if nargin < 1,  fig = figure('visible','off');  end
    tf = ~verLessThan('matlab','8.4.0');

catch
    tf = false;
end
if nargin < 1,  delete(fig);  end
end

author: marianux

repo: ecg-kit

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


text:
function v3 = r8vec_cross_product_3d ( v1, v2 )

%*****************************************************************************80
%
%% R8VEC_CROSS_PRODUCT_3D computes the cross product of two R8VEC's in 3D.
%
%  Discussion:
%
%    The cross product in 3D can be regarded as the determinant of the
%    symbolic matrix:
%
%          |  i  j  k |
%      det | x1 y1 z1 |
%          | x2 y2 z2 |
%
%      = ( y1 * z2 - z1 * y2 ) * i
%      + ( z1 * x2 - x1 * z2 ) * j
%      + ( x1 * y2 - y1 * x2 ) * k
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 October 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real V1(3), V2(3), the two vectors.
%
%    Output, real V3(3), the cross product vector.
%
  v3 = zeros ( 3, 1 );

  v3(1) = v1(2) * v2(3) - v1(3) * v2(2);
  v3(2) = v1(3) * v2(1) - v1(1) * v2(3);
  v3(3) = v1(1) * v2(2) - v1(2) * v2(1);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function a = fd1d_heat_implicit_matrix ( x_num, cfl )

%*****************************************************************************80
%
%% FD1D_HEAT_IMPLICIT_MATRIX: set the system matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2012
%
%  Author:
%
%    John Burkardt
%
  a = sparse ( [], [], [], x_num, x_num );

  a(1,1) = 1.0;

  for i = 2 : x_num - 1
    a(i,i-1) =           - cfl;
    a(i,i  ) = 1.0 + 2.0 * cfl;
    a(i,i+1) =           - cfl;
  end

  a(x_num,x_num) = 1.0;

  return
end


text:
function title = p20_title ( )

%*****************************************************************************80
%
%% P20_TITLE returns a title for problem 20.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'The Gregory and Karney Tridiagonal Matrix Function';

  return
end


text:
function [x_new,P_new,w_new]= gauss_prune(x,P,w,elim_threshold)
% gauss_prune Prune mixture components whose weights fall below a given
%   threshold
%
% TODO:
% ----
% * Generate documentation
%
% Credit:
% -------
% * This code has been taken from http://ba-tuong.vo-au.com/codes.html.
% * Minor modifications may have been applied to the original version.

idx= find( w > elim_threshold );
w_new= w(idx);
x_new= x(:,idx);
P_new= P(:,:,idx);

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function path=VMR2mrGrayNew(mmPerPix,VMRfile);
% function path=VMR2mrGrayNew(mmPerPix, VMRfile);
% Function to convert between Brainvoyager VMR file format and mrGray
% ARW 100300
% 040202: Note that BV flips L and R when it reads in anatomies.

if ~exist('mmPerPix','var')
   mmPerPix = [240/256 240/256 1.2];
   disp(['mmPerPix defaulting to [ ' num2str(mmPerPix,'%.4f ') ...
         	'].  I hope this is correct!']);
end

if ~exist('VMRfile','var')
    [VMRfile,VMRpath]=uigetfile('*.vmr');
end
VMRfile=[VMRpath,VMRfile];
disp(['Loading ',VMRfile]);

fid=fopen(VMRfile,'r');
header=fread(fid,3,'int16');
mainImg=fread(fid,'uchar');
fclose(fid);
mainImg=reshape(mainImg,[header(1),header(2),header(3)]);
disp('Flipping images...');

for thisIm=1:header(3)
   mainImg(:,:,thisIm)=rot90(fliplr(mainImg(:,:,thisIm)));
end
path = writeVolAnatShortHeader(mainImg, mmPerPix);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [h,hc,hl] = plotEllipse(cRow,cCol,ra,rb,phi,color,nPnts,lw,ls)
% Adds an ellipse to the current plot.
%
% USAGE
%  [h,hc,hl] = plotEllipse(cRow,cCol,ra,rb,phi,[color],[nPnts],[lw],[ls])
%
% INPUTS
%  cRow    - the row location of the center of the ellipse
%  cCol    - the column location of the center of the ellipse
%  ra      - semi-major axis radius length (in pixels) of the ellipse
%  rb      - semi-minor axis radius length (in pixels) of the ellipse
%  phi     - rotation angle (radians) of semimajor axis from x-axis
%  color   - ['b'] color for ellipse
%  nPnts   - [100] number of points used to draw each ellipse
%  lw      - [1] line width
%  ls      - ['-'] line style
%
% OUTPUT
%  h       - handle to ellipse
%  hc      - handle to ellipse center
%  hl      - handle to ellipse orient
%
% EXAMPLE
%  plotEllipse( 3, 2, 1, 5, pi/6, 'g');
%
% See also plotGaussEllipses
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.42
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

narginchk( 5, 9 );
if(nargin<6 || isempty(color)); color = 'b'; end
if(nargin<7 || isempty(nPnts)); nPnts = 100; end
if(nargin<8 || isempty(nPnts)); lw = 1; end
if(nargin<9 || isempty(ls)); ls = '-'; end

% plot ellipse (rotate a scaled circle):
ts = linspace(-pi,pi,nPnts+1);  cts = cos(ts); sts = sin(ts);
h = plot( ra*cts*cos(-phi) + rb*sts*sin(-phi) + cCol, ...
  rb*sts*cos(-phi) - ra*cts*sin(-phi) + cRow );

% plot center point and line indicating orientation
washeld = ishold; if(~washeld); hold('on'); end
hc = plot( cCol, cRow, 'k+' );
hl = plot( [cCol cCol+cos(-phi)*ra], [cRow cRow-sin(-phi)*ra] );
set([h hc hl],'Color',color,'LineWidth',lw,'LineStyle',ls);
if(~washeld); hold('off'); end;

end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function x = bezierCurve(p,t,tSpan)
% x = bezierCurve(p,t,tSpan)
%
% This function evaluates a bezier curve, defined by the set of control
% points p, at each value in t.
%
% INPUTS:
%   p = [nCurve x nPoint] = control points
%   t = [1 x nTime] = times at which to evaluate curve
%   tSpan = [1 x 2] = tSpan(1) <= t <= tSpan(2)
%
% OUTPUTS:
%   x = [nCurve x nTime] = bezier curve, evaluated at t
%
% NOTES:
%   It is not advisable to use this function for high-order polynomials.
%
% See also: RATIONALBEZIERCURVE, FITBEZIERCURVE
%

[nCurve, nPoint] = size(p);
nTime = length(t);

t = (t-tSpan(1))/diff(tSpan);

x = zeros(nCurve,nTime);
n = nPoint - 1;
for i=0:n
    tt = (t.^i).*(1-t).^(n-i);
    binom = nchoosek(n,i);
    x = x + binom*p(:,i+1)*tt;
end

end


text:
% MAIN  --  Quad-Rotor  --  Minimal-Throttle trajectory
%
% Find the minimal throttle-squared trajectory to move the quad-rotor from an
% arbitrary state to the origin.
%
% NOTES:
%   X = [x;y;z;p;r;w] = [x pos, y pos, z pos, pitch att, roll att, yaw att] = configuration
%  dX = [dx;dy;dz;dp;dr;dw] = [x vel, y vel, z vel, pitch rate, roll rate, yaw rate] = rate
% ddX = [ddx;ddy;ddz;ddp;ddr;ddw] = acceleration
%

clc; clear;

addpath ../../ ./utilities ./test

% Define environmental and plant model params
[p] = loadPlant_QuadRotor3d();

% Trajectory Parameters:
duration = 3;

% Initial State:
X0 = [1;0;0;0;0;0];   % initial configuration
dX0 = zeros(6,1);     % initial rates
z0 = [X0; dX0];  % initial state

XF = [0;0;0;0;0;0];   % final configuration
dXF = zeros(6,1);     % final rates
zF = [XF; dXF];  % final state

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                     Set up function handles                             %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

problem.func.dynamics = @(t,z,u)( dynQuadRotor3d(z,u,p) );
problem.func.pathObj = @(t,z,u)( sum(u.^2,1) );  % Throttle-squared cost function

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                     Set up problem bounds                               %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

problem.bounds.initialTime.low = 0;
problem.bounds.initialTime.upp = 0;
problem.bounds.finalTime.low = duration;
problem.bounds.finalTime.upp = duration;

problem.bounds.initialState.low = z0;
problem.bounds.initialState.upp = z0;
problem.bounds.finalState.low = zF;
problem.bounds.finalState.upp = zF;

problem.bounds.control.low = -p.uMax*[1;1;1;1];
problem.bounds.control.upp = p.uMax*[1;1;1;1];

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                    Initial guess at trajectory                          %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

problem.guess.time = [0,duration];
problem.guess.state = [z0, zeros(12,1)];
problem.guess.control = ones(4,2);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                         Solver options                                  %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

problem.options.nlpOpt = optimset(...
    'Display','iter',...
    'MaxFunEvals',1e5);

problem.options.method = 'trapezoid';
problem.options.trapezoid.nGrid = 16;

% problem.options.method = 'hermiteSimpson';
% problem.options.hermiteSimpson.nSegment = 30;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                            Solve!                                       %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%

soln = optimTraj(problem);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                        Display Solution                                 %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
% Plot the solution:
plotQuadRotor3d(soln)

repo: OptimTraj


text:
function h = nifti(varargin)
% Create a NIFTI-1 object
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id: nifti.m 1143 2008-02-07 19:33:33Z spm $

switch nargin
case 0,
    org = niftistruc;
    hdr = [];
    for i=1:length(org),
        hdr.(org(i).label) = feval(org(i).dtype.conv,org(i).def);
    end;
    h = struct('hdr',hdr,'dat',[],'extras',struct);
    h = class(h,'nifti');

case 1
    if ischar(varargin{1})
        if size(varargin{1},1)>1,
            h = nifti(cellstr(varargin{1}));
            return;
        end;
        fname  = deblank(varargin{1});
        vol    = read_hdr(fname);
        extras = read_extras(fname);

        if ~isfield(vol.hdr,'magic'),
            vol.hdr = mayo2nifti1(vol.hdr);

            % For SPM99 compatibility
            if isfield(extras,'M') && ~isfield(extras,'mat'),
                 extras.mat = extras.M;
                 if spm_flip_analyze_images,
                     extras.mat = diag([-1 1 1 1])*extras.mat;
                 end;
            end;

            % Over-ride sform if a .mat file exists
            if isfield(extras,'mat') && size(extras.mat,3)>=1,
                mat            = extras.mat(:,:,1);
                mat1           = mat*[eye(4,3) [1 1 1 1]'];
                vol.hdr.srow_x = mat1(1,:);
                vol.hdr.srow_y = mat1(2,:);
                vol.hdr.srow_z = mat1(3,:);
                vol.hdr.sform_code = 2;
                vol.hdr.qform_code = 2;
                vol.hdr = encode_qform0(mat,vol.hdr);
            end;
        end;

        if isfield(extras,'M'), extras = rmfield(extras,'M'); end;
        if isfield(extras,'mat') && size(extras.mat,3)<=1,
            extras = rmfield(extras,'mat');
        end;

        dim   = double(vol.hdr.dim);
        dim   = dim(2:(dim(1)+1));
        dt    = double(vol.hdr.datatype);
        offs  = double(vol.hdr.vox_offset);

        if ~vol.hdr.scl_slope && ~vol.hdr.scl_inter,
            vol.hdr.scl_slope = 1;
        end;
        slope = double(vol.hdr.scl_slope);
        inter = double(vol.hdr.scl_inter);

        dat   = file_array(vol.iname,dim,[dt,vol.be],offs,slope,inter);
        h     = struct('hdr',vol.hdr,'dat',dat,'extras',extras);
        h     = class(h,'nifti');

    elseif isstruct(varargin{1})
        h     = class(varargin{1},'nifti');

    elseif iscell(varargin{1})
        fnames = varargin{1};
        h(numel(fnames)) = struct('hdr',[],'dat',[],'extras',struct);
        h     = class(h,'nifti');
        for i=1:numel(fnames),
            h(i) = nifti(fnames{i});
        end;

    else
        error('Dont know what to do yet.');
    end;
otherwise
    error('Dont know what to do yet');
end;
return;

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function bias = gaussian_bias(Nt,L)
% GAUSSIAN_BIAS Computes bias of gaussian entropy estimates
%
%   -----
%   USAGE
%   -----
%
%       bias = gaussian_bias(nt,L)
%
%   where NT is the number of trials used for the estimation of entropy and
%   L is the dimensionality of the response. NT can either be an integer
%   scalar or a column arra of integers: in this case the bias is computed
%   for each of the values in NT.
%
%   See also ENTROPY, INFORMATION

%   Copyright (C) 2009 Cesare Magri
%   Version: 1.0.0

% -------
% LICENSE
% -------
% This software is distributed free under the condition that:
%
% 1. it shall not be incorporated in software that is subsequently sold;
%
% 2. the authorship of the software shall be acknowledged and the following
%    article shall be properly cited in any publication that uses results
%    generated by the software:
%
%      Magri C, Whittingstall K, Singh V, Logothetis NK, Panzeri S: A
%      toolbox for the fast information analysis of multiple-site LFP, EEG
%      and spike train recordings. BMC Neuroscience 2009 10(1):81;
%
% 3.  this notice shall remain in place in each source file.

Nt = Nt-1;

NtMat = Nt(:, ones(L,1));

pvec = 0:L-1;

NtMat = NtMat - pvec(ones(length(Nt),1), :);

bias = sum(psi(NtMat./2), 2) - L.*log(Nt(:,1)./2);
bias = bias ./ (2*log(2));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip

path: github-repos/MATLAB/fieldtrip-fieldtrip/fieldtrip-c2039be598a02d86b39aae76bfa7aaa720f9801c/external/ibtb/Bias Corrections/Bias Corrections for Guassian Method/gaussian_bias.m


text:
function gp_rec = gpmc_unpak(gp_rec, w, varargin)
%GPMC_UNPAK Set parameters of GP-MC structure
%
%  Description
%    GP_REC = GPMC_UNPAK(GP, W)
%    Given Gaussian process GP and hyperparameter weight matrix W,
%    create GP-MC structure with hyperparameter values filled from W.
%
%  See also
%    GPMC_PAK, GP_PAK, GP_UNPAK, GP_MC
%
% Copyright (c) 2014 Ville Tolvanen

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if ~isstruct(gp_rec)
  error('Input GP_REC must be either sampled GP (see gp_mc) or normal GP structure (see gp_set).');
end

ncf=length(gp_rec.cf);
nsamples=size(w,1);

% Set covariace function paramaters
ind=0;
if isfield(gp_rec, 'etr')
  gptmp=take_nth(gp_rec,1);
else
  gptmp=gp_rec;
end
for i1=1:ncf
  gpcf=gptmp.cf{i1};
  np=size(gpcf.fh.pak(gpcf),2);
  gpcftmp=gpcf;
  reccf=gpcftmp.fh.recappend(gpcftmp, gpcftmp);
  for i2=1:nsamples
    gpcftmp=gpcftmp.fh.unpak(gpcftmp, w(i2, ind+1:ind+np));
    reccf=gpcftmp.fh.recappend(reccf, i2, gpcftmp);
  end
  ind=ind+np;
  gp_rec.cf{i1}=reccf;
end
% Set likelihood function parameters
lik=gptmp.lik;
reclik=lik.fh.recappend(lik, lik);
if ind < size(w,2)
  np=size(lik.fh.pak(lik),2);
  for i2=1:nsamples
    lik=lik.fh.unpak(lik, w(i2, ind+1:ind+np));
    reclik=lik.fh.recappend(reclik, i2, lik);
  end
  gp_rec.lik=reclik;
end

if nargin==3
  gp_rec.latentValues=varargin{1};
end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function initialise(fa)
% Initialise file on disk
%
% This creates a file on disk with the appropriate size by explicitly
% writing data to prevent a sparse file.
%__________________________________________________________________________
% Copyright (C) 2013-2017 Wellcome Trust Centre for Neuroimaging

%
% $Id: initialise.m 7147 2017-08-03 14:07:01Z spm $

% first approach
% fa = subsasgn(fa, substruct('()',num2cell(size(fa))), 0);

% second approach
% fa = subsasgn(fa, substruct('()',repmat({':'},1,ndims(fa))), 0);

% third approach (problem if n > intmax('int32'))
% bs = 2^20;
% n  = prod(size(fa)); %#ok<PSIZE>
% fa = reshape(fa,n,1);
% for i=1:ceil(n/bs)
%    ii = ((((i-1)*bs)+1):min((i*bs),n))';
%    fa = subsasgn(fa,struct('type','()','subs',{{ii}}),zeros(numel(ii),1));
% end

d  = datatypes;
dt = find(cat(1,d.code)==fa.dtype);
if isempty(dt), error('Unknown datatype.'); end
d  = d(dt);
nbytes = d.nelem * d.size * prod(size(fa)); %#ok<PSIZE>
init(fa.fname, nbytes, struct('offset',fa.offset));

author: spm


text:
function y = sum( x, dim )

%   Disciplined convex/geometric programming information for SUM:
%      SUM(X) and SUM(X,DIM) is a vectorized version of addition. So
%      when SUM is used in a DCP or DGP, elements in each subvector must
%      satisfy the corresponding combination rules for addition (see
%      PLUS). For example, suppose that X looks like this:
%         X = [ convex concave affine  ;
%               affine concave concave ]
%      Then SUM(X,1) would be permittted, but SUM(X,2) would not,
%      because the top row contains the sum of convex and concave terms,
%      in violation of the DCP ruleset. For DGPs, addition rules dictate
%      that the elements of X must be log-convex or log-affine.

%
% Basic argument check
%

s = size( x );
switch nargin,
    case 0,
        error( 'Not enough input arguments.' );
    case 1,
        dim = [ find( s > 1 ), 1 ];
        dim = dim( 1 );
    case 2,
        if ~isnumeric( dim ) || dim <= 0 || dim ~= floor( dim ),
            error( 'Second argument must be a dimension.' );
        end
end

if dim > length( s ) || s( dim ) == 1,

    y = x;

elseif s( dim ) == 0,

    if ~any( s ), s = [1,1];
    else s( dim ) = 1; end
    y = cvx( s, sparse( 1, prod( s ) ) );

else

    p  = prod( s( 1 : dim - 1 ) );
    cc = 0 : prod( s ) - 1;
    cl = rem( cc, p );
    cr = floor( cc / ( p * s( dim ) ) );
    cc = cl + cr * p + 1;
    s( dim ) = 1;
    b = x.basis_;
    [ r, c, v ] = find( b );
    b = sparse( r, cc( c ), v, size( b, 1 ), prod( s ) );
    y = cvx( s, b );
    v = cvx_vexity( y );
    if any( isnan( v( : ) ) ),
        error( 'Disciplined convex programming error:\n   Illegal addition encountered (e.g., {convex} + {concave}).', 1 ); %#ok
    end

end

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
% example_radar.m
% ====================================================>
% The following seem to be descent parameters:
%   - ProbOfDetection = 0.8
%   - ProbOfDeath = 0.2
%   - VelocityErrVariance = 2
%   - ObsErrVariance = 50
%   - ProbOfConfirm = 0.9
%   - PHD.BirthScheme = {'Expansion', 5000}
%

%% Radar specific settings
RadarName = 'staddon';
switch(RadarName)
    case('staddon')
        % Surveillance region parameters
        V_bounds = [-8154.72944624983;... % X-min |
                    -212.289393440959;... % X-max | Surveillance region bounding
                    -7548.44272179096;... % Y-min | box coordinates (m)
                    4355.32645897434]';   % Y-max |
%         RadarCoords.lat = 50.33933333;
%         RadarCoords.lon = -4.12527778;
        RadarCoords.lat = 50.346069;
        RadarCoords.lon = -4.113670;
    case('longroom')
        % Surveillance region parameters
        V_bounds = [-8154.72944624983;... % X-min |
                    3000;...              % X-max | Surveillance region bounding
                    -7548.44272179096;... % Y-min | box coordinates (m)
                    4355.32645897434]';   % Y-max |
        RadarCoords.lat = 50.36286670714617;
        RadarCoords.lon = -4.156833300366998;
end

% Load dataset
load(strcat(RadarName,'.mat'));
N = size(DataList,2); % Simulation length

%% Focus on region of high clutter
lon_lim = [-4.195 -4.11];
lat_lim = [50.31 50.375];
% lon_lim = [-4.1538   -4.1152];
% lat_lim = [50.3097   50.3373];
% lon_lim = [-4.15 -4.13];
% lat_lim = [50.335 50.35];
% Convert LLA limits to NED
[y_lim, x_lim] = geodetic2ned(lat_lim,...
                              lon_lim, ...
                              0,...
                              RadarCoords.lat,...
                              RadarCoords.lon,...
                              0,...
                              referenceEllipsoid('wgs84'));
V_bounds = [x_lim(1);... % X-min |
            x_lim(2);... % X-max | Surveillance region bounding
            y_lim(1);... % Y-min | box coordinates (m)
            y_lim(2)]';  % Y-max |
for i = 1:N
    j = 1;
    while j <= size(DataList{i},2)
        if(DataList{i}(3,j)>lat_lim(2) || DataList{i}(3,j)<lat_lim(1) || DataList{i}(4,j)<lon_lim(1) || DataList{i}(4,j)>lon_lim(2))
            DataList{i}(:,j) = [];
            j = j -1;
        end
        j = j + 1;
    end
end

radar_meas_convert;
%% Plot & Recording settings
% Plot settings
ShowPlots = 1;              % Set to 0 to prevent showing any plots
ShowUpdate = 1;             % Set to 0 to skip showing update plots
ShowTrackInfo = 1;
NumPersistFrames = 50;

% Recording settings
Record = 1;                 % Set to (0|1) to turn video recording (off|on)
FrameRate = 0.5;            % Number of frames per second
VideoQuality = 100;         % Set to desired quality percentage
VideoPathName = strcat(RadarName,'_heavy_clutter_staddon.avi'); % Set to the desired path and name of produced recording

% Model parameter shortcuts
lambdaV = 10; % Expected number of clutter measurements over entire surveillance region
V = (abs(V_bounds(2)-V_bounds(1))*abs(V_bounds(4)-V_bounds(3))); % Total area of surveillance region
P_D = 0.7;    % Probability of detection
timestep_duration = duration(0,0,2);

%% Instantiation of necessary components

% Instantiate a Dynamic model (CV model with q = 1 m/s^2)
transition_model = ConstantVelocityX('VelocityErrVariance', 0.1,...
                                     'NumDims', 2,...
                                     'TimestepDuration', timestep_duration);
% Instantiate an Observation model (Variance of 50m^2 on each coordinate)
measurement_model = LinearGaussianX('NumMeasDims', 2,...
                                    'NumStateDims', 4,...
                                    'MeasurementErrVariance', 10^2,...
                                    'Mapping', [1 3]);
clutter_model = PoissonRateUniformPositionX('ClutterRate',lambdaV,...
                                            'Limits',[V_bounds(1:2);...
                                                      V_bounds(3:4)]);
detection_model = ConstantDetectionProbabilityX('DetectionProbability',P_D);

birth_model = DistributionBasedBirthModelX('Distribution', UniformDistributionX([V_bounds(1:2); ...
                                                                                [-5 5 ];...
                                                                                V_bounds(3:4);...
                                                                                [-5 5 ]]),...
                                           'BirthIntensity', 0.0001);
% Compile the State-Space model
model = StateSpaceModelX(transition_model, measurement_model,...
                        'Clutter',clutter_model,...
                        'Detection', detection_model,...
                        'Birth', birth_model);

%% Base Filter
obs_covar= measurement_model.covar();
PriorState = GaussianStateX(zeros(4,1), transition_model.covar() + blkdiag(obs_covar(1,1), 0, obs_covar(2,2),0));
base_filter = KalmanFilterX('Model', model, 'StatePrior', PriorState);

%% Data Associator
config.ClutterModel = clutter_model;
config.Clusterer = NaiveClustererX();
config.Gater = EllipsoidalGaterX(2,'GateLevel',10)';
config.DetectionModel = detection_model;
jpdaf = JointIntegratedProbabilisticDataAssocX(config);

%% Track Initiator

% Initiate Data Associator
config_phd.Model = model;
[priorParticles, priorWeights] = model.Birth.random(50000);
config_phd.StatePrior = ParticleStateX(priorParticles,10*priorWeights);
config_phd.BirthScheme = {'Expansion', 5000};
config_phd.SurvivalProbability = 0.995;

% Instantiate PHD filter
myphd = SMC_PHDFilterX(config_phd);

% Initiate Tag Generator
tag_gen = UuidTagGeneratorX();

% Prepare initiator parameters
config_ti.TagGenerator = tag_gen;
config_ti.InitFilter = base_filter;
config_ti.PhdFilter = myphd;
config_ti.ConfirmThreshold = 0.8;

% Create the track initiator
myti = PhdTrackInitiatorX(config_ti);

%% Track Deleter
config_td = struct('Fieldname', 'ExistenceProbability',...
                   'ReferenceValue', 0.1,...
                   'ReferenceOperand', 'lt');
mytd = FieldBasedDeleterX(config_td);

TrackList = [];

%% Create plot windows
if(ShowPlots)

    % Map plot
    figure('units','normalized','outerposition',[0 0 .5 1])
    ax(1) = gca;
    plot_google_map('Axis',ax(1),'APIKey','AIzaSyBXKujdtXRZiqya1soVS9pxBzYR4g7aGvM','Resize',3,'Scale',2,'MapType','satellite');
    axis(ax(1),[lon_lim(1) lon_lim(2) lat_lim(1) lat_lim(2)])
%     axis(ax(1),[-4.195 -4.11 50.31 50.375])

%     % PHD Intensity plot
%     figure('units','normalized','outerposition',[.5 0 .5 1])
%     ax(2) = gca;

    plots = [];
end

%% START OF SIMULATION
% ===================>
for k=2:N
    fprintf('Iteration = %d/%d\n================>\n',k,N);

    %% Extract DataList at time k
    MeasurementList = MeasurementScans(k);
    timestamp_km1 = MeasurementScans(k-1).Timestamp;
    timestamp_k = MeasurementList.Timestamp;
    dt = timestamp_k - timestamp_km1;
    transition_model.TimestepDuration = dt;
    fprintf('Timestamp = %s\n================>\n',timestamp_k);

    %% Process JPDAF
    jpdaf.MeasurementList = MeasurementList;
    jpdaf.TrackList = TrackList;
    jpdaf.predictTracks();
    jpdaf.associate(TrackList, MeasurementList);
    jpdaf.updateTracks();

    % Delete all plots (other than map)
    for i = 1:numel(plots)
        delete(plots(i))
    end
    plots = [];
    hold on;

%     if(MeasurementList.NumMeasurements>0)
%         % Convert measurements to LLA and plot them
%         [lat,lon,~] = ned2geodetic(MeasurementList.Vectors(2,:),...
%                                    MeasurementList.Vectors(1,:),...
%                                    0,...
%                                    RadarCoords.lat,...
%                                    RadarCoords.lon,...
%                                    0,...
%                                    referenceEllipsoid('wgs84'));
%         plots(end+1) = plot(ax(1), lon,lat,'y*','MarkerSize', 10);
%
%         for i=1:numel(lat)
%             plots(end+1) = text(ax(1), lon(:,i)+0.0001,lat(:,i)+0.00027,num2str(i),'FontSize',8,'Color','w');
%         end
% %             plot(ax(1), RadarCoords.lon,RadarCoords.lat,...
% %                                '-s','MarkerSize',20,...
% %                                'MarkerEdgeColor','red',...
% %                                'MarkerFaceColor',[1 .6 .6]);
%     end
%     pause(0.01);
     %% Perform Track initiation
    [TrackList] = myti.initiateTracks(jpdaf.TrackList, MeasurementList, jpdaf.AssocWeightsMatrix);

    %% Perform Track deletion
    TrackList = mytd.deleteTracks(TrackList);

    % Plot update step results
    if(ShowPlots && ShowUpdate)

        % Delete all plots (other than map)
        for i = 1:numel(plots)
            delete(plots(i))
        end
        plots = [];
        hold on;

        if(MeasurementList.NumMeasurements>0)
            % Convert measurements to LLA and plot them
            [lat,lon,~] = ned2geodetic(MeasurementList.Vectors(2,:),...
                                       MeasurementList.Vectors(1,:),...
                                       0,...
                                       RadarCoords.lat,...
                                       RadarCoords.lon,...
                                       0,...
                                       referenceEllipsoid('wgs84'));
            plots(end+1) = plot(ax(1), lon,lat,'y*','MarkerSize', 10);
%             plot(ax(1), RadarCoords.lon,RadarCoords.lat,...
%                                '-s','MarkerSize',20,...
%                                'MarkerEdgeColor','red',...
%                                'MarkerFaceColor',[1 .6 .6]);

        end

        % Plot all existing tracks
        for j=1:numel(TrackList)
            track = TrackList{j};
            % Convert track trajectory to LLA and plot it
            means = [track.Trajectory.Mean];
            [lat,lon,~] = ned2geodetic(means(3,:),...
                                       means(1,:),...
                                       0,...
                                       RadarCoords.lat,...
                                       RadarCoords.lon,...
                                       0,...
                                       referenceEllipsoid('wgs84'));
            traj_length = size(lon,2);
            if(traj_length>NumPersistFrames)
                start = traj_length-NumPersistFrames;
            else
                start = 1;
            end
            plots(end+1) = plot(ax(1), lon(:,start:end),lat(:,start:end),'-.w','LineWidth',2);
            plots(end+1) = plot(ax(1), lon(:,end),lat(:,end),'ws','MarkerSize',15);

            % Convert track velocity to LLA and plot it
            [lat_vel,lon_vel,~] = ned2geodetic(track.State.Mean(4,end),...
                                               track.State.Mean(2,end),...
                                               0,...
                                               lat(:,end),...
                                               lon(:,end),...
                                               0,...
                                               referenceEllipsoid('wgs84'));
            lat_vel = lat_vel-lat(:,end);
            lon_vel = lon_vel-lon(:,end);

            plots(end+1) = quiver(ax(1), lon(:,end),lat(:,end),20*lon_vel,20*lat_vel,'r','LineWidth',1.5);
            speed_kmph = sqrt(track.State.Mean(4,end)^2+track.State.Mean(2,end)^2)*3.6;
            speed_knot = speed_kmph/1.852;

            if(ShowTrackInfo)
                plots(end+1) = text(ax(1), lon(:,end)+0.001,lat(:,end)+0.00027,strcat("Sog:",num2str(speed_knot,2)," kt"),'FontSize',8,'Color','w');
                plots(end+1) = text(ax(1), lon(:,end)+0.001,lat(:,end)-0.00027,strcat("PoE:",num2str(track.ExistenceProbability*100,3)," %"),'FontSize',8,'Color','w');
            end
            % TODO: Convert heading and state covariance to LLA and plot them
            % [lat,lon,h] = ned2geodetic(North,East,0,50.346069,-4.113670,0,referenceEllipsoid('wgs84'));
            % h2 = plot_gaussian_ellipsoid(TrackList{j}.Filter.StateMean([1 3]), TrackList{j}.Filter.StateCovar([1 3],[1 3]),1,20,ax(1));
            % plots(end+1) = h2;
            % plots(end+1) = text(ax(1),TrackList{j}.Filter.StateMean(1)+20,TrackList{j}.Filter.StateMean(3)-5,int2str(TrackList{j}.TrackID));
            % plots(end+1) = text(ax(1),TrackList{j}.Filter.StateMean(1)+20,TrackList{j}.Filter.StateMean(3)-130,num2str(TrackList{j}.ProbOfExist,2));
        end

        % Add axis labels
        xlabel('Longitude')
        ylabel('Latitude')

        % Store video frame
        if(Record)
            F(k) = getframe(ax(1));
        end
    end

end

% Create video file and write to it
if(Record)
    %F = F(2:end);
    vidObj = VideoWriter(char(VideoPathName));
    vidObj.Quality = VideoQuality;
    vidObj.FrameRate = FrameRate;
    open(vidObj);
    writeVideo(vidObj, F);
    close(vidObj);
end


text:
function [x,w] = fnnls(XtX,Xty,tol)
%FNNLS	Non-negative least-squares.
%
% 	Adapted from NNLS of Mathworks, Inc.
%
%	x = fnnls(XtX,Xty) returns the vector X that solves x = pinv(XtX)*Xty
%	in a least squares sense, subject to x >= 0.
%	Differently stated it solves the problem min ||y - Xx|| if
%	XtX = X'*X and Xty = X'*y.
%
%	A default tolerance of TOL = MAX(SIZE(XtX)) * NORM(XtX,1) * EPS
%	is used for deciding when elements of x are less than zero.
%	This can be overridden with x = fnnls(XtX,Xty,TOL).
%
%	[x,w] = fnnls(XtX,Xty) also returns dual vector w where
%	w(i) < 0 where x(i) = 0 and w(i) = 0 where x(i) > 0.
%
%	See also NNLS and FNNLSb

%	L. Shure 5-8-87
%	Revised, 12-15-88,8-31-89 LS.
%	(Partly) Copyright (c) 1984-94 by The MathWorks, Inc.

%	Modified by R. Bro 5-7-96 according to
%       Bro R., de Jong S., Journal of Chemometrics, 1997, 11, 393-401
% 	Corresponds to the FNNLSa algorithm in the paper
%
%
%	Rasmus bro
%	Chemometrics Group, Food Technology
%	Dept. Dairy and Food Science
%	Royal Vet. & Agricultural
%	DK-1958 Frederiksberg C
%	Denmark
%	rb@kvl.dk
%	http://newton.foodsci.kvl.dk/rasmus.html

%  Reference:
%  Lawson and Hanson, "Solving Least Squares Problems", Prentice-Hall, 1974.

% Copyright (c) 1999, Rasmus Bro
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% initialize variables
if nargin < 3
    tol = 10*eps*norm(XtX,1)*max(size(XtX));
end
[m,n] = size(XtX);
P = zeros(1,n);
Z = 1:n;
x = P';
ZZ=Z;
w = Xty-XtX*x;

% set up iteration criterion
iter = 0;
itmax = 30*n;

% outer loop to put variables into set to hold positive coefficients
while any(Z) & any(w(ZZ) > tol)
    [wt,t] = max(w(ZZ));
    t = ZZ(t);
    P(1,t) = t;
    Z(t) = 0;
    PP = find(P);
    ZZ = find(Z);
    nzz = size(ZZ);
    z(PP')=(Xty(PP)'/XtX(PP,PP)');
    z(ZZ) = zeros(nzz(2),nzz(1))';
    z=z(:);
% inner loop to remove elements from the positive set which no longer belong

    while any((z(PP) <= tol)) & iter < itmax

        iter = iter + 1;
        QQ = find((z <= tol) & P');
        alpha = min(x(QQ)./(x(QQ) - z(QQ)));
        x = x + alpha*(z - x);
        ij = find(abs(x) < tol & P' ~= 0);
        Z(ij)=ij';
        P(ij)=zeros(1,max(size(ij)));
        PP = find(P);
        ZZ = find(Z);
        nzz = size(ZZ);
        z(PP)=(Xty(PP)'/XtX(PP,PP)');
        z(ZZ) = zeros(nzz(2),nzz(1));
        z=z(:);
    end
    x = z;
    w = Xty-XtX*x;
end

author: davidkun

repo: HyperSpectralToolbox


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
% ##2
%==============================================================================
%
% function [Tc,dT] = nnInter(T,omega,x,varargin);
%
% next neighbour interpolator for the data T given on a cell-cantered grid
% evaluated at x, see FAIR 3.2 p24.
%==============================================================================

function [Tc,dT] = nnInter(T,omega,x,varargin)

Tc = mfilename('fullpath'); dT = [];

if nargin == 0,

  testOneImgModel(mfilename);
  return;
elseif nargin == 1 && isempty(T),
  return;
end;

% flag for computing the derivative
doDerivative = (nargout>1);
for k=1:2:length(varargin), % overwrite default parameter
  eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

% get data size m, cell size h, dimension dim, and number n of interpolation points
dim = length(omega)/2;
m   = size(T);         if dim == 1, m = numel(T); end;
h   = (omega(2:2:end)-omega(1:2:end))./m;
n   = length(x)/dim;
dT  = sparse(n,dim*n);
x   = reshape(x,n,dim);
% map x from [h/2,omega-h/2] -> [1,m],
for i=1:dim, x(:,i) = (x(:,i)-omega(2*i-1))/h(i) + 0.5; end;

% round x
x = round(x);

Tc = zeros(n,1);                            % initialize output
Valid = @(j) (0<x(:,j) & x(:,j)<m(j)+1);    % determine indices of valid points
switch dim,
  case 1, valid = find( Valid(1) );
  case 2, valid = find( Valid(1) & Valid(2) );
  case 3, valid = find( Valid(1) & Valid(2) & Valid(3) );
end;

if isempty(valid),
  return;
end;

p = @(j) x(valid,j);
% increments for linearized ordering
i1 = 1; i2 = size(T,1); i3 = size(T,1)*size(T,2);
switch dim,
  case 1, p = p(1);
  case 2, p = p(1) + i2*(p(2) - 1);
  case 3, p = p(1) + i2*(p(2) - 1) + i3*(p(3) -1);
end;
Tc(valid) = T(p);
%==============================================================================


text:
function A = transfer2coeffs(H, freq, labelcmb, maxlag)

% TRANSFER2COEFFS converts a spectral transfer matrix into the time domain
% equivalent multivariate autoregressive coefficients up to a specified
% lag, starting from lag 1.

if nargin<3
  labelcmb = [];
end
if nargin<4
  maxlag = [];
end

% do a check on the input data
siz = size(H);
if numel(siz)==3 && siz(1)==siz(2)
  % assume chan_chan_freq
  isfull = true;
elseif numel(siz)==2
  % assume chancmb_freq
  isfull = false;
  %assert(~isempty(labelcmb), 'input data appears to be chancmb_freq, but labelcmb is missing');
else
  ft_error('dimensionality of input data is not supported');
end

dfreq = round(diff(freq)*1e5)./1e5; % allow for some numeric issues
if ~all(dfreq==dfreq(1))
  ft_error('FieldTrip:transfer2iis', 'frequency axis is not evenly spaced');
end

if freq(1)~=0
  ft_warning('FieldTrip:transfer2iis', 'when converting the transfer function to coefficients, the frequency axis should ideally start at 0, zero padding the spectral density');
  dfreq = mean(dfreq);
  npad  = freq(1)./dfreq;

  % update the freq axis and keep track of the frequency bins that are
  % expected in the output
  selfreq  = (1:numel(freq)) + npad;
  freq     = [(0:(npad-1))./dfreq freq];
  if isfull
    H = cat(3, zeros(siz(1),siz(2),npad), H);
  else
    H = cat(2, zeros(siz(1),npad), H);
  end
else
  selfreq  = 1:numel(freq);
end

% ensure H to be double precision
H = double(H);

% deal with the two different types of input
if isfull
  % check whether the last frequency bin is strictly real-valued.
  % if that's the case, then it is assumed to be the Nyquist frequency
  % and the two-sided spectral density will have an even number of
  % frequency bins. if not, in order to preserve hermitian symmetry,
  % the number of frequency bins needs to be odd.
  Hend = H(:,:,end);
  N    = numel(freq);
  m    = size(H,1);
  if all(imag(Hend(:))<abs(trace(Hend)./size(Hend,1)*1e-9))
    N2 = 2*(N-1);
  else
    N2 = 2*(N-1)+1;
  end

  % preallocate memory for efficiency
  Harr   = zeros(m,m,N2) + 1i.*zeros(m,m,N2);

  % the input cross-spectral density is assumed to be weighted with a
  % factor of 2 in all non-DC and Nyquist bins, therefore weight the
  % DC-bin with a factor of sqrt(2) to get a correct two-sided representation
  Harr(:,:,1) = H(:,:,1).*sqrt(2);
  for k = 2:N
    Harr(:,:,       k) = H(:,:,k);
    Harr(:,:,(N2+2)-k) = conj(H(:,:,k));
  end

  % the input cross-spectral density is assumed to be weighted with a
  % factor of 2 in all non-DC and Nyquist bins, therefore weight the
  % Nyquist bin with a factor of sqrt(2) to get a correct two-sided representation
  if mod(size(Harr,3),2)==0
    Harr(:,:,N) = Harr(:,:,N).*sqrt(2);
  end

  % invert the transfer matrix to get the fourier representation of the
  % coefficients, and add an identity matrix
  I = eye(siz(1));
  for k = 1:size(Harr,3)
    Harr(:,:,k) = I-inv(Harr(:,:,k));
  end

  % take the inverse fft to get the coefficients
  A = ifft(reshape(permute(Harr, [3 1 2]), N2, []), 'symmetric');
  A = A(2:end,:);
  A = ipermute(reshape(A, [N2-1 siz(1) siz(1)]), [3 1 2]);

  if ~isempty(maxlag)
    A = A(:,:,1:maxlag);
  end
else
  % check whether the last frequency bin is strictly real-valued.
  % if that's the case, then it is assumed to be the Nyquist frequency
  % and the two-sided spectral density will have an even number of
  % frequency bins. if not, in order to preserve hermitian symmetry,
  % the number of frequency bins needs to be odd.
  Hend = H(:,end);
  N    = numel(freq);
  m    = size(H,1);
  if all(imag(Hend(:))<max(abs(Hend))*1e-9)
    % the above heuristic may be a bit silly, FIXME
    N2 = 2*(N-1);
  else
    N2 = 2*(N-1)+1;
  end

  % preallocate memory for efficiency
  Harr   = zeros(m,N2) + 1i.*zeros(m,N2);

  % the input cross-spectral density is assumed to be weighted with a
  % factor of 2 in all non-DC and Nyquist bins, therefore weight the
  % DC-bin with a factor of sqrt(2) to get a correct two-sided representation
  Harr(:,1) = H(:,1).*sqrt(2);
  for k = 2:N
    Harr(:,       k) = H(:,k);
    Harr(:,(N2+2)-k) = conj(H(:,k));
  end

  % the input cross-spectral density is assumed to be weighted with a
  % factor of 2 in all non-DC and Nyquist bins, therefore weight the
  % Nyquist bin with a factor of sqrt(2) to get a correct two-sided representation
  if mod(size(Harr,3),2)==0
    Harr(:,N) = Harr(:,N).*sqrt(2);
  end

  % invert the transfer matrix to get the fourier representation of the
  % coefficients, and add an identity matrix
  %
  % this assumes Harr to be in the rows quadruplets of pairwise
  % decompositions, i.e. reshapable, without checking the labelcmb
  ncmb = size(Harr,1)./4;
  I = eye(2);
  for k = 1:N2
    Htmp = reshape(Harr(:,k), [2 2 ncmb]);
    Htmp = repmat(I, [1 1 ncmb]) - inv2x2(Htmp);
    Harr(:,k) = Htmp(:);
  end

  % take the inverse fft to get the coefficients
  A = ifft(permute(Harr, [2 1]), 'symmetric');
  A = A(2:end,:);
  A = ipermute(A, [2 1]);

  if ~isempty(maxlag)
    A = A(:,1:maxlag);
  end

end


text:
function D = spm_eeg_load(P)
% Load an M/EEG file in SPM format
% FORMAT D = spm_eeg_load(P)
%
% P        - filename of M/EEG file
% D        - MEEG object
%__________________________________________________________________________
%
% spm_eeg_load loads an M/EEG file using the SPM MEEG format. Importantly,
% the data array is memory-mapped and the struct is converted to MEEG object.
%__________________________________________________________________________
% Copyright (C) 2008-2012 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel
% $Id: spm_eeg_load.m 6102 2014-07-14 09:19:09Z vladimir $

%-Bypass if the input is already an MEEG object
%--------------------------------------------------------------------------
if nargin && isa(P, 'meeg')
    D = P;
    return;
end

%-Get filename
%--------------------------------------------------------------------------
if ~nargin
    [P, sts] = spm_select(1, 'mat', 'Select SPM M/EEG file');
    if ~sts, D = []; return; end
end

P = spm_file(spm_file(P, 'ext', '.mat'), 'cpath');

%-Load MAT file
%--------------------------------------------------------------------------
try
    load(P);
catch

    try
        load(spm_file(P,'filename'))

        warning('Ignoring path. The dataset was loaded from the current directory.');
    catch
        error('Cannot load file "%s".', P);
    end

end

%-Check whether there is a struct D
%--------------------------------------------------------------------------
if ~exist('D','var')
    error('File "%s" doesn''t contain SPM M/EEG data.', P);
end

%-Handle situations where the object has been directly saved in file
%--------------------------------------------------------------------------
if ~isa(D, 'struct')
    try
        D = struct(D);
    catch
        error('The file should contain an SPM M/EEG struct named D.');
    end
end

%-Save path and fname in structure
%--------------------------------------------------------------------------
D.path  = spm_file(P, 'path');
D.fname = spm_file(P, 'filename');

%-And return an MEEG object
%--------------------------------------------------------------------------
D = meeg(D);


text:
function [NeuronInformation, NeuronActiveMatrix] =GenerateNeuronInformation(Suite2P_file)

Suite2P_raw=load(Suite2P_file);
disp(['Suite2P result: ',Suite2P_file,' is loaded.']);

NeuronInformation=struct;
NeuronInformation.TotalROINumber=length(Suite2P_raw.stat(1,:));
NeuronInformation.FrameClockNumber=length(Suite2P_raw.F(1,:));
NeuronInformation.ROIStat=Suite2P_raw.stat;
NeuronInformation.CellPossibility=Suite2P_raw.iscell;
% NeuronInformation.TotalIsCellROINumber=length(Suite2P_raw.iscell(:,2)==1);
NeuronInformation.Options=Suite2P_raw.ops;
NeuronInformation.volumerate=NeuronInformation.Options.fs; %Hz
NeuronInformation.ImagingPlane=NeuronInformation.Options.nplanes;
NeuronInformation.ImagingChannel=NeuronInformation.Options.nchannels;

NeuronActiveMatrix=struct;
NeuronActiveMatrix.EventTrain=zeros(NeuronInformation.FrameClockNumber,2+NeuronInformation.TotalROINumber);
NeuronActiveMatrix.EventTrain(:,3:end)=Suite2P_raw.spks';
NeuronActiveMatrix.F_raw=zeros(NeuronInformation.FrameClockNumber,2+NeuronInformation.TotalROINumber);
NeuronActiveMatrix.F_raw(:,3:end)=Suite2P_raw.F';
NeuronActiveMatrix.F_neuropil=zeros(NeuronInformation.FrameClockNumber,2+NeuronInformation.TotalROINumber);
NeuronActiveMatrix.F_neuropil(:,3:end)=Suite2P_raw.Fneu';
NeuronActiveMatrix.EventTrain(:,1)=0:1/NeuronInformation.volumerate:(length(NeuronActiveMatrix.EventTrain(:,1))-1)/NeuronInformation.volumerate;
NeuronActiveMatrix.F_raw(:,1)=0:1/NeuronInformation.volumerate:(length(NeuronActiveMatrix.EventTrain(:,1))-1)/NeuronInformation.volumerate;
NeuronActiveMatrix.F_neuropil(:,1)=0:1/NeuronInformation.volumerate:(length(NeuronActiveMatrix.EventTrain(:,1))-1)/NeuronInformation.volumerate;

disp(['NeuronActiveMatrix for: ',Suite2P_file,' was generated.']);

end

author: kavli-ntnu


text:
function At = transpose(A)
At.func = A.transpose;
At.transpose = A.func;
At = class(At,'A_operator');

author: thomaskuestner


text:
function w = calculateKaiserwindow(N, beta)

denominator = calculateI0(beta);

val = 2 * [0 : (N - 1) / 2] / (N - 1) - 1;
w_half = calculateI0(beta * sqrt (1 - val .* val)) / denominator;

w = [w_half w_half(end - 1 : -1 : 1)].';


text:
function [new_data, lc, cs]= init_data(data,v,c)
% INPUT:
%data = data matrix
%c = column of the data matrix where class vector is given i.e. [6]
%v= columns of the data matrix where your data measurements are i.e. [1:5]
%
% OUTPUT
% new_data = your new data matrix ordered w.r.t. classes.
%lc = size of the class
%cs = starting points of the classes in your matrix new_data

data_v=data(:,v);
data_c=data(:,c);

data_v;
mins_v = min(data_v);
Ones = ones(size(data_v));
data_v = data_v+Ones*diag(abs(mins_v));
maxs_v = max(data_v);
data_v = data_v*diag(maxs_v.^(-1));

data = [data_v, data_c];

d_dim = size(data,1); % no of vectors
v_dim = size(data,2)-1; % dimension of vectors
data=sortrows(data, v_dim+1); % Rearrange according to classes

% rearrange according to class order 1,2,...
class_n=1;
temp = data(1, v_dim +1);
data(1, v_dim +1) = 1;

for i=2:d_dim
    if temp ~= data(i,v_dim +1);
       class_n=class_n+1;
       temp = data(i, v_dim +1);
    end
    data(i, v_dim +1) = class_n;
end

nc = class_n;
cs = []; % starting points of the classes
cs(1) = 1;

for j = 1 : nc % Beginning points and sizes of classes
    lc(j) = length(find(data(:, v_dim+1) == j)); %size of class 'j'.
    cs(j+1) = cs(j) + lc(j);
end
cs=cs(1:nc);
new_data = data;

author: Sable

repo: mcbench-benchmarks


text:
function a = circulant2 ( n )

%*****************************************************************************80
%
%% CIRCULANT2 returns the CIRCULANT2 matrix.
%
%  Formula:
%
%    K = 1 + mod ( J-I, N )
%    A(I,J) = K
%
%  Example:
%
%    N = 5
%
%    1  2  3  4  5
%    5  1  2  3  4
%    4  5  1  2  3
%    3  4  5  1  2
%    2  3  4  5  1
%
%  Properties:
%
%    A is generally not symmetric: A' /= A.
%
%    A is a circulant: each row is shifted once to get the next row.
%
%    A is Toeplitz: constant along diagonals.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    A commutes with any other circulant.
%
%    A is normal.
%
%    The inverse of A is a circulant matrix.
%
%    The eigenvector matrix is the Fourier matrix.
%
%    If R is an N-th root of unity, that is, R is a complex number such
%    that R**N = 1, then
%
%      Y = 1 + 2*R + 3*R**2 + ... + N*R**(N-1)
%
%    is an eigenvalue of A, with eigenvector
%
%     ( 1, R, R**2, ..., R**(N-1) )
%
%    and left eigenvector
%
%      ( R**(N-1), R**(N-2), ..., R**2, R, 1 ).
%
%    Although there are exactly N distinct roots of unity, the circulant
%    may have repeated eigenvalues, because of the behavior of the polynomial.
%    However, the matrix is guaranteed to have N linearly independent
%    eigenvectors.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Philip Davis,
%    Circulant Matrices,
%    John Wiley, 1979, QA188.D37.
%
%    Robert Gregory, David Karney,
%    A Collection of Matrices for Testing Computational Algorithms,
%    Wiley, New York, 1969, page 22,
%    LC: QA263.G68.
%
%    Morris Newman, John Todd,
%    The evaluation of matrix inversion programs,
%    Journal of the Society for Industrial and Applied Mathematics,
%    Volume 6, Number 4, pages 466-476, 1958.
%
%    Joan Westlake,
%    Test Matrix A24,
%    A Handbook of Numerical Matrix Inversion and Solution of Linear Equations,
%    John Wiley, 1968.
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  for i = 1 : n
    for j = 1 : n

      k = 1 + i4_modp ( j-i, n );
      a(i,j) = k;

    end
  end

  return
end


text:
function varargout = drawFaceNormals(varargin)
%DRAWFACENORMALS Draw normal vector of each face in a mesh.
%
%   drawFaceNormals(V, E, F)
%   Compute and draw the face normals of the mesh defined by vertices V,
%   edges E and faces F. See meshes3d for format of each argument.
%
%   H = drawFaceNormals(...)
%   Return handle array to the created objects.
%
%   Example
%   % draw face normals of a cube
%     [v, f] = createCubeOctahedron;
%     figure; drawMesh(v, f)
%     axis([-2 2 -2 2 -2 2]); axis equal; hold on;
%     drawFaceNormals(v, f)
%     view(3);
%
%   See also
%   meshes3d, drawMesh, drawVector3d, meshFaceNormals, meshFaceCentroids
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2010-10-06, using Matlab 7.9.0.529 (R2009b)
% Copyright 2010-2022 INRA - Cepia Software Platform

% extract vertices and faces
mesh = parseMeshData(varargin{:});

% compute vector data
c = meshFaceCentroids(mesh);
n = meshFaceNormals(mesh);

% display an arrow for each normal
h = quiver3(c(:,1), c(:,2), c(:,3), n(:,1), n(:,2), n(:,3));

% format output
if nargout > 0
    varargout{1} = h;
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function WriteOutFaceCheckersCNNbinary(locationTxt, faceCheckers)

    addpath('../PDM_helpers\');

    % use little-endian
    faceCheckerFile = fopen(locationTxt, 'w', 'l');

    views = numel(faceCheckers);

    % Type 0 - linear SVR, 1 - feed forward neural net, 2 - CNN, 3 - new
    % CNN
    fwrite(faceCheckerFile, 3, 'uint'); % 4 bytes

    % Number of face checkers
    fwrite(faceCheckerFile, views, 'uint'); % 4 bytes

    % Matrices representing view orientations
    for i=1:views
        % this indicates that we're writing a 3x1 double matrix
        writeMatrixBin(faceCheckerFile, faceCheckers(i).centres', 6);
    end

    for i = 1:views

        % The normalisation models
        % Mean of images
        writeMatrixBin(faceCheckerFile, faceCheckers(i).mean_ex, 6);

        % Standard deviation of images
        writeMatrixBin(faceCheckerFile, faceCheckers(i).std_ex, 6);

        cnn = faceCheckers(i).cnn;

        num_depth_layers = size(cnn.layers,2);

        % Get the number of layers
        fwrite(faceCheckerFile, num_depth_layers, 'uint'); % 4 bytes

        % For disambiguation between FC and conv layers
        res = vl_simplenn(cnn, single(faceCheckers(i).mask), [], []);

        for layers=1:num_depth_layers

            % write layer type: 0 - convolutional, 1 - max pooling (2x2 stride 2), 2 -
            % fully connected, 3 - relu, 4 - sigmoid
            if(cnn.layers{layers}.type == 'conv')

                % First check if it is an FC layer (they are represented
                % like that in matconvnet)
                if(numel(res(layers).x) == numel(cnn.layers{layers}.weights{1}(:,:,:,1)))
                    % This is the fully connected layer
                    fwrite(faceCheckerFile, 2, 'uint'); % 4 bytes

                    % the bias term
                    writeMatrixBin(faceCheckerFile, cnn.layers{layers}.weights{2}(:), 5);
                    % the weights

                    % Convert the filters to a matrix
                    weights_c = cnn.layers{layers}.weights{1};
                    size_w = size(weights_c);
                    weights = zeros(size_w(1)*size_w(2)*size_w(3), size_w(4));
                    weights(:) = weights_c;
                    writeMatrixBin(faceCheckerFile, weights, 5);
                else

                    % write the type (convolutional)
                    fwrite(faceCheckerFile, 0, 'uint'); % 4 bytes

                    num_in_map = size(cnn.layers{layers}.weights{1},3);

                    % write the number of input maps
                    fwrite(faceCheckerFile, num_in_map, 'uint'); % 4 bytes

                    num_out_kerns = size(cnn.layers{layers}.weights{1},4);

                    % write the number of kernels for each output map
                    fwrite(faceCheckerFile, num_out_kerns, 'uint'); % 4 bytes

                    % Write output map bias terms
                    for k2=1:num_out_kerns
                        fwrite(faceCheckerFile, cnn.layers{layers}.weights{2}(k2), 'float32'); % 4 bytes
                    end

                    for k=1:num_in_map
                        for k2=1:num_out_kerns
                            % Write out the kernel
                            W = squeeze(cnn.layers{layers}.weights{1}(:,:,k,k2));
                            writeMatrixBin(faceCheckerFile, W, 5);
                        end
                    end
                end
            elseif(cnn.layers{layers}.type == 'pool')
                fwrite(faceCheckerFile, 1, 'uint'); % 4 bytes, indicate max pooling layer, no params, assume (2x2 stride 2)
            elseif(cnn.layers{layers}.type == 'relu')
                fwrite(faceCheckerFile, 3, 'uint'); % 4 bytes, indicate relu layer, no params
            end
        end

        % Piecewise affine warp

        nPix = faceCheckers(i).nPix;
        minX = faceCheckers(i).minX;
        minY = faceCheckers(i).minY;

        destination = reshape(faceCheckers(i).destination, numel(faceCheckers(i).destination), 1);
        triangulation = faceCheckers(i).triangulation;
        triX = faceCheckers(i).triX;
        mask = faceCheckers(i).mask;
        alphas = faceCheckers(i).alphas;
        betas = faceCheckers(i).betas;

        fwrite(faceCheckerFile, nPix, 'uint'); % 4 bytes
        fwrite(faceCheckerFile, minX, 'float64'); % 8 bytes
        fwrite(faceCheckerFile, minY, 'float64'); % 8 bytes

        % Destination shape
        writeMatrixBin(faceCheckerFile, destination, 6);

        % Triangulation
        writeMatrixBin(faceCheckerFile, triangulation, 4);

        % Triangle map
        writeMatrixBin(faceCheckerFile, triX, 4);

        % Mask
        writeMatrixBin(faceCheckerFile, mask, 4);

        % Alphas
        writeMatrixBin(faceCheckerFile, alphas, 6);

        % Betas
        writeMatrixBin(faceCheckerFile, betas, 6);

    end

    fclose(faceCheckerFile);

end

author: TadasBaltrusaitis

repo: OpenFace


text:
%% Detection of ArUco Markers Demo
% Basic marker detection and pose estimation from single ArUco markers.
%
% Sources:
%
% * <https://docs.opencv.org/3.1.0/d5/dae/tutorial_aruco_detection.html>
% * <https://github.com/opencv/opencv_contrib/blob/3.1.0/modules/aruco/samples/detect_markers.cpp>
%

%% Parameters

% options
vidFile = '';              % Use video file instead of camera as input
markerLength = 0.1;        % Marker side length (in meters). Needed for correct scale in camera pose
dictionaryId = '6x6_250';  % Dictionary id
showRejected = false;      % Show rejected candidates too
estimatePose = false;      % Wheather to estimate pose or not
if estimatePose
    % Camera intrinsic parameters. Needed for camera pose
    load camera_parameters.mat -mat camMatrix distCoeffs
    %camMatrix = eye(3);
    %distCoeffs = zeros(1,5);
else
    camMatrix = [];
    distCoeffs = [];
end

% marker detector parameters
detectorParams = struct();
if false
    %detectorParams.nMarkers = 1024;
    detectorParams.adaptiveThreshWinSizeMin = 3;
    detectorParams.adaptiveThreshWinSizeMax = 23;
    detectorParams.adaptiveThreshWinSizeStep = 10;
    detectorParams.adaptiveThreshConstant = 7;
    detectorParams.minMarkerPerimeterRate = 0.03;
    detectorParams.maxMarkerPerimeterRate = 4.0;
    detectorParams.polygonalApproxAccuracyRate = 0.05;
    detectorParams.minCornerDistanceRate = 0.05;
    detectorParams.minDistanceToBorder = 3;
    detectorParams.minMarkerDistanceRate = 0.05;
    detectorParams.cornerRefinementMethod = 'None';
    detectorParams.cornerRefinementWinSize = 5;
    detectorParams.cornerRefinementMaxIterations = 30;
    detectorParams.cornerRefinementMinAccuracy = 0.1;
    detectorParams.markerBorderBits = 1;
    detectorParams.perspectiveRemovePixelPerCell = 8;
    detectorParams.perspectiveRemoveIgnoredMarginPerCell = 0.13;
    detectorParams.maxErroneousBitsInBorderRate = 0.04;
    detectorParams.minOtsuStdDev = 5.0;
    detectorParams.errorCorrectionRate = 0.6;
end
detectorParams.cornerRefinementMethod = 'Subpix';  % do corner refinement in markers

% dictionary
dictionary = {'Predefined', dictionaryId};

%% Input source
if ~isempty(vidFile) && exist(vidFile, 'file') == 2
    vid = cv.VideoCapture(vidFile);
    waitTime = 1;     % 1 sec
else
    vid = cv.VideoCapture(0);
    waitTime = 0.01;  % 10 msec
end
if ~vid.isOpened(), error('failed to initialize VideoCapture'); end

%% Main loop
totalTime = 0;
totalIterations = 0;
hImg = [];
while true
    % grab frame
    img = vid.read();
    if isempty(img), break; end

    tId = tic();

    % detect markers and estimate pose
    [corners, ids, rejected] = cv.detectMarkers(img, dictionary, ...
        'DetectorParameters',detectorParams);

    if estimatePose && ~isempty(ids)
        [rvecs, tvecs] = cv.estimatePoseSingleMarkers(corners, ...
            markerLength, camMatrix, distCoeffs);
    end

    % tic/toc
    currentTime = toc(tId);
    totalTime = totalTime + currentTime;
    totalIterations = totalIterations + 1;
    if mod(totalIterations, 30) == 0
        fprintf('Detection time = %f ms (Mean = %f ms)\n', ...
            1000*currentTime, 1000*totalTime/totalIterations);
    end

    % draw results
    if ~isempty(ids)
        img = cv.drawDetectedMarkers(img, corners, 'IDs',ids);

        if estimatePose
            for i=1:numel(ids)
                img = cv.drawAxis(img, camMatrix, distCoeffs, ...
                    rvecs{i}, tvecs{i}, markerLength * 0.5);
            end
        end
    end

    if showRejected && ~isempty(rejected)
        img = cv.drawDetectedMarkers(img, rejected, 'BorderColor',[255 0 100]);
    end

    if isempty(hImg)
        hImg = imshow(img);
    elseif ishghandle(hImg)
        set(hImg, 'CData',img);
    else
        break;
    end
    drawnow; pause(waitTime);
end
vid.release();

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [g]= spm_gx_adem_salience(x,v,a,P)
% returns the prediction for visual search (proprioception and vision)
% FORMAT [g]= spm_gx_adem_salience(x,v,a,P)
%
% x    - hidden states:
%   o(1) - oculomotor angle
%   o(2) - oculomotor angle
%
% v    - hidden causes
% P    - parameters
%
% g    - sensations:
%   g(1) - oculomotor angle (proprioception - x)
%   g(2) - oculomotor angle (proprioception - y)
%   g(3) - retinal input - channel 1
%   g(4) - retinal input - channel 2
%   g(5) - ...
%
% As for spm_dem_reach but with no visual target
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_gx_adem_salience.m 4580 2011-12-02 20:22:19Z karl $

% sensory input sampled from image
%--------------------------------------------------------------------------
global STIM

% retinotopic sampling
%--------------------------------------------------------------------------
s   = ADEM_sample_image(STIM.V,x,STIM.R);

% add proprioceptive angles in intrinsic (polar) coordinates
%--------------------------------------------------------------------------
g   = spm_vec(x,s);


text:
function [auc,sensitivity,specificity,accuracy] = predictionEstimationRF_LGG(X,Y,cat,nTrees,seed,testCost)
% -------------------------------------------------------------------------
% function [auc,sensitivity,specificity,accuracy] = predictionEstimation_RF(X,Y,cat,nTrees,seed,testCost)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function performs prediction estimation using out-of-bag estimates
% of random forests.
% -------------------------------------------------------------------------
% INPUTS:
% 1. X: Array of size [nPatient X nFeatures] in "Table" format.
%       --> Ex: See dataStruct.mat created in organizeRF_LGG.m
% 2. Y: Column vector of size [nPatient X 1] specifying the outcome status
%       (1 or 0) for all instances.
% 3. cat: Vecor of size [nFeatures X 1] specifying if each of the features
%         in 'X' are either continuous (0) or categorical (1).
% 4. nTrees: Number of decision-trees in the random forests.
%            --> Ex: 500
% 5. seed: Numerical number to use as seed for random-forest construction.
%          --> Ex: 54288
% 6. testCost: Vector of cost values (emphasis on instances of the minority
%              class) to be tested in the random forest construction. The
%              one providing the highest performance on out-of-bag
%              estimations will eventually be chosen.
% -------------------------------------------------------------------------
% OUTPUTS:
% 1. auc: Area under the receiving operating characteristic curve on
%         out-of-bag estimates.
% 2. sensitivity: Sensitivity of prediction on out-of-bag estimates.
% 3. specificity: Specificity of prediction on out-of-bag estimates.
% 4. accuracy: Accuracy of prediction on out-of-bag estimates.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2017
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015-2017  Martin Vallieres
%
%    This package is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

cost = findBestCost(nTrees,X,Y,cat,testCost,seed);
rng(seed), RF = TreeBagger(nTrees,X,Y,'OOBPrediction','on','Cost',[0,1/cost;1,0],'CategoricalPredictors',cat);
[prob] = oobPredict_table(RF,X); % Probability that y = 1 for all out-of-bag estimates
[auc,sensitivity,specificity,accuracy] = calcPerformMetrics(prob,Y,0.5);

end

author: mvallieres

repo: radiomics


text:
%% Magnetic Anisotropy
%
% by Dr. Marco Witte, Salzgitter Mannesmann Forschung, March 2020.
%
% The magnetocrystalline anisotropy energy Ea is necessary to rotate all
% magnetic moments towards the direction of an external field (saturation).
% This texture dependent propterty is of special interest for electrical
% steels.
%
% This example script is based on the followin literature:
% * Wu et al. (2017): Journal of Magnetism and Magnetic Materials 444:
% 211–217.
% * Bunge, H.-J. (2013): Texture analysis in materials science:
% mathematical methods.
% * Bunge, H. J. (1989): Texture and magnetic properties. Textures and
% microstructures, 11.
% * Szpunar, J. (1989): Texture, Stress, and Microstructure 11, 2-4:
% 93-105.
% * Landgraf et al.(2003): Journal of Magnetism and Magnetic Materials
% 254–255: 364–66.

%% Define Parameters

% alloy content of Si and Al in wt.%
x_Si = 3;
x_Al = 1;

% magnetic anisotropy constant
K1 = 4.77 - 0.21256*x_Si - 0.03816*x_Al;

% magnetic saturation J_S depends on alloy composition
J_S = 2.162 - 0.043*x_Si - 0.0625*x_Al;

% crystal symmetry
CS = crystalSymmetry('m-3m');
SS = specimenSymmetry('-1');

% define an odf - here with Goss orientation, as desired for grain oriented
% electrical steel
odf = unimodalODF(orientation('Euler',0,45*degree,0,CS,SS),CS,SS);

% simulate random orientations from the ODF
oris = discreteSample(odf,10000);

%% Calculate anisotropy energy for different external magnetic field directions in sheet plane
%
% Due to the rotation of an electrical motor the field direction is
% changing constantly inside the sheet plane.

% rotation from 0 - 90 degree
rot_mag = 0:10:90;

% initialize variable for storage
E_a = zeros(length(rot_mag),length(oris));

% loop over all rotations
for j = 1:length(rot_mag)

  %rotate orientations
  rot = rotation('axis',zvector,'angle',rot_mag(j)*degree);
  ori_rot = rot*oris;

  %Determine uvw of orientations (makes loop necessary)
  ori_uvw = inv(ori_rot)*Miller(1,0,0,crystalSymmetry('1'));

  %direction cosines with direction <100> of easy magnetization
  cos_val_1 = cos(angle(Miller(1,0,0,CS),ori_uvw,'noSymmetry'));
  cos_val_2 = cos(angle(Miller(0,1,0,CS),ori_uvw,'noSymmetry'));
  cos_val_3 = cos(angle(Miller(0,0,1,CS),ori_uvw,'noSymmetry'));

  %calculate magnetic anisotropy energy
  E_a(j,:) = K1*(cos_val_1.^2.*cos_val_2.^2 + ...
    cos_val_2.^2.*cos_val_3.^2 + cos_val_1.^2.*cos_val_3.^2);

end

% calculate magentization J_50 at H = 5000 A/m (J_50 is suposed to depend
% only on stexture
J_50 = J_S*(1-0.19*mean(E_a,2));

%plot results
figure
plot(rot_mag,mean(E_a,2))
ylabel('magnetic anisotropy energy [10^4J/m^3]')
xlabel('angle of external field to RD [°]')
grid on
title(['Mean magnetic anisotropy energy = ' num2str(mean(mean(E_a)),'%.2f') ' 10^4J/m^3.'],'fontsize',14)

figure
plot(rot_mag,mean(J_50,2))
ylabel('magnetic polarization [T]')
xlabel('angle of external field to RD [°]')
grid on
title(['Mean magnetic anisotropy energy = ' num2str(mean(mean(J_50)),'%.2f') ' T.'],'fontsize',14)

author: mtex-toolbox

repo: mtex


text:
function h = plotPub(X, Y, N, opt)
% function h = plotPub(X, Y, N, opt)
% This function plots X{i} vs Y{i}, changes the properties of the generated
% figure and exports it as a publication quality image file. The resolution
% of the image can be chosen by the user. Supported image formats are EPS,
% PDF, PNG, JPEG and TIFF. The figure properties are specified by the
% options structure 'opt'.
%
% Parameters:
% X:   cell array of x coordinates
% Y:   cell array of y coordinates
% N:   number of plots to be created. N <= length(X)
% opt: options structure:
%   BoxDim:       vector [width, height]: size of the axes box in inches; default: [6, 2.5]
%   ShowBox:      'on' = show or 'off' = hide bounding box; default: 'on'
%   FontName:     string: font name; default: 'Arial'
%   FontSize:     integer; default: 26
%   LineWidth:    vector [width1, width2, ..]: element i changes the property of i-th dataset; default: 2
%   LineStyle:    cell array {'style1', 'style2', ..}: element i changes the property of i-th dataset; default: '-'
%   Markers:      cell array {'marker1', 'marker2', ..}: element i changes the property of i-th dataset; default: 'None'
%   MarkerSpacing:vector [space1, space2, ..]: element i changes the property of i-th dataset; default: 0
%   Colors:       3xN matrix, [red, green, blue] where N is the number of datasets.
%   AxisColor:    [red, green, blue]; color of the axis lines; default: black
%   AxisLineWidth:Witdth of the axis lines; default: 2
%   XLabel:       X axis label
%   YLabel:       Y axis label
%   ZLabel:       Z axis label
%   Title:        Plot title
%   XTick:        [tick1, tick2, ..]: major ticks for X axis.
%   YTick:        [tick1, tick2, ..]: major ticks for Y axis.
%   ZTick:        [tick1, tick2, ..]: major ticks for Z axis.
%   XMinorTick:   'on' or 'off': show X minor tick?
%   YMinorTick:   'on' or 'off': show Y minor tick?
%   ZMinorTick:   'on' or 'off': show Z minor tick?
%   TickDir:      tick direction: 'in' or 'out'; default: 'in'
%   TickLength:   tick length; default: [0.02, 0.02]
%   XLim:         [min, max]: X axis limit.
%   YLim:         [min, max]: Y axis limit.
%   ZLim:         [min, max]: Z axis limit.
%   XScale:       'linear' or 'log': X axis scale.
%   YScale:       'linear' or 'log': Y axis scale.
%   ZScale:       'linear' or 'log': Z axis scale.
%   XGrid:        'on' or 'off': show grid in X axis?
%   YGrid:        'on' or 'off': show grid in Y axis?
%   ZGrid:        'on' or 'off': show grid in Z axis?
%   XDir:         'in' or 'out': X axis tick direction
%   YDir:         'in' or 'out': Y axis tick direction
%   ZDir:         'in' or 'out': Z axis tick direction
%   Legend:       {'legend1','legend2',...}
%   LegendBox:    bounding box of legend: 'on'/'off'; default: 'off'
%   LegendBoxColor: color of the bounding box of legend; default: 'none'
%   LegendTextColor: color of the legend text; default: [0,0,0]
%   LegendEdgeColor: color of the legend edges; default: [0,0,0]
%   LegendLoc:    'NorthEast', ..., 'SouthWest': legend location
%   Resolution:   Resolution (dpi) for bitmapped file. Default:600.
%   HoldLines:    true/false. true == only modify axes settings, do not touch plot lines/surfaces. Default false.
%   FileName:     Save? Give a file name.
%
% hfig:         Figure handle (optional). Default: current figure.
%
%
% Written by: K M Masum Habib (http://masumhabib.com)
% Copyright (c) K M Masum Habib 2012-2015.
%
% Distributed under the BSD License.
%
% Version: 2.2
%

if nargin ~= 4
    fprintf('Usage: function h = plotPub(X, Y, N, opt)');
end

% create figure window
hfig = figure;

% plot the functions and set style
for ii=1:N

    hold on;
    hp{ii} = plot(X{ii},Y{ii});
end
hold off;

% Set plot properties
setPlotProp(opt, hfig);

h = hfig;

author: masumhabib

repo: PlotPub

save_path: github-repos/MATLAB/masumhabib-PlotPub