text:
%% panel for updating neuron parameters
panel_neuron = uipanel(...
'Parent',main_fig,...
'FontUnits',get(0,'defaultuipanelFontUnits'),...
'Units','characters',...
'Title','Neuron',...
'Position',[5 29 67 11],...
'BackgroundColor',get(0,'defaultuipanelBackgroundColor'),...
'BackgroundColorMode',get(0,'defaultuipanelBackgroundColorMode'),...
'ParentMode','manual',...
'Tag','data_panel',...
'FontSize',18,...
'FontWeight','bold');

% neuron size
text_gsiz = uicontrol(...
'Parent',panel_neuron,...
'FontUnits',get(0,'defaultuicontrolFontUnits'),...
'Units','characters',...
'String','Size',...
'Style','pushbutton',...
'Position',[2 7 10 1.75],...
'BackgroundColor',get(0,'defaultuicontrolBackgroundColor'),...
'BackgroundColorMode',get(0,'defaultuicontrolBackgroundColorMode'),...
'Callback','%automatic',...
'Children',[],...
'ParentMode','manual',...
'Tag','data_dir',...
'FontSize',14, 'fontweight', 'bold');

edit_gsiz = uicontrol(...
'Parent',panel_neuron,...
'FontUnits',get(0,'defaultuicontrolFontUnits'),...
'Units','characters',...
'String', 15,...
'Style','edit',...
'Position',[13 7 5 1.75],...
'BackgroundColor',get(0,'defaultuicontrolBackgroundColor'),...
'BackgroundColorMode',get(0,'defaultuicontrolBackgroundColorMode'),...
'Callback','%automatic',...
'Children',[],...
'ParentMode','manual',...
'Tag','data_dir',...
'FontSize',14, 'fontweight', 'bold');

% gaussian width
text_gsig = uicontrol(...
'Parent',panel_neuron,...
'FontUnits',get(0,'defaultuicontrolFontUnits'),...
'Units','characters',...
'String','sig',...
'Style','pushbutton',...
'Position',[20 7 10  1.75],...
'BackgroundColor',get(0,'defaultuicontrolBackgroundColor'),...
'BackgroundColorMode',get(0,'defaultuicontrolBackgroundColorMode'),...
'Callback','%automatic',...
'Children',[],...
'ParentMode','manual',...
'Tag','data_dir',...
'FontSize',14, 'fontweight', 'bold');

edit_gsig = uicontrol(...
'Parent',panel_neuron,...
'FontUnits',get(0,'defaultuicontrolFontUnits'),...
'Units','characters',...
'String', 4,...
'Style','edit',...
'Position',[31 7 5 1.75],...
'BackgroundColor',get(0,'defaultuicontrolBackgroundColor'),...
'BackgroundColorMode',get(0,'defaultuicontrolBackgroundColorMode'),...
'Callback','%automatic',...
'Children',[],...
'ParentMode','manual',...
'Tag','data_dir',...
'FontSize',14, 'fontweight', 'bold', 'callback', 'edit_gsig_callback;');

%% searching method
text_method = uicontrol(...
'Parent',panel_neuron,...
'FontUnits',get(0,'defaultuicontrolFontUnits'),...
'Units','characters',...
'String','Method',...
'Style','text',...
'Position',[2 4 12 1.75],...
'BackgroundColor',get(0,'defaultuicontrolBackgroundColor'),...
'BackgroundColorMode',get(0,'defaultuicontrolBackgroundColorMode'),...
'Callback','%automatic',...
'Children',[],...
'ParentMode','manual',...
'Tag','data_dir',...
'FontSize',14, 'fontweight', 'bold');

radio_ellipse = uicontrol(...
'Parent',panel_neuron,...
'FontUnits',get(0,'defaultuicontrolFontUnits'),...
'Units','characters',...
'String','ellipse',...
'Style','radiobutton',...
'value', 1, ...
'Position',[15 5 12 1.75],...
'BackgroundColor',get(0,'defaultuicontrolBackgroundColor'),...
'BackgroundColorMode',get(0,'defaultuicontrolBackgroundColorMode'),...
'Callback','%automatic',...
'Children',[],...
'ParentMode','manual',...
'Tag','data_dir',...
'FontSize',14, ...
'callback', 'radio_ellipse_callback;');

text_ellipse = uicontrol(...
'Parent',panel_neuron,...
'FontUnits',get(0,'defaultuicontrolFontUnits'),...
'Units','characters',...
'String','ratio',...
'Style','pushbutton',...
'Position',[30 5 10 1.75],...
'BackgroundColor',get(0,'defaultuicontrolBackgroundColor'),...
'BackgroundColorMode',get(0,'defaultuicontrolBackgroundColorMode'),...
'Callback','%automatic',...
'Children',[],...
'ParentMode','manual',...
'Tag','data_dir',...
'FontSize',14);

edit_ellipse = uicontrol(...
'Parent',panel_neuron,...
'FontUnits',get(0,'defaultuicontrolFontUnits'),...
'Units','characters',...
'String',3,...
'Style','edit',...
'Position',[41 5 5 1.75],...
'BackgroundColor',get(0,'defaultuicontrolBackgroundColor'),...
'BackgroundColorMode',get(0,'defaultuicontrolBackgroundColorMode'),...
'Callback','%automatic',...
'Children',[],...
'ParentMode','manual',...
'Tag','data_dir',...
'FontSize',14);

 radio_dilate = uicontrol(...
'Parent',panel_neuron,...
'FontUnits',get(0,'defaultuicontrolFontUnits'),...
'Units','characters',...
'String','dilate',...
'Style','radiobutton',...
'value', 0, ...
'Position',[15 3 12 1.75],...
'BackgroundColor',get(0,'defaultuicontrolBackgroundColor'),...
'BackgroundColorMode',get(0,'defaultuicontrolBackgroundColorMode'),...
'Callback','%automatic',...
'Children',[],...
'ParentMode','manual',...
'Tag','data_dir',...
'FontSize',14,...
'callback', 'radio_dilate_callback;');

text_dilate = uicontrol(...
'Parent',panel_neuron,...
'FontUnits',get(0,'defaultuicontrolFontUnits'),...
'Units','characters',...
'String','radius',...
'Style','pushbutton',...
'Position',[30 3 10 1.75],...
'BackgroundColor',get(0,'defaultuicontrolBackgroundColor'),...
'BackgroundColorMode',get(0,'defaultuicontrolBackgroundColorMode'),...
'Callback','%automatic',...
'Children',[],...
'ParentMode','manual',...
'Tag','data_dir',...
'FontSize',14);

edit_dilate = uicontrol(...
'Parent',panel_neuron,...
'FontUnits',get(0,'defaultuicontrolFontUnits'),...
'Units','characters',...
'String',3,...
'Style','edit',...
'Position',[41 3 5 1.75],...
'BackgroundColor',get(0,'defaultuicontrolBackgroundColor'),...
'BackgroundColorMode',get(0,'defaultuicontrolBackgroundColorMode'),...
'Callback','%automatic',...
'Children',[],...
'ParentMode','manual',...
'Tag','data_dir',...
'FontSize',14);

author: zhoupc


text:
function areasmooth = LowPassFilterArea(area,filterorder,maxfreq)

f = fdesign.lowpass('N,F3db',filterorder,maxfreq);
h = design(f,'butter');
h.PersistentMemory = true;

h.filter(fliplr(area));
areasmooth = h.filter(area);

author: kristinbranson

repo: JAABA


text:
% eeg_emptyset() - Initialize an EEG dataset structure with default values.
%
% Usage:
%   >> EEG = eeg_emptyset();
%
% Outputs:
%   EEG    - empty dataset structure with default values.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: eeglab()

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

% 01-25-02 reformated help & license -ad

function EEG = eeg_emptyset();

EEG.setname     = '';
EEG.filename    = '';
EEG.filepath    = '';
EEG.subject     = '';
EEG.group       = '';
EEG.condition   = '';
EEG.session     = [];
EEG.comments    = '';
EEG.nbchan      = 0;
EEG.trials      = 0;
EEG.pnts        = 0;
EEG.srate       = 1;
EEG.xmin        = 0;
EEG.xmax        = 0;
EEG.times       = [];
EEG.data        = [];
EEG.icaact      = [];
EEG.icawinv     = [];
EEG.icasphere   = [];
EEG.icaweights  = [];
EEG.icachansind = [];
EEG.chanlocs    = [];
EEG.urchanlocs  = [];
EEG.chaninfo    = [];
EEG.ref         = [];
EEG.event       = [];
EEG.urevent     = [];
EEG.eventdescription = {};
EEG.epoch       = [];
EEG.epochdescription = {};
EEG.reject      = [];
EEG.stats       = [];
EEG.specdata    = [];
EEG.specicaact  = [];
EEG.splinefile  = '';
EEG.icasplinefile = '';
EEG.dipfit      = [];
EEG.history     = '';
EEG.saved       = 'no';
EEG.etc         = [];

%EEG.reject.threshold  = [1 0.8 0.85];
%EEG.reject.icareject  = [];
%EEG.reject.compreject = [];
%EEG.reject.gcompreject= [];
%EEG.reject.comptrial  = [];
%EEG.reject.sigreject  = [];
%EEG.reject.elecreject = [];

%EEG.stats.kurta      = [];
%EEG.stats.kurtr      = [];
%EEG.stats.kurtd      = [];
%EEG.stats.eegentropy = [];
%EEG.stats.eegkurt    = [];
%EEG.stats.eegkurtg   = [];
%EEG.stats.entropy    = [];
%EEG.stats.kurtc      = [];
%EEG.stats.kurtt      = [];
%EEG.stats.entropyc   = [];

return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function string=getProgressStr(obj,verbosity)
% give description of text outcome
%
% Inputs:
%   obj             MOxUnitTestOutcome instance.
%   verbosity       integer in range 0 to 2
%
% Output:
%   string          String containing a line (if verbosity==2); or a single
%                   character ('.','F','E','s'; if verbosity==1); or the
%                   empty string

    outcome=getOutcomeStr(obj, verbosity);

    if verbosity>=2
        test_node=getTest(obj);
        string=sprintf('%10s in % 7.2f s: %s\n', ...
                            outcome, getDuration(obj), ...
                            str(test_node));
    elseif verbosity>=1
        string=outcome;
    else
        string='';
    end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [R,current_eps]=DT_USTM(snrdB,T,L,Mt,Mr,epsilon,prec,filename)
%
% Function to compute the DT lower bound for a Rayleigh block-fading
% channel with no CSI at transmitter and receiver. The bound assumes that
% USTM is chosen as input distribution
% snrdB: SNR in dB
% Mt: number of transmit antennas
% Mr: number of receive antennas (Mt<=Mr);
% T: size of coherence interval
% L: number of independent coherence intervals; L*T is the blocklength
% epsilon: maximal block error probability
% prec: it controls the number of samples for the Monte Carlo simulation; Note nsamples=2^prec; One should have nsamples>> 100 x 1/epsilon
% filename: data file where the samples of the information density are saved for possible future refinements.
%
% The outputs of the program are
% R: R^*(T*L,current_eps,snrdB)
% current_eps: actual upper bound on the maximal error probability
%
%-------------------------------------------------------------------
%                       SET-UP PARAMETERS
%-------------------------------------------------------------------

SAVE=0; % if this flag is active the samples of the information density are saved for possibe further refinemes
MAT=0; % save with .mat extension

K = 2^prec; % number of monte carlo simulations (it should be at least 100 x 1/epsilon)
rho = 10.^(snrdB/10); % SNR in linear scale

%-------------------------------------------------------------------
%                       MONTE CARLO SIMULATION
%-------------------------------------------------------------------
I = zeros(K,1); %allocate for the montecarlo runs
%-------------------------------------------------------------------
%                       CONSTANTS
%-------------------------------------------------------------------
rho_tilde = T*rho/Mt;
D= [sqrt(1+rho_tilde)*eye(Mt), zeros(Mt,T-Mt);
    zeros(T-Mt,Mt), eye(T-Mt)]; % D matrix (covariance matrix of equivalent noise)

