text:
function [f,t]=enframe(x,win,inc)
%ENFRAME split signal up into (overlapping) frames: one per row. [F,T]=(X,WIN,INC)
%
%	F = ENFRAME(X,LEN) splits the vector X(:) up into
%	frames. Each frame is of length LEN and occupies
%	one row of the output matrix. The last few frames of X
%	will be ignored if its length is not divisible by LEN.
%	It is an error if X is shorter than LEN.
%
%	F = ENFRAME(X,LEN,INC) has frames beginning at increments of INC
%	The centre of frame I is X((I-1)*INC+(LEN+1)/2) for I=1,2,...
%	The number of frames is fix((length(X)-LEN+INC)/INC)
%
%	F = ENFRAME(X,WINDOW) or ENFRAME(X,WINDOW,INC) multiplies
%	each frame by WINDOW(:)
%
%   The second output argument, T, gives the time in samples at the centre
%   of each frame. T=i corresponds to the time of sample X(i).
%

nx=length(x(:));
nwin=length(win);
if (nwin == 1)
   len = win;
else
   len = nwin;
end
if (nargin < 3)
   inc = len;
end
nf = fix((nx-len+inc)/inc);
f=zeros(nf,len);
indf= inc*(0:(nf-1)).';
inds = (1:len);
f(:) = x(indf(:,ones(1,len))+inds(ones(nf,1),:));
if (nwin > 1)
    w = win(:)';
    f = f .* w(ones(nf,1),:);
end
if nargout>1
    t=(1+len)/2+indf;
end

author: bastamon


text:
%% Lag extraction demo script
% this script runs an EEGLAB plugin used for lag extraction on single trial data.
%
% $Id: lagextraction_demo.m 4 2009-08-15 21:10:35Z gramfort $

%% Load data

USE_SYNTH = true;
USE_SYNTH = false;

clear options
close all

if USE_SYNTH

    randn('seed',0);

    clear EEG

    EEG.pnts = 500;
    EEG.srate = 1000; % 0.5 sec. recording
    EEG.trials = 100;
    EEG.nbchan = 1;
    EEG.comments = '';
    EEG.chanlocs(1).labels = 'Cz';
    EEG.chanlocs(1).X = 0;
    EEG.chanlocs(1).Y = 0;
    EEG.chanlocs(1).Z = 0;
    EEG.filename = 'synth';
    EEG.filepath = cd;

    EEG.times = linspace(0,EEG.pnts / EEG.srate,EEG.pnts);

    snr = 1.5;
    clear generate_options
    generate_options.snr = snr;

    if 1
        [EEG,ref,latencies] = generate_synth_eeg_data(EEG,generate_options);
    else
        clear generate_options

        generate_options.std_latency = 0.03;
        generate_options.mean_frequency = 3;
        generate_options.mean_latency = 0.3;

        % generate_options.mean_latency = 0.2;
        % generate_options.std_latency = 0.0;
        % generate_options.std_strech = 0.1;
        % generate_options.std_strech = 0.2;

        [EEG,ref,latencies] = generate_synth(EEG,generate_options);
    end
    EEG.times = linspace(0,1000 * EEG.pnts / EEG.srate,EEG.pnts);
    disp(['SNR : ',num2str(snr)]);

    channel = 1;
    time_win = [EEG.times(1) EEG.times(end)]; % (ms)
    bad_trials = []; % set bad trials
    use_ica = false;

else

    load('data/oddball3-num1-512Hz-chan10.set','-mat');

    channel = 1;
    time_win = [150 500]; % (ms)
    bad_trials = []; % set bad trials
    use_ica = false;

    latencies = eeg_getepochevent( EEG, {'200'}, [], 'latency');
    [tmp,motor_order] = sort(latencies,'descend');

    data = squeeze(EEG.data(1,:,:))';

    if  0
        % Hack to see result on motor response
        options.order = motor_order;
    end

end

%% Set parameters

options.distance = 'corr2';
% options.distance = 'corr'; % with corr

options.verbose = true; % show figures
options.disp_log = false; % print log messages

options.show_pca = true;
options.show_pca = false;

options.sigma = 0.01:0.01:0.2;

if 0
    % one alpha : no crossvalidation
    options.alpha = [0.1];
else
    % multiple alphas : crossvalidation (longer)
    options.alpha = [0.001,0.01,0.1];
end

%% Run extraction

[EEG, com, order, lags, event_type, E_lags] = pop_extractlag( EEG , use_ica, channel, time_win, options);

%% View results

if ~USE_SYNTH
    if options.verbose
        figure,
        imagesc(data(motor_order,:))
        title('Reordering with motor response')
        figure,
        imagesc(data(order,:))
        title('Reordering with Laplacian')
    end
end

if USE_SYNTH & options.verbose
    smart_figure('evoked_potentials');
    hold on
    plot(EEG.times,ref,'g','LineWidth',3);
    hold off
end

smart_figure('erp_data_ordered_lags');
if USE_SYNTH
    savefig(['lags_synth_alpha_',num2str(options.alpha)],22,{'pdf'})
else
    savefig(['lags_oddball_alpha_',num2str(options.alpha)],22,{'pdf'})
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function sF = vertcat(varargin)
% overloads vertcat

sF = cat(1,varargin{:});

end

author: mtex-toolbox

repo: mtex

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


text:
function featS = PyradWrapper(scan3M, mask3M, voxelSizeV, paramFilePath)

% This is a wrapper function to create a .nrrd file from CERR and call
% pyradiomics.
% Requires specifying path to Pyradiomics
%
% RKP, 03/22/2018
% AI, 06/10/2020

%Get path to pyradiomics wrapper
CERRPath = getCERRPath;
CERRPathSlashes = strfind(getCERRPath,filesep);
topLevelCERRDir = CERRPath(1:CERRPathSlashes(end-1));
pyradiomicsWrapperPath = fullfile(topLevelCERRDir,'Unit_Testing','tests_for_cerr', 'pyFeatureExtraction.py');

%Add python module to system path & iImport
pyModule = 'pyFeatureExtraction';

P = py.sys.path;
currentPath = pwd;
cd(fullfile(topLevelCERRDir,'Unit_Testing','tests_for_cerr'));

try
    if count(P,pyradiomicsWrapperPath) == 0
        %Insert paths to scipy and pyradiomics.
        insert(P,int64(0),'C:\Miniconda3\Lib\site-packages\radiomics'); %Replace with path to pyradiomic pkg
        P = py.sys.path;
        insert(P,int64(0),'C:\Miniconda3\Lib\site-packages\scipy'); %Replace with path to scipy pkg
        P = py.sys.path;
        insert(P,int64(0),pyradiomicsWrapperPath);
    end
    py.importlib.import_module(pyModule);
catch e
    error('Python module %s could not be imported. %s',pyModule,e.message);
end

%Write scan & mask to NRRD format
fprintf('\nWriting scan and mask to NRRD format...\n');
cd(currentPath);

originV = [0,0,0];
encoding = 'raw';

mask3M = uint16(mask3M);
mask3M = permute(mask3M, [2 1 3]);
%mask3M = flip(mask3M,3);

scan3M = permute(scan3M, [2 1 3]);
%scan3M = flip(scan3M,3);

scanFilename = strcat(tempdir,'scan.nrrd');
scanRes = nrrdWriter(scanFilename, scan3M, voxelSizeV, originV, encoding);

maskFilename = strcat(tempdir, 'mask.nrrd');
maskRes = nrrdWriter(maskFilename, mask3M, voxelSizeV, originV, encoding);

%Call feature extractor
try

    pyFeatDict = py.pyFeatureExtraction.extract(scanFilename, maskFilename,...
        paramFilePath, tempdir);

    %Convert python dictionary to matlab struct
    featS = struct(pyFeatDict);

catch e
    error('Feature extraction failed with message %s',e.message)
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [Eft, Varft, lpyt, Eyt, Varyt] = gpep_pred(gp, x, y, varargin)
%GPEP_PRED  Predictions with Gaussian Process EP approximation
%
%  Description
%    [EFT, VARFT] = GPEP_PRED(GP, X, Y, XT, OPTIONS)
%    takes a GP structure together with matrix X of training
%    inputs and vector Y of training targets, and evaluates the
%    predictive distribution at test inputs XT. Returns a posterior
%    mean EFT and variance VARFT of latent variables.
%
%    [EFT, VARFT, LPYT] = GPEP_PRED(GP, X, Y, XT, 'yt', YT, OPTIONS)
%    returns also logarithm of the predictive density LPYT of the
%    observations YT at test input locations XT. This can be used
%    for example in the cross-validation. Here Y has to be a vector.
%
%    [EFT, VARFT, LPYT, EYT, VARYT] = GPEP_PRED(GP, X, Y, XT, OPTIONS)
%    returns also the posterior predictive mean EYT and variance VARYT.
%
%    [EF, VARF, LPY, EY, VARY] = GPEP_PRED(GP, X, Y, OPTIONS)
%    evaluates the predictive distribution at training inputs X
%    and logarithm of the predictive density LPY of the training
%    observations Y.
%
%    OPTIONS is optional parameter-value pair
%      predcf - an index vector telling which covariance functions are
%               used for prediction. Default is all (1:gpcfn).
%               See additional information below.
%      tstind - a vector/cell array defining, which rows of X belong
%               to which training block in *IC type sparse models.
%               Default is []. In case of PIC, a cell array
%               containing index vectors specifying the blocking
%               structure for test data. IN FIC and CS+FIC a
%               vector of length n that points out the test inputs
%               that are also in the training set (if none, set
%               TSTIND = [])
%      yt     - optional observed yt in test points (see below)
%      z      - optional observed quantity in triplet (x_i,y_i,z_i)
%               Some likelihoods may use this. For example, in case of
%               Poisson likelihood we have z_i=E_i, that is, expected value
%               for ith case.
%      zt     - optional observed quantity in triplet (xt_i,yt_i,zt_i)
%               Some likelihoods may use this. For example, in case of
%               Poisson likelihood we have z_i=E_i, that is, the expected
%               value for the ith case.
%      fcorr  - Method used for latent marginal posterior corrections.
%               Default is 'off'. For EP possible method is 'fact'.
%               If method is 'on', 'fact' is used for EP.
%
%    NOTE! In case of FIC and PIC sparse approximation the
%    prediction for only some PREDCF covariance functions is just
%    an approximation since the covariance functions are coupled in
%    the approximation and are not strictly speaking additive
%    anymore.
%
%    For example, if you use covariance such as K = K1 + K2 your
%    predictions Eft1 = gpep_pred(GP, X, Y, X, 'predcf', 1) and
%    Eft2 = gpep_pred(gp, x, y, x, 'predcf', 2) should sum up to
%    Eft = gpep_pred(gp, x, y, x). That is Eft = Eft1 + Eft2. With
%    FULL model this is true but with FIC and PIC this is true only
%    approximately. That is Eft \approx Eft1 + Eft2.
%
%    With CS+FIC the predictions are exact if the PREDCF covariance
%    functions are all in the FIC part or if they are CS
%    covariances.
%
%    NOTE! When making predictions with a subset of covariance
%    functions with FIC approximation the predictive variance can
%    in some cases be ill-behaved i.e. negative or unrealistically
%    small. This may happen because of the approximative nature of
%    the prediction.
%
%  See also
%    GPEP_E, GPEP_G, GP_PRED, DEMO_SPATIAL, DEMO_CLASSIFIC

