text:
function [x,f,funEvals] = minConF_BC(funObj,x,LB,UB,options)
% function [x,f] = minConF_BC(funObj,x,LB,UB,options)
%
% Function for using Two-Metric Projection to solve problems of the form:
%   min funObj(x)
%   s.t. LB_i <= x_i <= UB_i
%
%   @funObj(x): function to minimize (returns gradient as second argument)
%
%   options:
%       verbose: level of verbosity (0: no output, 1: final, 2: iter (default), 3:
%       debug)
%       optTol: tolerance used to check for progress (default: 1e-7)
%       maxIter: maximum number of calls to funObj (default: 250)
%       numDiff: compute derivatives numerically (0: use user-supplied
%       derivatives (default), 1: use finite differences, 2: use complex
%       differentials)
%       method: 'sd', 'lbfgs', 'newton'

nVars = length(x);

% Set Parameters
if nargin < 5
    options = [];
end
[verbose,numDiff,optTol,maxIter,suffDec,interp,method,corrections,damped] = ...
    myProcessOptions(...
    options,'verbose',3,'numDiff',0,'optTol',1e-6,'maxIter',500,'suffDec',1e-4,...
    'interp',1,'method','lbfgs','corrections',100,'damped',0);

% Output Log
if verbose >= 3
    fprintf('%10s %10s %15s %15s %15s\n','Iteration','FunEvals','Step Length','Function Val','Opt Cond');
end

% Make objective function (if using numerical derivatives)
funEvalMultiplier = 1;
if numDiff
    if numDiff == 2
        useComplex = 1;
    else
        useComplex = 0;
    end
    funObj = @(x)autoGrad(x,useComplex,funObj);
    funEvalMultiplier = nVars+1-useComplex;
end

% Evaluate Initial Point
x = projectBounds(x,LB,UB);
if strcmp(method,'newton')
    [f,g,H] = funObj(x);
    secondOrder = 1;
else
    [f,g] = funObj(x);
    secondOrder = 0;
end
funEvals = 1;

% Compute Working Set
working = ones(nVars,1);
working((x < LB+optTol*2) & g >= 0) = 0;
working((x > UB-optTol*2) & g <= 0) = 0;
working = find(working);

% Check Optimality
if isempty(working)
    if verbose >= 1
        fprintf('All variables are at their bound and no further progress is possible at initial point\n');
    end
    return;
elseif norm(g(working)) <= optTol
    if verbose >=1
        fprintf('All working variables satisfy optimality condition at initial point\n');
    end
    return;
end

if verbose >= 3
    switch method
        case 'sd'
            fprintf('Steepest Descent\n');
        case 'lbfgs'
            fprintf('L-BFGS\n');
        case 'bfgs'
            fprintf('BFGS\n');
        case 'newton'
            fprintf('Newton\n');
    end
end

i = 1;
while funEvals <= maxIter

    % Compute Step Direction
    d = zeros(nVars,1);
    switch(method)
        case 'sd'
            d(working) = -g(working);
        case 'lbfgs'
            if i == 1
                d(working) = -g(working);
                old_dirs = zeros(nVars,0);
                old_stps = zeros(nVars,0);
                Hdiag = 1;
            else
                if damped
                    [old_dirs,old_stps,Hdiag] = dampedUpdate(g-g_old,x-x_old,corrections,verbose==3,old_dirs,old_stps,Hdiag);
                else
                    [old_dirs,old_stps,Hdiag] = lbfgsUpdate(g-g_old,x-x_old,corrections,verbose==3,old_dirs,old_stps,Hdiag);
                end
                    curvSat = sum(old_dirs(working,:).*old_stps(working,:)) > 1e-10;
               d(working) = lbfgs(-g(working),old_dirs(working,curvSat),old_stps(working,curvSat),Hdiag);
            end
            g_old = g;
            x_old = x;
        case 'bfgs'
            if i == 1
                d(working) = -g(working);
                B = eye(nVars);
            else
                y = g-g_old;
                s = x-x_old;

                ys = y'*s;

                if i == 2
                    if ys > 1e-10
                        B = ((y'*y)/(y'*s))*eye(nVars);
                    end
                end
                if ys > 1e-10
                    B = B + (y*y')/(y'*s) - (B*s*s'*B)/(s'*B*s);
                else
                    if verbose == 2
                        fprintf('Skipping Update\n');
                    end
                end
                d(working) = -B(working,working)\g(working);
            end
            g_old = g;
            x_old = x;

        case 'newton'
            [R,posDef] = chol(H(working,working));

            if posDef == 0
                d(working) = -R\(R'\g(working));
            else
                if verbose == 3
                    fprintf('Adjusting Hessian\n');
                end
                H(working,working) = H(working,working) + eye(length(working)) * max(0,1e-12 - min(real(eig(H(working,working)))));
                d(working) = -H(working,working)\g(working);
            end
        otherwise
            fprintf('Unrecognized Method: %s\n',method);
            break;
    end

    % Check that Progress can be made along the direction
    f_old = f;
    gtd = g'*d;
    if gtd > -optTol
        if verbose >= 2
            fprintf('Directional Derivative below optTol\n');
        end
        break;
    end

    % Select Initial Guess to step length
    if i == 1 && ~secondOrder
        t = min(1,1/sum(abs(g(working))));
    else
        t = 1;
    end

    % Evaluate the Objective and Projected Gradient at the Initial Step
    x_new = projectBounds(x+t*d,LB,UB);
    if secondOrder
        [f_new,g_new,H] = funObj(x_new);
    else
        [f_new,g_new] = funObj(x_new);
    end
    funEvals = funEvals+1;

    % Backtracking Line Search
    lineSearchIters = 1;
    while f_new > f + suffDec*g'*(x_new-x) || ~isLegal(f_new)
        temp = t;
        if interp == 0 || ~isLegal(f_new) || ~isLegal(g_new)
            if verbose == 3
                fprintf('Halving Step Size\n');
            end
            t = .5*t;
        else
            if verbose == 3
                fprintf('Cubic Backtracking\n');
            end
            t = polyinterp([0 f gtd; t f_new g_new'*d]);
        end

        % Adjust if change is too small
        if t < temp*1e-3
            if verbose == 3
                fprintf('Interpolated value too small, Adjusting\n');
            end
            t = temp*1e-3;
        elseif t > temp*0.6
            if verbose == 3
                fprintf('Interpolated value too large, Adjusting\n');
            end
            t = temp*0.6;
        end

        % Check whether step has become too small
        if sum(abs(t*d)) < optTol
            if verbose == 3
                fprintf('Line Search failed\n');
            end
            t = 0;
            f_new = f;
            g_new = g;
            break;
        end

        % Evaluate New Point
        x_new = projectBounds(x+t*d,LB,UB);
        [f_new,g_new] = funObj(x_new);
        funEvals = funEvals+1;
        lineSearchIters = lineSearchIters+1;

    end

    % Take Step
    x = x_new;
    f = f_new;
    g = g_new;

    % Compute Working Set
    working = ones(nVars,1);
    working((x < LB+optTol*2) & g >= 0) = 0;
    working((x > UB-optTol*2) & g <= 0) = 0;
    working = find(working);

    % Output Log
    if verbose >= 2
        fprintf('%10d %10d %15.5e %15.5e %15.5e\n',i,funEvals*funEvalMultiplier,t,f,sum(abs(g(working))));
    end

    % Check Optimality
    if isempty(working)
        if verbose >= 1
            fprintf('All variables are at their bound and no further progress is possible\n');
        end
        break;
    elseif norm(g(working)) <= optTol
        if verbose >=1
            fprintf('All working variables satisfy optimality condition\n');
        end
        break;
    end

    % Check for lack of progress
    if sum(abs(t*d)) < optTol
        if verbose >= 1
            fprintf('Step size below optTol\n');
        end
        break;
    end

    if abs(f-f_old) < optTol
        if verbose >= 1
            fprintf('Function value changing by less than optTol\n');
        end
        break;
    end

    if funEvals*funEvalMultiplier > maxIter
        if verbose >= 1
            fprintf('Function Evaluations exceeds maxIter\n');
        end
        break;
    end

    % If necessary, compute Hessian
    if secondOrder && lineSearchIters > 1
        [f_new,g_new,H] = funObj(x);
    end

    i = i + 1;
end
end

function [x] = projectBounds(x,LB,UB)
x(x < LB) = LB(x < LB);
x(x > UB) = UB(x > UB);
end

author: ziyuw

repo: rembo

save_path: github-repos/MATLAB/ziyuw-rembo


text:
% test for denoising

name = 'lena';
name = 'peppers';
name = 'boat';
name = 'lena';
name = 'polygons_blurred';
name = 'barb';

n = 128*2;
M = load_image(name);
eta = 0.12;
c = size(M)/2;
% c = [120 200]; % lena hat
M0 = rescale( M(c(1)-n/2+1:c(1)+n/2,c(2)-n/2+1:c(2)+n/2 ), eta, 1-eta );

% add some noise
sigma = 0.04 * max( M0(:) );
randn('state',1234);    % to have reproductible results
M = M0 + randn(n)*sigma;

options.sigma = sigma;
options.repres2 = 'daub3';
options.parent = 1;
M1 = perform_blsgsm_denoising(M, options);

pnoisy = psnr(M,M0);
pwav = psnr(M0,M1);

display_image_layout({M0 M M1}, ...
    {'Original' sprintf('Noisy,psnr=%.2f', pnoisy) sprintf('Denoised,psnr=%.2f', pwav)  });

author: gpeyre

repo: matlab-toolboxes


text:
function [posr, trir, texr] = refine(pos, tri, method, varargin)

% REFINE a 3D surface that is described by a triangulation
%
% Use as
%   [pos, tri]          = refine(pos, tri)
%   [pos, tri]          = refine(pos, tri, 'banks')
%   [pos, tri, texture] = refine(pos, tri, 'banks', texture)
%   [pos, tri]          = refine(pos, tri, 'updown', numtri)
%
% If no method is specified, the default is to refine the mesh globally by bisecting
% each edge according to the algorithm described in Banks, 1983.
%
% The Banks method allows the specification of a subset of triangles to be refined
% according to Banks' algorithm. Adjacent triangles will be gracefully dealt with.
%
% The alternative 'updown' method refines the mesh a couple of times
% using Banks' algorithm, followed by a downsampling using the REDUCEPATCH
% function.
%
% If the textures of the vertices are specified, the textures for the new
% vertices are computed
%
% The Banks method is a memory efficient implementation which remembers the
% previously inserted vertices. The refinement algorithm executes in linear
% time with the number of triangles. It is mentioned in
% http://www.cs.rpi.edu/~flaherje/pdf/fea8.pdf, which also contains the original
% reference.

% Copyright (C) 2002-2014, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%

if nargin<3
  method = 'banks';
end

texture = [];
numtri  = [];

if nargin>3
  switch lower(method)
    case 'banks'
      texture = varargin{1};
    case 'updown'
      numtri = varargin{1};
  end % switch
end

switch lower(method)
  case 'banks'
    if ~isempty(texture)
      npnt   = size(pos,1);
      ntri   = size(tri,1);
      ntex   = size(texture,1);

      assert(ntex==npnt, 'invalid size of texture');

      insert = spalloc(3*npnt,3*npnt,3*ntri);
      trir  = zeros(4*ntri,3);      % allocate memory for the new triangles
      posr  = zeros(npnt+3*ntri,3); % allocate memory for the maximum number of new vertices
      texr  = zeros(ntex+3*ntri,2);
      posr(1:npnt,:) = pos;         % insert the original vertices
      texr(1:ntex,:) = texture;
      current = npnt;

      for i=1:ntri

        if ~insert(tri(i,1),tri(i,2))
          current = current + 1;
          posr(current,:) = (pos(tri(i,1),:) + pos(tri(i,2),:))/2;
          texr(current,:) = (texture(tri(i,1),:) + texture(tri(i,2),:))/2;
          insert(tri(i,1),tri(i,2)) = current;
          insert(tri(i,2),tri(i,1)) = current;
          v12 = current;
        else
          v12 = insert(tri(i,1),tri(i,2));
        end

        if ~insert(tri(i,2),tri(i,3))
          current = current + 1;
          posr(current,:) = (pos(tri(i,2),:) + pos(tri(i,3),:))/2;
          texr(current,:) = (texture(tri(i,2),:) + texture(tri(i,3),:))/2;
          insert(tri(i,2),tri(i,3)) = current;
          insert(tri(i,3),tri(i,2)) = current;
          v23 = current;
        else
          v23 = insert(tri(i,2),tri(i,3));
        end

        if ~insert(tri(i,3),tri(i,1))
          current = current + 1;
          posr(current,:) = (pos(tri(i,3),:) + pos(tri(i,1),:))/2;
          texr(current,:) = (texture(tri(i,3),:) + texture(tri(i,1),:))/2;
          insert(tri(i,3),tri(i,1)) = current;
          insert(tri(i,1),tri(i,3)) = current;
          v31 = current;
        else
          v31 = insert(tri(i,3),tri(i,1));
        end

        % add the 4 new triangles with the correct indices
        trir(4*(i-1)+1, :) = [tri(i,1) v12 v31];
        trir(4*(i-1)+2, :) = [tri(i,2) v23 v12];
        trir(4*(i-1)+3, :) = [tri(i,3) v31 v23];
        trir(4*(i-1)+4, :) = [v12 v23 v31];

      end
      posr = posr(1:current, :);
      texr = texr(1:current, :);

    else
      % there is no texture

      npnt   = size(pos,1);
      ntri   = size(tri,1);
      insert = spalloc(3*npnt,3*npnt,3*ntri);

      trir  = zeros(4*ntri,3);      % allocate memory for the new triangles
      posr  = zeros(npnt+3*ntri,3); % allocate memory for the maximum number of new vertices
      posr(1:npnt,:) = pos;         % insert the original vertices
      current = npnt;

      for i=1:ntri

        if ~insert(tri(i,1),tri(i,2))
          current = current + 1;
          posr(current,:) = (pos(tri(i,1),:) + pos(tri(i,2),:))/2;
          insert(tri(i,1),tri(i,2)) = current;
          insert(tri(i,2),tri(i,1)) = current;
          v12 = current;
        else
          v12 = insert(tri(i,1),tri(i,2));
        end

        if ~insert(tri(i,2),tri(i,3))
          current = current + 1;
          posr(current,:) = (pos(tri(i,2),:) + pos(tri(i,3),:))/2;
          insert(tri(i,2),tri(i,3)) = current;
          insert(tri(i,3),tri(i,2)) = current;
          v23 = current;
        else
          v23 = insert(tri(i,2),tri(i,3));
        end

        if ~insert(tri(i,3),tri(i,1))
          current = current + 1;
          posr(current,:) = (pos(tri(i,3),:) + pos(tri(i,1),:))/2;
          insert(tri(i,3),tri(i,1)) = current;
          insert(tri(i,1),tri(i,3)) = current;
          v31 = current;
        else
          v31 = insert(tri(i,3),tri(i,1));
        end

        % add the 4 new triangles with the correct indices
        trir(4*(i-1)+1, :) = [tri(i,1) v12 v31];
        trir(4*(i-1)+2, :) = [tri(i,2) v23 v12];
        trir(4*(i-1)+3, :) = [tri(i,3) v31 v23];
        trir(4*(i-1)+4, :) = [v12 v23 v31];
      end
      % remove the space for the vertices that was not used
      posr = posr(1:current, :);
    end

  case 'updown'
    ntri = size(tri,1);
    while ntri<numtri
      % increase the number of triangles by a factor of 4
      [pos, tri] = refine(pos, tri, 'banks');
      ntri = size(tri,1);
    end
    % reduce number of triangles using MATLAB function
    [trir, posr] = reducepatch(tri, pos, numtri);

  otherwise
    ft_error('unsupported method "%s"', method);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ score ] = bliinds2_score( image )
%BLIINDS2_SCORE Summary of this function goes here
%   Detailed explanation goes here

    import bliinds2.*

    features = bliinds2_feature_extraction(image);
    score = bliinds_prediction(features);

end

author: dsoellinger


text:
%CONVERTTO  Converts an array to another data type with optional scaling
%
%     dst = cv.convertTo(src)
%     dst = cv.convertTo(src, 'OptionName', optionValue, ...)
%
% ## Input
% * __src__ Input matrix.
%
% ## Output
% * __dst__ Output image of the same size as `src`, and the specified type.
%
% ## Options
% * __RType__ desired output matrix type (`uint8`, `int32`, `double`, etc.),
%   or rather the depth since the number of channels are the same as the input
%   has; if it is negative, the output matrix will have the same type as the
%   input. Default -1
% * __Alpha__ optional scale factor. default 1.0
% * __Beta__ optional delta added to the scaled values. default 0.0
%
% The method converts source pixel values to the target data type.
% Saturation is applied at the end to avoid possible overflows:
%
%     dst = cast(src*alpha + beta, RType);
%
% See also: cv.copyTo, cv.normalize, cast, typecast, im2double, im2uint8,
%  mat2gray
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [model] = bp_finetune(model)
% Discriminative finetuning of CDBN.
% Input: generatively pretrained model.
% Output: discriminative CNN model.

rng('shuffle');
kernels;
debug = 0;

data_list = read_data_list(model.data_path, model.classnames, ....
    model.volume_size + 2 * model.pad_size, 'train', debug);

param = [];
param.epochs = 100;
param.lr = 0.1;
param.weight_decay = 5*10^-4;
param.momentum = 0.9;
param.batch_size = 32;
param.snapshot_iter = 10;
param.snapshot_name = 'bp_finetune_iter';
param.test_iter = 5;
batch_size = param.batch_size;

fprintf('Begin discriminative funetuning the CDBN\n');
fprintf('lr = %f, wd = %d, momentum = %f\n', param.lr, param.weight_decay, param.momentum);

% prepare data and label
[new_list, label] = balance_data(data_list, batch_size);
n = length(new_list);
batch_num = n / batch_size;
assert(batch_num == floor(batch_num));

% prepare model: replace the topmost layer
numLayer = model.numLayer;
numClass = model.classes;
model.layers{numLayer}.layerSize = [numClass, 1];
model.layers{numLayer}.w = rand([prod(model.layers{numLayer-1}.layerSize), prod(model.layers{numLayer}.layerSize)], 'single');
model.layers{numLayer}.w = (model.layers{numLayer}.w - 0.5) * 2 * sqrt( 6 / (prod(model.layers{numLayer}.layerSize) + prod(model.layers{numLayer-1}.layerSize)));
model.layers{numLayer}.c = zeros([1, model.layers{numLayer}.layerSize], 'single');
for l = 2 : numLayer
    model.layers{l} = rmfield(model.layers{l},'b');
    model.layers{l}.grdw = zeros(size(model.layers{l}.w), 'single');
    model.layers{l}.grdc = zeros(size(model.layers{l}.c), 'single');
    model.layers{l}.histw = zeros(size(model.layers{l}.w), 'single');
    model.layers{l}.histc = zeros(size(model.layers{l}.c), 'single');
end

% start training
for iter = 1 : param.epochs
    loss_all = 0;
    shuffle_index = randperm(n);
    for b = 1 : batch_num
        batch_index = shuffle_index((b-1)*batch_size + 1 : b * batch_size);
        batch = read_batch(model, new_list(batch_index), false);
        batch_label = label(batch_index,:);
        [model, activation] = bp_forward(model, batch);
        [model, loss] = bp_backward(model, activation, batch_label);
        loss_all = loss_all + loss;
        model = bp_update(model, param);
    end
    loss_all = loss_all / batch_num;
    fprintf('iteration: %d, loss: %f\n', iter, loss_all);

    if mod(iter, param.snapshot_iter) == 0
        fprintf('snapshoting to %s_%d\n', param.snapshot_name, iter);
        snapshot_name = sprintf('%s_%d', param.snapshot_name, iter);
        save(snapshot_name, 'model');
    end

    if mod(iter, param.test_iter) == 0
        test_loss = bp_test(model);
        fprintf('test loss: %f\n', test_loss);
    end
end

for l = 2 : numLayer
    model.layers{l} = rmfield(model.layers{l},'grdw');
    model.layers{l} = rmfield(model.layers{l},'grdc');
    model.layers{l} = rmfield(model.layers{l},'histw');
    model.layers{l} = rmfield(model.layers{l},'histc');
end

save('bp_finetuned_model', 'model');

author: zhirongw


text:
function [ next_joint ] = neighbour_joint_list(  )

next_joint = [2 3 9 10 4 5 8 9 13 13 10 11 14];

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
function [obX,obY,obXM,obYM,lp] = mrTransByButton(numofanats,obXM,obYM,lp,curInplane,transDelta,direction)
%NAME:    [obX,obY,obXM,obYM,lp] = mrTransByButton(numofanats,obXM,obYM,lp,curInplane,direction)
%AUTHOR:  Poirson
%DATE:	  08.04.96
%PURPOSE: One of a set of routines that allows user
%         to set and select a set of oblique planes in saggital slice.
%	  The routines are mrTransInplanes.m, mrRotInplanes.m,
%	  mrClipInplanes.m, mrSelInplane.m, mrSetupInplanes.m
%         mrSpreadInplanes.m
%HISTORY: Started with mrGetOblPlane from G. Boynton 4/6/96
%NOTES:

global sagwin
figure(sagwin)

% One more pair of points for the perpendicular line
nPtPairs = numofanats + 1;

% See if the user has already been working on some inPlanes
if size(obXM,1) == 0
	disp('You must first create a candidate set of Inplanes');
	disp('Choose Set_Up_Inplanes first');
	return
end

switch direction
  case 1
    obXM = obXM - transDelta;
  case 2
    obXM = obXM + transDelta;
  case 3
    obYM = obYM - transDelta;
  case 4
    obYM = obYM + transDelta;
end

if (curInplane ~= 0)
	obX = obXM(curInplane,:);
	obY = obYM(curInplane,:);
else
	obX = [0,0];
	obY = [0,0];
end

return

for i=1:length(lp)
  delete(lp(i));
end

for i=1:nPtPairs
  if i == curInplane
    lp(i)=line(obXM(i,:),obYM(i,:),'Color','r');
    % Set obX and obY to new values.
  else
    lp(i)=line(obXM(i,:),obYM(i,:),'Color','b');
  end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ n_data, n, a, x, fx ] = gegenbauer_poly_values ( n_data )

%*****************************************************************************80
%
%% GEGENBAUER_POLY_VALUES returns some values of the Gegenbauer polynomials.
%
%  Discussion:
%
%    The Gegenbauer polynomials are also known as the "spherical
%    polynomials" or "ultraspherical polynomials".
%
%    In Mathematica, the function can be evaluated by:
%
%      GegenbauerC[n,m,x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, the order parameter of the function.
%
%    Output, real A, the real parameter of the function.
%
%    Output, real X, the argument of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 38;

  a_vec = [ ...
      0.5E+00, ...
      0.5E+00, ...
      0.5E+00, ...
      0.5E+00, ...
      0.5E+00, ...
      0.5E+00, ...
      0.5E+00, ...
      0.5E+00, ...
      0.5E+00, ...
      0.5E+00, ...
      0.5E+00, ...
      0.0E+00, ...
      1.0E+00, ...
      2.0E+00, ...
      3.0E+00, ...
      4.0E+00, ...
      5.0E+00, ...
      6.0E+00, ...
      7.0E+00, ...
      8.0E+00, ...
      9.0E+00, ...
     10.0E+00, ...
      3.0E+00, ...
      3.0E+00, ...
      3.0E+00, ...
      3.0E+00, ...
      3.0E+00, ...
      3.0E+00, ...
      3.0E+00, ...
      3.0E+00, ...
      3.0E+00, ...
      3.0E+00, ...
      3.0E+00, ...
      3.0E+00, ...
      3.0E+00, ...
      3.0E+00, ...
      3.0E+00, ...
      3.0E+00 ];

  fx_vec = [ ...
       1.0000000000E+00, ...
       0.2000000000E+00, ...
      -0.4400000000E+00, ...
      -0.2800000000E+00, ...
       0.2320000000E+00, ...
       0.3075200000E+00, ...
      -0.0805760000E+00, ...
      -0.2935168000E+00, ...
      -0.0395648000E+00, ...
       0.2459712000E+00, ...
       0.1290720256E+00, ...
       0.0000000000E+00, ...
      -0.3600000000E+00, ...
      -0.0800000000E+00, ...
       0.8400000000E+00, ...
       2.4000000000E+00, ...
       4.6000000000E+00, ...
       7.4400000000E+00, ...
      10.9200000000E+00, ...
      15.0400000000E+00, ...
      19.8000000000E+00, ...
      25.2000000000E+00, ...
      -9.0000000000E+00, ...
      -0.1612800000E+00, ...
      -6.6729600000E+00, ...
      -8.3750400000E+00, ...
      -5.5267200000E+00, ...
       0.0000000000E+00, ...
       5.5267200000E+00, ...
       8.3750400000E+00, ...
       6.6729600000E+00, ...
       0.1612800000E+00, ...
      -9.0000000000E+00, ...
     -15.4252800000E+00, ...
      -9.6969600000E+00, ...
      22.4409600000E+00, ...
     100.8892800000E+00, ...
     252.0000000000E+00 ];

  n_vec = [ ...
     0,  1,  2, ...
     3,  4,  5, ...
     6,  7,  8, ...
     9, 10,  2, ...
     2,  2,  2, ...
     2,  2,  2, ...
     2,  2,  2, ...
     2,  5,  5, ...
     5,  5,  5, ...
     5,  5,  5, ...
     5,  5,  5, ...
     5,  5,  5, ...
     5,  5 ];

  x_vec = [ ...
      0.20E+00, ...
      0.20E+00, ...
      0.20E+00, ...
      0.20E+00, ...
      0.20E+00, ...
      0.20E+00, ...
      0.20E+00, ...
      0.20E+00, ...
      0.20E+00, ...
      0.20E+00, ...
      0.20E+00, ...
      0.40E+00, ...
      0.40E+00, ...
      0.40E+00, ...
      0.40E+00, ...
      0.40E+00, ...
      0.40E+00, ...
      0.40E+00, ...
      0.40E+00, ...
      0.40E+00, ...
      0.40E+00, ...
      0.40E+00, ...
     -0.50E+00, ...
     -0.40E+00, ...
     -0.30E+00, ...
     -0.20E+00, ...
     -0.10E+00, ...
      0.00E+00, ...
      0.10E+00, ...
      0.20E+00, ...
      0.30E+00, ...
      0.40E+00, ...
      0.50E+00, ...
      0.60E+00, ...
      0.70E+00, ...
      0.80E+00, ...
      0.90E+00, ...
      1.00E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    a = 0.0;
    x = 0.0;
    fx = 0.0;
  else
    n = n_vec(n_data);
    a = a_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [A,F]=pisar(xt,sin_num)

%   PISAR Pisarenko Harmonic Decomposition
%
%   PISAR(XT,SIN_NUM) Subplots the input signal (Xt) and the
%   signal generated by the Pisarenko Harmonic Decomposition
%   algorithm for a sum of SIN_NUM sinusoid waves.
%
%   [A F]=PISAR(XT,SIN_NUM) Generates SIN_NUM sinusoids, where the A
%   column-vector represents the amplitudes and the column-vector F
%   represents the normalized frequencies of those sinusoids.
%
%   Considerations:
%      XT must be a column vector;
%      For best results --  SIN_NUM << length of XT
%
%   Plese send any comments to palafox@ieee.org
%   			       Luis E. Palafox Maestre
%

[N t]=size(xt);
rxx=xcorr(xt,'biased');
rxx=rxx(N:(2*sin_num)+N);

%Frequencies estimation
Rxx=toeplitz(rxx);
ev=eig(Rxx);
[S i]=min(ev);
[V D]=eig(Rxx);
a=V(:,i);
rts=roots(a);
w_est=[];
for i=1:sin_num
   w_est(i)= abs(angle(rts(2*i)));
end
F=(w_est/(2*pi))';

%Amplitudes estimation
mcos=[];
for n=1:sin_num
   vcos=[];
   for i=1:sin_num
      vcos=[vcos cos(n*w_est(i))];
   end
   mcos=[mcos; vcos];
end
rxx=rxx(2:sin_num+1);
rxx=2*rxx;
A=inv(mcos)*rxx;
A=A.^(1/2);

if nargout==0,
  xe=[];
  for n=1:N
    xe(n)=0;
    for i=1:sin_num
      xe(n)= xe(n)+(A(i)*cos(w_est(i)*n));
    end
  end
  f=figure;
  subplot(2,1,1);
  plot(xt);
  title('Input Signal')
  xlabel('n');
  ylabel('x(n)');
  subplot(2,1,2);
  plot(xe);
  title('Estimated Signal');
  xlabel('n');
  ylabel('x(n)');
end

author: Sable

repo: mcbench-benchmarks


text:
%%***************************************************************************
%% blkeig: compute eigenvalue decomposition of a cell array
%%         whose contents are square matrices or the diagonal
%%         of a diagonal matrix.
%%
%% [d,V] = blkeig(blk,X);
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%***************************************************************************

    function [d,V] = blkeig(blk,X);

    spdensity = 0.5;

    if ~iscell(X);
       if strcmp(blk{1},'s');
          blktmp = blk{2};
          if (length(blktmp) == 1);
             if (nargout == 1);
                d = eig(full(X));
             elseif (nargout == 2);
                [V,d] = eig(full(X));
                d = diag(d);
             end
          else
             if (nargout == 2);
                V = sparse(length(X),length(X));
             end
	     d = zeros(sum(blktmp),1);
             xx = mexsvec(blk,X,0);
             blktmp2 = blktmp.*(blktmp+1)/2;
             s2 = [0, cumsum(blktmp2)];
             blksub{1,1} = 's'; blksub{1,2} = 0;
             s = [0, cumsum(blktmp)];
             for i = 1:length(blktmp)
                pos = [s(i)+1 : s(i+1)];
                blksub{2} = blktmp(i);
                Xsub = mexsmat(blksub,xx(s2(i)+1:s2(i+1)),0);
                if (nargout == 1);
                   lam = eig(Xsub);
                elseif (nargout == 2);
                   [evec,lam] = eig(Xsub);
                   lam = diag(lam);
                   V(pos,pos) = sparse(evec);
                end
                d(pos,1) = lam;
             end
          end
          n2 = sum(blktmp.*blktmp);
          if (nargout == 2);
             if (nnz(V) <= spdensity*n2);
                V = sparse(V);
             else
                V = full(V);
             end
          end
       elseif strcmp(blk{1},'l');
          if (nargout == 2);
             V = ones(size(X)); d = X;
          elseif (nargout == 1);
             d = X;
          end
       end
    else
       if (nargout == 2);
          V = cell(size(X));  d = cell(size(X));
          for p = 1:size(blk,1);
	     [d{p},V{p}] = blkeig(blk(p,:),X{p});
          end
       elseif (nargout == 1);
          d = cell(size(X));
          for p = 1:size(blk,1);
	     d{p} = blkeig(blk(p,:),X{p});
          end
       end
    end
%%***************************************************************************

repo: REST


text:
p560

% at zero pose
t = fkine(p560, qz)
t
q = ikine560(p560, t)
q
fkine(p560, q)
ikine560(p560, t, 'r')
ikine560(p560, t, 'rn')

ikine(p560, t)

% at nominal pose
qn
t = fkine(p560, qn)
t
%q = ikine560(p560, t)
q
fkine(p560, q)
ikine(p560, t, [0, 0.7, 3, 0, 0.7, 0])

% along trajectory
[q,qd,qdd] = jtraj(qz, qr, 20)
fkine(p560, q)

t1 = fkine(p560, qz)
t2 = fkine(p560, qr)
traj = ctraj(t1, t2, 5)
ikine(p560, traj)

author: petercorke

repo: robotics-toolbox-matlab


text:
function pneopt = mpopt2pneopt(mpopt, model, alg)
%MPOPT2PNEOPT   Create/modify PNES_MASTER options struct from MPOPT.
%
%   PNEOPT = MPOPT2PNEOPT(MPOPT, MODEL)
%   PNEOPT = MPOPT2PNEOPT(MPOPT, MODEL, ALG)
%
%   Uses a MATPOWER options struct, MPOPT, to create or modify an
%   PNES_MASTER options struct.
%
%   Inputs (default values in parentheses):
%       MPOPT : MATPOWER options struct
%       MODEL ('PNE') : (optional) one of the following model types, required
%               for selection of solver in case ALG is 'DEFAULT' (solver
%               precedence for each model type list in parentheses):
%           'PNE' - nonlinear equation with all continuous variables
%                   (NEWTON, FSOLVE)
%       ALG ('DEFAULT') : (optional) any valid value of OPT.alg for
%               PNES_MASTER.
%
%   Output:
%       PNEOPT : an options struct for use by PNES_MASTER and friends
%
%   See PNES_MASTER, MPOPTION.

%   MP-Opt-Model
%   Copyright (c) 2015-2021, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

%% set default args
if nargin < 3
    alg = 'DEFAULT';
    if nargin < 2
        model = '';
    end
end

%% use pf.alg to set nleqs_opt for corrector
switch mpopt.pf.alg
    case 'FSOLVE'
        nleqs_alg = 'FSOLVE';
    otherwise
        nleqs_alg = 'NEWTON';
end

%% create PNES_MASTER options struct
pneopt = struct( ...
    'alg',              alg, ...
    'verbose',          mpopt.verbose, ...
    'nleqs_opt',        struct( ...
        'alg',              nleqs_alg, ...
        'tol',              mpopt.pf.tol    ), ...
    'stop_at',          mpopt.cpf.stop_at, ...
    'solve_base',       1, ...
    'parameterization', mpopt.cpf.parameterization, ...
    'step',             mpopt.cpf.step, ...
    'step_min',         mpopt.cpf.step_min, ...
    'step_max',         mpopt.cpf.step_max, ...
    'adapt_step',       mpopt.cpf.adapt_step, ...
    'adapt_step_damping',   mpopt.cpf.adapt_step_damping, ...
    'adapt_step_tol',   mpopt.cpf.adapt_step_tol, ...
    'target_lam_tol',   mpopt.cpf.target_lam_tol, ...
    'nose_tol',         mpopt.cpf.nose_tol, ...
    'output_fcn',       @(varargin)cpf_pne_output_fcn(obj, ad, varargin{:}), ...
    'plot',             struct( ...
        'level',            mpopt.cpf.plot.level, ...
        'idx',              mpopt.cpf.plot.bus, ...
        'idx_default',      @()cpf_plot_idx_default(obj, dm, ad), ...
        'yname',            'V', ...
        'yfcn',             @(v_,idx)cpf_plot_yfcn(obj, dm, ad, v_, idx), ...
        'title',            'Voltage at Bus %d', ...
        'title2',           'Voltage at Multiple Buses', ...
        'ylabel',           'Voltage Magnitude', ...
        'legend',           'Bus %d'    ) ...
);

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
%%*********************************************************************
%% sortA: sort columns of At{p} in ascending order according to the
%%        number of nonzero elements.
%%
%% [At,C,b,X0,Z0,permA,permZ] = sortA(blk,At,C,b,X0,Z0);
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*********************************************************************

function [At,C,X0,Z0,permA,permZ] = sortA(blk,At,C,b,X0,Z0)

global spdensity smallblkdim
%%
if isempty(spdensity); spdensity = 0.4; end
if isempty(smallblkdim); smallblkdim = 50; end
%%
numblk = size(blk,1);
m = length(b);
nnzA = zeros(numblk,m);
permA = kron(ones(numblk,1),1:m);
permZ = cell(size(blk,1),1);
%%
for p=1:size(blk,1)
    pblk = blk(p,:);
    n = sum(pblk{2});
    % numblk = length(pblk{2});
    if strcmp(pblk{1},'s') && (max(pblk{2}) > smallblkdim)
        % n2 = sum(pblk{2}.*pblk{2});
        n22 = sum(pblk{2}.*(pblk{2}+1))/2;
        m1 = size(At{p,1},2);
        if (length(pblk{2}) == 1)
            tmp = abs(C{p}) + abs(Z0{p});
            if  (~isempty(At{p,1}))
                tmp = tmp + smat(blk(p,:),abs(At{p,1})*ones(m1,1),1);
            end
            if (nnz(tmp) < spdensity*n22);
                per = symamd(tmp);
                invper = zeros(n,1); invper(per) = 1:n;
                permZ{p} = invper;
                if (~isempty(At{p,1}))
                    isspAt = issparse(At{p,1});
                    for k = 1:m1
                        Ak = smat(pblk,At{p,1}(:,k),1);
                        At{p,1}(:,k) = svec(pblk,Ak(per,per),isspAt);
                    end
                end
                C{p}  = C{p}(per,per);
                Z0{p} = Z0{p}(per,per);
                X0{p} = X0{p}(per,per);
            else
                per = [];
            end
            if (length(pblk) > 2) && (~isempty(per))
                % m2 = length(pblk{3});
                P = spconvert([(1:n)', per', ones(n,1)]);
                At{p,2} = P*At{p,2};
            end
        end
        if ~isempty(At{p,1}) && (mexnnz(At{p,1}) < m*n22/2)
            for k = 1:m1
                Ak = At{p,1}(:,k);
                nnzA(p,k) = length(find(abs(Ak) > eps));
            end
            [dummy,permAp] = sort(nnzA(p,1:m1)); %#ok
            At{p,1}  = At{p,1}(:,permAp);
            permA(p,1:m1) = permAp;
        end
    elseif strcmp(pblk{1},'q') || strcmp(pblk{1},'l') || strcmp(pblk{1},'u');
        if ~issparse(At{p,1});
            At{p,1} = sparse(At{p,1});
        end
    end
end
%%*********************************************************************

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function delta = correct_elevation(delta)
%CORRECT_ELEVATION ensures elevation angle between -pi/2 and pi/2
%
%   Usage: delta = correct_elevation(delta)
%
%   Input parameters:
%       delta     - elevation / rad. Can be a single value or a matrix.
%
%   Output paramteres:
%       delta     - angle between -pi/2 and +pi/2 / rad
%
%   See also: correct_azimuth, get_ir

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 1;
nargmax = 1;
narginchk(nargmin,nargmax);

%% ===== Computation ====================================================
% Ensure -pi <= delta <= pi
delta = correct_azimuth(delta);
% Ensure -pi/2 <= delta <= pi/2
delta(delta<-pi/2) = -delta(delta<-pi/2) - pi;
delta(delta>pi/2) = -delta(delta>pi/2) + pi;

author: sfstoolbox

repo: sfs-matlab


text:
function imageData = whitebalance(imageData)
% WHITEBALANCE forces the average image color to be gray.
% Copyright 2013 The MathWorks, Inc.

% Find the average values for each channel.
avg_rgb = mean(mean(imageData));

% Find the average gray value and compute the scaling
% factor.
factors = max(mean(avg_rgb), 128)./avg_rgb;

% Adjust the image to the new gray value.
imageData(:,:,1) = uint8(imageData(:,:,1)*factors(1));
imageData(:,:,2) = uint8(imageData(:,:,2)*factors(2));
imageData(:,:,3) = uint8(imageData(:,:,3)*factors(3));


text:
function z = proj_cone(z,c)
    z = z + proj_dual_cone(-z,c);
end

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%% SPACEXPLOSION Analyze the SpaceX explosion that occurred on 1 Sep 2016
% All supporting functions are in the infrasoundGT directory (or GISMO)
%

%% setup
clear all
close all
matfile = '/Users/glennthompson/Dropbox/Rockets/analysis/20160901_SpaceXplosion/explosion2.mat';
if exist(matfile,'file')
    load(matfile);
else

    % waveform data parameters
    %ds = datasource('antelope', '/raid/data/rockets/dbspacexplosion');
    ds = datasource('antelope', '/Users/glennthompson/Dropbox/Rockets/db/20160901_explosion');
    snum=datenum(2016,9,1,13,0,0);
    enum = snum + 1/24;
    scnl = scnlobject('BCHH', '*', 'FL');

    % Geographical coordinates
    % %%%%%%%%%%%%% SCAFFOLD - could load these from Antelope or Excel
    lat = [28.574182 28.573894 28.574004 28.574013 28.574013 28.574013];
    lon = [-80.572410 -80.572352 -80.572561 -80.572360  -80.572360 -80.572360];
    source.lat = 28.562106; % SLC40 - SpaceX launch complex
    source.lon = -80.57718;
    % Wind tower data - could read this from Excel instead
    relativeHumidity = 92; % percent from NASA weather tower data
    temperatureF = 80; % 80 Fahrenheit according to weather tower data from NASA
    wind_direction_from = 150; % degrees - this is the direction FROM according to NASA, see Lisa Huddleston email of Oct 10th
    wind_direction = mod(wind_direction_from + 180, 360);
    wind_speed_knots = 10; % knots
    wind_speed = wind_speed_knots * 0.514444; % m/s
    % rmpath(genpath('/raid/apps/src/GISMO'))

    %% compute speed of sound based on temperature & rel. humidity
    temperatureC = fahrenheit2celsius(temperatureF);
    speed_of_sound = computeSpeedOfSound(temperatureC, relativeHumidity);
    disp(sprintf('speed of sound at %.1f Celsius and %f percent relative humidity is %.1f',temperatureC, relativeHumidity, speed_of_sound));

    %% load waveform data
    disp('Loading waveform data...')
    w=waveform(ds,scnl,snum,enum);

    save(matfile);
end

make_figures = true
%addpath('infrasoundGT')
figureOutDirectory = '20160901_results';
mkdir('.',figureOutDirectory);

save(matfile)

%% plot raw waveform data
if make_figures
    figure
    plot_panels(w);
    outfile = sprintf('%s/waveforms_raw.png',figureOutDirectory);
    feval('print', '-dpng', outfile);
    close
end

%% compute predicted travel times for infrasound waves based on GPS coords & wind
%% also add lat, lon, distance and bacaz fields to waveform objects
disp('Predicting travel times based on GPS coordinates and wind vector...')
fprintf('\n_______________________________________________\n');
fprintf('PREDICTED TRAVEL TIME BASED ON:\n');
fprintf('  sound speed    %.1fm/s\n', speed_of_sound);
fprintf('  wind speed     %.1fm/s\n', wind_speed);
fprintf('  wind direction %.1f degrees\n', wind_direction);
fprintf('------\t--------\t-----------\t----------\n');
fprintf('Channel\tDistance\tBackAzimuth\tTravelTime\n');
fprintf('------\t--------\t-----------\t----------\n');
for c=1:length(lat)
    [arclen(c), backaz(c)] = distance(lat(c), lon(c), source.lat, source.lon, 'degrees');
    arclen(c) = deg2km(arclen(c))*1000;
    effective_speed = speed_of_sound + wind_speed * cos(deg2rad( (180+backaz(c)) - wind_direction) );
    predicted_traveltime_seconds(c) = arclen(c)/effective_speed;
    fprintf('%s\t%.1fm\t\t%.1f degrees\t%.3fs\n',get(w(c),'channel'), arclen(c), backaz(c), predicted_traveltime_seconds(c));
    w(c) = addfield(w(c), 'lat', lat(c));
    w(c) = addfield(w(c), 'lon', lon(c));
    w(c) = addfield(w(c), 'distance', arclen(c));
    w(c) = addfield(w(c), 'backaz', backaz(c));

end
fprintf('_______________________________________________\n');
fprintf('Program name: %s\n',mfilename('fullpath'))
save spacexplosion.mat

%% plot array map & compute eastings and northings
if make_figures
    disp('Plotting array map')
    close all
    deg2m = deg2km(1) * 1000;
    cols = 'rwbggg';
    for c=1:length(lat)
        chan = get(w(c),'channel');
        easting(c) = distance(lat(c), lon(c), lat(c), source.lon) * deg2m;
        northing(c) = distance(lat(c), lon(c), source.lat, lon(c)) * deg2m;
        plot(easting(c),northing(c),'o','MarkerFaceColor',cols(c),'MarkerSize',10)
        hold on
        quiver(easting(c),northing(c),-easting(c)/100,-northing(c)/100,0); % /100 just gives arrow length
        text(easting(c)+1,northing(c),chan(1:3));
    end
    grid on
    quiver(440,1325,wind_speed*sin(deg2rad(wind_direction)), wind_speed*cos(deg2rad(wind_direction)) ,0,'k');
    text(440,1325,'wind')
    hold off
    title('Beach House array position relative to SLC40');
    xlabel('metres east');
    ylabel('metres north');
    axis equal;
    outfile = sprintf('%s/arraymap.png',figureOutDirectory);
    feval('print', '-dpng', outfile);
    close
    save spacexplosion.mat
end

%% Load arrivals
disp('Loading arrivals...')
arrivals=Arrival.retrieve('antelope', '/raid/data/rockets/dbspacexplosion');
save spacexplosion.mat

%% Subset out X1 arrivals
disp('Subsetting arrivals...')
arrivals = arrivals.subset('iphase', 'X1');
save spacexplosion.mat

%% Associate events
disp('Associating arrivals into events...')
maxTimeDiff = 1; % seconds
eventOn = false;
eventNumber = 0;
for c=2:numel(arrivals.daynumber)
    if arrivals.daynumber(c-1) + maxTimeDiff/86400 > arrivals.daynumber(c)
        if ~eventOn % start new event
            eventOn = true;
            eventNumber = eventNumber + 1;
            infrasoundEvent(eventNumber).FirstArrivalTime = arrivals.daynumber(c-1);
        else % event already in progress
        end
        infrasoundEvent(eventNumber).LastArrivalTime = arrivals.daynumber(c);
    else
        if eventOn % write out last event
            eventOn = false;
        end
    end
end
numEvents = numel(infrasoundEvent);
save spacexplosion.mat

%% Filter waveform data to center around zero
disp('Filtering waveform data...')
wfilt = detrend(w);
f=filterobject('h',[10],3);
wfilt=filtfilt(f,wfilt);
save spacexplosion.mat

%% Segment event waveforms
pretrigger = 1;
posttrigger = 1;
wevent = segment_event_waveforms(wfilt, infrasoundEvent, pretrigger, posttrigger);
save spacexplosion.mat

%% Plot infrasound events
if make_figures
    plot_events(wevent, 'waveforms_infrasoundEvent', figureOutDirectory);
end

%% correlate
% loop through infrasound channels
% take a 0.3-second snippet starting 0.1s before FirstArrivalTime, till 0.2s
% after it
% correlate this against the whole wevent for each infrasound
% trace
% this should result in a correlation matrix
% from this record the time lag matrix for each infrasound channel against each other
% infrasound channel
disp('CORRELATION ...')
disp('_______________')
%infrasoundEvent = xcorr3C(wevent, infrasoundEvent, make_figures, figureOutDirectory, pretrigger);
infrasoundEvent = xcorr3C(wevent, infrasoundEvent, false, figureOutDirectory, pretrigger);

% for eventNumber=1:numEvents
%     fprintf('- processing event %d of %d\n', eventNumber, numEvents);
%     haystacks = wevent{eventNumber};
%     infrasoundEvent(eventNumber).maxCorr = eye(3);
%     infrasoundEvent(eventNumber).secsDiff = eye(3);
%     precorrtime = 0.1; % NEEDLE seconds of data to add before first arrival
%     postcorrtime = 0.2; % NEEDLE seconds of data to add after first arrival
%     for chanNumber=1:3
%         needle = extract(haystacks(chanNumber), 'time', infrasoundEvent(eventNumber).FirstArrivalTime-precorrtime/86400, infrasoundEvent(eventNumber).FirstArrivalTime+postcorrtime/86400);
%         for haystackNum = 1:3
%             fprintf('  - looking for needle %d in haystack %d\n', chanNumber, haystackNum);
%             haystack = haystacks(haystackNum);
%             [acor,lag] = xcorr(get(needle,'data'),get(haystack,'data'));
%             [m,I] = max(abs(acor));
%             infrasoundEvent(eventNumber).maxCorr(chanNumber,haystackNum) = m;
%             infrasoundEvent(eventNumber).secsDiff(chanNumber,haystackNum) = lag(I)/get(haystack,'freq') + pretrigger - precorrtime;
%             if make_figures
%                 figure;
%                 subplot(3,1,1),plot(haystack);
%                 subplot(3,1,2),plot(needle);
%                 subplot(3,1,3),plot(lag,acor);
%                 outfile = sprintf('figs/xcorr_infrasoundEvent%03d_%d_%d.png',eventNumber,chanNumber,haystackNum);
%                 feval('print', '-dpng', outfile);
%                 close
%             end
%         end
%     end
% end
save spacexplosion.mat

%% Construct a master infrasound event, from all the individual ones
disp('Constructing master infrasound event from individual event statistics')
masterEvent.FirstArrivalTime = infrasoundEvent(1).FirstArrivalTime;
masterEvent.LastArrivalTime = infrasoundEvent(1).LastArrivalTime;

% find the mean xcorr time lag difference for non-identical infrasound components - it should be close to zero if only one event in time window
% e.g. needle 1 and haystack 2 should have same magnitude but opposite sign
% time delay to needle 2 and haystack 1 if there is only one clear N wave
% in the haystacks
disp('- finding events with a mean time lag difference of close to 0 - these are the events we can use')
indexes = find(abs([infrasoundEvent.meanSecsDiff]) < 0.01); % these are events with probably only one event in wevent time window
fprintf('- found %d events we can use\n', numel(indexes));
disp('- event indexes to use');
disp(indexes)

masterEvent.secsDiff = zeros(3,3);
masterEvent.stdSecsDiff = zeros(3,3);
for row=1:3
    for column=1:3
        a = [];
        for eventNumber = indexes
            thisEvent = infrasoundEvent(eventNumber);
            a = [a thisEvent.secsDiff(row, column)];
        end

        % eliminate any events which have a difference from the mean of
        % greater than the standard deviation
        diffa = abs( (a-mean(a))/std(a) );

        % now set the mean and std for the master event
        masterEvent.secsDiff(row, column) = mean(a(diffa<1.0));
        masterEvent.stdSecsDiff(row, column) = std(a(diffa<1.0));
    end
end
disp('  - mean:');
disp(masterEvent.secsDiff)
disp('  - std:')
disp(masterEvent.stdSecsDiff)
disp('  - fractional std:')
disp(masterEvent.stdSecsDiff ./ masterEvent.secsDiff)

%% compute sound speed based on GPS coordinates and master event differential travel times
disp('- Estimating sound speed for each component pair using GPS coordinates')
clear speed
speed = zeros(3,3)*NaN;
for row=1:3
    for column = 1:3
        if row ~= column
            radialDistanceDifference = ( get(w(row),'distance') - get(w(column),'distance') );
            timeDifference = masterEvent.secsDiff(row, column);
            s = radialDistanceDifference / timeDifference;
            disp(sprintf('row %d column %d distance difference %.1f time difference %.4f speed %.1fm/s',row,column,radialDistanceDifference,timeDifference,s));
            speed(row,column) = s;
        end
    end
end
speed
meanspeed = nanmean(nanmean(abs(speed)));
stdspeed =  nanstd(nanstd(abs(speed)));
fprintf('- mean sound speed %.1f, std sound speed %.1f\n', meanspeed, stdspeed);
save spacexplosion.mat

%% beamforming
disp('Beamforming to estimate backazimuth of source from array')
[bestbackaz,bestsoundspeed,distanceDiff,speedMatrix] = beamform2(easting(1:3), northing(1:3), masterEvent.secsDiff, 199.0, 348.6);
distanceDiff
speedMatrix
sourceDist = get(w(2),'distance');
[beamformingsourcelat, beamformingsourcelon] = reckon(lat(2), lon(2), km2deg(sourceDist/1000), bestbackaz);
distFromSLC40 = deg2km(distance(beamformingsourcelat, beamformingsourcelon, source.lat, source.lon)) * 1000;
disp(sprintf('- source location estimated to be lat = %.4f lon = %.4f, which is %1.fm from true position',beamformingsourcelat, beamformingsourcelon,distFromSLC40));

%%
if make_figures
    outfile = sprintf('%s/beamforming.png',figureOutDirectory);
    feval('print', '-dpng', outfile);
    close
end
save spacexplosion.mat

%% estimate travel times from cross-correlation derived differential travel times and actual source location, wind speed, and predicted sound speed
[minimumPredictedTravelTime,index] =  min(predicted_traveltime_seconds);
for component=1:3
    traveltime_secs(component) = minimumPredictedTravelTime + mean([masterEvent.secsDiff(component,index) -masterEvent.secsDiff(index,component)]);
end
traveltime_secs(4) = minimumPredictedTravelTime + (get(w(4),'distance') - get(w(index),'distance')) / bestsoundspeed;
traveltime_secs(5) = traveltime_secs(4);
traveltime_secs(6) = traveltime_secs(4);
save spacexplosion.mat

%% Shift waveforms based on travel times
wshift = wfilt;
for c=1:numel(wshift)
     starttime = get(w(c), 'start');
     newstarttime = starttime - traveltime_secs(c)/86400;
     disp(sprintf('moving channel %d from %s to %s\n', c, datestr(starttime, 'HH:MM:SS.FFF'), datestr(newstarttime, 'HH:MM:SS.FFF')));
     wshift(c)=set(wshift(c),'start', newstarttime);
     disp(sprintf('%s\n',datestr(get(wshift(c),'start'))));
end
save spacexplosion.mat

%% Segment time shifted event waveforms
disp('Segmenting traveltime-corrected event waveforms...')
arrivalTimeCorrection = minimumPredictedTravelTime;
preplot = 0.15;
postplot = 0.15;
weventshift = segment_event_waveforms(wshift, infrasoundEvent, preplot, postplot, arrivalTimeCorrection);
save spacexplosion.mat

%% Plot shifted events
if make_figures
    plot_events(weventshift, 'waveforms_infrasoundEvent_shifted', figureOutDirectory)
end

%% Pick events
close all
for eventNumber=1:numEvents

    % plot waveforms for this event
    w2=weventshift{eventNumber};
    fh=plot_panels(w2);
    ah=get(fh,'Children');
    set(fh, 'Position', [0 0 1600 1000]);
    infrasoundEvent(eventNumber).maxAmp=zeros(1,6);
    infrasoundEvent(eventNumber).minAmp=zeros(1,6);
    infrasoundEvent(eventNumber).maxTime=zeros(1,6);
    infrasoundEvent(eventNumber).minTime=zeros(1,6);

    for chanNum=1:6

        % scan over this event on this channel and find the greatest
        % max/min difference in 0.1s
        y = get(w2(chanNum),'data');
        fs = get(w2(chanNum),'freq');
        numSamples = length(y);
        windowSize = round(fs/25);
        maxA = 0;
        for startSamp = 1:windowSize:round(numSamples*0.7)-windowSize
            samples = startSamp:startSamp+windowSize-1;
            [maxy, maxindex] = max(y(samples));
            [miny, minindex] = min(y(samples));
            if (maxy-miny) > maxA
                maxSecs = ((maxindex+samples(1)-1)/fs);
                minSecs = ((minindex+samples(1)-1)/fs);
                maxA = maxy-miny;

                infrasoundEvent(eventNumber).maxTime(chanNum) = tstart + maxSecs/86400;
                infrasoundEvent(eventNumber).minTime(chanNum) = tstart + minSecs/86400;
                infrasoundEvent(eventNumber).maxAmp(chanNum) = maxy;
                infrasoundEvent(eventNumber).minAmp(chanNum) = miny;
            end
        end

        % plot max & min
        axisnum = 8 - chanNum;
        axes(ah(axisnum));
        hold on
        plot(ah(axisnum),maxSecs, infrasoundEvent(eventNumber).maxAmp(chanNum), 'g*');
        plot(ah(axisnum),minSecs, infrasoundEvent(eventNumber).minAmp(chanNum), 'r*');

    end
    feval('print', '-dpng', sprintf('%s/picked_event_%03d',figureOutDirectory, eventNumber) );
    close
end
save spacexplosion.mat

%% Compute relative calibrations of infrasound sensors
clc
tolerance = 0.02; % times have to be within 2 hundreds of a second
component_pairs = [ [1,3]; [1,2]; [2,3] ];

for component_pair_num = 1:length(component_pairs)
    chanNum = nan(2,1);
    chanNum(1) = component_pairs(component_pair_num, 1);
    chanNum(2) = component_pairs(component_pair_num, 2);
    maxA = nan(2, numEvents);
    minA = nan(2, numEvents);
    for eventNum = 1:numEvents
        ev = infrasoundEvent(eventNum);
        if std(ev.maxTime) < tolerance/86400 & std(ev.minTime) < tolerance/86400
            for c=1:2
                maxA(c, eventNum) = ev.maxAmp(chanNum(c));
                minA(c, eventNum) = ev.minAmp(chanNum(c));
            end
        end
    end
    avA = (maxA - minA)/2;

    % first compute raw ratio, mean and std
    ratio = avA(1,:)./avA(2,:);
    m=nanmean(ratio);
    s=nanstd(ratio);

    % now recompute after throwing out any measurement that has an error larger
    % than the standard deviation
    ratio2 = ratio(ratio>m-s & ratio<m+s);
    m2=nanmean(ratio2);
    s2=nanstd(ratio2);
    distratio = get(w(chanNum(2)),'distance') / get(w(chanNum(1)), 'distance')

    fprintf('relative amplitude of signals on sensors %d and %d: %f (+/- %f)\n',chanNum(1), chanNum(2), m2, s2);
    fprintf('distance ratio = %f\n', distratio );

    relAmp(component_pair_num) = m2/distratio;
    fprintf('relative amplitude of signals corrected for distance is %f\n', m2/distratio);

end

%% compute new calibration constants assume sensor 3 has correct calib
oldcalib = get(w,'calib');
newcalib(1) = oldcalib(1) / relAmp(1);
newcalib(2) = oldcalib(2) / relAmp(3);
newcalib(3) = get(w(3),'calib');
for chanNum=1:3
    fprintf('Using channel 3 as a reference, channel 1 calib old = %f, new = %f\n',oldcalib(chanNum),newcalib(chanNum))
end

%% compute reduced pressures
for eventNum = 1:numEvents
    for chanNum = 1:3
        ev = infrasoundEvent(eventNum);
        pchange = ev.maxAmp(chanNum) - ev.maxAmp(chanNum);
        distanceInKm = get(w(chanNum),'distance') / 1000.0;
        infrasoundEvent(eventNum).pchange(chanNum) = pchange;
        infrasoundEvent(eventNum).preduced(chanNum) = pchange / 2 * distanceInKm;
    end
end
% add origin time - time of earliest arrival corrected for travel time

%%

SCAFFOLD

%% write out events
fout = fopen(sprintf('%s/eventlist.csv',figureOutDirectory),'w');
for eventNum = 1:numEvents
    ev = infrasoundEvent(eventNum);
    % event number
    fprintf(fout, '%d,', eventNum);
    % origin time
    fprintf(fout, '%s,', datestr(ev.originTime));
    % pressure change (from stack)
    fprintf(fout, '%.1f,', nanmean(ev.pchange));
    % reduced pressure (from stack)
    fprintf(fout, '%.1f,', nanmean(ev.preduced));
    % energy (from stack)
    fprintf(fout, '%.1f\n', nanmean(ev.preduced)*some_correction_factor);
    % by component
        % arrival time (not time corrected)
        % maxA
        % minA
end
fclose(fout)

%% Particle motion analysis
thisw = wevent{1};
t = threecomp(thisw([6 5 4])',backaz(1));
tr = t.rotate()
tr2 = tr.particlemotion();
tr2.plotpm()
tr2.plot3()

author: geoscience-community-codes

repo: GISMO


text:
function [e1,e2,E1_l,E2_l,E1_s1,E2_s2] = pluckerEndpoints(L,s1,s2)

% PLUCKERENDPOINTS  Plucker line and abscissas to endpoints conversion.
%   [E1,E2] = PLUCKERENDPOINTS(L,S1,S2) are the endpoints of the Plucker
%   line L at abscissas S1 and S2.
%
%   [E1,E2,E1_l,E2_l,E1_s1,E2_s2] = PLUCKERENDPOINTS(L,S1,S2) returns the
%   Jacobians wrt the line L and the abscissas S1 and S2.
%
%   See also LS2E, LS2SEG.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

v = L(4:6);

if nargout == 1

    vn = normvec(v);
    p0 = pluckerOrigin(L);
    e1 = p0 + s1*vn;
    e2 = p0 + s2*vn;

else % Jac

    [vn,VN_v] = normvec(v,0);
    [p0,P0_l] = pluckerOrigin(L);

    P0_n  = P0_l(:,1:3);
    P0_v  = P0_l(:,4:6);

    e1    = p0 + s1*vn;
    E1_s1 = vn;

    E1_n  = P0_n;
    E1_v  = P0_v + s1*VN_v;

    E1_l  = [E1_n E1_v];

    e2    = p0 + s2*vn;
    E2_s2 = vn;

    E2_n  = P0_n;
    E2_v  = P0_v + s2*VN_v;

    E2_l  = [E2_n E2_v];

end

return

%% jac

syms n1 n2 n3 v1 v2 v3 s1 s2 real
L = [n1;n2;n3;v1;v2;v3];

[e1,e2,E1_l,E2_l,E1_s1,E2_s2] = pluckerEndpoints(L,s1,s2);

simplify(E1_l - jacobian(e1,L))
simplify(E2_l - jacobian(e2,L))
simplify(E1_s1 - jacobian(e1,s1))
simplify(E2_s2 - jacobian(e2,s2))

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

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function SetDefaultValue(position, argName, defaultValue)
%SETDEFAULTVALUE Initialise a missing or empty value in the caller.
%
% SETDEFAULTVALUE(POSITION, ARGNAME, DEFAULTVALUE) checks to see if the
% argument named ARGNAME in position POSITION of the caller function is
% missing or empty, and if so, assigns it the value DEFAULTVALUE.
%
% Example:
% function x = TheCaller(x)
% SetDefaultValue(1, 'x', 10);
% end
% TheCaller()    % 10
% TheCaller([])  % 10
% TheCaller(99)  % 99

% $Author: rcotton $  $Date: 2010/10/01 13:58:35 $ $Revision: 1.1 $
% Copyright: Health and Safety Laboratory 2010

if evalin('caller', 'nargin') < position || ...
      isempty(evalin('caller', argName))
   assignin('caller', argName, defaultValue);
end
end


text:
clear all, close all, clc
n = 1000;
X = zeros(n,n);
X(n/4:3*n/4,n/4:3*n/4) = 1;
[U,S,V] = svd(X);
imagesc(X), hold on;
semilogy(diag(S),'-o','color',cm(1,:)), hold on, grid on

nAngles = 12;  % sweep through 12 angles, from 0:4:44
Xrot = X;
for j=2:nAngles
    Y = imrotate(X,(j-1)*4,'bicubic'); % rotate (j-1)*4
    startind = floor((size(Y,1)-n)/2);
    Xrot1 = Y(startind:startind+n-1, startind:startind+n-1);
    Xrot2 = Xrot1 - Xrot1(1,1);
    Xrot2 = Xrot2/max(Xrot2(:));
    Xrot(Xrot2>.5) = j;

    [U,S,V] = svd(Xrot1);
    subplot(1,2,1), imagesc(Xrot), colormap([0 0 0; cm])
    subplot(1,2,2), semilogy(diag(S),'-o','color',cm(j,:))
end

author: dynamicslab


text:
% zica() - Z-transform of ICA activations; useful for studying component SNR
%
% Usage: >> [zact,basesd,maz,mazc,mazf] = zica(activations,frames,baseframes)
%
% Inputs:
%   activations - activations matrix produced by runica()
%   frames      - frames per epoch {0|default ->  length(activations)}
%   baseframes  - vector of frames in z-defining baseline period {default frames}
%
% Outputs:
%   zact        - activations z-scaled and reorded in reverse order of max abs
%   basesd      - standard deviations in each activation row (reverse ordered)
%   maz         - maximum absolute z-value for each activation row (rev ordered)
%   mazc        - component indices of the reverse-sorted max abs z-values
%                 (this is the act -> zact reordering)
%   mazf        - frame indices of the max abs z-values (reverse ordered)
%
% Author: Scott Makeig, SCCN/INC/UCSD, La Jolla, 2-25-98
%
% See also: runica()

function [zact,basesd,maxabsz,maxc,maxabszf] = zica(activations,frames,baseframes)

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

% 3-2-98 added frames variable -sm
% 1-25-01 put revsort subfunction in the core of the zica program -ad
% 01-25-02 reformated help & license, added link -ad

if nargin<1
   help zica
   return
end

[chans,framestot] = size(activations);

if nargin < 3
   baseframes = 0;
end
if nargin < 2
   frames = 0;
end
if frames == 0
   frames = framestot;
end
if baseframes == 0
   baseframes = 1:frames
end
epochs = floor(framestot/frames);
if frames*epochs ~= framestot
   fprintf('zica(): indicated frames does not divide data length.\n');
   return
end

if length(baseframes) < 3
  fprintf('\n  zica() - too few baseframes (%d).\n',length(baseframes));
  help zica
  return
end

if min(baseframes) < 1 | max(baseframes) > frames
  fprintf('\n  zica() - baseframes out of range.\n');
  help zica
  return
end

baselength = length(baseframes);
baseact = zeros(epochs*baselength,chans);
for e=1:epochs
 baseact((e-1)*baselength+1:e*baselength,:) = ...
                        matsel(activations,frames,baseframes,0,e)';
end
basesd = sqrt(covary(baseact));
zact = activations./(basesd'*ones(1,framestot));
[maxabsz,maxabszf] = sort(abs(zact'));
maxabsz  = maxabsz(frames,:);
maxabszf = maxabszf(frames,:);

%
% reorder outputs in reverse order of max abs z
[maxabsz,maxc] = revsort(maxabsz);
zact = zact(maxc,:);
basesd = basesd(maxc);
maxabszf = maxabszf(maxc);

% revsort  - reverse sort columns (biggest 1st, ...)
function [out,i] = revsort(in)

if size(in,1) == 1
   in = in'; % make column vector
end

[out,i] = sort(in);
out = out(size(in,1):-1:1,:);
  i = i(size(in,1):-1:1,:);
return;

author: PatternRecognition

repo: OpenBMI


text:
function [stats,st] = SStat_CoxExt(sID_ext,X_ext,d_ext,t_ext,e)
% [stats,st] = SStat_CoxExt(sID_ext,X_ext,d_ext,t_ext,e)
%
% Parameter estimation for the extended Cox model. This function uses as
% input the output of SStat_X_ext but with X_ext having new time-dependent
% columns added by the user (eg. the product of a column of X_ext with
% t_ext. Survival time ties are handled using Efron's method.
%
% Input
% sID_ext: Extended subjects' IDs (nx1).
% X_ext: Extended design matrix (nxp).
% d_ext: Extended censorship status vector (nx1).
% t_ext: Extended survival time vector (nx1).
% e: Convergence epsilon (gradient's norm). Default 10^-3;
%
% Output
% stats.Bhat: Estimated vector of the population regression parameters.
% stats.CovBhat: Estimated covariance matrix of the population regression
% parameters.
% stats.llh: Values of the maximum log-likelihood across the optimization
% process.
% st: Termination state (1 for convergence and 0 otherwise).
%
% $Revision: 1.1 $  $Date: 2015/01/06 17:04:02 $
% Original Author: Jorge Luis Bernal Rusiel
% CVS Revision Info:
%    $Author: mreuter $
%    $Date: 2015/01/06 17:04:02 $
%    $Revision: 1.1 $
% References: Kleinbaum, D.G., Klein, M., 2005. Survival analysis. A self-
% learning approach, second edition. New York: Springer..
%
if nargin < 4
    error('Too few inputs');
elseif nargin < 5
    e = 0.001;
end;
tic;
[n,p] = size(X_ext);
if (length(sID_ext)~=n) || (length(d_ext)~=n) || (length(t_ext)~=n)
    error(['The design matrix X_ext, censorship status vector d_ext, time'...
        ' vector t_ext and subject ID vector sID_ext must all have the same'...
        ' number of rows.']);
end;
%indices of unique failure times in ft_ix (last index when ties happen)
st_ix = find(d_ext==1);
[~,ft_ix] = unique(t_ext(st_ix),'last');
ft_ix = st_ix(ft_ix);
%Starting values
Bhat = zeros(p,1);

%% Iterations
nit = 50;
gnorm = e+1;
it = 1;
display('Starting Newton-Raphson iterations');
while (gnorm>e) && (it<=nit)
    gr = SStat_Gradient(X_ext,t_ext,Bhat,ft_ix);
    He = SStat_Hessian(X_ext,t_ext,Bhat,ft_ix);
    if (cond(He) < 1e+10)
        invHe = He\eye(p);
    else
        [Vtemp,Dtemp] = eig(He);
        invHe = Vtemp*diag(1./max(diag(Dtemp),1e-5))*Vtemp';
    end
    Bhat = Bhat - invHe*gr;
    %log-likelihood
    llh = SStat_Likelihood(X_ext,t_ext,Bhat,ft_ix);
    display(['Likelihood at iteration ' num2str(it) ' : ' num2str(llh)]);
    gnorm = norm(gr);
    display(['Gradient norm: ' num2str(gnorm)]);
    it = it+1;
end;
%% Termination
z_sc = Bhat./sqrt(diag(-invHe));
pval = 2*(1-normcdf(abs(z_sc),0,1));
stats = struct('Bhat',Bhat,'zscore',z_sc,'pval',pval,'CovBhat',-invHe,'llh',llh);
if (gnorm<=e)
    st = 1;
else
    st = 0;
    display(['Algorithm does not converge after ' num2str(nit)...
        ' iterations!!!']);
end;
et = toc;
display(['Total elapsed time is ' num2str(et) ' seconds']);
end

%% Likelihood, Gradient and Hessian

function llk = SStat_Likelihood(X_ext,t_ext,Bhat,ft_ix)
%
% Log-likelihood value.
%
% Input
% X_ext: Extended design matrix.
% t_ext: Extended survival time vector.
% Bhat: Estimated vector of the population regression parameters.
% ft_ix: Failure time indices in X_ext (last index if any tie).
%
% Output
% llk: Log-likelihood value.
%
llk = 0;
nft = length(ft_ix);
for j=1:nft
    term = 0;
    ties = ft_ix(t_ext(ft_ix(j))==t_ext(ft_ix));
    nties = length(ties)-1;
    lpr = X_ext(ties,:)*Bhat;
    aux1 = sum(exp(lpr));
    aux2 = sum(exp(X_ext(t_ext(ft_ix(j))==t_ext,:)*Bhat));
    for l=0:nties
        term = term + log(aux2-l*aux1/(nties+1));
    end;
    llk = llk + sum(lpr) - term;
end;
end

function gr = SStat_Gradient(X_ext,t_ext,Bhat,ft_ix)
%
% Gradient vector for the log-likelihood.
%
% Input
% X_ext: Extended design matrix.
% t_ext: Extended survival time vector.
% Bhat: Estimated vector of the population regression parameters.
% ft_ix: Failure time indices in X_ext (last index if any tie).
%
% Output
% gr: Gradient vector.
%
p = size(X_ext,2);
gr = zeros(p,1);
nft = length(ft_ix);
for j=1:nft
    ties = ft_ix(t_ext(ft_ix(j))==t_ext(ft_ix));
    nties = length(ties)-1;
    riskset = t_ext(ft_ix(j))==t_ext;
    term1 = sum(X_ext(ties,:),1);
    term2 = exp(X_ext(riskset,:)*Bhat);
    term3 = exp(X_ext(ties,:)*Bhat);
    term4 = term2'*X_ext(riskset,:);
    term5 = term3'*X_ext(ties,:);
    term = 0;
    for l=0:nties
        term = term + (term4-l*term5/(nties+1))/(sum(term2)-l*sum(term3)/(nties+1));
    end;
    gr = gr + (term1 - term)';
end;
end

function He = SStat_Hessian(X_ext,t_ext,Bhat,ft_ix)
%
% Hessian matrix for the log-likelihood.
%
% Input
% X_ext: Extended design matrix.
% t_ext: Extended survival time vector.
% Bhat: Estimated vector of the population regression parameters.
% ft_ix: Failure time indices in X_ext (last index if any tie).
%
% Output
% He: Hessian matrix.
%
p = size(X_ext,2);
He = zeros(p,p);
nft = length(ft_ix);
for j=1:nft
    ties = ft_ix(t_ext(ft_ix(j))==t_ext(ft_ix));
    nties = length(ties)-1;
    rsk_ix = find(t_ext(ft_ix(j))==t_ext);
    m = length(rsk_ix);
    term1 = 0;
    for i=1:m
        term1 = term1 + exp(X_ext(rsk_ix(i),:)*Bhat)*X_ext(rsk_ix(i),:)'*X_ext(rsk_ix(i),:);
    end;
    term2 = 0;
    for i=1:nties
        term2 = term2 + exp(X_ext(ties(i),:)*Bhat)*X_ext(ties(i),:)'*X_ext(ties(i),:);
    end;
    term3 = sum(exp(X_ext(rsk_ix,:)*Bhat));
    term4 = sum(exp(X_ext(ties,:)*Bhat));
    term5 = exp(X_ext(rsk_ix,:)*Bhat)'*X_ext(rsk_ix,:);
    term6 = exp(X_ext(ties,:)*Bhat)'*X_ext(ties,:);
    term = 0;
    for l=0:nties
        Z = term5 - l*term6/(nties+1);
        phi = term3 - l*term4/(nties+1);
        term = term + (term1-l*term2/(nties+1))/phi - Z'*Z/(phi*phi);
    end;
    He = He - term;
end;
end


text:
function [pf,f]=v_lpcar2pf(ar,np)
%V_LPCAR2PF Convert AR coefs to power spectrum PF=(AR,NP)
%
%  Inputs: ar(nf,n)     AR coefficients, one frame per row
%          np           Size of output spectrum is np+1 [n]
%
% Outputs: pf(nf,np+1)  Power spectrum from DC to Nyquist
%          f(1,np+1)    Normalized frequencies (0 to 0.5)
%
% For high speed make np equal to a power of 2

%      Copyright (C) Mike Brookes 1998-2014
%      Version: $Id: v_lpcar2pf.m 10865 2018-09-21 17:22:45Z 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);
if nargin<2
    if nargout
        np=p1-1;
    else
        np=128;
    end
end
pf=abs(v_rfft(ar.',2*np).').^(-2);
f=(0:np)/(2*np);
if ~nargout
        plot(f,pf.');
        xlabel('Normalized frequency f/f_s');
        ylabel('LPC Power Spectrum');
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function [CM]=fourthOrderMat(C)

ind_C_all=1:numel(C);
[I,J,K,L]=ind2sub(size(C),ind_C_all);

%Create the 9x9 matrix indices
p=3*(I-1)+K;
q=3*(J-1)+L;

%Treat posible symbolic class
switch class(C)
    case 'sym'
        CM=sym(zeros(9,9));
    otherwise
        CM=zeros(9,9);
end

%Set values
[ind_pq]=sub2ind(size(CM),p,q);
CM(ind_pq(:))=C(:);

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function Us = get_PCproj(S1, row, col, wPCA, maskMaxChans)

[nT, nChan] = size(S1);
dt = -21 + [1:size(wPCA,1)];
inds = repmat(row', numel(dt), 1) + repmat(dt', 1, numel(row));

clips = reshape(S1(inds, :), numel(dt), numel(row), nChan);

mask = repmat([1:nChan], [numel(row) 1]) - repmat(col, 1, nChan);
Mask(1,:,:) = abs(mask)<maskMaxChans;

clips = bsxfun(@times, clips , Mask);

Us = wPCA' * reshape(clips, numel(dt), []);
Us = reshape(Us, size(wPCA,2), numel(row), nChan);

Us = permute(Us, [3 2 1]);

author: cortex-lab

repo: KiloSort

save_path: github-repos/MATLAB/cortex-lab-KiloSort