lambda=1+rho_tilde;
lambda1=1/lambda;
lambda2=rho_tilde*lambda1;
c2 = logComplexGammaRatio(Mt,Mr,T); %gamma constant
c1 = Mt*(T-Mt)*log(lambda2); % SNR constant
const = c1+c2;

noise_norm=sqrt(.5);

%-------------------------------------------------------------------
%                       MONTE CARLO
%-------------------------------------------------------------------

for k = 1:K %do K montecarlo runs
        i_L = 0;
        Z = (randn(T,Mr,L)+1i*randn(T,Mr,L))*noise_norm;
        for l = 1:L %Create each realization
            Sigma = svd(D*Z(:,:,l)).^2;
            if (Mt==1)
                M = createMalt(Mt,Mr,T,Sigma,lambda2); %create  matrix
                logdetM=log(det(M))+Sigma(1)*lambda2-(T-Mr)*log(Sigma(1));
                partial_sum=sum(Sigma)-logdetM;
            elseif Mr >= Mt
                M = createM_rx_larger(Mt,Mr,T,Sigma,lambda2);
                logdetM=logdet(M);
                logdetSigma = (T-Mr)*sum(log(Sigma));% compute logdet(Sigma^(T-M))
                partial_sum=lambda1*sum(Sigma) - logdetM +logdetSigma;
            else
                logdetM = createM_tx_larger(Mt, Mr, T, Sigma, lambda2); %already in log domain
                logdetSigma = (T-Mr)*sum(log(Sigma));% compute logdet(Sigma^(T-M))
                const2 = sum(gammaln(T-((Mt+1):T)+1)) - sum(gammaln(1:(T-Mr)));
                partial_sum = sum(Sigma) -logdetM + logdetSigma - const2;
            end

            vanderTerm = det(vander(Sigma)); %get the determinant of the vandermode matrix
            TraceZ=real(trace(Z(:,:,l)'*Z(:,:,l)));

            i_temp = const- TraceZ  +partial_sum + log(vanderTerm);  %Information density for time l (i(x_l, y_l)
            i_L = i_L + i_temp; %add it to the total i_L
        end

        I(k) =  i_L; %put all computations on a pile to compute the average later

end

if (SAVE==1)
  if (MAT==1)
    save(filename,'I')
  else
    save(filename,'I','-ascii','-append')
  end
end

%---------------------------------------
%   START SEARCHING FOR THE RATE
%---------------------------------------

I=sort(I);

Kcurrent=length(I); % redefine K to account for append

current_prec=floor(log2(Kcurrent)); % actual precision

K=2^(current_prec); % round off K to avoid search errors

step=K/2;
index=step;

onevec=ones(K,1);

while(step>1),

   th=I(index);

   current_eps=sum(exp(-max(0,I-th)))/K;

   step=step/2;

   if current_eps> epsilon,

       index=index-step;

   else

       index=index+step;

   end

end

current_eps=sum(exp(-max(0,I-I(index))))/K;

if (I(index)>500),

  R=I(index)/(L*T*log(2));

else

  R=log2(exp(I(index))+1)/(L*T); % factor 2 removed from DT bound to account for max error probability

end

end

function val = logComplexGammaRatio(Mt,Mr,T)
    k=T-Mt+1:1:T;
    val = -sum(gammaln(k));
    r=1:1:Mt;
    val=val+sum(gammaln(r));
end

function M = createMalt(Mt,Mr,T,Sigma,lambda)
M = nan(Mr,Mr); %(l is the row, k is the column)
for l = 1:Mr
    for k=1:Mt
        M(l,k)= exp((Mt-k)*log(Sigma(l)) + log(gammainc(lambda*Sigma(l),T+k-Mt-Mr)) + lambda*(Sigma(l)-Sigma(1)) + (T-Mr)*log(Sigma(1)/Sigma(l)));
    end

    for k=Mt+1:1:Mr
        M(l,k)= Sigma(l)^(Mr-k);
    end
end
end

function M = createM_rx_larger(Mt,Mr,T,Sigma,lambda)
M = nan(Mr,Mr); %(l is the row, k is the column)
P=Mr;

for l = 1:Mr
    for k=1:Mt
        M(l,k)=(Sigma(l)^(Mt-k))*gammainc(lambda*Sigma(l),T+k-Mt-P);
    end

    for k=Mt+1:Mr
        M(l,k)=exp((T-k)*log(Sigma(l)) -Sigma(l)*lambda); % case Mt<Mr
    end
end

end

function C = createM_tx_larger(M,N,T,Sigma,lambda)
inf_flag = 0;
A = zeros(max(M,N), max(M,N));

for i = 1:M
    for j = 1:M
        if j <= N
            if Sigma(j) == 0 && M-i == 0
                A(i,j) = exp(lambda*Sigma(j) + log(gammainc(Sigma(j)*lambda, T-N-M+i-1)));
            else
                if (T-N-M+i-1) > 0
                    A(i,j) = exp((M-i)*log(Sigma(j))  + lambda*Sigma(j) + log(gammainc(Sigma(j)*lambda, T-N-M+i-1)));
                else
                    A(i,j) = exp((M-i)*log(Sigma(j))  + lambda*Sigma(j));
                end
            end
            if isinf(A(i,j)) == 1
                inf_flag = 1;
                break;
            end
            %A(i,j) = exp((M-i)*log(lambda(j))  + gammaln(T-M) +  log(gammainc(lambda(j)*p, T+i-2*M)));
        else
            A(i,j) = exp((T-j-(M-i))*log(lambda) + sum(log(T-j-(0:(M-i-1)))));
            %A(i,j) = exp((T-j-(M-i))*log(p)  + sum(log(T-j-(0:(M-i-1)))));
        end
    end
    if inf_flag == 1
        break;
    end
end
if inf_flag == 1
    for i = 1:M
        for j = 1:M
            if j <= N
                if (T-N-M+i-1) > 0
                    A(i,j) = exp((M-i)*log(Sigma(j))  +  log(gammainc(Sigma(j)*lambda, (T-N-M+i-1))));
                else
                    A(i,j) = exp((M-i)*log(Sigma(j)));
                end
            else
                A(i,j) = exp((T-j-(M-i))*log(lambda)  + sum(log(T-j-(0:(M-i-1)))));
            end
        end
    end
end

if inf_flag == 0
    C = logdet(A);
else
    C = logdet(A) + lambda*sum(Sigma);
end
end

function v = logdet(A, op)

assert(isfloat(A) && ndims(A) == 2 && size(A,1) == size(A,2), ...
    'logdet:invalidarg', ...
    'A should be a square matrix of double or single class.');

if nargin < 2
    use_chol = 0;
else
    assert(strcmpi(op, 'chol'), ...
        'logdet:invalidarg', ...
        'The second argument can only be a string ''chol'' if it is specified.');
    use_chol = 1;
end

%% computation

if use_chol
    v = 2 * sum(log(diag(chol(A))));
else
    [L, U, P] = lu(A);
    du = diag(U);
    c = det(P) * prod(sign(du));
    v = log(c) + sum(log(abs(du)));
end
end

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function box_behnken_test ( )

%*****************************************************************************80
%
%% BOX_BEHNKEN_TEST tests the BOX_BEHNKEN library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 October 2008
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BOX_BEHNKEN_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the BOX_BEHNKEN library.\n' );

  box_behnken_test01 ( );
  box_behnken_test02 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BOX_BEHNKEN_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function newdata=struct2jdata(data,varargin)
%
% newdata=struct2jdata(data,opt,...)
%
% convert a JData object (in the form of a struct array) into an array
%
% authors:Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
%
% input:
%      data: a struct array. If data contains JData keywords in the first
%            level children, these fields are parsed and regrouped into a
%            data object (arrays, trees, graphs etc) based on JData
%            specification. The JData keywords are
%               "_ArrayType_", "_ArraySize_", "_ArrayData_"
%               "_ArrayIsSparse_", "_ArrayIsComplex_"
%      opt: (optional) a list of 'Param',value pairs for additional options
%           The supported options include
%               'Recursive', if set to 1, will apply the conversion to
%                            every child; 0 to disable
%
% output:
%      newdata: the covnerted data if the input data does contain a JData
%               structure; otherwise, the same as the input.
%
% examples:
%      obj=struct('_ArrayType_','double','_ArraySize_',[2 3],
%                 '_ArrayIsSparse_',1 ,'_ArrayData_',null);
%      ubjdata=struct2jdata(obj);
%
% license:
%     BSD License, see LICENSE_BSD.txt files for details
%
% -- this function is part of JSONLab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

fn=fieldnames(data);
newdata=data;
len=length(data);
if(jsonopt('Recursive',0,varargin{:})==1)
  for i=1:length(fn) % depth-first
    for j=1:len
        if(isstruct(getfield(data(j),fn{i})))
            newdata(j)=setfield(newdata(j),fn{i},jstruct2array(getfield(data(j),fn{i})));
        end
    end
  end
end
if strncmpi('x0x5F_ArrayType_',fn,16) && strncmpi('x0x5F_ArrayData_',fn,16)
  newdata=cell(len,1);
  for j=1:len
    ndata=cast(data(j).x0x5F_ArrayData_,data(j).x0x5F_ArrayType_);
    iscpx=0;
    if strncmpi('x0x5F_ArrayIsComplex_',fn,21)
        if(data(j).x0x5F_ArrayIsComplex_)
           iscpx=1;
        end
    end
    if strncmpi('x0x5F_ArrayIsSparse_',fn,20)
        if(data(j).x0x5F_ArrayIsSparse_)
            if strncmpi('x0x5F_ArraySize_',fn,16)
                dim=double(data(j).x0x5F_ArraySize_);
                if(iscpx && size(ndata,2)==4-any(dim==1))
                    ndata(:,end-1)=complex(ndata(:,end-1),ndata(:,end));
                end
                if isempty(ndata)
                    % All-zeros sparse
                    ndata=sparse(dim(1),prod(dim(2:end)));
                elseif dim(1)==1
                    % Sparse row vector
                    ndata=sparse(1,ndata(:,1),ndata(:,2),dim(1),prod(dim(2:end)));
                elseif dim(2)==1
                    % Sparse column vector
                    ndata=sparse(ndata(:,1),1,ndata(:,2),dim(1),prod(dim(2:end)));
                else
                    % Generic sparse array.
                    ndata=sparse(ndata(:,1),ndata(:,2),ndata(:,3),dim(1),prod(dim(2:end)));
                end
            else
                if(iscpx && size(ndata,2)==4)
                    ndata(:,3)=complex(ndata(:,3),ndata(:,4));
                end
                ndata=sparse(ndata(:,1),ndata(:,2),ndata(:,3));
            end
        end
    elseif strncmpi('x0x5F_ArraySize_',fn,16)
        if(iscpx && size(ndata,2)==2)
             ndata=complex(ndata(:,1),ndata(:,2));
        end
        ndata=reshape(ndata(:),data(j).x0x5F_ArraySize_);
    end
    newdata{j}=ndata;
  end
  if(len==1)
      newdata=newdata{1};
  end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
classdef ShapePrinter < CompositeResultsPrinter

    methods (Access = public)

        function obj = ShapePrinter(d)
            obj.init(d);
        end

    end

    methods (Access = protected)

        function storeFieldsToPrint(obj,d)
            for iF = 1:numel(d.fieldToPrint)
                d.fields = d.fieldToPrint{iF}.value;
                obj.printers{iF}.storeFieldsToPrint(d)
            end
        end

        function createHeadPrinter(obj,d,dh)
            for iP = 1:numel(obj.printers)
              p = obj.printers{iP};
              hasGauss(iP) = p.getHasGaussData;
            end
            hasGauss = find(hasGauss);
            firstPrinterWithGauss = obj.printers{hasGauss(1)};
            p = firstPrinterWithGauss;
            p.createHeadPrinter(d,dh);
            h = p.getHeadPrinter();
            obj.headPrinter = h;
        end

        function createPrinters(obj,d)
            for iF = 1:numel(d.fieldToPrint)
                type = d.fieldToPrint{iF}.type;
                d.name = d.fieldToPrint{iF}.name;
                obj.printers{iF} = ResultsPrinter.create(type,d);
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function test22(nmat)
%TEST22 test pos.def and indef. matrices
% Example:
% test22(nmat)
%
% if nmat <= 0, just test problematic matrices
% See also cholmod_test

% Copyright 2006-2007, Timothy A. Davis, University of Florida

fprintf ('=================================================================\n');
fprintf ('test22: test pos.def and indef. matrices\n') ;

index = UFget ;

[ignore f] = sort (index.nrows) ;

if (nargin > 0)
    problematic = (nmat <= 0) ;
    if (problematic)
	% Matrices for which MATLAB and CHOLMOD differ, for which the error
	% is high, or other issues arose during debugging
	fprintf ('testing matrices for which MATLAB and CHOLMOD differ\n') ;
	f = [
	    186 % HB/jgl011
	    109 % HB/curtis54
	    793 % Qaplib/lp_nug05
	    607 % LPnetlib/lp_brandy
	    707 % LPnetlib/lp_bore3d
	    231 % HB/plskz362
	    794 % Qaplib/lp_nug06
	    673 % LPnetlib/lp_scorpion
	   1156 % Sandia/oscil_dcop_45 (also 1157:1168)
	    795 % Qaplib/lp_nug07
	    796 % Qaplib/lp_nug08
	    260 % HB/well1033
	    261 % HB/well1850
	    230 % HB/plsk1919
	    649 % LPnetlib/lp_pds_02
	    660 % LPnetlib/lp_qap12
	    609 % LPnetlib/lp_cre_a
	    619 % LPnetlib/lp_dfl001
	    661 % LPnetlib/lp_qap15
	    650 % LPnetlib/lp_pds_06
	    379 % Cote/vibrobox
	    638 % LPnetlib/lp_ken_11
	    799 % Qaplib/lp_nug20
	    ]' ;
    else
	nmat = max (0,nmat) ;
	nmat = min (nmat, length (f)) ;
	f = f (1:nmat) ;
    end
end

skip = [
    811, ...		% Simon/appu, which is a random matrix
    937:939, ...	% large ND/ problems
    1157:1168 ...	% duplicates
    799			% rather large: Qaplib/lp_nug20
    ] ;

tlimit = 0.1 ;
fprintf ('test22: chol and chol2 are repeated so each take >= %g sec\n',tlimit);

klimit = 1 ;

% warmup, for more accurate timing
[R,p] = chol (sparse (1)) ;						    %#ok
[R,p] = chol2 (sparse (1)) ;						    %#ok
clear R p

for i = f

    if (any (i == skip))
	continue ;
    end

    Problem = UFget (i) ;
    A = Problem.A ;
    [m n] = size (A) ;
    fprintf ('\n================== %4d: Problem: %s  m: %d n: %d nnz: %d', ...
	i, Problem.name, m, n, nnz (A)) ;
    clear Problem ;

    try	% create a symmetric version of the matrix
	if (m == n)
	    if (nnz (A-A') > 0)
		A = A+A' ;
	    end
	else
	    A = A*A' ;
	end
    catch
	fprintf ('skip\n') ;
	continue
    end

    fprintf (' %d\n', nnz (A)) ;

    p = amd2 (A) ;
    A = A (p,p) ;
    anorm = norm (A,1) ;

    % Run each code for at least 'tlimit' seconds

    % MATLAB
    k = 0 ;
    t1 = 0 ;
    while (t1 < tlimit & k < klimit)					    %#ok
	tic ;
	[R1,p1] = chol (A) ;
	t = toc ;
	t1 = t1 + t ;
	k = k + 1 ;
    end
    t1 = t1 / k ;

    % CHOLMOD
    k = 0 ;
    t2 = 0 ;
    while (t2 < tlimit & k < klimit)					    %#ok
	tic ;
	[R2,p2] = chol2 (A) ;
	t = toc ;
	t2 = t2 + t ;
	k = k + 1 ;
    end
    t2 = t2 / k ;

    if (klimit == 1)
	rmin = full (min (abs (diag (R2)))) ;
	rmax = full (max (abs (diag (R2)))) ;
	if (p2 ~= 0 | isnan (rmin) | isnan (rmax) | rmax == 0)		    %#ok
	    rcond = 0 ;
	else
	    rcond = rmin / rmax ;
	end
	fprintf ('rcond: %30.20e\n', rcond) ;
    end

    if (p1 == 1)
	% MATLAB does not follow its own definitions.  If p is 1, then R is
	% supposed to be 0-by-n, not 0-by-0.  CHOLMOD fixes this bug.
	% Here, A is m-by-m
	R1 = sparse (0,m) ;
    end

    kerr = 0 ;
    if (p1 ~= p2)
	% MATLAB and CHOLMOD don't agree.  See if both are correct,
	% because differences in roundoff errors can make one go
	% a little farther than the other.

	% if p1 is zero, it means MATLAB was fully successful
	k1 = p1 ;
	if (k1 == 0)
	    k1 = n ;
	end

	% if p2 is zero, it means CHOLMOD was fully successful
	k2 = p2 ;
	if (k2 == 0)
	    k2 = n ;
	end

	if (k1 > k2)
	    % MATLAB went further than CHOLMOD. This is OK if MATLAB found
	    % a small entry where CHOLMOD stopped.
	    k = k2 ;
	    kerr = R1 (k,k) ;
	    % now reduce R1 in size, to compare with R2
	    R1 = R1 (1:k2-1,:) ;
	else
	    % CHOLMOD went further than MATLAB. This is OK if CHOLMOD found
	    % a small entry where MATLAB stopped.
	    k = k1 ;
	    kerr = R2 (k,k) ;
	    % now reduce R2 in size, to compare with R1
	    R2 = R2 (1:k1-1,:) ;
	end
    end

    err = norm (R1-R2,1) / max (anorm,1) ;
    fprintf ('p: %6d %6d MATLAB: %10.4f CHOLMOD: %10.4f speedup %6.2f err:', ...
	p1, p2, t1, t2, t1/t2) ;

    if (err == 0)
	fprintf ('      0') ;
    else
	fprintf (' %6.0e', err) ;
    end

    if (kerr == 0)
	fprintf ('      0\n') ;
    else
	fprintf (' %6.0e\n', kerr) ;
    end

%    if (err > 1e-6)
%	error ('!') ;
%    end

% pause
    clear R1 R2 p1 p2 p A

end

fprintf ('test22: all tests passed\n') ;


text:
% SCRIPT TEST FOR THE KINEMATIC PROBLEM FOR SERIAL ROBOTS

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
%
% qinv =
%
%          0         0         0         0    3.1416    3.1416    3.1416    3.1416
%    -0.0412   -0.0412    1.6597    1.6597   -1.5798   -1.5798   -0.3550   -0.3550
%     0.0720    0.0720    3.0696    3.0696   -0.3825   -0.3825   -2.7591   -2.7591
%     0.0000   -3.1416    3.1416   -0.0000    0.0000    3.1416    0.0000   -3.1416
%    -0.0308    0.0308   -1.5539    1.5539   -1.1793    1.1793   -0.0275    0.0275
%    -0.0000    3.1416   -3.1416         0    3.1416   -0.0000    3.1416   -0.0000

function test_kinematics_jaco

close all

fprintf('\nTHE DEMO PRESENTS THE DIRECT AND INVERSE KINEMATIC PROBLEM')

%load robot parameters. You can try different robots%
robot=load_robot('KINOVA', 'JACO');
%adjust 3D view as desired
adjust_view(robot)

q =[pi/2 -pi/2 0 0 0 0];
q0 = q + [0.1 0.1 0.1 0.1 0.1 0.1];

T = directkinematic(robot, q)

%Call the inversekinematic for this robot. All the possible solutions are
%stored at qinv. At least, one of the possible solutions should match q
qinv = inversekinematic(robot, T, q0)

T2 = directkinematic(robot, qinv)

qinv(:,1)
q

T-T2

function test_movement(robot, qinv, T)

q0 = [0, 0, 0, 0, 0, 0];
n = 5
for i=1:size(qinv,2)
   qq = [];
   qi = qinv(:,i);
   for j=1:6
        v = linspace(q0(j), qi(j), n);
        qq = [qq; v];
   end
   animate(robot, qq)
end

function test_solutions(robot, qinv, T)
n_solutions = 8;

fprintf('\nNOW WE CAN REPRESENT THE DIFFERENT SOLUTIONS TO ACHIEVE THE SAME POSITION AND ORIENTATION\n')
fprintf('\nNote that some solutions may not be feasible since some joints may be out of range.\n')
correct=zeros(1,n_solutions);
%check that all of them are possible solutions!
for i=1:size(qinv,2)

    Ti = directkinematic(robot, qinv(:,i)) %Ti is constant for the different solutions

    % Note that all the solutions may not be feasible. Some of the joints may
    % be out of range. You can test this situation with test_joints
    test_joints(robot, qinv(:,i));

    %now draw the robot to see the solution
    drawrobot3d(robot, qinv(:,i))

    pause(0.5);

    k=sum(sum((T-Ti).^2));
    if k < 0.01 % a simple threshold to find differences in the solution
        correct(1,i)= 1;
    else
        correct(1,i)= 0; %uncorrect solution
        fprintf('\nERROR: One of the solutions seems to be uncorrect. Sum of errors: %f', i, k);
    end
end

fprintf('\n************** RESULTS **************')

%Display a message if any of the solutions is not correct
if sum(correct)==n_solutions
    fprintf('\nTEST 1--> OK: Every solution in qinv yields the same position/orientation T');
else
    fprintf('\nTEST 1--> ERROR: One or more of the solutions seem to be uncorrect.');
end

fprintf('\n************** ****** **************\n')

repo: ARTE


text:
function prob_test077 ( )

%*****************************************************************************80
%
%% TEST077 tests GAMMA_CDF, GAMMA_PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  a = 1.0;
  b = 1.5;
  c = 3.0;

  check = gamma_check ( a, b, c );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST077 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST077\n' );
  fprintf ( 1, '  For the Gamma PDF:\n' );
  fprintf ( 1, '  GAMMA_CDF evaluates the CDF.\n' );
  fprintf ( 1, '  GAMMA_PDF evaluates the PDF.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =             %14f\n', a );
  fprintf ( 1, '  PDF parameter B =             %14f\n', b );
  fprintf ( 1, '  PDF parameter C =             %14f\n', c );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  X  PDF   CDF\n' );
  fprintf ( 1, '\n' );

  for i = 0 : 10

    x = i / 5.0;

    cdf = gamma_cdf ( x, a, b, c );

    pdf = gamma_pdf ( x, a, b, c );

    fprintf ( 1, '  %12f  %12f  %12f\n', x, pdf, cdf );

  end

  return
end


text:
function [ ap, det, inert ] = chpdi ( ap, n, ipvt, job )

%*****************************************************************************80
%
%% CHPDI: determinant, inertia and inverse of a complex hermitian matrix.
%
%  Discussion:
%
%    The routine uses the factors from CHPFA.
%
%    The matrix is stored in packed form.
%
%    A division by zero will occur if the inverse is requested and CHPCO has
%    set RCOND == 0.0 or CHPFA has set INFO ~= 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, complex AP(N*(N+1)/2); the factored matrix
%    from CHPFA.
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer IPVT(N), the pivot vector from CHPFA.
%
%    Input, integer JOB, has the decimal expansion ABC where:
%    if C ~= 0, the inverse is computed,
%    if B ~= 0, the determinant is computed,
%    if A ~= 0, the inertia is computed.
%    For example, JOB = 111 gives all three.
%
%    Output, complex AP(N*(N+1)/2); if the inverse was requested, then
%    the upper triangle of the inverse of the original matrix, stored in packed
%    form.  The columns of the upper triangle are stored sequentially in a
%    one-dimensional array.
%
%    Output, real DET(2), if requested, the determinant of the original matrix.
%    Determinant = DET(1) * 10.0**DET(2) with 1.0 <= abs ( DET(1) ) < 10.0
%    or DET(1) = 0.0.
%
%    Output, integer INERT(3), if requested, the inertia of the original matrix.
%    INERT(1) = number of positive eigenvalues.
%    INERT(2) = number of negative eigenvalues.
%    INERT(3) = number of zero eigenvalues.
%
  det = [];
  inert = [];

  noinv = mod ( job, 10 ) == 0;
  nodet = floor ( mod ( job, 100 ) / 10 ) == 0;
  noert = floor ( mod ( job, 1000 ) / 100 ) == 0;

  if ( ~nodet | ~noert )

    if ( ~noert )
      inert(1:3) = 0;
    end

    if ( ~nodet )
      det(1) = 1.0;
      det(2) = 0.0;
    end

    t = 0.0;
    ik = 0;

    for k = 1 : n

      kk = ik + k;
      d = real ( ap(kk) );
%
%  Check if 1 by 1
%
      if ( ipvt(k) <= 0 )
%
%  2 by 2 block
%  Use DET (D  S; S  C)  =  ( D / T * C - T ) * T, T = abs ( S )
%  to avoid underflow/overflow troubles.
%  Take two passes through scaling.  Use T for flag.
%
        if ( t == 0.0 )
          ikp1 = ik + k;
          kkp1 = ikp1 + k;
          t = abs ( ap(kkp1) );
          d = ( d / t ) * real ( ap(kkp1+1) ) - t;
        else
          d = t;
          t = 0.0;
        end

      end

      if ( ~noert )

        if ( 0.0 < d )
          inert(1) = inert(1) + 1;
        elseif ( d < 0.0 )
          inert(2) = inert(2) + 1;
        elseif ( d == 0.0 )
          inert(3) = inert(3) + 1;
        end

      end

      if ( ~nodet )

        det(1) = det(1) * d;

        if ( det(1) ~= 0.0 )

          while ( abs ( det(1) ) < 1.0 )
            det(1) = det(1) * 10.0;
            det(2) = det(2) - 1.0;
          end

          while ( 10.0 <= abs ( det(1) ) )
            det(1) = det(1) / 10.0;
            det(2) = det(2) + 1.0;
          end

        end

      end

      ik = ik + k;

    end

  end
%
%  Compute inverse(A).
%
  if ( ~noinv )

    k = 1;
    ik = 0;

    while ( k <= n )

      km1 = k - 1;
      kk = ik + k;
      ikp1 = ik + k;
      kkp1 = ikp1 + k;
%
%  1 by 1
%
      if ( 0 <= ipvt(k) )

        ap(kk) = 1.0 / real ( ap(kk) );

        if ( 1 <= km1 )

          work(1:km1) = ap(ik+1:ik+km1);

          ij = 0;
          for j = 1 : km1
            jk = ik + j;
            ap(jk) = conj ( ap(ij+1:ij+j) ) * tranpose ( work(1:j) );
            ap(ik+1:ik+j-1) = ap(ik+1:ik+j-1) + work(j) * ap(ij+1:ij+j-1);
            ij = ij + j;
          end

          ap(kk) = ap(kk) + real ( conj ( work(1:km1) ) * transpose ( ap(ik+1:ik+km1) ) );

        end

        kstep = 1;
%
%  2 by 2
%
      else

        t = abs ( ap(kkp1) );
        ak = real ( ap(kk) ) / t;
        akp1 = real ( ap(kkp1+1) ) / t;
        akkp1 = ap(kkp1) / t;
        d = t * ( ak * akp1 - 1.0 );
        ap(kk) = akp1 / d;
        ap(kkp1+1) = ak / d;
        ap(kkp1) = -akkp1 / d;

        if ( 1 <= km1 )

          work(1:km1) = ap(ikp1+1:ikp1+km1);

          ij = 0;
          for j = 1 : km1
            jkp1 = ikp1 + j;
            ap(jkp1) = conj ( ap(ij+1:ij+j) ) * transpose ( work(1:j) );
            ap(ikp1+1:ikp1+j-1) = ap(ikp1+1:ikp1+j-1) ...
              + work(j) * ap(ij+1:ij+j-1);
            ij = ij + j;
          end

          ap(kkp1+1) = ap(kkp1+1) ...
            + real ( conj ( work(1:km1) ) * transpose ( ap(ikp1+1) ) );

          ap(kkp1) = ap(kkp1) ...
            + conj ( ap(ik+1:ik+km1) ) * transpose ( ap(ikp1+1:ikp1+km1) );

          work(1:km1) = ap(ik+1:ik+km1);

          ij = 0;

          for j = 1 : km1
            jk = ik + j;
            ap(jk) = conj ( ap(ij+1:ij+j) ) * transpose ( work(1:j) );
            ap(ik+1:ik+j-1) = ap(ik+1:ik+j-1) + work(j) * ap(ij+1:ij+j-1);
            ij = ij + j;
          end

          ap(kk) = ap(kk) ...
            + real ( conj ( work(1:km1) ) * transpose ( ap(ik+1:ik+km1) ) );

        end

        kstep = 2;

      end
%
%  Swap
%
      ks = abs ( ipvt(k) );

      if ( ks ~= k )

        iks = ( ks * ( ks - 1 ) ) / 2;

        temp             = ap(iks+1:iks+ks);
        ap(iks+1:iks+ks) = ap(ik+1:ik+ks);
        ap(ik+1:ik+ks)   = temp;

        ksj = ik + ks;

        for jb = ks : k

          j = k + ks - jb;
          jk = ik + j;

          temp    = conj ( ap(jk) );
          ap(jk)  = conj ( ap(ksj) );
          ap(ksj) = temp;

          ksj = ksj - ( j - 1 );

        end

        if ( kstep ~= 1 )

          kskp1 = ikp1 + ks;

          temp      = ap(kskp1);
          ap(kskp1) = ap(kkp1);
          ap(kkp1)  = temp;

        end

      end

      ik = ik + k;

      if ( kstep == 2 )
        ik = ik + k + 1;
      end

      k = k + kstep;

    end

  end

  return
end


text:
function images = U_prepare_image(im,IMAGE_MEAN,cnn_name,opt)
% ------------------------------------------------------------------------
if ~exist('IMAGE_MEAN','var') || isempty(IMAGE_MEAN)
    switch cnn_name
    case {'alex','nin','gnet'} % alexnet, nin, inception
        d = load('data/ilsvrc_2012_mean.mat');
        IMAGE_MEAN = d.image_mean;
    case 'vgg16'
        IMAGE_MEAN = repmat(reshape([103.939, 116.779, 123.68],[1,1,3]),[256 256]);
    end
end
if ~exist('opt','var');opt=0;end

IMAGE_DIM = 256;
switch opt
    case -1
        % visual: denormalize
        switch cnn_name
            case 'alex';CROPPED_DIM = 227;
            case {'vgg16','nin','gnet'};CROPPED_DIM = 224;
        end
        indices = [0 IMAGE_DIM-CROPPED_DIM] + 1;
        center = floor(indices(2) / 2)+1;
        images = bsxfun(@plus,permute(im(:,:,[3 2 1],:),[2,1,3,4]) , IMAGE_MEAN(center:center+CROPPED_DIM-1,center:center+CROPPED_DIM-1,3:-1:1));
    case -2
        % visual: normalize
        switch cnn_name
            case 'alex';CROPPED_DIM = 227;
            case {'vgg16','nin','gnet'};CROPPED_DIM = 224;
        end
        indices = [0 IMAGE_DIM-CROPPED_DIM] + 1;
        center = floor(indices(2) / 2)+1;
        images = bsxfun(@minus,im(:,:,[3 2 1],:) , IMAGE_MEAN(center:center+CROPPED_DIM-1,center:center+CROPPED_DIM-1,:));
        images = permute(images,[2,1,3,4]);
    case {0,1,5,6}
        switch cnn_name
            case 'alex'
                % alexnet
                CROPPED_DIM = 227;
                % resize to fixed input size
                im = single(im);
                im = imresize(im, [IMAGE_DIM IMAGE_DIM], 'bilinear');
                % permute from RGB to BGR (IMAGE_MEAN is already BGR)
                im = im(:,:,[3 2 1]) - IMAGE_MEAN;

                % oversample (4 corners, center, and their x-axis flips)
                images = zeros(CROPPED_DIM, CROPPED_DIM, 3, 10, 'single');
                indices = [0 IMAGE_DIM-CROPPED_DIM] + 1;
                curr = 1;
                for i = indices
                    for j = indices
                        images(:, :, :, curr) = ...
                            permute(im(i:i+CROPPED_DIM-1, j:j+CROPPED_DIM-1, :), [2 1 3]);
                        images(:, :, :, curr+5) = images(end:-1:1, :, :, curr);
                        curr = curr + 1;
                    end
                end
                center = floor(indices(2) / 2)+1;
                images(:,:,:,5) = ...
                    permute(im(center:center+CROPPED_DIM-1,center:center+CROPPED_DIM-1,:), ...
                    [2 1 3]);
                images(:,:,:,10) = images(end:-1:1, :, :, curr);
            case 'vgg16'
                CROPPED_DIM = 224;

                % resize to fixed input size
                im = single(im);

                if size(im, 1) < size(im, 2)
                    im = imresize(im, [IMAGE_DIM NaN]);
                else
                    im = imresize(im, [NaN IMAGE_DIM]);
                end

                % RGB -> BGR
                im = im(:, :, [3 2 1]);

                % oversample (4 corners, center, and their x-axis flips)
                images = zeros(CROPPED_DIM, CROPPED_DIM, 3, 10, 'single');

                indices_y = [0 size(im,1)-CROPPED_DIM] + 1;
                indices_x = [0 size(im,2)-CROPPED_DIM] + 1;
                center_y = floor(indices_y(2) / 2)+1;
                center_x = floor(indices_x(2) / 2)+1;

                % hack the mean image
                IMAGE_MEAN =  repmat(reshape([103.939, 116.779, 123.68],[1,1,3]),[CROPPED_DIM CROPPED_DIM]);
                curr = 1;
                for i = indices_y
                    for j = indices_x
                        images(:, :, :, curr) = ...
                            permute(im(i:i+CROPPED_DIM-1, j:j+CROPPED_DIM-1, :)-IMAGE_MEAN, [2 1 3]);
                        images(:, :, :, curr+5) = images(end:-1:1, :, :, curr);
                        curr = curr + 1;
                    end
                end
                images(:,:,:,5) = ...
                    permute(im(center_y:center_y+CROPPED_DIM-1,center_x:center_x+CROPPED_DIM-1,:)-IMAGE_MEAN, ...
                    [2 1 3]);
                images(:,:,:,10) = images(end:-1:1, :, :, curr);
        end

        if opt==1
            images = images(:,:,:,5);
        end

end

author: donglaiw

repo: mNeuron

save_path: github-repos/MATLAB/donglaiw-mNeuron


text:
% model_setup_ui
% Run this to specify a design for optimization with GA2
%
% The script creates two struct variables, mspec and conditions,
% which are then used as inputs to ga2.m
%
% mspec and conditions contain fields, accessed with mspec.fieldname,
% that contain information about the design.
% Numbers in these variables describe timing info in TRs! (except the runtime and TR fields, which are in s)
%
% This function is good for:
%   1 - creating a design matrix with combinations of convolved, unconvolved, and shape-free elements
%   2 - creating designs AT THE RESOLUTION OF THE TR
%   3 - creating designs for optimization of temporal jitter and estimation of parts of a single trial
%   4 - creating designs with completely independent conditions (can overlap, etc.)
%
% This function is NOT as good for:
%   1 - optimizing the ordering of stimulus presentation for multiple trial types
%   2 - creating designs that have 1 or 2 events in some conditions.  confuses the program.
%   3 - creating designs in which trial types are mutually exclusive (present A or B at time t, but not both)

% scanlength, tr

% event, epoch, block
% convolve
% regular / random

% each condition is a trial type
% each trial type can be
%   event
%   epoch
%   multi-part trial

% each condition has these fields:
%   onsets      vector, for fixed onsets, times in s
%               [onset range; trial length range] for random onsets, in s
%               0 is first second of run
%
%   parts       epochs within a trial, for multipart trials
%               for each part, specify onsets (vector or range)
%               onsets can be an integer, for fixed onsets rel to trial or prev part.

% for multi-part
%   no of parts
%

clear mspec, clear conditions

conditions(1).onsets = 1;
i = 1;

disp(['Model setup user-interface: for specifying a design matrix.'])
disp(['by Tor Wager, version a1'])
disp(['___________________________________________________________'])

mspec.TR = input('Enter TR in s ');
mspec.runtime = input('Enter run length in s ');
mspec.hrf = spm_hrf(mspec.TR); mspec.hrf = mspec.hrf ./ max(mspec.hrf);
mspec.numframes = mspec.runtime ./ mspec.TR;

tmp = input('All trials at least N TRs apart? (Enter N or return to skip.) ');
if ~isempty(tmp), mspec.toverlap = [1 tmp];,end

disp(['Condition setup:'])
disp(['Press return when asked for onsets for a condition to exit.'])
disp(['___________________________________________________________'])

while ~isempty(conditions(i).onsets)

    conditions(i).onsets = input(['Onsets for condition ' num2str(i) ' in s ']);
    conditions(i).onsets = conditions(i).onsets ./ mspec.TR;

    if isempty(conditions(i).onsets), break, end

    conditions(i).parts = input(['Parts for condition ' num2str(i) ': enter 1 or num parts ']);
    % if this is a multi-part trial, parts > 1

        for j = 1:conditions(i).parts

            rel_to = input(['Onsets for condition ' num2str(i) ' part ' num2str(j) ' relative to: 1 for start of trial, 2 for prev. part ']);
            if rel_to == 1, str = ' start of trial ';, elseif rel_to == 2, str = ' previous part ';, else, error('Choose 1 or 2'), end

            if rel_to == 2 & j == 1, error('Cannot specify onset of 1st trial part relative to previous part.'),end
            conditions(i).subcond(j).rel_to = rel_to;
            conditions(i).subcond(j).onsets = input(['Onsets in s for condition ' num2str(i) ' part ' num2str(j) ' from' str]);
            if length(conditions(i).subcond(j).onsets) == 1,
                conditions(i).subcond(j).onsets = repmat(conditions(i).subcond(j).onsets,1,length(conditions(i).onsets));
            end
            conditions(i).subcond(j).onsets = conditions(i).subcond(j).onsets ./ mspec.TR;

            conditions(i).subcond(j).stimlength = input(['Stimulation length in TRs for condition ' num2str(i) ' part ' num2str(j) ' ']);

            conditions(i).subcond(j).convolve = input(['Convolve condition ' num2str(i) ' part ' num2str(j) ' with canonical HRF? 1/0 ']);
            if ~conditions(i).subcond(j).convolve
                conditions(i).subcond(j).hrfest = input(['Estimate shape-free HRF for condition ' num2str(i) ' part ' num2str(j) ' ? 1 or TP to estimate ']);
            else
                conditions(i).subcond(j).hrfest = [];
            end

        end

    if conditions(i).parts == 1, conditions(i).stimlength = conditions(i).subcond(1).stimlength;, end

    i = i + 1;
    conditions(i).onsets = 1;
end

conditions = conditions(1:end-1);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%UI Controller class to create an EventDetection algorithm from a UI
classdef EventDetectionConfigurator < handle
    properties (Access = public)
        uiElements EventDetectionConfiguratorUIElements;
        algorithmConfigurator;
        preprocessingConfigurator;
        delegate;
    end

    properties (Access = private)
        eventDetectionAlgorithms;
    end

    methods (Access = public)
        function obj = EventDetectionConfigurator(signals, preprocessingSignalAlgorithms,...
                eventDetectionAlgorithms, uiElements,delegate)

            obj.uiElements = uiElements;
            obj.eventDetectionAlgorithms = eventDetectionAlgorithms;

            if nargin > 4
                obj.delegate = delegate;
            end

            obj.preprocessingConfigurator = PreprocessingConfigurator(...
                signals,...
                preprocessingSignalAlgorithms,...
                uiElements.preprocessingSignalsList,...
                uiElements.preprocessingAlgorithmsList,...
                uiElements.preprocessingAlgorithmsVariablesTable,...
                obj);

            obj.algorithmConfigurator = AlgorithmConfigurator(...
                eventDetectionAlgorithms,...
                uiElements.eventDetectionList,...
                uiElements.eventDetectionVariablesTable,...
                obj);
        end

        function algorithm = createEventDetectorWithUIParameters(obj)
            preprocessingAlgorithm = obj.createPreprocessingAlgorithmWithUIParameters();
            if isempty(preprocessingAlgorithm)
                algorithm = [];
            else
                eventDetector = obj.algorithmConfigurator.createAlgorithmWithUIParameters();
                if isa(eventDetector,'NoOp')
                    algorithm = [];
                else
                    preprocessingAlgorithm.lastAlgorithms{1}.addNextAlgorithm(eventDetector);
                    algorithm = preprocessingAlgorithm;
                end
            end
        end

        function algorithm = createPreprocessingAlgorithmWithUIParameters(obj)
            algorithm = obj.preprocessingConfigurator.createPreprocessingAlgorithmWithUIParameters();
        end

        function setSignals(obj,signals)
            obj.preprocessingConfigurator.setSignals(signals);
        end

        %% handles
        function handlePreprocessingAlgorithmChanged(obj,preprocessingAlgorithm,~)
            outputType = preprocessingAlgorithm.outputPort;
            obj.algorithmConfigurator.algorithms = ...
                Palette.FilterAlgorithmsToInputType(obj.eventDetectionAlgorithms,outputType);

            obj.algorithmConfigurator.reloadUI();
        end

        function handleAlgorithmChanged(obj,algorithm,~)
            if ~isempty(obj.delegate)
                obj.delegate.handleAlgorithmChanged(algorithm,obj);
            end
        end

    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
function model = utl_complete_model(model,func,bundle)
% Internal. Complete the definition of a predicitive model with some meta-data.
% utl_complete_model(Model, PredictionFunction, Bundle)
%
% This function ensures that model structs created by BCI paradigm plugins have a uniform format
% and contain the set of fields needed by downstream processing functions (e.g., the offline and
% online prediction framework). This function relieves the BCI paradigms from specifying each and
% every of those fields manually, by automatically appending fields that can be deduced from the
% data or other parts of the model.
%
% This function also performs some early validation of the model structure to catch possible errors
% early on rather than deep into a subsequent processing run.
%
% In:
%   Model : a predictive model that will be finalized by this function; may have the following fields:
%           .tracking.filter_graph : cell array of symbolic filter expressions, one per inlet of the
%                                    prediction function; can be used to filter streams or bundles
%
%                                    Can also be passed either as an EEGLAB data set (with .tracking
%                                    field) or as a stream bundle; in this case, the processing that
%                                    has been applied to those data sets will be used as a template
%                                    to lay out the filter graph.
%
%                                    If not passed, will be filled in by the Bundle argument.
%
%           .tracking.prediction_function : function that takes a window of the output of filter graph's
%                                           output and computes a prediction from it; If not passed,
%                                           will be filled in based on the PredictionFunction argument
%
%           .tracking.prediction_channels : cell array of chanlocs structures expected by the prediction
%                                           function; has one element for each inlet of the prediction
%                                           function.
%
%                                           If not passed, will be filled in based on the Bundle argument.
%
%
%   PredictionFunction : handle of the prediction function (supplied by the framework code that calls
%                        utl_complete_model)
%
%   Bundle : a stream bundle fully preprocessed by the filter graph, used to extract some
%            meta-data (supplied by framework code)
%
% Out:
%   Model : a completed model struct; the following fields are ensured to exist:
%           .tracking                --> tracking information about the model for offline/online processing
%            .filter_graph           --> a cell array of (symbolic) filter expressions, denoting the
%                                        filter operations that yield each stream inlet that is
%                                        accepted by the model's prediction function. Each expression
%                                        may contain leaf nodes of the type @rawdata, which are free
%                                        inlets of the filter graph (tagged with the channels expected
%                                        at that inlet)
%            .prediction_function    --> a function that takes a stream bundle with one stream for
%                                        each entry of the filter graph; produces a prediction for each
%                                        epoch/segment in the data (therefore works both online and offline)
%            .prediction_channels    --> a cell array of chanlocs structs that defines that channels
%                                        expected by the prediction function at each of its inlets
%                                        (this is redundant additional information)
%
% See also:
%   bci_train
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2011-08-28
dp;

% input validation
if ~isstruct(model) || ~isscalar(model)
    error('The given model must be a 1x1 struct, but was: %s',hlp_tostring(model,10000)); end
if exist('bundle','var')
    if ~isstruct(bundle) || ~isscalar(bundle)
        error('The given stream bundle must be a 1x1 struct, but was: %s',hlp_tostring(bundle,10000)); end
    if ~isfield(bundle,'streams') && isfield(bundle,'tracking') && isfield(bundle.tracking,'online_expression')
        bundle.streams = {bundle}; end
    for s=1:length(bundle.streams)
        utl_check_fields(bundle.streams{s},{'tracking','chanlocs'},'signal','signal'); end
end
if exist('func','var') && ~(isa(func,'function_handle') || isa(func,'char'))
    error('The given PredictionFunction argument must either be a function handle or a string.'); end

% the .tracking field holds the data that is of interest to the offline/online evaluation system
% it is allowed for the paradigm to fill in the following field manually, in order to yield custom behavior
if ~isfield(model,'tracking')
    model.tracking = struct(); end

% add the filter_graph, if necessary
if ~isfield(model.tracking,'filter_graph')
    if ~exist('bundle','var')
        error('Cannot automatically deduce the .tracking.filter_graph field for this model, because the data on which its state shall be initialized is not known or available;\nthis most likely the case if your model was calibrated on a dataset collection and the framework could not decide which dataset in the collection to initialize the filters on.\nYour can set the field directly in your paradigm''s calibrate() function.'); end
    for s=1:length(bundle.streams)
        model.tracking.filter_graph{s} = bundle.streams{s}.tracking.online_expression; end
else
    % a filter graph was passed in; check format
    if all(isfield(model.tracking.filter_graph,{'tracking','chanlocs'}))
        % ... as a data set; override the bundle with it
        bundle = struct('streams',{{model.tracking.filter_graph}});
        % set the filter graph accordingly
        model.tracking.filter_graph = {bundle.streams{1}.tracking.online_expression};
    elseif isfield(model.tracking.filter_graph,'streams')
        % ... as a stream bundle: override the bundle with it
        bundle = model.tracking.filter_graph;
        % and initialize the filter graph accordingly
        model.tracking.filter_graph = {};
        for s=1:length(bundle.streams)
            model.tracking.filter_graph{s} = bundle.streams{s}.tracking.online_expression; end
    elseif all(isfield(model.tracking.filter_graph,{'head','parts'}))
        % ... as an abstract expression
        model.tracking.filter_graph = {model.tracking.filter_graph};
    elseif iscell(model.tracking.filter_graph)
        % ... in the proper format (cell array)
        if isempty(model.tracking.filter_graph)
            error('The given model''s .tracking.filter_graph is empty.'); end
        for c=1:length(model.tracking.filter_graph)
            if ~all(isfield(model.tracking.filter_graph{c},{'head','parts'}))
                error('The given model''s .tracking.filter_graph is not a valid filter graph (at least one of its cells is lacking the required fields .head and/or .parts): %s',hlp_tostring(model.tracking.filter_graph)); end
        end
    else
        error('The given model''s .tracking.filter_graph is not a valid filter graph: %s',hlp_tostring(model.tracking.filter_graph));
    end
end

% add the prediction function, if necessary
if ~isfield(model.tracking,'prediction_function')
    if ~exist('func','var')
        error('If the given model is lacking a .tracking.prediction_function field, the prediction function should be supplied as second input to utl_complete_model.'); end
    model.tracking.prediction_function = func;
end

% validate the model's prediction function
if ischar(model.tracking.prediction_function)
    if strncmp(model.tracking.prediction_function,'Paradigm',8)
        try
            % Paradigm class reference: try to instantiate
            instance = eval(model.tracking.prediction_function); %#ok<NASGU>
        catch e
            error('Failed to instantiate the paradigm referred to by the model''s .tracking.prediction_function with error: %s',hlp_handleerror(e));
        end
    else
        if ~exist(model.tracking.prediction_function,'file')
            error('The function referred to by the model''s .tracking.prediction_function field does not exist: %s',model.tracking.prediction_function); end
    end
else
    tmp = char(model.tracking.prediction_function);
    if tmp(1) ~= '@' && ~exist(tmp,'file')
        error('The function referred to by the model''s .tracking.prediction_function field does not exist: %s',char(model.tracking.prediction_function)); end
end

% add and validate the channel locations (these are technically redundant, as they can be
% derived from the preprocessing chain) (deduced from the bundle)
try
    if ~isfield(model.tracking,'prediction_channels')
        if ~exist('bundle','var')
            error('Cannot automatically assign the .tracking.prediction_channels field because the data set from which these channels shall be determined could not be deduced by the framework; consider setting this field manually in the calibrate() function of your paradigm.'); end
        for s=1:length(bundle.streams)
            if ~((isstruct(bundle.streams{s}.chanlocs) && isfield(bundle.streams{s}.chanlocs,'labels')) || iscellstr(bundle.streams{s}.chanlocs))
                error('The .chanlocs fields in one of the streams of the stream bundle is malformed; needs to be a chanlocs struct (or cell array of channel labels), but was: %s',s,hlp_tostring(bundle.streams{s}.chanlocs)); end
            if iscellstr(bundle.streams{s}.chanlocs)
                bundle.streams{s}.chanlocs = struct('labels',bundle.streams{s}.chanlocs); end
            model.tracking.prediction_channels{s} = bundle.streams{s}.chanlocs;
        end
    else
        for s=1:length(model.tracking.prediction_channels)
            if ~((isstruct(model.tracking.prediction_channels{s}) && isfield(model.tracking.prediction_channels{s},'labels')) || iscellstr(model.tracking.prediction_channels{s}))
                error('The given model''s .tracking.prediction_channels{k} field is malformed for k=%i; needs to be a chanlocs struct or cell array of channel labels, but was: %s',s,hlp_tostring(model.tracking.prediction_channels{s})); end
        end
    end
catch e %#ok<CTCH>
    fprintf('Could not initialize prediction channels: %s\n',hlp_handleerror(e));
end

if length(model.tracking.filter_graph) ~= length(model.tracking.prediction_channels)
    error('The given filter graph has a different number of elements than the given prediction channels; they should have the same length.'); end

% finally add also a timestamp (so that we can sort them by creation time)
model.timestamp = now;

repo: REST


text:
function boolVal=linTargetIsUntrackable(targetParams,PD,lambda,AbsTol,RelTol,maxIter)
%%LINTARGETISUNTRACKABLE Assuming that an asymptotic inoovation covariance
%           is available or that a target follows a particular specified
%           linear dynamic model (e.g. position, velocity, and
%           acceleration, and it could have correlation terms like a Gauss-
%           Markov model) with known process noise covariance and Cartesian
%           measurement covariance, and assuming that a probability of
%           detection is known, determine whether a target is untrackable.
%           A target is deemed untrackable if the cost function used in a
%           multiple hypothesis tracker (MHT)  is always greater than or
%           equal to  the null hypothesis than for any track hypothesis
%           when given asymptotic track accuracy parameters. If a track is
%           untrackable, even an MHT of infinite length would not be able
%           to track the target. This assumes Cartesian measurements. The
%           use of range rate, complex target amplitude, micro-Doppler and
%           other information could potentially render untrackable targets
%           trackable.
%
%INPUTS: targetParams This can be one of two types of inputs. If an
%          asymptotic innovation covariance matrix is known, then this is
%          the innovation covariance matrix. Otherwise, this is a structure
%          holding parameters for a linear system so that an asymptotic
%          innovation covariance matrix can be computed. The members of the
%          structure are:
%           H The zDim X xDim measurement matrix such that H*x+w is the
%             measurement, where x is the state and w is zero-mean Gaussian
%             noise with covariance matrix R.
%           F The xDim X xDim state transition matrix The state at
%             discrete-time k+1 is modeled as F times the state at time k
%             plus zero- mean Gaussian process noise with covariance matrix
%             Q.
%           R The zDim X zDim measurement covariance matrix.
%           Q The xDim X xDim process noise covariance matrix. This can
%             not be a singular matrix.
%       PD The optional detection probability of the target at each scan.
%          PD should not be 1.
%   lambda The false alarm density in Cartesian coordinates. For example,
%          this might be 1e-8 false alarms per unit volume in cubic meters.
%   RelTol The maximum relative error tolerance allowed when computing the
%          asymptotic predictive covariance matrix. This is a positive
%          scalar. If omitted or an empty matrix is passed, the default
%          value of 1e-13 is used. This value is only used if targetParams
%          is a structure.
%   AbsTol The absolute error tolerance allowed, a positive scalar. When
%          computing the asymptotic predictive covariance matrix. This is a
%          positive scalar. If omitted or an empty matrix is passed, the
%          default value of 1e-10 is used. This value is only used if
%          targetParams is a structure.
%  maxIter An optional integer specifying the maximum number of iterations
%          to use when cmputing the asymptotic predictive covariance
%          matrix. By default, if omitted, this is 5000. This value is only
%          used if targetParams is a structure.
%
%As discussed in Chapter 7.5 of [1], track oriented MHTs declare target
%existence based on finite sums of log-likelihood ratios. These ratios are
%typically made into dimensionless score functions as in Chapter 7.5 of [1]
%and in [2]. For the probability that a measurement z is from a particular
%target, the contribution to the score function is f(z)*PD/lambda, where
%f(z) is the likelihood of the measurement conditioned on the predicted
%state estimate. The null hypothesis receives a score of 1-PD. If the null
%hypothesis is always greater than or equal to the track association
%hypothesis then a target is untrackable.
%
%Given a Cartesian measurement model (z=H*x+w where x is the target state,
%w is measurement noise with covariance matrix R and H is the measurement
%matrix) and assuming a linear dynamic model (xpredicted=F*xPrior+v where v
%is noise with covariance matrix Q) where the state transition matrix F and
%the process noise covariance matrix Q are constant (implying a uniform
%revisit rate) as well as having a known detection probability of PD,  one
%determine an asymptotic lower bound on the estimation accuracy of a
%predicted state measurement using the RiccatiPredNoClutter function. Given
%PPred, a lower bound on the state precision covariance, the score function
%contribution for a target-originated measurement is maximum when the
%predicted measurement equals the measurement value. From Chapter 7.5.3 of
%[1], for a linear measurement model, this is
%PD/(lambda*sqrt(det(2*pi*SPred))) where SPred is the innovation covariance
%matrix, which as given in the derivation of the Kalman filter in Chapter 5
%of [3] is H*PPred*H'+R. Thus, this function determines whether this
%maximum score for a track-originated measurement is greater than the
%alternative (1-PD). If not, then the target is untrackable.
%
%The trackability bound here can be considered something of a "stealth"
%bound, because for a fixed false alarm level in a system, a target can
%become untrackable either by lowing its detection probability sufficiently
%(becoming "stealth") and/or by increasing its maneuverability, which,
%here, is reflected in increasign the process noise covariance.
%
%Note that if the innovation covariance matrix is passed for targetParams,
%this function just uses its size and determiniant.
%
%EXAMPLE:
% T=1;
% targetParams=[];
% targetParams.F=FPolyKal(T,6,1);
% q0=1;
% targetParams.Q=QPolyKal(T,6,1,q0);
% targetParams.R=diag([10;10;10]);
% targetParams.H=[1,0,0,0,0,0;
%    0,1,0,0,0,0;
%    0,0,1,0,0,0];
% PD=0.5;
% lambda=1e-8;
% boolVal=linTargetIsUntrackable(targetParams,PD,lambda)
% %The above example is trackable (boolVal=false). However, if the
% %detection probability decreases/ the false alarm rate increases, then
% %one gets an untrackable target. Consider:
% PD=0.1;
% lambda=1e-6;
% boolVal1=linTargetIsUntrackable(targetParams,PD,lambda)
%
%REFERENCES:
%[1] Y. Bar-Shalom, P. K. Willett, and X. Tian, Tracking and Data Fusion.
%    Storrs, CT: YBS Publishing, 2011.
%[2] Y. Bar-Shalom, S. S. Blackman, and R. J. Fitzgerald, "Dimensionless
%    score function for multiple hypothesis tracking," IEEE Transactions on
%    Aerospace and Electronic Systems, vol. 43, no. 1, pp. 392-400, Jan.
%    2007.
%[3] Y. Bar-Shalom, X. R. Li, and T. Kirubarajan, Estimation with
%    Applications to Tracking and Navigation. New York: John Wiley and
%    Sons, Inc, 2001.
%
%January 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(AbsTol))
    AbsTol=1e-13;
end

if(nargin<5||isempty(RelTol))
    RelTol=1e-10;
end

if(nargin<6||isempty(maxIter))
   maxIter=5000;
end

if(isstruct(targetParams))
    H=targetParams.H;
    F=targetParams.F;
    R=targetParams.R;
    Q=targetParams.Q;
    m=size(H,1);

    %Asymptotic predicted state covariance matrix with given PD assuming
    %completely correct association hypotheses.
    PPred=RiccatiPredNoClutter(H,F,R,Q,PD,AbsTol,RelTol,maxIter);

    %Asymptotic innovation covariance matrix.
    detSPred=det(H*PPred*H'+R);
else
    S=targetParams;
    m=size(S,1);
    detSPred=det(S);
end

maxIsTargetScore=PD/(lambda*sqrt((2*pi)^m*detSPred));
missedDetectScore=1-PD;

boolVal=maxIsTargetScore<=missedDetectScore;
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function p = predictOneVsAll(all_theta, X)
%PREDICT Predict the label for a trained one-vs-all classifier. The labels
%are in the range 1..K, where K = size(all_theta, 1).
%  p = PREDICTONEVSALL(all_theta, X) will return a vector of predictions
%  for each example in the matrix X. Note that X contains the examples in
%  rows. all_theta is a matrix where the i-th row is a trained logistic
%  regression theta vector for the i-th class. You should set p to a vector
%  of values from 1..K (e.g., p = [1; 3; 1; 2] predicts classes 1, 3, 1, 2
%  for 4 examples)

m = size(X, 1);
num_labels = size(all_theta, 1);

% You need to return the following variables correctly
p = zeros(size(X, 1), 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
%               your learned logistic regression parameters (one-vs-all).
%               You should set p to a vector of predictions (from 1 to
%               num_labels).
%
% Hint: This code can be done all vectorized using the max function.
%       In particular, the max function can also return the index of the
%       max element, for more information see 'help max'. If your examples
%       are in rows, then, you can use max(A, [], 2) to obtain the max
%       for each row.
%

all_predictions = all_theta * X';
[M, I] = max(all_predictions);
p = I(:);

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

end

author: zlotus


text:
function test_issue1985

% WALLTIME 00:10:00
% MEM 2gb

% DEPENDENCY filter_with_correction fir_filterdcpadded

nchan = 1;
fsample = 1000;
nsample = 2*fsample;

dat = zeros(nchan, nsample);
dat(1,fsample) = 1;

%%

fhp = 0.3;
flp = 30;
fbp = [fhp flp];

order = 30;
type = 'firws';

% use defaulls where possible
filt1 = ft_preproc_bandpassfilter(dat, fsample, fbp, order, type, 'onepass-zerophase');
filt2 = ft_preproc_bandpassfilter(dat, fsample, fbp, order, type, 'onepass-reverse-zerophase');
filt3 = ft_preproc_bandpassfilter(dat, fsample, fbp, order, type, 'onepass-minphase');

%%

figure; hold on
plot(dat * 0.1)
plot(filt1+0.01)
plot(filt2+0.02)
plot(filt3+0.03)
legend


text:
% Fig. 6.13   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%

clear all;
%close all;
clf

num=10;
den=[1 1 0];
w=logspace(-2,1,100);
[m,p]=bode(num,den,w);
wa=[.01 1.2];
ma=[1000 8.3333];
loglog(w,m,w,10*ones(1,100),'--',wa,ma);
xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Fig. 6.13 Determination of Kv from the Bode plot');
bodegrid;

author: Sable

repo: mcbench-benchmarks


text:
function y = cvx_s_sparse( m, n, i, j )
%CVX_S_SPARSE Matrices with a fixed sparsity pattern.

if nargin < 3,
    error( 'Sparsity structure missing.' );
elseif ~isnumeric( i ) || ~isnumeric( j ),
    error( 'Sparsity arguments must be vectors of nonnegative integers.' );
end
i = i(:);
j = j(:);
if any( i <= 0 ) || any( j <= 0 ) || any( i ~= floor( i ) ) || any( j ~= floor( j ) ),
    error( 'Sparsity arguments must be vectors nonnegative integers.' );
elseif length( i ) ~= 1 && length( j ) ~= 1 && length( i ) ~= length( j ),
    error( 'Sparsity arguments have incompatible size.' );
elseif any( i > m ) || any( j > n ),
    error( 'One or more indices are out of range.' );
end
nz = max( length( i ), length( j ) );
y = min( sparse( 1 : nz, i + ( j - 1 ) * m, 1, nz, m * n ), 1 );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.


text:
function res = set_partition(varargin)
% Obtain a subset of a data set, using an index set.
% Signal = set_partition(Signal, IndexSet, EpochBounds)
%
% A partitioning function (as required by the generic functions utl_crossval, utl_searchmodel,
% utl_nested_crossval) that takes a data set (either epoched or continuous), and index set (either
% of the raw samples, the target events, or of the epochs), and returns a reduced (subset) data set.
%
% In:
%   Signal   : Continuous or epoched input data set.
%
%   IndexSet : Either [] or a vector of indices;
%               * if [], the function returns the cardinality (highest allowed index) of the index set,
%               * if a vector of indices, the function returns the data set reduced to those indices;
%                 if the indices refer to target events, but the data is continuous, then the
%                 partition function tries to respect the implied epoch bounds when partitioning
%
%   EpochBounds : Only required when partitioning a continuous data set based on target events;
%                 this is done with the help of implied epoch bounds (as in set_makepos)
%
% Out:
%   Result : Reduced data set or index set cardinality.
%
% Examples:
%   % for an epoched data set or a continuous set with target markers, get the number of trials
%   numtrials = set_partition(eeg,[])
%
%   % for an epoched data set or a continuous set with target markers, retain only trials 20:40
%   eeg = set_partition(eeg,20:40)
%
%   % for a continuous data set without target markers, get the number of samples
%   numsamples = set_partition(eeg,[])
%
%   % for a continuous data set without target markers, retain only samples 1:10000 and 50000:100000
%   eeg = set_partition(eeg,[1:10000 50000:100000])
%
% See also:
%   set_makepos, set_selepos, set_selinterval, utl_crossval
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-01
dp;

% set_partition_version<1.01> -- for the cache

if ~exp_beginfun('editing', 'memoize',0, 'fingerprint_create','passthrough') return; end

declare_properties('independent_channels',true,'independent_trials',false);

arg_define(varargin, ...
    arg_norep({'signal','Signal'}), ...
    arg({'idxset','IndexSet'},[],uint32([1 1000000]), 'Index set. This is either [] (indicating that the function should return the index set cardinality), or a vector of indices (indicating that the function shall partition the data set accordingly).','shape','row'), ...
    arg({'epoch_bounds','EpochBounds'},[],[], 'Epoch bounds. Only required when partitioning on a continuous data set based on target events; this is performed using the epoch bounds (as in set_makepos).','shape','row'));

% input validation
utl_check_fields(signal,{'data','srate','event'},'signal','signal'); %#ok<NODEF>
is_continuous = (~isfield(signal,'epoch') || isempty(signal.epoch)) && (isempty(signal.data) || size(signal.data,3) == 1);

if isempty(idxset) %#ok<NODEF>
    % --- calc index set size ---

    if ~is_continuous
        % epoched data set: return the number of epochs
        res = size(signal.data,3);
    else
        % continuous data set
        if ~isempty(signal.event) && isfield(signal.event,'target')
            targetmask = ~cellfun('isempty',{signal.event.target});
            if any(targetmask)
                % return the number of target events
                res = nnz(targetmask);
                exp_endfun; return;
            end
        end
        % otherwise return the number of samples
        res = size(signal.data,2);
    end

else
    % --- partition the data ---

    if ~is_continuous
        % set is epoched: strictly partition in terms of epochs
        res = exp_eval(set_selepos(signal,idxset));
    else
        % set is continuous
        if ~isempty(signal.event) && isfield(signal.event,'target')
            targetmask = ~cellfun('isempty',{signal.event.target});
            if any(targetmask)
                % ... and it has target events: strictly partition in terms of target events
                % this partitioning is perfectly consistent with what would happen if the data were epoched
                % (assuming that target markers with boundary-crossing epochs were already removed in
                %  set_targetmarkers, and assuming that the given epoch bounds are an upper bound of those
                %  used in set_makepos)
                if isempty(epoch_bounds) %#ok<NODEF>
                    disp_once('WARNING: Attempting to partition a continuous data set based on target events, but no EpochBounds argument is present. Assuming 1-sample epochs.');
                    ival = 0;
                else
                    % make sure that the epoch bounds include the event itself (and add a bit more slack
                    % to not lose it later even under resampling and rounding errors...)
                    epoch_bounds(1) = min(epoch_bounds(1),-0.25);
                    epoch_bounds(2) = max(epoch_bounds(2),+0.25);
                    ival = round(epoch_bounds(1)*signal.srate) : round(epoch_bounds(2)*signal.srate);
                end

                % check for sorting issues
                if ~issorted(idxset)
                    warn_once('Continuous-time data cannot reasonably be used with randomized partitioning; consider using chronological/blockwise partitioning! The index set will be sorted in the following.');
                    idxset = sort(idxset);
                end

                % build a logical mask that marks all retained intervals, starting with all false
                select_mask = false(1,size(signal.data,2));
                deselect_mask = false(1,size(signal.data,2));

                % find all target event latencies and the selected latencies
                if ~isfield(signal.event,'latency')
                    error('The given signal is lacking a .event.latency field.'); end
                latency_numels = cellfun('prodofsize',{signal.event.latency});
                if any(latency_numels == 0)
                    error('The signal contains events with an empty .latency field; this is not permitted.'); end
                if any(latency_numels ~= 1)
                    error('The signal contains events with a .latency value that is not a scalar; this is not permitted.'); end

                targlats = round([signal.event(targetmask).latency]);
                selected = targlats(idxset);

                % find the deselected latencies
                deselmask = true(size(targlats));
                deselmask(idxset) = false;
                deselected = targlats(deselmask);

                % add to the selected & deselected masks
                select_mask(min(length(select_mask),max(1,vec(bsxfun(@plus,ival,selected'))))) = true;
                deselect_mask(min(length(select_mask),max(1,vec(bsxfun(@plus,ival,deselected'))))) = true;

                % warn if the ranges overlap
                if any(select_mask & deselect_mask)
                    disp_once('Note: the epochs for some selected markers are overlapping with epochs of deselected markers. Make sure that you are using long enough safety margins if perfoming cross-validation.'); end

                % identify gaps in the select_mask
                gaps = find(diff([false select_mask false]));
                gaps = reshape(gaps(2:end-1),2,[])';
                gaps(:,2) = gaps(:,2)-1;

                % find all gaps that do not contain a deselected sample and fuse those...
                for g=1:size(gaps,1)
                    gap = gaps(g,1):gaps(g,2);
                    if ~any(deselect_mask(gap))
                        select_mask(gap) = true; end
                end

                % disable the now de-selected target events (settings their target flag back to empty)
                % (note that partitioning is strictly in terms of target markers)
                targetinds = find(targetmask);
                [signal.event(targetinds(deselmask)).target] = deal([]);

                % finally, re-extract the intervals and perform the actual partitioning
                retain_intervals = reshape(find(diff([false select_mask false])),2,[])';
                retain_intervals(:,2) = retain_intervals(:,2)-1;
                res = exp_eval(set_selinterval(signal,retain_intervals,'samples'));

                exp_endfun; return;
            end
        end

        % otherwise partition the raw data based on sample intervals

        % create a logical mask from the index set
        mask = false(1,size(signal.data,2)); mask(idxset) = true;
        % find the inclusive, 1-based intervals of non-zeros
        retain_intervals = reshape(find(diff([false mask false])),2,[])';
        retain_intervals(:,2) = retain_intervals(:,2)-1;
        res = exp_eval(set_selinterval(signal,retain_intervals,'samples'));
    end
end

exp_endfun;


text:
classdef S2FunHandle < S2Fun
% a class represeneting a function on the sphere

properties
  fun
  antipodal = false
end

methods
  function S2F = S2FunHandle(fun)
    S2F.fun = fun;
  end

  function f = eval(S2F,v)
    f = S2F.fun(v);
  end

end

methods (Static = true)

  function S2F = Kachanov(lambda)

    S2F = S2FunHandle(@(v) fun(v,lambda));

    function values = fun(v,lambda)

      phi = v.theta;
      values =  ((lambda.^2 + 1) * exp(-lambda * phi) + ...
        lambda*exp((-lambda*pi)/2))./(2*pi);

      values = values(:);

    end

  end

end

end

author: mtex-toolbox

repo: mtex

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


text:
 function ir_display_struct(ob, varargin)
%function ir_display_struct(ob, varargin)
%|
%| display a structure concisely
%|
%| option
%|	'pre'		prefix
%|
%| Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if streq(ob, 'test'), ir_display_struct_test, return, end

arg.pre = ' ';
arg = vararg_pair(arg, varargin);

%pn = jf_protected_names;
% display(struct(pn))
% tmp = pn.struct_recurse(ob);

if ~isstruct(ob)
	fail('%s not a struct', class(ob))
end

printf([inputname(1) ':'])
%printf('pre = "%s"', arg.pre)
ir_display_struct_do(ob, arg.pre)

% ir_display_struct_do()
function ir_display_struct_do(ob, pre)

fnames = fieldnames(ob);

nmax = 1; % find max length field name
for ii=1:length(fnames)
	nmax = max(nmax, length(fnames{ii}));
end
%nmax = min(nmax, 40);

format = sprintf('%d', nmax);
format = [pre '%' format 's: %s'];

for ii=1:length(fnames)
	fname = fnames{ii};
	tmp = ob.(fname);

	if isempty(tmp)
		str = sprintf('[] (%s)', class(tmp));
	else
		switch class(tmp)
		case 'function_handle'
			str = func2str(tmp);
			if (str(1) ~= '@')
				str = ['@' str];
			end
		case 'char'
			str = ['"' tmp '"'];
		case 'struct'
			str = '(struct)';
		case 'cell'
			str = ['{' num2str(size(tmp)) '}'];
		case {'double', 'single', 'logical'}
			if numel(tmp) == 1
				str = num2str(tmp);
			else
				str = ['[' num2str(size(tmp)) ']'];
				if numel(tmp) < 5
					str = [str ' : ' num2str(tmp(:)')];
				end
			end
			str = [str sprintf(' (%s)', class(tmp))];
		otherwise
			str = sprintf('? (%s)', class(tmp));
		end
	end

	printf(format, fname, str)
	if isstruct(tmp)
		white = sprintf('%d', nmax+2); % 1 more because of ':'
		white = [pre '%' white 's'];
		white = sprintf(white, ' ');
		ir_display_struct_do(tmp, white);
	end
end

function ir_display_struct_test
s.a = 'a';
s.b = 2;
s.c = {1,2,3};
s.d = ones(3);
s.e = 1+2i;
s.f1 = @(x) x+1;
s.f2 = @ir_display_struct;
s.g = complex(1);
t.s = 'sub';
u.v = 'subsub';
t.u = u;
s.t = t;
ir_display_struct(s)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
clear all
close all

wing_tp

wing_lmi

wing_sim


text:
function varargout = main(varargin)
% MAIN M-file for main.fig
%      MAIN, by itself, creates a new MAIN or raises the existing
%      singleton*.
%
%      H = MAIN returns the handle to a new MAIN or the handle to
%      the existing singleton*.
%
%      MAIN('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in MAIN.M with the given input arguments.
%
%      MAIN('Property','Value',...) creates a new MAIN or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before main_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to main_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help main

%   Copyright 1984-2010 The MathWorks, Inc.
% Last Modified by GUIDE v2.5 30-Sep-2005 11:05:33

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @main_OpeningFcn, ...
                   'gui_OutputFcn',  @main_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before main is made visible.
function main_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to main (see VARARGIN)

% Choose default command line output for main
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes main wait for user response (see UIRESUME)
% uiwait(handles.figure1);

setappdata(0  , 'hMainGui'    , gcf);
setappdata(gcf,   'thresh'    , 121);
setappdata(gcf, 'fhUpdateAxes', @updateAxes);

% --- Outputs from this function are returned to the command line.
function varargout = main_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

hMainGui = getappdata(0, 'hMainGui');

fileName = uigetfile('*.jpg');

setappdata(hMainGui, 'fileName', fileName);

updateAxes

function updateAxes

hMainGui = getappdata(0, 'hMainGui');
fileName = getappdata(hMainGui, 'fileName');
thresh   = getappdata(hMainGui, 'thresh');

I = imread(fileName);
I = (I > thresh);

hAxes = findobj(hMainGui, 'type', 'axes');
axes(hAxes)
imshow(I)

% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

threshGui;


text:
function [X,Y,U] = CollectData(SimPar,IC1,IC2,f1,f2,umin,umax,Ntraj,SimLength)
% Collect data for Burgers equation
disp('Starting data collection ...')

% random input
Ubig= rand(2,SimLength,Ntraj)* (umax-umin)+umin;

% Transition mapping of the controleld dynamical system
% control is a linear combination of u1 and u2
f = @(x,u)(BurgerSolver(x,u(1)*f1+u(2)*f2,SimPar));

% initialize
X = []; Y = []; U=[];

% loop pver trajectories
for i = 1:Ntraj
    xx = [];
    % Intial state is a random convex combination of IniitialConditions
    b = rand;
    a = [b,1-b];
    xx =b*IC1 + (1-b)*IC2;
    tic
    fprintf('Trajectory %d out of %d \n',i,Ntraj)
    % loop over each time step
    for j = 1:SimLength
        xx = [xx f(xx(:,end),Ubig(:,j,i))];
        U  = [U,Ubig(:,j,i)];
        % if the solution diverges, go to the next trajectory
        if ~isempty(find(isnan(xx(:,end)),1))
            break
        end

    end
    toc
    % Store
    X = [X xx(:,1:end-1)];
    Y = [Y xx(:,2:end)];

end

save('BurgersTrajectoryData','X','Y','U','SimLength','Ntraj')
end

author: arbabiha


text:
% EX_LAPLACE_ISO_PLATE: solve the Poisson problem in one quarter of a plate with a hole, discretized with NURBS (isoparametric approach).

% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = 'geo_plate_with_hole.txt';

% Type of boundary conditions for each side of the domain
problem_data.nmnn_sides   = [3 4];
problem_data.drchlt_sides = [1 2];

% Physical parameters
problem_data.c_diff  = @(x, y) ones(size(x));

% Source and boundary terms
problem_data.f = @(x, y) zeros(size(x));
problem_data.g = @test_plate_mixed_bc_g_nmnn;
problem_data.h = @(x, y, ind)  exp(x).*sin(y);

% Exact solution (optional)
problem_data.uex     = @(x, y) exp(x).*sin (y);
problem_data.graduex = @(x, y) cat (1, ...
                       reshape (exp(x).*sin(y), [1, size(x)]), ...
                       reshape (exp(x).*cos(y), [1, size(x)]));

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.degree     = [3 3];       % Degree of the splines
method_data.regularity = [2 2];       % Regularity of the splines
method_data.nsub       = [8 8];       % Number of subdivisions
method_data.nquad      = [4 4];       % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
[geometry, msh, space, u] = solve_laplace_iso (problem_data, method_data);

% 4) POST-PROCESSING
% 4.1) EXPORT TO PARAVIEW

output_file = 'Plate_NRB_Deg3_Reg2_Sub8';

vtk_pts = {linspace(0, 1, 21), linspace(0, 1, 21)};
fprintf ('The result is saved in the file %s \n \n', output_file);
sp_to_vtk (u, space, geometry, vtk_pts, output_file, 'u')

% 4.2) PLOT IN MATLAB. COMPARISON WITH THE EXACT SOLUTION

[eu, F] = sp_eval (u, space, geometry, vtk_pts);
[X, Y]  = deal (squeeze(F(1,:,:)), squeeze(F(2,:,:)));
subplot (1,2,1)
surf (X, Y, eu)
title ('Numerical solution'), axis tight
subplot (1,2,2)
surf (X, Y, problem_data.uex (X,Y))
title ('Exact solution'), axis tight

% Display errors of the computed solution in the L2 and H1 norm
[error_h1, error_l2] = ...
           sp_h1_error (space, msh, u, problem_data.uex, problem_data.graduex)

%!demo
%! ex_laplace_iso_plate

%!test
%! problem_data.geo_name = 'geo_plate_with_hole.txt';
%! problem_data.nmnn_sides   = [3 4];
%! problem_data.drchlt_sides = [1 2];
%! problem_data.c_diff  = @(x, y) ones(size(x));
%! problem_data.f = @(x, y) zeros(size(x));
%! problem_data.g = @test_plate_mixed_bc_g_nmnn;
%! problem_data.h = @(x, y, ind)  exp(x).*sin(y);
%! problem_data.uex     = @(x, y) exp(x).*sin (y);
%! problem_data.graduex = @(x, y) cat (1, ...
%!                        reshape (exp(x).*sin(y), [1, size(x)]), ...
%!                        reshape (exp(x).*cos(y), [1, size(x)]));
%! method_data.degree     = [3 3];       % Degree of the splines
%! method_data.regularity = [2 2];       % Regularity of the splines
%! method_data.nsub       = [8 8];       % Number of subdivisions
%! method_data.nquad      = [4 4];       % Points for the Gaussian quadrature rule
%! [geometry, msh, space, u] = solve_laplace_iso (problem_data, method_data);
%! [error_h1, error_l2] = sp_h1_error (space, msh, u, problem_data.uex, problem_data.graduex);
%! assert (msh.nel, 128)
%! assert (space.ndof, 231)
%! assert (error_h1, 6.16178098362953e-04, 1e-16)
%! assert (error_l2, 4.44336164898602e-05, 1e-16)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


function price = PROJ_ForwardStarting(N, alph, r, q, T1, T2, S_0, call, rnCHF1, rnCHF2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% About: Pricing Function for Forward Starting Options using PROJ method (uses cubic B-splines)
% Models Supported: Levy Processes, including jump diffusions and Black-Scholes model
% Returns: price of contract
% Author: Justin Lars Kirkby
%
% ----------------------
% Contract/Model Params
% ----------------------
% S_0 = initial stock price (e.g. 100)
% W   = strike  (e.g. 100)
% r   = interest rate (e.g. 0.05)
% q   = dividend yield (e.g. 0.05)
% T1  =  time to maturity (in years), e.g. T1=1
% T2  = T - T1, how much time remains in contract after the forward start date (choose 0 < T2 < T1)
% call  = 1 for call (else put)
% rnCHF1 = risk netural characteristic function of process up to T1 (function handle with single argument)
% rnCHF2 = risk netural characteristic function of process with T2 = T - T1 remaining time to maturity after forward start date
% ----------------------
% Numerical (PROJ) Params
% ----------------------
% alph  = grid with is 2*alph
% N     = budget: resolution = 2*alph/(N-1), where support is of length 2*alph
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

T = T1 + T2;

dx = 2*alph/(N-1); a = 1/dx;
dw = 2*pi/(N*dx);

omega = dw*(1:N-1);  %We calcuate coefficient of w=0 explicitly
nbar = N/2;
xmin = (1 - N/2)*dx;

%%%% Cubic Spline
b0    = 1208/2520; b1 = 1191/2520; b2 = 120/2520; b3 = 1/2520;
grand = @(w)rnCHF2(w).*(sin(w/(2*a))./w).^4./(b0 + b1*cos(w/a) +b2*cos(2*w/a) +b3*cos(3*w/a));
beta  = real(fft([1/(32*a^4) exp(-1i*xmin*omega).*feval(grand,omega)]));

%FIND Value of E[(1 - exp(X_tau))^+]
G = zeros(1,N);

G(nbar +1) = 1*(1/24 - 1/20*exp(dx)*(exp(-7/4*dx)/54 + exp(-1.5*dx)/18 + exp(-1.25*dx)/2 + 7*exp(-dx)/27));

G(nbar )   =  1*(.5 -.05*(28/27 + exp(-7/4*dx)/54 + exp(-1.5*dx)/18 + exp(-1.25*dx)/2 + 14*exp(-dx)/27 ...
            + 121/54*exp(-.75*dx) + 23/18*exp(-.5*dx) + 235/54*exp(-.25*dx)));

G(nbar -1) = 1*( 23/24 - exp(-dx)/90*( (28 + 7*exp(-dx))/3 ...
            + ( 14*exp(dx) + exp(-7/4*dx) + 242*cosh(.75*dx) + 470*cosh(.25*dx))/12 ...
            +.25*(exp(-1.5*dx) + 9*exp(-1.25*dx) + 46*cosh(.5*dx))) );

vartheta_star = 1/90*( 14/3*(2+cosh(dx)) ...
                + .5*(cosh(1.5*dx) + 9*cosh(1.25*dx) +23*cosh(.5*dx))...
                +  1/6*(cosh(7/4*dx) + 121*cosh(.75*dx) +235*cosh(.25*dx)));

G(1: nbar -2) = 1 - 1*exp(xmin +dx*(0:nbar-3))*vartheta_star;
Cons = 32*a^4;

Vbar2 = Cons/N*G(1,1:(nbar +1))*(beta(1,1:(nbar +1))');

%%% Find Second Expansion
grand = @(w)rnCHF1(w).*(sin(w/(2*a))./w).^4./(b0 + b1*cos(w/a) +b2*cos(2*w/a) +b3*cos(3*w/a));
beta  = real(fft([1/(32*a^4) exp(-1i*xmin*omega).*feval(grand,omega)]));

G(1:N) = exp(xmin + (0:(N-1))*dx);
Vbar1 = Cons/N*G*beta';
Val_put = exp(-r*T)*Vbar2*Vbar1*S_0*vartheta_star;

%%%%  PUT CALL PARITY/PRICING FORMULA
if call ==1
    Val_Proj = S_0*(exp(-q*T)-exp(-r*T2)*exp(-q*T1)) + Val_put;  %check this formula
else
    Val_Proj = Val_put;
end

price = Val_Proj;

end


text:
% Data File LAN
% Free linear vibrations of a system
% with one degree of freedom
s     = 1; % degree of freedom
L     = '1/2*(a*qt1^2 - c*q1^2)'; % Lagrangian
QN{1} = '-b*qt1'; % generalized non potential force
qj0   = 'q0';     % initial coordinate
qtj0  = 'qt0';    % initial velocity
Tend  = 20;       % upper bound of integration
eps   = 1e-10;    % desirable accuracy
np    = 3;        % number of parameters
P{1}  = 'a';      % generalized coefficient of inertia
P{2}  = 'b';      % generalized coefficient of resistance
P{3}  = 'c';      % generalized coefficient of stiffness