% Copyright (c) 2007-2010 Jarno Vanhatalo
% Copyright (c) 2010      Heikki Peura
% Copyright (c) 2011      Pasi Jylänki
% Copyright (c) 2012 Aki Vehtari

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

  ip=inputParser;
  ip.FunctionName = 'GPEP_PRED';
  ip.addRequired('gp', @isstruct);
  ip.addRequired('x', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addRequired('y', @(x) ~isempty(x) && isreal(x) && all(isfinite(x(:))))
  ip.addOptional('xt', [], @(x) isempty(x) || (isreal(x) && all(isfinite(x(:)))))
  ip.addParamValue('yt', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('z', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('zt', [], @(x) isreal(x) && all(isfinite(x(:))))
  ip.addParamValue('predcf', [], @(x) isempty(x) || ...
                   isvector(x) && isreal(x) && all(isfinite(x)&x>0))
  ip.addParamValue('tstind', [], @(x) isempty(x) || iscell(x) ||...
                   (isvector(x) && isreal(x) && all(isfinite(x)&x>0)))
  ip.addParamValue('fcorr', 'off', @(x) ismember(x, {'off', ...
                   'cm2', 'fact', 'on','lr'}))
  if numel(varargin)==0 || isnumeric(varargin{1})
    % inputParser should handle this, but it doesn't
    ip.parse(gp, x, y, varargin{:});
  else
    ip.parse(gp, x, y, [], varargin{:});
  end
  xt=ip.Results.xt;
  yt=ip.Results.yt;
  z=ip.Results.z;
  zt=ip.Results.zt;
  predcf=ip.Results.predcf;
  tstind=ip.Results.tstind;
  fcorr=ip.Results.fcorr;
  if isempty(xt)
    xt=x;
    if isempty(tstind)
      if iscell(gp)
        gptype=gp{1}.type;
      else
        gptype=gp.type;
      end
      switch gptype
        case {'FULL' 'VAR' 'DTC' 'SOR'}
          tstind = [];
        case {'FIC' 'CS+FIC'}
          tstind = 1:size(x,1);
        case 'PIC'
          if iscell(gp)
            tstind = gp{1}.tr_index;
          else
            tstind = gp.tr_index;
          end
      end
    end
    if isempty(yt)
      yt=y;
    end
    if isempty(zt)
      zt=z;
    end
  end

  if isfield(gp, 'monotonic') && gp.monotonic
      [gp,x,y,z,xt,zt] = gp.fh.setUpDataForMonotonic(gp,x,y,z,xt,zt);
  end

  [tn, tnin] = size(x);
  [n, nout] = size(y);

  if isfield(gp.lik, 'nondiagW')
    switch gp.type
      % ============================================================
      % FULL
      % ============================================================
      case 'FULL'
        %[e, edata, eprior, tautilde, nutilde, BKnu, B, cholP, invPBKnu]= gpep_e(gp_pak(gp), gp, x, y, 'z', z);
        [e, edata, eprior, p]= gpep_e(gp_pak(gp), gp, x, y, 'z', z);
        if isnan(e)
            Eft=NaN; Varft=NaN; lpyt=NaN; Eyt=NaN; Varyt=NaN;
            return
        end
        [nutilde, BKnu, B, cholP, invPBKnu]=deal(p.nutilde, p.BKnu, p.B, p.cholP, p.invPBKnu);

        if isfield(gp, 'comp_cf')  % own covariance for each ouput component
          multicf = true;
          if length(gp.comp_cf) ~= nout
            error('GPEP_PRED: the number of component vectors in gp.comp_cf must be the same as number of outputs.')
          end
          if ~isempty(predcf)
            if ~iscell(predcf) || length(predcf)~=nout
              error(['GPEP_PRED: if own covariance for each output component is used,'...
                     'predcf has to be cell array and contain nout (vector) elements.   '])
            end
          else
            predcf = gp.comp_cf;
          end
        else
          multicf = false;
          for i1=1:nout
            predcf2{i1} = predcf;
          end
          predcf=predcf2;
        end

        ntest=size(xt,1);
        % covariances between the training and test latents
        Kt = zeros(ntest,n,nout);
        if multicf
          for i1=1:nout
            Kt(:,:,i1) = gp_cov(gp,xt,x,predcf{i1});
          end
        else
          for i1=1:nout
            Kt(:,:,i1) = gp_cov(gp,xt,x,predcf{i1});
          end
        end

        % full ep with non-diagonal site covariances
        zz=zeros(n*nout,1);
        for k1=1:nout
          zz((1:n)+(k1-1)*n)=BKnu(:,k1)-B(:,:,k1)*invPBKnu;
        end

        %- posterior predictive mean
        Eft=zeros(ntest*nout,1);
        for z1=1:nout
          Eft((1:ntest)+(z1-1)*ntest)=Kt(:,:,z1)*(nutilde(:,z1)-zz((1:n)+(z1-1)*n));
        end

        if nargout > 1
          % posterior predictive covariance
          Covf=zeros(nout, nout, ntest);

          invcholPBKt=zeros(n,ntest,nout);
          for k1=1:nout
            invcholPBKt(:,:,k1)=cholP\(B(:,:,k1)*Kt(:,:,k1)');
          end

          %- update posterior covariance
          for k1=1:nout
            % covariances for the test latents
            kstarstar = gp_trvar(gp,xt,predcf{i1});

            Covf(k1,k1,:)=kstarstar-sum(Kt(:,:,k1)'.*(B(:,:,k1)*Kt(:,:,k1)'))'+sum(invcholPBKt(:,:,k1).*invcholPBKt(:,:,k1))';
            for j1=(k1+1):nout
              Covf(k1,j1,:)=sum(invcholPBKt(:,:,k1).*invcholPBKt(:,:,j1));
              Covf(j1,k1,:)=Covf(k1,j1,:);
            end
          end
          Varft=Covf;

        end

        % ============================================================
        % FIC
        % ============================================================
      case 'FIC'        % Predictions with FIC sparse approximation for GP
                        % ============================================================
                        % PIC
                        % ============================================================
      case {'PIC' 'PIC_BLOCK'}        % Predictions with PIC sparse approximation for GP
                                      % ============================================================
                                      % CS+FIC
                                      % ============================================================
      case 'CS+FIC'        % Predictions with CS+FIC sparse approximation for GP
    end

  else % isfield(gp.lik, 'nondiagW')
    switch gp.type
      % ============================================================
      % FULL
      % ============================================================
      case 'FULL'        % Predictions with FULL GP model
        %[e, edata, eprior, tautilde, nutilde, L] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
        if isfield(gp.lik, 'int_likparam')

          [e, edata, eprior, p] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
          if isnan(e)
              Eft=NaN; Varft=NaN; lpyt=NaN; Eyt=NaN; Varyt=NaN;
              return
          end
          [tautildee, nutildee, L, L2] = deal(p.tautilde, p.nutilde, p.L, p.La2);

          tautilde=tautildee(:,1);
          nutilde=nutildee(:,1);
          if isfield(gp.lik,'int_likparam') && gp.lik.int_likparam && ~gp.lik.inputparam
             % Give q(theta) to likelihood function to integrate ovet
            zt=[p.mf2 L2'*L2];
          end
          if isfield(gp.lik, 'int_magnitude') && gp.lik.int_magnitude && ~gp.lik.inputmagnitude
            zt=[zt p.mf3 p.La3'*p.La3];
          end
          if (isfield(gp.lik, 'int_likparam') && gp.lik.inputparam) || ...
             (isfield(gp.lik, 'int_magnitude') && gp.lik.inputmagnitude) ...
              || (isfield(gp.lik, 'int_likparam') && isfield(gp, 'comp_cf'))
            [K,C]=gp_trcov(gp,x,gp.comp_cf{1});
            kstarstar = gp_trvar(gp, xt, gp.comp_cf{1});
            K_nf=gp_cov(gp,xt,x,gp.comp_cf{1});
          else
            [K, C]=gp_trcov(gp,x);
            kstarstar = gp_trvar(gp, xt, predcf);
            K_nf=gp_cov(gp,xt,x,predcf);
          end
          ntest=size(xt,1);
          [n,nin] = size(x);

          if size(tautildee,2)==1 && all(tautilde > 0) && ~isequal(gp.latent_opt.optim_method, 'robust-EP')
            % This is the usual case where likelihood is log concave
            % for example, Poisson and probit
            sqrttautilde = sqrt(tautilde(:,1));
            Stildesqroot = sparse(1:n, 1:n, sqrttautilde, n, n);

            if ~isfield(gp,'meanf')
              if issparse(L) % If compact support covariance functions are used
                             % the covariance matrix will be sparse
                zz=Stildesqroot*ldlsolve(L,Stildesqroot*(C*nutilde));
              else
                zz=Stildesqroot*(L'\(L\(Stildesqroot*(C*nutilde))));
              end
              Eft=K_nf*(nutilde-zz);    % The mean, zero mean GP
            else
              zz = Stildesqroot*(L'\(L\(Stildesqroot*(C))));

              Eft_zm=K_nf*(nutilde-zz*nutilde); % The mean, zero mean GP
              Ks = eye(size(zz)) - zz;           % inv(K + S^-1)*S^-1
              Ksy = Ks*nutilde;
              [RB RAR] = mean_predf(gp,x,xt,K_nf',Ks,Ksy,'EP',Stildesqroot.^2);

              Eft = Eft_zm + RB;            % The mean
            end

            % Compute variance
            if nargout > 1
              if issparse(L)
                V = ldlsolve(L, Stildesqroot*K_nf');
                Varft = kstarstar - sum(K_nf.*(Stildesqroot*V)',2);
              else
                V = (L\Stildesqroot)*K_nf';
                Varft = kstarstar - sum(V.^2)';
              end
              if isfield(gp,'meanf')
                Varft = Varft + RAR;
              end
            end
          else
            % We might end up here if the likelihood is not log concave
            % For example Student-t likelihood.

            %{
            zz=tautilde.*(L'*(L*nutilde));
            Eft=K_nf*(nutilde-zz);

            if nargout > 1
              S = diag(tautilde);
              V = K_nf*S*L';
              Varft = kstarstar - sum((K_nf*S).*K_nf,2) + sum(V.^2,2);
            end
            %}

            % An alternative implementation for avoiding negative variances
            [Eft,V]=pred_var(tautilde,K,K_nf,nutilde);
            Varft=kstarstar-V;

          end
          if isfield(gp.lik, 'int_likparam') && gp.lik.int_likparam && gp.lik.inputparam
            tautilde=tautildee(:,2);
            nutilde=nutildee(:,2);
            [K, C]=gp_trcov(gp,x, gp.comp_cf{2});
            kstarstar = gp_trvar(gp, xt, gp.comp_cf{2});
            K_nf=gp_cov(gp,xt,x,gp.comp_cf{2});

            [Eft(:,2),V]=pred_var(tautilde,K,K_nf,nutilde);
            Varft(:,2)=kstarstar-V;
          end
          if isfield(gp.lik, 'int_magnitude') && gp.lik.int_magnitude && gp.lik.inputmagnitude
            tautilde=tautildee(:,end);
            nutilde=nutildee(:,end);
            [K, C]=gp_trcov(gp,x, gp.comp_cf{end});
            kstarstar = gp_trvar(gp, xt, gp.comp_cf{end});
            K_nf=gp_cov(gp,xt,x,gp.comp_cf{end});

            [Eft(:,end+1),V]=pred_var(tautilde,K,K_nf,nutilde);
            Varft(:,end+1)=kstarstar-V;
          end

        else % isfield(gp.lik, 'int_likparam')

          [e, edata, eprior, p] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
          if isnan(e)
              Eft=NaN; Varft=NaN; lpyt=NaN; Eyt=NaN; Varyt=NaN;
              return
          end
          [tautilde, nutilde, L] = deal(p.tautilde, p.nutilde, p.L);

          [K, C]=gp_trcov(gp,x);
          kstarstar = gp_trvar(gp, xt, predcf);
          ntest=size(xt,1);
          K_nf=gp_cov(gp,xt,x,predcf);
          [n,nin] = size(x);

          if all(tautilde > 0) ...
                  && ~(isequal(gp.latent_opt.optim_method, 'robust-EP') )
            % This is the usual case where likelihood is log concave
            % for example, Poisson and probit
            sqrttautilde = sqrt(tautilde);
            nstt=length(sqrttautilde);
            Stildesqroot = sparse(1:nstt, 1:nstt, sqrttautilde, nstt,  nstt);

            if ~isfield(gp,'meanf')
              if issparse(L)          % If compact support covariance functions a  re used
                                      % the covariance matrix will be sparse
                zz=Stildesqroot*ldlsolve(L,Stildesqroot*(C*nutilde));
              else
                zz=Stildesqroot*(L'\(L\(Stildesqroot*(C*nutilde))));
              end

              Eft=K_nf*(nutilde-zz);    % The mean, zero mean GP
            else
              zz = Stildesqroot*(L'\(L\(Stildesqroot*(C))));

              Eft_zm=K_nf*(nutilde-zz*nutilde); % The mean, zero mean GP
              Ks = eye(size(zz)) - zz;           % inv(K + S^-1)*S^-1
              Ksy = Ks*nutilde;
              [RB RAR] = mean_predf(gp,x,xt,K_nf',Ks,Ksy,'EP',Stildesqroot.^2);

              Eft = Eft_zm + RB;            % The mean
            end

            % Compute variance
            if nargout > 1
              if issparse(L)
                V = ldlsolve(L, Stildesqroot*K_nf');
                Varft = kstarstar - sum(K_nf.*(Stildesqroot*V)',2);
              else
                V = (L\Stildesqroot)*K_nf';
                Varft = kstarstar - sum(V.^2)';
              end
              if isfield(gp,'meanf')
                Varft = Varft + RAR;
              end
            end
          else
            % We might end up here if the likelihood is not log concave
            % For example Student-t likelihood.

            %{
            zz=tautilde.*(L'*(L*nutilde));
            Eft=K_nf*(nutilde-zz);

            if nargout > 1
              S = diag(tautilde);
              V = K_nf*S*L';
              Varft = kstarstar - sum((K_nf*S).*K_nf,2) + sum(V.^2,2);
            end
            %}

            % An alternative implementation for avoiding negative variances
            [Eft,V]=pred_var(tautilde,K,K_nf,nutilde);
            Varft=kstarstar-V;

          end

        end % isfield(gp.lik, 'int_likparam')
        % ============================================================
        % FIC
        % ============================================================
      case 'FIC'        % Predictions with FIC sparse approximation for GP
                        %[e, edata, eprior, tautilde, nutilde, L, La, b] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
        [e, edata, eprior, p] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
        if isnan(e)
            Eft=NaN; Varft=NaN; lpyt=NaN; Eyt=NaN; Varyt=NaN;
            return
        end
        [tautilde, nutilde, L, La, b] = deal(p.tautilde, p.nutilde, p.L, p.La2, p.b);

        % Here tstind = 1 if the prediction is made for the training set
        if nargin > 6
          if ~isempty(tstind) && length(tstind) ~= size(x,1)
            error('tstind (if provided) has to be of same length as x.')
          end
        else
          tstind = [];
        end

        u = gp.X_u;
        m = size(u,1);

        K_fu = gp_cov(gp,x,u,predcf);          % f x u
        K_nu=gp_cov(gp,xt,u,predcf);
        K_uu = gp_trcov(gp,u,predcf);          % u x u, noiseless covariance K_uu
        K_uu = (K_uu+K_uu')./2;                % ensure the symmetry of K_uu

        kstarstar=gp_trvar(gp,xt,predcf);

        if all(tautilde > 0) && ~isequal(gp.latent_opt.optim_method, 'robust-EP')

          % From this on evaluate the prediction
          % See Snelson and Ghahramani (2007) for details
          %        p=iLaKfu*(A\(iLaKfu'*mutilde));
          p = b';

          ntest=size(xt,1);

          Eft = K_nu*(K_uu\(K_fu'*p));

          % if the prediction is made for training set, evaluate Lav also for prediction points
          if ~isempty(tstind)
            [Kv_ff, Cv_ff] = gp_trvar(gp, xt(tstind,:), predcf);
            Luu = chol(K_uu)';
            B=Luu\(K_fu');
            Qv_ff=sum(B.^2)';
            Lav = Kv_ff-Qv_ff;
            Eft(tstind) = Eft(tstind) + Lav.*p;
          end

          % Compute variance
          if nargout > 1
            %Varft(i1,1)=kstarstar(i1) - (sum(Knf(i1,:).^2./La') - sum((Knf(i1,:)*L).^2));
            Luu = chol(K_uu)';
            B=Luu\(K_fu');
            B2=Luu\(K_nu');
            Varft = kstarstar - sum(B2'.*(B*(repmat(La,1,m).\B')*B2)',2)  + sum((K_nu*(K_uu\(K_fu'*L))).^2, 2);

            % if the prediction is made for training set, evaluate Lav also for prediction points
            if ~isempty(tstind)
              Varft(tstind) = Varft(tstind) - 2.*sum( B2(:,tstind)'.*(repmat((La.\Lav),1,m).*B'),2) ...
                  + 2.*sum( B2(:,tstind)'*(B*L).*(repmat(Lav,1,m).*L), 2)  ...
                  - Lav./La.*Lav + sum((repmat(Lav,1,m).*L).^2,2);
            end
          end

        else
          % Robust-EP
          [Eft,V]=pred_var2(tautilde,nutilde,L,K_uu,K_fu,b,K_nu);
          Varft=kstarstar-V;

        end

        % ============================================================
        % PIC
        % ============================================================
      case {'PIC' 'PIC_BLOCK'}        % Predictions with PIC sparse approximation for GP
                                      % Calculate some help matrices
        u = gp.X_u;
        ind = gp.tr_index;
        %[e, edata, eprior, tautilde, nutilde, L, La, b] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
        [e, edata, eprior, p] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
        if isnan(e)
            Eft=NaN; Varft=NaN; lpyt=NaN; Eyt=NaN; Varyt=NaN;
            return
        end
        [L, La, b] = deal(p.L, p.La2, p.b);

        K_fu = gp_cov(gp, x, u, predcf);         % f x u
        K_nu = gp_cov(gp, xt, u, predcf);         % n x u
        K_uu = gp_trcov(gp, u, predcf);    % u x u, noiseles covariance K_uu

        % From this on evaluate the prediction
        % See Snelson and Ghahramani (2007) for details
        %        p=iLaKfu*(A\(iLaKfu'*mutilde));
        p = b';

        iKuuKuf = K_uu\K_fu';

        w_bu=zeros(length(xt),length(u));
        w_n=zeros(length(xt),1);
        for i=1:length(ind)
          w_bu(tstind{i},:) = repmat((iKuuKuf(:,ind{i})*p(ind{i},:))', length(tstind{i}),1);
          K_nf = gp_cov(gp, xt(tstind{i},:), x(ind{i},:), predcf);              % n x u
          w_n(tstind{i},:) = K_nf*p(ind{i},:);
        end

        Eft = K_nu*(iKuuKuf*p) - sum(K_nu.*w_bu,2) + w_n;

        % Compute variance
        if nargout > 1
          kstarstar = gp_trvar(gp, xt, predcf);
          KnfL = K_nu*(iKuuKuf*L);
          Varft = zeros(length(xt),1);
          for i=1:length(ind)
            v_n = gp_cov(gp, xt(tstind{i},:), x(ind{i},:), predcf);              % n x u
            v_bu = K_nu(tstind{i},:)*iKuuKuf(:,ind{i});
            KnfLa = K_nu*(iKuuKuf(:,ind{i})/chol(La{i}));
            KnfLa(tstind{i},:) = KnfLa(tstind{i},:) - (v_bu + v_n)/chol(La{i});
            Varft = Varft + sum((KnfLa).^2,2);
            KnfL(tstind{i},:) = KnfL(tstind{i},:) - v_bu*L(ind{i},:) + v_n*L(ind{i},:);
          end
          Varft = kstarstar - (Varft - sum((KnfL).^2,2));

        end
        % ============================================================
        % CS+FIC
        % ============================================================
      case 'CS+FIC'        % Predictions with CS+FIC sparse approximation for GP
                           % Here tstind = 1 if the prediction is made for the training set
        if nargin > 6
          if ~isempty(tstind) && length(tstind) ~= size(x,1)
            error('tstind (if provided) has to be of same length as x.')
          end
        else
          tstind = [];
        end

        u = gp.X_u;
        m = length(u);
        n = size(x,1);
        n2 = size(xt,1);

        %[e, edata, eprior, tautilde, nutilde, L, La, b] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
        [e, edata, eprior, p] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
        if isnan(e)
            Eft=NaN; Varft=NaN; lpyt=NaN; Eyt=NaN; Varyt=NaN;
            return
        end
        [L, La, b] = deal(p.L, p.La2, p.b);

        % Indexes to all non-compact support and compact support covariances.
        cf1 = [];
        cf2 = [];
        % Indexes to non-CS and CS covariances, which are used for predictions
        predcf1 = [];
        predcf2 = [];

        ncf = length(gp.cf);
        % Loop through all covariance functions
        for i = 1:ncf
          % Non-CS covariances
          if ~isfield(gp.cf{i},'cs')
            cf1 = [cf1 i];
            % If used for prediction
            if ~isempty(find(predcf==i))
              predcf1 = [predcf1 i];
            end
            % CS-covariances
          else
            cf2 = [cf2 i];
            % If used for prediction
            if ~isempty(find(predcf==i))
              predcf2 = [predcf2 i];
            end
          end
        end
        if isempty(predcf1) && isempty(predcf2)
          predcf1 = cf1;
          predcf2 = cf2;
        end

        % Determine the types of the covariance functions used
        % in making the prediction.
        if ~isempty(predcf1) && isempty(predcf2)       % Only non-CS covariances
          ptype = 1;
          predcf2 = cf2;
        elseif isempty(predcf1) && ~isempty(predcf2)   % Only CS covariances
          ptype = 2;
          predcf1 = cf1;
        else                                           % Both non-CS and CS covariances
          ptype = 3;
        end

        K_fu = gp_cov(gp,x,u,predcf1);   % f x u
        K_uu = gp_trcov(gp,u,predcf1);     % u x u, noiseles covariance K_uu
        K_uu = (K_uu+K_uu')./2;     % ensure the symmetry of K_uu
        K_nu=gp_cov(gp,xt,u,predcf1);

        Kcs_nf = gp_cov(gp, xt, x, predcf2);

        p = b';
        ntest=size(xt,1);

        % Calculate the predictive mean according to the type of
        % covariance functions used for making the prediction
        if ptype == 1
          Eft = K_nu*(K_uu\(K_fu'*p));
        elseif ptype == 2
          Eft = Kcs_nf*p;
        else
          Eft = K_nu*(K_uu\(K_fu'*p)) + Kcs_nf*p;
        end

        % evaluate also Lav if the prediction is made for training set
        if ~isempty(tstind)
          [Kv_ff, Cv_ff] = gp_trvar(gp, xt(tstind,:), predcf1);
          Luu = chol(K_uu)';
          B=Luu\(K_fu');
          Qv_ff=sum(B.^2)';
          Lav = Kv_ff-Qv_ff;
        end

        % Add also Lav if the prediction is made for training set
        % and non-CS covariance function is used for prediction
        if ~isempty(tstind) && (ptype == 1 || ptype == 3)
          Eft(tstind) = Eft(tstind) + Lav.*p;
        end

        % Evaluate the variance
        if nargout > 1
          Knn_v = gp_trvar(gp,xt,predcf);
          Luu = chol(K_uu)';
          B=Luu\(K_fu');
          B2=Luu\(K_nu');
          p = amd(La);
          iLaKfu = La\K_fu;
          % Calculate the predictive variance according to the type
          % covariance functions used for making the prediction
          if ptype == 1 || ptype == 3
            % FIC part of the covariance
            Varft = Knn_v - sum(B2'.*(B*(La\B')*B2)',2) + sum((K_nu*(K_uu\(K_fu'*L))).^2, 2);
            % Add Lav2 if the prediction is made for the training set
            if  ~isempty(tstind)
              % Non-CS covariance
              if ptype == 1
                Kcs_nf = sparse(tstind,1:n,Lav,n2,n);
                % Non-CS and CS covariances
              else
                Kcs_nf = Kcs_nf + sparse(tstind,1:n,Lav,n2,n);
              end
              % Add Lav2 inside Kcs_nf
              Varft = Varft - sum((Kcs_nf(:,p)/chol(La(p,p))).^2,2) + sum((Kcs_nf*L).^2, 2) ...
                      - 2.*sum((Kcs_nf*iLaKfu).*(K_uu\K_nu')',2) + 2.*sum((Kcs_nf*L).*(L'*K_fu*(K_uu\K_nu'))' ,2);
              % In case of both non-CS and CS prediction covariances add
              % only Kcs_nf if the prediction is not done for the training set
            elseif ptype == 3
              Varft = Varft - sum((Kcs_nf(:,p)/chol(La(p,p))).^2,2) + sum((Kcs_nf*L).^2, 2) ...
                      - 2.*sum((Kcs_nf*iLaKfu).*(K_uu\K_nu')',2) + 2.*sum((Kcs_nf*L).*(L'*K_fu*(K_uu\K_nu'))' ,2);
            end
            % Prediction with only CS covariance
          elseif ptype == 2
            Varft = Knn_v - sum((Kcs_nf(:,p)/chol(La(p,p))).^2,2) + sum((Kcs_nf*L).^2, 2) ;
          end
        end
        % ============================================================
        % DTC/(VAR)
        % ============================================================
      case {'DTC' 'VAR' 'SOR'}        % Predictions with DTC or variational sparse approximation for GP
                                      %[e, edata, eprior, tautilde, nutilde, L, La, b] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
        [e, edata, eprior, p] = gpep_e(gp_pak(gp), gp, x, y, 'z', z);
        if isnan(e)
            Eft=NaN; Varft=NaN; lpyt=NaN; Eyt=NaN; Varyt=NaN;
            return
        end
        [L, La, b] = deal(p.L, p.La2, p.b);

        % Here tstind = 1 if the prediction is made for the training set
        if nargin > 6
          if ~isempty(tstind) && length(tstind) ~= size(x,1)
            error('tstind (if provided) has to be of same length as x.')
          end
        else
          tstind = [];
        end

        u = gp.X_u;
        m = size(u,1);

        K_fu = gp_cov(gp,x,u,predcf);         % f x u
        K_nu=gp_cov(gp,xt,u,predcf);
        K_uu = gp_trcov(gp,u,predcf);          % u x u, noiseles covariance K_uu
        K_uu = (K_uu+K_uu')./2;          % ensure the symmetry of K_uu

        kstarstar=gp_trvar(gp,xt,predcf);

        % From this on evaluate the prediction
        p = b';

        ntest=size(xt,1);

        Eft = K_nu*(K_uu\(K_fu'*p));

        % if the prediction is made for training set, evaluate Lav also for prediction points
        if ~isempty(tstind)
          [Kv_ff, Cv_ff] = gp_trvar(gp, xt(tstind,:), predcf);
          Luu = chol(K_uu)';
          B=Luu\(K_fu');
          Qv_ff=sum(B.^2)';
          Lav = Kv_ff-Cv_ff;
          Eft(tstind) = Eft(tstind);% + Lav.*p;
        end

        if nargout > 1
          % Compute variances of predictions
          %Varft(i1,1)=kstarstar(i1) - (sum(Knf(i1,:).^2./La') - sum((Knf(i1,:)*L).^2));
          Luu = chol(K_uu)';
          B=Luu\(K_fu');
          B2=Luu\(K_nu');

          Varft = sum(B2'.*(B*(repmat(La,1,m).\B')*B2)',2)  + sum((K_nu*(K_uu\(K_fu'*L))).^2, 2);
          switch gp.type
            case {'VAR' 'DTC'}
              Varft = kstarstar - Varft;
            case 'SOR'
              Varft = sum(B2.^2,1)' - Varft;
          end
        end
    end
  end
  if ~isequal(fcorr, 'off')
    % Do marginal corrections for samples
    [pc_predm, fvecm] = gp_predcm(gp, x, y, xt, 'z', z, 'ind', 1:size(xt,1), 'fcorr', fcorr);
    for i=1:size(xt,1)
      % Remove NaNs and zeros
      pc_pred=pc_predm(:,i);
      dii=isnan(pc_pred)|pc_pred==0;
      pc_pred(dii)=[];
      fvec=fvecm(:,i);
      fvec(dii)=[];
      % Compute mean correction
      Eft(i) = trapz(fvec.*(pc_pred./sum(pc_pred)));
    end
   end

  % ============================================================
  % Evaluate also the predictive mean and variance of new observation(s)
  % ============================================================
  if ~isequal(fcorr, 'off')
    if nargout == 3
      if isempty(yt)
        lpyt=[];
      else
        lpyt = gp.lik.fh.predy(gp.lik, fvecm', pc_predm', yt, zt);
      end
    elseif nargout > 3
      [lpyt, Eyt, Varyt] = gp.lik.fh.predy(gp.lik, fvecm', pc_predm', yt, zt);
    end
  else
    if nargout == 3
      if isempty(yt)
        lpyt=[];
      else
        lpyt = gp.lik.fh.predy(gp.lik, Eft, Varft, yt, zt);
      end
    elseif nargout > 3
      [lpyt, Eyt, Varyt] = gp.lik.fh.predy(gp.lik, Eft, Varft, yt, zt);
    end
  end
end

function [m,S]=pred_var(tau_q,K,A,b)

% helper function for determining
%
% m = A * inv( K+ inv(diag(tau_q)) ) * inv(diag(tau_q)) *b
% S = diag( A * inv( K+ inv(diag(tau_q)) ) * A)
%
% when the site variances tau_q may be negative
%

  ii1=find(tau_q>0); n1=length(ii1); W1=sqrt(tau_q(ii1));
  ii2=find(tau_q<0); n2=length(ii2); W2=sqrt(abs(tau_q(ii2)));

  m=A*b;
  b=K*b;
  S=zeros(size(A,1),1);
  u=0;
  U=0;
  L1=[];
  if ~isempty(ii1)
    % Cholesky decomposition for the positive sites
    L1=(W1*W1').*K(ii1,ii1);
    L1(1:n1+1:end)=L1(1:n1+1:end)+1;
    L1=chol(L1);

    U = bsxfun(@times,A(:,ii1),W1')/L1;
    u = L1'\(W1.*b(ii1));

    m = m-U*u;
    S = S+sum(U.^2,2);
  end

  if ~isempty(ii2)
    % Cholesky decomposition for the negative sites
    V=bsxfun(@times,K(ii2,ii1),W1')/L1;
    if isempty(V)
        V=0;
    end
    L2=(W2*W2').*(V*V'-K(ii2,ii2));
    L2(1:n2+1:end)=L2(1:n2+1:end)+1;

    [L2,pd]=chol(L2);
    if pd==0
      U = bsxfun(@minus, bsxfun(@times,A(:,ii2),W2')/L2,U*(bsxfun(@times,V,W2)'/L2));
      u = L2'\(W2.*b(ii2)) -L2'\(bsxfun(@times,V,W2)*u);

      m = m+U*u;
      S = S-sum(U.^2,2);
    else
      fprintf('Posterior covariance is negative definite.\n')
    end
  end

end

function [m_q,S_q]=pred_var2(tautilde,nutilde,L,K_uu,K_fu,D,K_nu)

% function for determining the parameters of the q-distribution
% when site variances tau_q may be negative
%
% q(f) = N(f|0,K)*exp( -0.5*f'*diag(tau_q)*f + nu_q'*f )/Z_q = N(f|m_q,S_q)
%
% S_q = inv(inv(K)+diag(tau_q)) where K is sparse approximation for prior
%       covariance
% m_q = S_q*nu_q;
%
% det(eye(n)+K*diag(tau_q))) = det(L1)^2 * det(L2)^2
% where L1 and L2 are upper triangular
%
% see Expectation consistent approximate inference (Opper & Winther, 2005)

  n=length(nutilde);

  U = K_fu;
  S = 1+tautilde.*D;
  B = tautilde./S;
  BUiL = bsxfun(@times, B, U)/L';
  % iKS = diag(B) - BUiL*BUiL';

  Ktnu = D.*nutilde + U*(K_uu\(U'*nutilde));
  m_q = nutilde - B.*Ktnu + BUiL*(BUiL'*Ktnu);
  kstar = K_nu*(K_uu\K_fu');
  m_q = kstar*m_q;

  S_q = sum(bsxfun(@times,B',kstar.^2),2) - sum((kstar*BUiL).^2,2);
  % S_q = kstar*iKS*kstar';

end

author: gpstuff-dev

repo: gpstuff

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


text:
%% Ex. 13 Making a quick plot
x = [0:0.1:20];
y = sin(x);
plot(x,y)

% remaks : Remarks: This only serves as a very quick example of what Matlab can do in making
% plots.The first line is equivalent to x = [0 0.1 0.2 0.3 ... 19.8 19.9 20]. It
% assigns the content of x which is an array of 201 elements. The "0:0.1:20" means the
% 201 numbers are evenly spaced. They start from 0 and end at 20 with an increment of
% 0.1. The second line gives the content of the new array, y, as
 % y = [sin(x(1)) sin(x(2)) sin(x(3)) ... sin(x(200)) sin(x(201))] ,
% or
% y = [sin(0) sin(0.1) sin(0.2) ... sin(19.9) sin(20)] .
% The 3rd line makes a plot of y vs. x.

author: TheAlgorithms

repo: MATLAB-Octave


text:
function showmovi()
    %  this is the file showmovi.m. It displays a movie with the variable
    %  name 'm'.
    %   Stefan Wiemer 11/94
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    report_this_filefun();

    % This is the info window text
    %
    ttlStrmov='The Movie Window                                ';
    hlpStr1mov= ...
        ['                                                '
        ' This window displays a movie, a series of      '
        ' equivally spaced time cuts. Dispalyed are      '
        ' z-values of the selected function (e.g AS(t))  '
        ' in map view. The colorbar scaling is the same  '
        ' for eah frame, the maximum and minimum are the '
        ' overll maximum for all frames.                 '
        ' Menu options:                                  '
        '                                                '
        ' Circle: select the ni closest earthquakes to a '
        '       point selected with the mouse.           '
        ' Play:  Plays the movie. First the movie is     '
        '       loaded into the memory, then it is played'
        '       n-times (depending on the setting of the '
        '       <# of run> input box)                    '
        ' Speed: This is the number of frames per second '
        '       displayed. Your computer may not be able '
        '       to display the movie in high speed.       '
        ' Forward one frame (>) : Displays the next frame'
        ' Backward one frame: Displays the previous frame'
        ' Colormap: Seelect one of the colormaps in the  '
        '       pulldownmenu.                            '
        '                                                '];

    % find out if figure exists
    mov=findobj('Type','Figure','-and','Name','Movie Window');

    % Set up the Movie window Enviroment
    %
    %if isempty(mov)
    mov =   figure_w_normalized_uicontrolunits( ...
        'Name','Movie Window',...
        'NumberTitle','off', ...
        'NextPlot','new', ...
        'Visible','on', ...
        'Position',fs_m);

    %end % if exist

    figure(mov);
    set(gca,'NextPlot','add')
    mov = gcf;
    whitebg(color_fbg);
    h = gcf;
    speed = 0.5;
    clf reset
    drawnow;
    set(gcf,'resize','off');
    set(gca,'pos',rect1)

    axis off;
    cur_color = 'jet';
    colormap(jet);
    m0 = uicontrol('style','text','unit','norm','pos',[.06 .5 .8 .1]);
    set(m0,'string','Please Wait .... loading Data');
    set(m0,'background',color_fbg);
    drawnow;

    delete(m0);
    [a1,b1] = size(m);

    m1 = 'Loading movie to the Graphics Server ... please wait';

    movie(m(:,1),1,30);

    rect = [0.15 0.90 0.50 0.05 ];
    axes('position',rect)
    pco2 = pcolor([ZG.minc:0.1:ZG.maxc ; ZG.minc:0.1:ZG.maxc]);
    set(gca,'visible','on')
    h4 = gca;
    set(h4,'YTick',[-10 10])
    set(h4,'XTick',[-1000 1000])
    set(h4,'FontWeight','bold')
    shading flat

    rect = [0.15 0.90 0.50 0.05 ];
    pco5 = axes('position',rect);
    h5 = gca;
    axis([ ZG.minc ZG.maxc 0 1  ])
    set(h5,'YTick',[-10 10])
    set(h5,'FontWeight','bold')

    %rect = rect1;
    rect = [0.05 0.20 0.82 0.73 ];

    pco6 = axes('position',rect1);
    axis([min(gx) max(gx) min(gy)  max(gy) ])
    set(gca,'NextPlot','add')
    h6 = gca;
    hmo = gca;
    set(h6,'FontWeight','bold')

    movie(m(:,1),1,30);

    cs = uicontrol('style','popupmenu','string','HSV|Hot|Cool|Pink|Bone|Gary|Jet');
    set(cs,'unit','norm','pos',[.763 .05 .202 .05],...
        'callback',@callbackfun_001);
    h2 = uicontrol('style','text','unit','norm','string','Colormap');
    set(h2,'pos',[.763 .1 .2 .05],'background',color_fbg);

    frame_slide = uicontrol('style','slider','max',b1,'min',1,'uni','norm');
    set(frame_slide,'value',1,'pos',[.750 .25 .04 .5], 'callback',@callbackfun_002);
    frame = uicontrol('style','edit','value',10,'string',num2str(i), 'Callback', @callbackfun_003);
    flabel = uicontrol('style','text','units','norm','pos',[.55 .13 .2 .05]);
    set(flabel,'string','Speed','background',color_fbg);
    set(frame,'units','norm','pos',[.60 .07 .1 .05],'min',0.1,'max',30);

    uicontrol('style','text','units','norm','pos',[.80 .90 .20 .05],...
        'String','Forward 1 ','background',color_fbg);
    uicontrol('style','text','units','norm','pos',[.80 .80 .20 .05],...
        'String','Backward 1 ','background',color_fbg);

    next = uicontrol('style','pushbutton','unit','norm','pos',[0.75 .90 .04 .05]);
    set(next,'string','>','ForeGroundColor','k');
    set(next, 'callback',@callbackfun_004);
    bac = uicontrol('style','pushbutton','unit','norm','pos',[0.75 .80 .04 .05]);
    set(bac,'string','<');
    set(bac, 'callback',@callbackfun_005);
    time = uicontrol('style','edit','value',3,'string',num2str(3), 'Callback', @callbackfun_006);
    set(time,'units','norm','pos',[.23 .07 .1 .05],'min',1,'max',1000);
    tlabel = uicontrol('style','text','units','norm','pos',[.18 .13 .2 .05]);
    set(tlabel,'string','# of runs','background',color_fbg);
    start = uicontrol('style','pushbutton','unit','norm','pos',[.06 .05 .15 .1]);
    set(start,'interruptible','yes','string','Play');
    set(start, 'callback',@callbackfun_007);
    mc = 'close(gcf)';

    circ = uicontrol('units','norm','pos',[.40 .10 .10 .1],'style','pushbutton');
    set(circ,'string','Circle', 'callback',@callbackfun_008);

    set_ni = uicontrol('style','edit','value',100,'string',num2str(100));

    set(set_ni, 'Callback', @callbackfun_009);

    set(set_ni,'units','norm','pos',[.40 .02 .15 .05],'min',10,'max',10000);
    nilabel = uicontrol('style','text','units','norm','pos',[.36 .02 .04 .05]);
    set(nilabel,'string','ni:','background',color_fbg);

    uicontrol('Units','normal',...
        'Position',[.0 .95 .08 .06],'String','Close ',...
        'callback',@callbackfun_010)

    uicontrol('Units','normal',...
        'Position',[.0 .85 .08 .06],'String','Info ',...
        'callback',@callbackfun_011)

    uicontrol('Units','normal',...
        'Position',[.0 .75 .08 .06],'String','Refresh ',...
        'callback',@callbackfun_012);

    %si = signatur('ZMAP','',[0.01 0.02]);
    %set(si,'Color','k')

    axes(h6);
    watchoff(mess)
    watchoff(mov)

    function callbackfun_001(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        v = get(cs,'value');
        cur_color=lower(v);
        colormap(cur_color);
        movie(m(:,i),1,1);
    end

    function callbackfun_002(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        i=(get(frame_slide,'Value'));
        movie(m(:,i),1,1);
    end

    function callbackfun_003(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        speed=str2double(frame.String);
    end

    function callbackfun_004(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        i=i+1;
        if i > b1
            i=b1;
        end
        movie(m(:,i),1,1);
        frame_slide.Value=i;
    end

    function callbackfun_005(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        i=i-1;
        if i < 1
            i=1;
        end
        movie(m(:,i),1,1);
        frame_slide.Value=i;
    end

    function callbackfun_006(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        set(time,'Value',str2double(get(time,'string')));
    end

    function callbackfun_007(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        disp(m1);
        movie(m,fix(get(time,'value')*15/size(m,2)),speed);
    end

    function callbackfun_008(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        circmo;
    end

    function callbackfun_009(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        ni=str2double(set_ni.String);
        set_ni.String=num2str(ni);
    end

    function callbackfun_010(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        clear m;
        close(mov);

    end

    function callbackfun_011(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        zmaphelp(ttlStrmov,hlpStr1mov);
    end

    function callbackfun_012(mysrc,myevt)

        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        close;
        showmovi;
    end

end

author: CelsoReyes


text:
function [average, weight] = CrossWeightAverage(DATA,WEIGHT)
% =========================================================================
% Compute the weighted average of DATA, computing the weights form the
% matrix WEIGHT. Note that if there are NaNs in DATA the weights are
% rescaled. WEIGHTS CANNOT contain NaN
% =========================================================================
% [average, weight] = CrossWeightAverage(DATA,WEIGHT)
% -------------------------------------------------------------------------
% INPUT
%	- DATA: panel of time series DATA(T,N) of T observation x N series
%       (NaN are accepted and weights rescaled)
%   - WEIGHT: time series of FLOWS or WEIGHTS (TxN). You can also
%       provide a (1xN) matrix in this case the code would automatically
%       compute fixed weights
% -------------------------------------------------------------------------
% OUTPUT
%	- average : weighted average vector (T x 1)
%	- weight  : matrix of weights (T x N)
% =========================================================================
% Example 1: foreign variables as in GVAR (example with fixed weights)
% x = [1 NaN  3 ;
%      2   3  4 ;
%      3   4  5 ;
%      4   5  6];
% w = [0  .4 .6];
% [average, weight] = CrossWeightAverage(x,w)
% -------------------------------------------------------------------------
% Example 2: weighted average (with time-varying weights)
% x = [ 1 NaN   3 ;
%       2   3   4 ;
%       3   4   5 ;
%       4   5   6];
% w = [.4  .2  .4 ;
%      .3  .4  .3 ;
%      .4  .3  .3 ;
%      .5  .3  .2];
% [average, weight] = CrossWeightAverage(x,w)
% =========================================================================
% Ambrogio Cesa Bianchi, March 2015

%% Preliminaries
% Check inputs
[r1, c1] = size(DATA);
[r2, c2] = size(WEIGHT);
if r2==1 && c2==c1
    % Create a matrix of fixed weights
    aux = [];
    for ii=1:r1
        aux = [aux; WEIGHT];
    end
    WEIGHT = aux;
end

% Re-Check inputs
[r1, c1] = size(DATA);
[r2, c2] = size(WEIGHT);
if r1~=r2 || c1~=c2
    disp(' ')
    disp('The matrices of data and weights are not conformable.')
    return
end

%% Cross section weihgted average
[nobs, nvars] = size(DATA);
nans = isnan(DATA);
WEIGHT(nans) = NaN;

average(1:nobs,1) = NaN;
weight(1:nobs,1:nvars) = NaN;

for ii=1:nobs
    aux_weight = WEIGHT(ii,:)./nansum(WEIGHT(ii,:));
    weight(ii,:) = aux_weight;
    aux_weight(nans(ii,:)) = 0;
    DATA(ii,nans(ii,:)) = 0;
    average(ii,1) = DATA(ii,:) * aux_weight';
end

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function w = load_seisan(request)
   %LOAD_SEISAN loads a waveform from SEISAN files
   % combineWaves isn't currently used!

   % Glenn Thompson 2016/05/25 based on load_miniseed
   % request.combineWaves is ignored
   w=[];
   if isstruct(request)
      [thisSource, chanInfo, startTimes, endTimes, ~] = unpackDataRequest(request);
      for i=1:numel(chanInfo)
         for j=1:numel(startTimes)
            thisfilename = getfilename(thisSource,chanInfo(i),startTimes(j))
            thisw = load_seisan_file(thisfilename{1}); %, startTimes(j), endTimes(j));
            w = [w thisw];
         end
      end
      w = w(:);

   else
      %request should be a filename
      thisFilename = request;
      if exist(thisFilename, 'file')
        w = load_seisan_file(thisFilename);
      else
          w = waveform();
          warning(sprintf('File %s does not exist',thisFilename));
      end
   end
end

function w = load_seisan_file(fn)%dataRequest)
   % Read a v7.0 or later SEISAN file
   w  = waveform; %initialize, in case nothing works
   MACHINEFORMAT = 'l';
   %disp('Little Endian');
   fid = fopen(fn,'r',MACHINEFORMAT);
   bytesToRead = fread(fid,1,'uint32');
   if bytesToRead ~= 80
      fclose(fid);
      MACHINEFORMAT = 'b';
      %disp('Switching to Big Endian');
      fid = fopen(fn,'r',MACHINEFORMAT);
      bytesToRead = fread(fid,1,'uint32');
      if bytesToRead ~= 80

         warning('Waveform:load_seisan:invalidFileFormat',['File does not appear to be a SEISAN file,]'...
            '[ or is an old PC version.\n This program is currently]',...
            '[ incapable of opening PC SEISAN files older than V7.0']);
         fclose(fid);
         return
      end
   end

   fseek(fid,0,-1);

   nsta = 0;
   for i=1:333;
      bytesToRead = fread(fid,1,'uint32');
      if bytesToRead == 80
         st(1:80) = fread(fid,80,'char');
         if i== 1
            h = parseSeisanHeader(st);
            nsta = h.stationCount;
         end

         if i >= 3
            for j=1:3
               thisbit = st(((j-1) * 26 + 1):(j*26));
               ch = parseSeisamChannelHeader(thisbit);
            end
         end
         bytesRead = fread(fid,1,'uint32');

      else
         break
      end
   end

   for j = 1: nsta
      dt = fread(fid,1040,'char');
      mydata = parseSeisamData(dt);
      w(j) = waveform;
      w(j) = set(w(j),'station',mydata.station,'channel',mydata.channel,...
         'start',getStartTime(mydata),'freq',mydata.sampleRate);

      nbytesread = fread(fid,1,'uint32'); %bytecount
      bytestoread = fread(fid,1,'uint32'); %bytecount

      rawdata = fread(fid,bytestoread / 4,'int32');
      w(j) = set(w(j),'data',rawdata);

      nbytesread = fread(fid,1,'uint32'); %bytecount
      nbytesToread2 = fread(fid,1,'uint32'); %bytecount
   end
   fclose(fid);
end

function h = parseSeisanHeader(headerBlock)

   hFieldName = {'network','stationCount','centuryCode','year','doy','month','day','hr','min','sec','totalTime'};
   hFieldpos = {2:30,31:33,34,35:36,38:40,42:43,45:46,48:49,51:52,54:59,61:69};
   hFieldType = ['sddddddddff']; %s=string, d=integer, f=float

   for n=1:numel(hFieldName)
      switch hFieldType(n)
         case 's'
            h.(hFieldName{n}) = deblank(char(headerBlock([hFieldpos{n}])));
         case 'd'
            h.(hFieldName{n}) = str2double(char(headerBlock([hFieldpos{n}]')));
         case 'f'
            h.(hFieldName{n}) = str2double(char(headerBlock([hFieldpos{n}]')));
      end
   end
end
function ch = parseSeisamChannelHeader(channelBlock)

   chFieldName= {'station','comp','lastStaLetter','startRelToEvent','staDataInterval'};
   chFieldpos = {2:5,6:9,10,11:17,19:26};
   chFieldType = ['sssff']; %s=string, d=integer, f=float

   for n=1:numel(chFieldName)
      switch chFieldType(n)
         case 's'
            ch.(chFieldName{n}) = deblank(char(channelBlock([chFieldpos{n}])));
         case 'd'
            ch.(chFieldName{n}) = str2double(char(channelBlock([chFieldpos{n}]')));
         case 'f'
            ch.(chFieldName{n}) = str2double(char(channelBlock([chFieldpos{n}]')));
      end
   end
end
function d = parseSeisamData(dataBlock)

   dFieldName= {'station','channel','centuryCode','year','doy',...
      'month','day','hr','min','timingIndicator','second','sampleRate',...
      'nsamp','lat','lon','elev','gainFactor','howManyBits'};
   dFieldpos = {1:5,6:9,10,11:12,14:16,...
      18:19,21:22,24:25,27:28,29,30:35,37:43,...
      45:50,52:59,61:69,71:75,76,77};
   dFieldType = ['ssdddddddsffdffdsd']; %s=string, d=integer, f=float

   for n=1:numel(dFieldName)
      switch dFieldType(n)
         case 's'
            temp = char(dataBlock([dFieldpos{n}]));
            d.(dFieldName{n}) = deblank(temp(:)');
         case 'd'
            d.(dFieldName{n}) = str2double(char(dataBlock([dFieldpos{n}]')));
         case 'f'
            d.(dFieldName{n}) = str2double(char(dataBlock([dFieldpos{n}]')));
      end
   end
end

function st = getStartTime(myData)

   switch myData.centuryCode
      case 0
         baseyear = 1900;
      case 1
         baseyear = 2000;
      otherwise
         baseyear = 0;
   end
   st = datenum(myData.year + baseyear,...
      myData.month,...
      myData.day,...
      myData.hr,...
      myData.min,...
      myData.second);

   return
end

author: geoscience-community-codes

repo: GISMO


text:
function value = tfn ( h, a )

%*****************************************************************************80
%
%% TFN calculates the T function of Owen.
%
%  Discussion:
%
%    Owen's T function is useful for computation of the bivariate normal
%    distribution and the distribution of a skewed normal distribution.
%
%    Although it was originally formulated in terms of the bivariate
%    normal function, the function can be defined more directly as
%
%      T(H,A) = 1 / ( 2 * pi ) *
%        Integral ( 0 <= X <= A ) e^( -H^2 * (1+X^2) / 2) / (1+X^2) dX
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 August 2006
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    D B Owen,
%    Tables for computing the bivariate normal distribution,
%    Annals of Mathematical Statistics,
%    Volume 27, pages 1075-1090, 1956.
%
%    J C Young and C E Minder,
%    Algorithm AS 76,
%    An Algorithm Useful in Calculating Non-Central T and
%    Bivariate Normal Distributions,
%    Applied Statistics,
%    Volume 23, Number 3, 1974, pages 455-457.
%
%  Parameters:
%
%    Input, real H, A, the arguments of the T function.
%
%    Output, real VALUE, the value of the T function.
%
  ngauss = 10;

  two_pi_inverse = 0.1591549430918953;
  tv1 = 1.0E-35;
  tv2 = 15.0;
  tv3 = 15.0;
  tv4 = 1.0E-05;
  weight = [ ...
    0.666713443086881375935688098933E-01, ...
    0.149451349150580593145776339658E+00, ...
    0.219086362515982043995534934228E+00, ...
    0.269266719309996355091226921569E+00, ...
    0.295524224714752870173892994651E+00, ...
    0.295524224714752870173892994651E+00, ...
    0.269266719309996355091226921569E+00, ...
    0.219086362515982043995534934228E+00, ...
    0.149451349150580593145776339658E+00, ...
    0.666713443086881375935688098933E-01];
  xtab = [
   -0.973906528517171720077964012084E+00, ...
   -0.865063366688984510732096688423E+00, ...
   -0.679409568299024406234327365115E+00, ...
   -0.433395394129247190799265943166E+00, ...
   -0.148874338981631210884826001130E+00, ...
    0.148874338981631210884826001130E+00, ...
    0.433395394129247190799265943166E+00, ...
    0.679409568299024406234327365115E+00, ...
    0.865063366688984510732096688423E+00, ...
    0.973906528517171720077964012084E+00 ];
%
%  Test for H near zero.
%
  if ( abs ( h ) < tv1 )
    value = atan ( a ) * two_pi_inverse;
%
%  Test for large values of abs(H).
%
  elseif ( tv2 < abs ( h ) )
    value = 0.0;
%
%  Test for A near zero.
%
  elseif ( abs ( a ) < tv1 )
    value = 0.0;
%
%  Test whether abs(A) is so large that it must be truncated.
%  If so, the truncated value of A is H2.
%
  else

    hs = - 0.5 * h * h;
    h2 = a;
    as = a * a;
%
%  Computation of truncation point by Newton iteration.
%
    if ( tv3 <= log ( 1.0 + as ) - hs * as )

      h1 = 0.5 * a;
      as = 0.25 * as;

      while ( 1 )

        rt = as + 1.0;
        h2 = h1 + ( hs * as + tv3 - log ( rt ) ) ...
          / ( 2.0 * h1 * ( 1.0 / rt - hs ) );
        as = h2 * h2;

        if ( abs ( h2 - h1 ) < tv4 )
          break
        end

        h1 = h2;

      end

    end
%
%  Gaussian quadrature on the interval [0,H2].
%
    rt = 0.0;
    for i = 1 : ngauss
      x = 0.5 * h2 * ( xtab(i) + 1.0 );
      rt = rt + weight(i) * exp ( hs * ( 1.0 + x * x ) ) / ( 1.0 + x * x );
    end

    value = rt * ( 0.5 * h2 ) * two_pi_inverse;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% This script is used to evluation 3D face reconstruction using NME metric,
% The script contains the process of 3D face mesh reconstruction from the
% 3DMM parameters, and the process of Iterative Closest Point (ICP).
clear all;
close all;
addpath(genpath('./visualize'));
addpath(genpath('./3D_ICP-master'));
aa = textread('./keypoints.txt', '%s');
bb = dir('./results/AFLW-2000-3D_grdth/*.mat');
load('Model_Expression.mat');
load('Model_Shape.mat');
load('BaseSample.mat');

std_size = 120;
vertex_mean = reshape(mu_shape, 3, length(mu_shape)/3);
mu_id = mu_shape + mu_exp;

base_ind = keypoints;
base_ind1 = [(3 * base_ind - 2); (3 * base_ind -1); (3 * base_ind)];
base_ind1 = base_ind1(:);
mu_base = mu_id(base_ind1);
w_base = w_shape(base_ind1,:);
w_exp_base = w_exp(base_ind1,:);

pts = zeros(204, 1);
for i = 1:204
    xx = aa(i);
    xx = xx{1};
    pts(i) = str2num(xx(1:end-1));
end
pts = reshape(pts, 3, 68) + 1;
index = pts(3,:)/3;

LL = 40;
nme_list_2DASSL = zeros(2, LL);
nme_list_3DDFA = zeros(2, LL);
smpNum = 50;
for ii = 1:LL
    ii
    imgName = bb(ii).name;
    img = imread(strcat('./results/AFLW-2000-3D_grdth/', strrep(imgName,'mat', 'jpg')));

    [height, width, nChannels] = size(img);
    info = load(strcat('./results/AFLW-2000-3D_grdth/', imgName));
    pt3d_68 = info.pt3d_68;
    Pose_Para = info.Pose_Para;
    pts1 = pt3d_68;

    % Reconstruct 3D point of base point
    project_base_point = pt3d_68;
    bbox = [min(project_base_point(1,:)), min(project_base_point(2,:)), max(project_base_point(1,:)), max(project_base_point(2,:))];

    center = [(bbox(1)+bbox(3))/2, (bbox(2)+bbox(4))/2];
    radius = max(bbox(3)-bbox(1), bbox(4)-bbox(2)) / 2;
    bbox = [center(1) - radius, center(2) - radius, center(1) + radius, center(2) + radius];

    bbox = double(bbox);
    widthb = vertex_mean(1,keypoints(17)) - vertex_mean(1,keypoints(1));
    heightb = vertex_mean(2, keypoints(20)) - vertex_mean(2, keypoints(9));

    mean_x = vertex_mean(1, keypoints(31));
    mean_y = vertex_mean(2, keypoints(31));

    f0 = ((bbox(3) - bbox(1)) / widthb + (bbox(4) - bbox(2)) / heightb) / 2;
    t3d0(1) = (bbox(3) + bbox(1))/2 - f0*mean_x;
    temp = height + 1 - (bbox(2) + bbox(4))/2;
    t3d0(2) = temp - f0 * mean_y;
    t3d0(3) = 0;

    Pose_Para0 = [Pose_Para(1:3), t3d0(1), t3d0(2), t3d0(3), f0];
%    Shape_Para0 = zeros(size(Shape_Para));
%    Exp_Para0 = zeros(size(Exp_Para));
    Shape_Para = info.Shape_Para;
    Tex_Para = info.Tex_Para;
    Pose_Para = info.Pose_Para;

    [phi, gamma, theta, t3d, f] = ParaMap_Pose(Pose_Para);
    R = RotationMatrix(phi, gamma, theta);

    alpha_shape = Shape_Para;
    alpha_tex = Tex_Para;
    alpha_exp = Exp_Para;
    express = w_exp * alpha_exp; express = reshape(express, 3, length(express)/3);
    shape = mu_shape + w_shape * alpha_shape; shape = reshape(shape, 3, length(shape)/3);
    tex = mu_tex + w_tex * alpha_tex; tex = reshape(tex, 3, length(tex)/3);
    vertex = shape + express;
    grdVertex = f * R * vertex + repmat(t3d, 1, size(vertex, 2));
%     grdVertex(1,:) = grdVertex(1,:) - min(grdVertex(1,:));
%     grdVertex(2,:) = grdVertex(2,:) - min(grdVertex(2,:));
    grdVertex(3,:) = grdVertex(3,:) - min(grdVertex(3,:));

%     figure
%     imshow(uint8(img), [])
%     hold on

    minx = min(grdVertex(1,:)); miny = min(grdVertex(2,:));
    maxx = max(grdVertex(1,:)); maxy = max(grdVertex(2,:));
    llen = sqrt((maxx-minx)*(maxy-miny));

    grdVertex(3,:,:) = grdVertex(3,:,:) - min(grdVertex(3,:,:));
    grdVertex(2,:,:) = 1-grdVertex(2,:,:)+450;
%     pcshow(grdVertex')
%     view(2)

    vertex1 = load(strcat('./results/2DASL_results/', imgName));
    vertex1 = vertex1.vertex;
    vertex1(3,:) =  vertex1(3,:) - min(vertex1(3,:)) + 120;
%     vertex1(1,:) =  vertex1(1,:) - min(vertex1(1,:));
%     vertex1(2,:) =  vertex1(2,:) - min(vertex1(2,:));
    vertex1(3,:) =  vertex1(3,:) - min(vertex1(3,:));
%     pcshow(vertex1')
%     view(2)

    grdVertex_ori = grdVertex(:,1:smpNum:53215);
    vertex1 = vertex1(:,1:smpNum:53215);
%     plot_3d_2(grdVertex_ori', vertex1', -90); % ��ʾ����ǰ�����㼯
    [grdVertex1, vertex1] = icp_process_xgtu(grdVertex_ori', vertex1');
%     plot_3d_2(grdVertex1, vertex1, -90);
    dis1 = (grdVertex1' - vertex1');
    dis1 = sqrt(sum(dis1.^2));
    dis1 = mean(dis1)/llen;

    vertex2 = load(strcat('./results/3DDFA_results/', imgName));
    vertex2 = vertex2.vertex;
    vertex2(3,:,:) =  vertex2(3,:,:) - min(vertex2(3,:,:)) + 120;
    vertex2(3,:) =  vertex2(3,:) - min(vertex2(3,:));
%     pcshow(vertex2')
%     view(2)

    vertex2 = vertex2(:,1:smpNum:53215)
%     plot_3d_2(grdVertex_ori', vertex2', -90); % ��ʾ����ǰ�����㼯
    [grdVertex2, vertex2] = icp_process_xgtu(grdVertex_ori', vertex2');
%     plot_3d_2(grdVertex2, vertex2, -90);
    dis2 = (grdVertex2' - vertex2');
    dis2 = sqrt(sum(dis2.^2));
    dis2 = mean(dis2)/llen;

    nme_list_2DASSL(1, ii) = dis1;
    nme_list_3DDFA(1, ii) = dis2;
end

% save('nme_list', 'nme_list');
dis_2DASSL = nme_list_2DASSL(1,:);
dis_3DDFA = nme_list_3DDFA(1,:);

[s_dis_2DASSL, index] = sort(dis_2DASSL);
[s_dis_3DDFA, index] = sort(dis_3DDFA);

s_dis_2DASSL = s_dis_2DASSL(1:LL);
s_dis_3DDFA = s_dis_3DDFA(1:LL);

x_len = 0:1:length(s_dis_2DASSL);
x_len = x_len(1:length(s_dis_2DASSL));
plot(s_dis_2DASSL*100, x_len, 'g', 'linewidth',2);
hold on
plot(s_dis_3DDFA*100, x_len, 'r', 'linewidth',2);

axis([0 5 0 LL])
set(gca,'XLim',[0 5]);%
set(gca,'YLim',[0 LL]);%

grid on
grid minor

ax = gca;
ax.GridColor = [0 .5 .5];
ax.GridLineStyle = '--';
ax.GridAlpha = 0.5;
ax.Layer = 'top';

h = legend('2DASSL', '3DDFA', 'Location','southeast')
%set(h,'Orientation','horizon', 'Fontsize',12)
set(h,'Fontsize',12)

xlabel('NME normalized by bounding box', 'fontsize', 12)
ylabel('Number of images', 'fontsize',12)

author: XgTu


text:
% Fig. 5.14   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
% script to initiate the RLTOOL like in Example 5.7
% click on Compensator Editor tab
% then right click in Dynamics Box
% add pole, set value to -12. That will produce Fig. 5.14 and the locus
%    in Fig. 5.11
% drag the pole from -12 to -9.   That will produce the RL in Fig. 5.13
% drag the pole from -9 to -4.   That will produce the RL in Fig. 5.12
% Play with the pole location more for your amusement
% If you had added the compensator pole and zero through the compensartor editor,
%   you could have dragged either one around to see how it affected the RL

clf
numL=[1 1];
denL=[1 0 0];
sysL=tf(numL,denL);
rltool(sysL)

author: Sable

repo: mcbench-benchmarks


text:
function arx=lpcbwexp(ar,bw)
%LPCBWEXP expand formant bandwidths of LPC filter ARX=(AR,BW)
%minimum bandwidth will be BW*fs where fs is the sampling frequency
%the radius of each pole will be multiplied by R=exp(-BW*pi)
% To set the maximum pole radius to R use BW=-log(R)/PI.

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: lpcbwexp.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p1]=size(ar);
k=exp(-pi*(0:p1-1)*bw);
arx=ar.*k(ones(nf,1),:);

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function rect = pos2rect(obj_center, obj_size, win_size)
%POS2RECT Get rectangle [x,y,w,h] from obj_center [cx,cy] and size [w,h]
% Parameters:
%   obj_center Rectangle center location [cx, cy]
%   obj_size   Rectangle dimensions [w,h]
%   win_size   (optional) If [width, height] are given, the rectangle will
%              stay within the boundaries [1, 1, width, height]
  rect = [round(obj_center - obj_size./2), obj_size];
  if exist('win_size','var')
    if rect(1) < 1
      corr = abs(rect(1)) + 1;
      rect(1) = 1;
      rect(3) = rect(3) - corr;
    end
    if rect(2) < 1
      corr = abs(rect(2)) + 1;
      rect(2) = 1;
      rect(4) = rect(4) - corr;
    end
    if rect(1) + rect(3) > win_size(1)
      rect(3) = win_size(1) - rect(1);
    end
    if rect(2) + rect(4) > win_size(2)
      rect(4) = win_size(2) - rect(2);
    end
  end
end

author: flyers

repo: drone-tracking


text:
function varargout=quadratic_over_affine_expanded(varargin)
%quadratic_over_affine (internal)

switch class(varargin{1})
    case 'double'
        varargout{1} = varargin{1}'*varargin{1}/varargin{2};

    case 'char' % YALMIP send 'graph' when it wants the epigraph or hypograph
        switch varargin{1}
            case 'graph'
                % Description using epigraphs
                t = varargin{2};
                q = varargin{3};
                y = varargin{4};
                varargout{1} = [cone([2*q;t-y],t+y)];
                varargout{2} = CreateBasicOperator('convex','positive','graph');
                varargout{3} = [q;y];

            case {'exact','integer','callback'}

                t = varargin{2};
                q = varargin{3};
                y = varargin{4};
                varargout{1} = [];
                varargout{2} = CreateBasicOperator('positive','callback');
                varargout{3} = [q;y];

            otherwise
                error([upper(mfilename) ' called with weird argument']);
        end
    otherwise
        error([upper(mfilename) ' called with weird argument']);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [ gd,filtertype ] = gsp_jtv_design_can_dual(g,filtertype)
%GSP_JTV_DESIGN_CAN_DUAL This function returns the canonical dual of the time-vertex filterbank g
%   Usage:  [gd,filtertype] = gsp_jtv_design_can_dual( g,filtertype );
%
%   Inputs parameters:
%       g          : cell array of time-vertex filters
%       filtertype : Filter domain (ts,js,ts-array,js-array)
%
%   Ouputs parameters:
%       g          : cell array of canonical dual time-vertex filters
%       filtertype : Filter domain (ts,js)
%
%   This function returns the canonical dual of the time-vertex filterbank g
%

% Author: Francesco Grassi
% Date:   September 2016

Nf = size(g,1);
gd = cell(Nf,1);

for n = 1:Nf
    gd{n} = @(x,t) can_dual(g,filtertype,n,x,t);
end

switch filtertype
    case {'ts','ts-array'}
        filtertype = 'ts';
    case {'js','js-array'}
        filtertype = 'js';
end

end

function sol = can_dual(g,ft,n,x,t)

if ~isvector(x);x=x(:,1);end
if ~isvector(t);t=t(1,:);end

sol = gsp_jtv_evaluate_can_dual( g,ft,x,t,n );

end

repo: gspbox


text:
function b = r8bb_vxm ( n1, n2, ml, mu, a, x )

%*****************************************************************************80
%
%% R8BB_VXM multiplies a vector by a R8BB matrix.
%
%  Discussion:
%
%    The R8BB storage format is for a border banded matrix.  Such a
%    matrix has the logical form:
%
%      A1 | A2
%      ---+---
%      A3 | A4
%
%    with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
%    are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
%    respectively.
%
%    A should be defined as a vector.  The user must then store
%    the entries of the four blocks of the matrix into the vector A.
%    Each block is stored by columns.
%
%    A1, the banded portion of the matrix, is stored in
%    the first (2*ML+MU+1)*N1 entries of A, using standard LINPACK
%    general band format.  The reason for the factor of 2 in front of
%    ML is to allocate space that may be required if pivoting occurs.
%
%    The following formulas should be used to determine how to store
%    the entry corresponding to row I and column J in the original matrix:
%
%    Entries of A1:
%
%      1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
%
%      Store the I, J entry into location
%      (I-J+ML+MU+1)+(J-1)*(2*ML+MU+1).
%
%    Entries of A2:
%
%      1 <= I <= N1, N1+1 <= J <= N1+N2.
%
%      Store the I, J entry into location
%      (2*ML+MU+1)*N1+(J-N1-1)*N1+I.
%
%    Entries of A3:
%
%      N1+1 <= I <= N1+N2, 1 <= J <= N1.
%
%      Store the I, J entry into location
%      (2*ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
%
%    Entries of A4:
%
%      N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
%
%      Store the I, J entry into location
%      (2*ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
%      (same formula used for A3).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N1, N2, the order of the banded and dense blocks
%    N1 and N2 must be nonnegative, and at least one must be positive.
%
%    Input, integer ML, MU, the lower and upper bandwidths.
%    ML and MU must be nonnegative and no greater than N1-1.
%
%    Input, real A((2*ML+MU+1)*N1 + 2*N1*N2 + N2*N2), the R8BB matrix.
%
%    Input, real X(N1+N2), the vector to multiply A.
%
%    Output, real B(N1+N2), the product X times A.
%

%
%  Initialize B.
%
  b(1:n1+n2) = 0.0E+00;
%
%  Multiply by A1.
%
  for j = 1 : n1
    ilo = max ( 1, j - mu - ml );
    ihi = min ( n1, j + ml );
    ij = (j-1) * (2*ml+mu+1) - j + ml + mu + 1;
    for i = ilo : ihi
      b(j) = b(j) + x(i) * a(ij+i);
    end
  end
%
%  Multiply by A2.
%
  for j = n1+1 : n1+n2
    ij = (2*ml+mu+1)*n1+(j-n1-1)*n1;
    for i = 1 : n1
      b(j) = b(j) + x(i) * a(ij+i);
    end
  end
%
%  Multiply by A3 and A4.
%
  for j = 1 : n1+n2
    ij = (2*ml+mu+1)*n1+n1*n2+(j-1)*n2-n1;
    for i = n1+1 : n1+n2
      b(j) = b(j) + x(i) * a(ij+i);
    end
  end

  return
end


text:
function [yt Xt Xbart]=stvoltmat(Y,X,n,T)

% recover yt
yt=reshape(Y',[n,1,T]);

% recover Xbart
Xt=cell(T,1);
Xbart=cell(T,1);
for ii=1:T
Xt{ii,1}=X(ii,:);
Xbart{ii,1}=kron(speye(n),X(ii,:));
end

author: european-central-bank

repo: BEAR-toolbox


text:
function [pE,pC] = spm_L_priors(dipfit,pE,pC)
% prior moments for the lead-field parameters of ERP models
% FORMAT [pE,pC] = spm_L_priors(dipfit)
%
% dipfit    - forward model structure:
%
%    dipfit.type     - 'ECD', 'LFP' or 'IMG'
%    dipfit.symmetry - distance (mm) for symmetry constraints (ECD)
%    dipfit.location - allow changes in source location       (ECD)
%    dipfit.Lpos     - x,y,z source positions (mm)            (ECD)
%    dipfit.Nm       - number of modes                        (IMG)
%    dipfit.Ns       - number of sources
%    dipfit.Nc       - number of channels
%
% pE - prior expectation
% pC - prior covariance
%
% adds spatial parameters
%--------------------------------------------------------------------------
%    pE.Lpos - position                    - ECD
%    pE.L    - orientation                 - ECD
%              coefficients of local modes - Imaging
%              gain of electrodes          - LFP
%    pE.J    - contributing states (length(J) = number of states per source
%
%__________________________________________________________________________
%
% David O, Friston KJ (2003) A neural mass model for MEG/EEG: coupling and
% neuronal dynamics. NeuroImage 20: 1743-1755
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_L_priors.m 7409 2018-08-27 11:39:00Z bernadette $

% defaults
%--------------------------------------------------------------------------
try, model    = dipfit.model;    catch, model    = 'LFP'; end
try, type     = dipfit.type;     catch, type     = 'LFP'; end
try, location = dipfit.location; catch, location = 0;     end
try, pC;                         catch, pC       = [];    end

% number of sources
%--------------------------------------------------------------------------
try
    n = dipfit.Ns;
    m = dipfit.Nc;
catch
    n = dipfit;
    m = n;
end

% location priors (4 mm)
%--------------------------------------------------------------------------
if location, V = 2^2; else, V = 0; end

% parameters for electromagnetic forward model
%==========================================================================
switch type

    case{'ECD'} % mean           and variance
        %------------------------------------------------------------------
        pE.Lpos = dipfit.Lpos;   pC.Lpos = ones(3,n)*V;    % positions
        pE.L    = zeros(3,n);    pC.L    = ones(3,n)*64;   % orientations

        Sc = find(~cellfun(@isempty,dipfit.silent_source)); % silence sources for CSD
        if(Sc); pC.L(:,Sc) = pC.L(:,Sc)*0; end

    case{'IMG'}
        %------------------------------------------------------------------
        m       = dipfit.Nm;                               % number modes
        pE.Lpos = sparse(3,0);   pC.Lpos = sparse(3,0);    % positions
        pE.L    = zeros(m,n);    pC.L    = ones(m,n)*64;   % modes

        Sc = find(~cellfun(@isempty,dipfit.silent_source)); % silence sources for CSD
        if(Sc); pC.L(:,Sc) = pC.L(:,Sc)*0; end

    case{'LFP'}
        %------------------------------------------------------------------
        pE.Lpos = sparse(3,0);   pC.Lpos = sparse(3,0);    % positions
        pE.L    = ones(1,m);     pC.L    = ones(1,m)*64;   % gains

    otherwise
        warndlg('Unknown spatial model')

end

% contributing states (encoded in J)
%==========================================================================
if ischar(model), mod.source = model; model = mod; end
pE.J = {};
pC.J = {};

for i = 1:numel(model)

    switch upper(model(i).source)

        case{'ERP','SEP'}
            %--------------------------------------------------------------
            pE.J{end + 1} = sparse(1,9,1,1,9);               % 9 states
            pC.J{end + 1} = sparse(1,[1 7],1/32,1,9);

        case{'CMC','TFM'}
            %--------------------------------------------------------------
            pE.J{end + 1} = sparse(1,3,1,1,8);               % 8 states
            pC.J{end + 1} = sparse(1,[1 7],1/32,1,8);

        case{'LFP'}
            %--------------------------------------------------------------
            pE.J{end + 1} = sparse(1,9,1,1,13);              % 13 states
            pC.J{end + 1} = sparse(1,[1 7],1/32,1,13);

        case{'NMM'}
            %--------------------------------------------------------------
            pE.J{end + 1} = sparse(1,3,1,1,9);               % 9 states
            pC.J{end + 1} = sparse(1,[1,2],1/32,1,9);

        case{'NMDA'}
            %--------------------------------------------------------------
            pE.J{end + 1} = sparse(1,3,1,1,12);              % 12 states
            pC.J{end + 1} = sparse(1,[1,2],1/32,1,12);

        case{'CMM'}
            %--------------------------------------------------------------
            pE.J{end + 1} = sparse(1,2,1,1,12);              % 12 states
            pC.J{end + 1} = sparse(1,[3,4],1/32,1,12);

        case{'CMM_NMDA'}
            %--------------------------------------------------------------
            pE.J{end + 1} = sparse(1,2,1,1,16);              % 12 states
            pC.J{end + 1} = sparse(1,[3,4],1/32,1,16);

        case{'MFM'}
            %--------------------------------------------------------------
            pE.J{end + 1} = sparse(1,3,1,1,36);              % 36 (9 + 27)
            pC.J{end + 1} = sparse(1,[1,2],1/32,1,36);       % states

        case{'DEM','NFM'}
            %--------------------------------------------------------------
            pE.J{end + 1} = [];                              % null
            pC.J{end + 1} = [];

        case{'BGT'}
            %--------------------------------------------------------------
            %assuming data are from STN
            pE.J{end + 1} = sparse(1,5,1,1,10);               % 10 states
            pC.J{end + 1} = sparse(1,10);

        case{'MMC'}
            %--------------------------------------------------------------
            pE.J{end + 1} = sparse(1,[1 3 7],[.2 .2 .6],1,8);   % 8 states
            pC.J{end + 1} = sparse(1,8);

        case{'NULL'}
            %--------------------------------------------------------------
            nx   = size(pE.A,1)/m;
            pE.J{end + 1} = ones(1,nx);                      % nx states
            pC.J{end + 1} = ones(1,nx);

        otherwise
            warndlg('Unknown neural model')

    end

    % Cardinal sources
    %----------------------------------------------------------------------
    if isfield(model(i),'J')
        if numel(model(i).J)
            pE.J{i} = spm_zeros(pE.J{i});
            pE.J{i}(model(i).J) = 1;
        end
    end

    %  subsidiary (free) sources
    %----------------------------------------------------------------------
    if isfield(model(i),'K')
        if numel(model(i).K)
            pC.J{i} = spm_zeros(pE.J{i});
            pC.J{i}(model(i).K) = 1/32;
        end
    end

end

% replace J with a vector if there is only one sort of model
%--------------------------------------------------------------------------
if numel(pE.J) == 1
    pE.J = pE.J{:};
    pC.J = pC.J{:};
end

author: spm


text:
function [lat, lon, gam, k] = tranmerc_inv(lat0, lon0, x, y, ellipsoid)
%TRANMERC_INV  Inverse transverse Mercator projection
%
%   [LAT, LON] = TRANMERC_INV(LAT0, LON0, X, Y)
%   [LAT, LON, GAM, K] = TRANMERC_INV(LAT0, LON0, X, Y, ELLIPSOID)
%
%   performs the inverse transverse Mercator projection of points (X,Y) to
%   (LAT,LON) using (LAT0,LON0) as the center of projection.  These input
%   arguments can be scalars or arrays of equal size.  The ELLIPSOID vector
%   is of the form [a, e], where a is the equatorial radius in meters, e is
%   the eccentricity.  If ellipsoid is omitted, the WGS84 ellipsoid (more
%   precisely, the value returned by DEFAULTELLIPSOID) is used.  GEODPROJ
%   defines the projection and gives the restrictions on the allowed ranges
%   of the arguments.  The forward projection is given by TRANMERC_FWD.
%
%   GAM and K give metric properties of the projection at (LAT,LON); GAM is
%   the meridian convergence at the point and K is the scale.
%
%   LAT0, LON0, LAT, LON, GAM are in degrees.  The projected coordinates X,
%   Y are in meters (more precisely the units used for the equatorial
%   radius).  K is dimensionless.
%
%   This implementation of the projection is based on the series method
%   described in
%
%     C. F. F. Karney, Transverse Mercator with an accuracy of a few
%     nanometers, J. Geodesy 85(8), 475-485 (Aug. 2011);
%     Addenda: http://geographiclib.sf.net/tm-addenda.html
%
%   This extends the series given by Krueger (1912) to sixth order in the
%   flattening.  This is a substantially better series than that used by
%   the MATLAB mapping toolbox.  In particular the errors in the projection
%   are less than 5 nanometers withing 3900 km of the central meridian (and
%   less than 1 mm within 7600 km of the central meridian).  The mapping
%   can be continued accurately over the poles to the opposite meridian.
%
%   This routine depends on the MATLAB File Exchange package "Geodesics on
%   an ellipsoid of revolution":
%
%     http://www.mathworks.com/matlabcentral/fileexchange/39108
%
%   See also GEODPROJ, TRANMERC_FWD, GEODRECKON, DEFAULTELLIPSOID.

% Copyright (c) Charles Karney (2012) <charles@karney.com>.
%
% This file was distributed with GeographicLib 1.29.

  if nargin < 4, error('Too few input arguments'), end
  if nargin < 5, ellipsoid = defaultellipsoid; end
  try
    Z = lat0 + lon0 + x + y;
    Z = zeros(size(Z));
  catch err
    error('lat0, lon0, x, y have incompatible sizes')
  end
  if length(ellipsoid(:)) ~= 2
    error('ellipsoid must be a vector of size 2')
  end

  degree = pi/180;
  maxpow = 6;

  a = ellipsoid(1);
  f = ecc2flat(ellipsoid(2));
  e2 = f * (2 - f);
  e2m = 1 - e2;
  cc = sqrt(e2m) * exp(e2 * atanhee(1, e2));
  n = f / (2 -f);
  bet = betf(n);
  b1 = (1 - f) * (A1m1f(n) + 1);
  a1 = b1 * a;

  if isscalar(lat0) && lat0 == 0
    y0 = 0;
  else
    [sbet0, cbet0] = SinCosNorm((1-f) * sind(lat0), cosd(lat0));
    y0 = a1 * (atan2(sbet0, cbet0) + ...
               SinCosSeries(true, sbet0, cbet0, C1f(n)));
  end
  y = y + y0;

  xi = y / a1;
  eta = x / a1;
  xisign = 1 - 2 * (xi < 0 );
  etasign = 1 - 2 * (eta < 0 );
  xi = xi .* xisign;
  eta = eta .* etasign;
  backside = xi > pi/2;
  xi(backside) = pi - xi(backside);

  c0 = cos(2 * xi); ch0 = cosh(2 * eta);
  s0 = sin(2 * xi); sh0 = sinh(2 * eta);
  ar = 2 * c0 .* ch0; ai = -2 * s0 .* sh0;
  j = maxpow;
  xip0 = Z; yr0 = Z;
  if mod(j, 2)
    xip0 = xip0 + bet(j);
    yr0 = yr0 - 2 * maxpow * bet(j);
    j = j - 1;
  end
  xip1 = Z; etap0 = Z; etap1 = Z;
  yi0 = Z; yr1 = Z; yi1 = Z;
  for j = j : -2 : 1
    xip1  = ar .* xip0 - ai .* etap0 - xip1 - bet(j);
    etap1 = ai .* xip0 + ar .* etap0 - etap1;
    yr1 = ar .* yr0 - ai .* yi0 - yr1 - 2 * j * bet(j);
    yi1 = ai .* yr0 + ar .* yi0 - yi1;
    xip0  = ar .* xip1 - ai .* etap1 - xip0 - bet(j-1);
    etap0 = ai .* xip1 + ar .* etap1 - etap0;
    yr0 = ar .* yr1 - ai .* yi1 - yr0 - 2 * (j-1) * bet(j-1);
    yi0 = ai .* yr1 + ar .* yi1 - yi0;
  end
  ar = ar/2; ai = ai/2;
  yr1 = 1 - yr1 + ar .* yr0 - ai .* yi0;
  yi1 =   - yi1 + ai .* yr0 + ar .* yi0;
  ar = s0 .* ch0; ai = c0 .* sh0;
  xip  = xi  + ar .* xip0 - ai .* etap0;
  etap = eta + ai .* xip0 + ar .* etap0;
  gam = atan2(yi1, yr1);
  k = b1 ./ hypot(yr1, yi1);
  s = sinh(etap);
  c = max(0, cos(xip));
  r = hypot(s, c);
  lam = atan2(s, c);
  taup = sin(xip)./r;
  tau = tauf(taup, e2);
  phi = atan(tau);
  gam = gam + atan(tan(xip) .* tanh(etap));
  c = r ~= 0;
  k(c) = k(c) .* sqrt(e2m + e2 * cos(phi(c)).^2) .* ...
         hypot(1, tau(c)) .* r(c);
  c = ~c;
  if any(c)
    phi(c) = pi/2;
    lam(c) = 0;
    k(c) = k(c) * cc;
  end
  lat = phi / degree .* xisign;
  lon = lam / degree;
  lon(backside) = 180 - lon(backside);
  lon = lon .* etasign;
  lon = AngNormalize(lon + AngNormalize(lon0));
  gam = gam/degree;
  gam(backside) = 180 - gam(backside);
  gam = gam .* xisign .* etasign;
end

function bet = betf(n)
  bet = zeros(1,6);
  nx = n^2;
  bet(1) = n*(n*(n*(n*(n*(384796*n-382725)-6720)+932400)-1612800)+ ...
              1209600)/2419200;
  bet(2) = nx*(n*(n*((1695744-1118711*n)*n-1174656)+258048)+80640)/ ...
           3870720;
  nx = nx * n;
  bet(3) = nx*(n*(n*(22276*n-16929)-15984)+12852)/362880;
  nx = nx * n;
  bet(4) = nx*((-830251*n-158400)*n+197865)/7257600;
  nx = nx * n;
  bet(5) = (453717-435388*n)*nx/15966720;
  nx = nx * n;
  bet(6) = 20648693*nx/638668800;
end

function tau = tauf(taup, e2)
  overflow = 1/eps^2;
  tol = 0.1 * sqrt(eps);
  numit = 5;
  e2m = 1 - e2;
  tau = taup / e2m;
  stol = tol * max(1, abs(taup));
  g = ~(abs(taup) < overflow);
  tau(g) = taup(g);
  g = ~g;
  for i = 1 : numit
    if ~any(g), break, end
    tau1 = hypot(1, tau);
    sig = sinh(e2 * atanhee( tau ./ tau1, e2 ) );
    taupa = hypot(1, sig) .* tau - sig .* tau1;
    dtau = (taup - taupa) .* (1 + e2m .* tau.^2) ./ ...
           (e2m * tau1 .* hypot(1, taupa));
    tau(g) = tau(g) + dtau(g);
    g = g & abs(dtau) >= stol;
  end
end


text:
% Makefile
% run
% >> make
% in MATLAB for compile the c++ files.

mex icpCpp.cpp
mex kdtree.cc


text:
function fj=CF_SVj(xt,vt,tau,mu,a,uj,bj,rho,sig,phi)
%--------------------------------------------------------------------------
%PURPOSE: implements the CF fj of Heston. Uses Heston's notations.
%--------------------------------------------------------------------------

xj = bj - rho.*sig.*phi.*i;
dj = sqrt( xj.^2 - (sig.^2).*( 2.*uj.*phi.*i - phi.^2 ) );
gj = ( xj+dj )./( xj-dj );
D  = ( xj+dj )./(sig.^2).* ( 1-exp(dj.*tau) )./( 1-gj.*exp(dj.*tau)  ) ;
xx = ( 1-gj.*exp(dj.*tau) )./( 1-gj );
C  = mu.*phi.*i.*tau + a./( sig.^2 ) .* ( (xj+dj) .* tau - 2.*log(xx) );
fj = exp( C + D.*vt + i.*phi.*xt );


text:
function test_HyperTangents()

xVec = 0:0.001:0.5;
sF=1;
coeff=40;

for i=1:length(xVec)

    x = xVec(i);

    if x < 0.1

        uX_Tar(i) = sF*( tanh(coeff*(0.05-x) ) );

    elseif x <= 0.2

        uX_Tar(i) = sF*( tanh(coeff*(x-0.15) ) );

    elseif x < 0.3

        uX_Tar(i) = sF*( tanh(coeff*(0.25-x) ) );

    elseif x<0.40

        uX_Tar(i) = sF*( tanh(coeff*(x-0.35) ) );

    else
        uX_Tar(i) = sF*( tanh(coeff*(0.45-x) ) );
    end

end

plot(xVec,uX_Tar,'.'); hold on;

author: nickabattista


text:
function im = imread_list(pathstr, file, TOGRAY)
% IMREAD_LIST  Read list of image files into volume.
%
% IMREAD_LIST reads a list of image files, converts them to grayscale if
% desired, and collates them to form a 3D or 4D image volume.
%
% IM = IMREAD_LIST(PATHSTR, FILE)
%
%   PATHSTR is a string with the path to the files. All images must be the
%   same size and have the same number of channels.
%
%   FILE is an array of structs obtained with dir().
%
%   IM is the resulting image volume with indices:
%
%     IM(rows, columns, slice, channel)
%
% IM = IMREAD_LIST(..., TOGRAY)
%
%   TOGRAY is a boolean to convert colour images to grayscale. By default,
%   TOGRAY='false'.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2014 University of Oxford
% Version: 0.1.1
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see
% <http://www.gnu.org/licenses/>.

% check arguments
narginchk(2, 3);
nargoutchk(0, 1);

% defaults
if (nargin < 3 || isempty(TOGRAY))
    TOGRAY = false;
end

% number of files
N = length(file);

% no files to read
if (N == 0)
    im = [];
    return
end

% get number of rows and columns in the blockface volume, and whether
% images are colour. The easiest way is to just load the first image,
% instead of using imfinfo and then having to convert e.g. BitDepth to a
% Matlab class
im0 = imread([pathstr filesep file(1).name]);

% number of channels
nchan = size(im0, 3);

% initialize output volume
if (nchan == 1 || TOGRAY)
    im = zeros(size(im0, 1), size(im0, 2), N, class(im0));
else
    im = zeros(size(im0, 1), size(im0, 2), N, nchan, class(im0));
end
im(:, :, 1, :) = im0;

% loop image files
for I = 2:N

    % load image from file
    if (nchan == 1 || TOGRAY)
        im(:, :, I) = rgb2gray(imread([pathstr filesep file(I).name]));
    else
        im(:, :, I, :) = imread([pathstr filesep file(I).name]);
    end

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function data = parseDataInputs(data, pref)
%PARSEDATAINPUTS   Parse inputs from the DATA structure and assign defaults.
%   DATA = PARSEDATAINPUTS(DATA, PREF) parses inputs and options passed to the
%   constructor via its DATA argument and returns a new DATA structure with the
%   results, including defaults for any missing values.  The preferences in
%   PREF are used as necessary to determine the results.
%
%   The DATA argument to the constructor is a structure that holds all
%   additional parameters required for construction aside from the OP (the
%   function handle that is sampled or numeric data that is used to build the
%   TRIGTECH) and the constructor preferences.

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

if ( ~isfield(data, 'vscale') || isempty(data.vscale) )
    data.vscale = 0;
end

if ( ~isfield(data, 'hscale') || isempty(data.hscale) )
    data.hscale = 1;
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function C = OCP_C(u,x,p,parIdx) %#codegen
    C = OCP_GEN_C(u,x,p,parIdx);
end

author: deng-haoyang

repo: ParNMPC


text:
classdef Controller < handle
%% CONTROLLER is a template for all user-defined controllers.
%
% controller Properties:
%   type            - Name of the controller
%   state_estimate  - Current estimate of the robot's state
%   prev_ticks      - Previous tick count on the left and right wheels
%
% controller Methods:
%   execute         - Executes the control loop (one iteration).
%   update_odometry - Approximates the location of the robot.
%   uni_to_diff     - Converts from unicycle to differential drive model.
%   deg2rad         - Converts from degrees to radians.

    properties
        %% PROPERTIES Controller memory banks

        type                % Name of the controller
    end % properties

    methods
        %% METHODS

        function obj = Controller(type)
        %% CONTROLLER Constructor
        %   obj = controller() is the default constructor that sets the
        %   name of the controller to 'basic', the state estimate to 0, and
        %   the previous tick count to 0.

            % initialize memory banks
            obj.type = type;
        end
    end

    methods

        function outputs = execute(obj, robot, state_estimate, inputs)
        %% EXECUTE Executes the control loop (one iteration).
        %   [...] = obj.execute(robot, ...) is called by the supervisor once every
        %   iteration. This function should compute the control for one
        %   time step only.
        %
        %   This function accepts a variable number of input and output
        %   arguments that allow the user to define any set of inputs and
        %   outputs for the control loop.
            outputs = [];
        end

        function reset(obj)

        end

    end % methods

end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
function [F,xnames,CRLB, Fall]=SimCRLB(obj,Prot,xvalues,sigma,vars)
% Protocol design for qMR: Optimize the stability of fitting parameters
% toward gaussian noise.
% Use the Cramer-Rao Lower bound for objective function: Wikipedia
% [F,xnames,CRLB]=SimCRLB(obj,Prot,xvalues,sigma)
% https://en.wikipedia.org/wiki/Cramer-Rao_bound
% Based on: Alexander, D.C., 2008. A general framework for experiment design in diffusion MRI and its application in measuring direct tissue-microstructure features. Magn. Reson. Med. 60, 439?448.
%
% [F,xnames,CRLB]=SimCRLB(obj,Prot,xvalues,sigma)
% Outputs:
%   F           minimum COV per variable

if 0%~isprop(obj,'fx'), variables = 1:length(obj.xnames);
else
    variables=find(~obj.fx);
end
F=zeros(size(xvalues,1),length(variables));

for ix=1:size(xvalues,1)

    CRLB = (SimFisherMatrix(obj,Prot,xvalues(ix,:),variables,sigma)+eps)^(-1);
    F(ix,:) = diag(CRLB)'./xvalues(ix).^2;

end
Fall = F(:);
F = mean(F(:));
xnames=obj.xnames(variables);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab