text:
%kcknn_printsig 'Print a K-Nearest Neighbors Signature '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros cknn_printsig.pane file
%
% Parameters:
% InputFile: isig 'Input Signature', required: 'input signature file'
% MultiChoice: omode 'Output Mode', default: 1: 'output mode toggle'
%    Choices are:
%   1: 'Text'
%   2: 'LaTeX'
%   3: 'HTML'
%   4: 'CSV'
% OutputFile: o 'Output', optional: 'Resulting output data object'
%
% Example: o = kcknn_printsig(isig, {'isig','';'omode',1;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% cknn_printsig - Print a K-Nearest Neighbors Signature
%
%  DESCRIPTION
% This operator prints out the contents of a K-Nearest Neighbors Signature file created with "cknn_signature". It expects a valid K-Nearest Neighbors signature in the parameter [-isig] so its input must be either the output of a \fIcknn_signature\fP or the output of \fIcsigappend\fP which can be used to append signatures of variable size. Its output can be formatted so it can be used in reports, etc. and will be redirected to kstdout unless a file name is passed with the parameter [-o].
% The output of this program can be in four formats, depending on the value of the parameter [-omode] :
% The Text (ASCII) format option ([-omode 1]) formats the output so it can be viewed with simple text editors or in kstdout.
% The LaTeX format option ([-omode 2]) formats the output so it can be included in a LaTeX file with the \\input command.
% The HTML format option ([-omode 3]) formats the output in a HTML table so it can be read into a HTML file.
% The CSV format option ([-omode 4]) formats the output in a Comma-delimited format so it can be read by spreadsheets, etc.
% One of the problems with the generation of LaTeX tables is that a single table is generated, and if there are too many signatures in the input the table will be too big to be displayed in a single page. In this case the user must edit the .tex file and break the table manually.
%
%
%
%  EXAMPLES
% All examples for the Classify toolbox are listed on the Classify Toolbox Manual. For an example of this routine, please see the example workspace Classify:workspaces:KNN-Signatures.
%
%  "SEE ALSO"
% cknn_signature, csigappend
%
%  RESTRICTIONS
% Expect the input file to be a valid signature file for the K-Nearest Neighbors classifier.
%
%  REFERENCES
% All references for the Classify toolbox are listed on the Classify Toolbox Manual.
%
%  COPYRIGHT
% Copyright (C) 1997 Rafael Santos. Khoros (C) Khoral Research, Inc.
%

function varargout = kcknn_printsig(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kcknn_printsig(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'isig', '__input';'omode', 1;'o', '__output'};
maxval={0,0,1};
minval={0,0,1};
istoggle=[0,0,1];
was_set=istoggle * 0;
paramtype={'InputFile','MultiChoice','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'cknn_printsig"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function value = not1d ( x1, x2, x3 )

%*****************************************************************************80
%
%% NOT1D differentiates a factor for serendipity basis functions.
%
%  Discussion:
%
%    not1(x1,x2,x3) evaluates at the point x1, the basis factor that
%    is 0 at x2 and 1 at x3:
%
%      not1(x1,x2,x3) = ( x1 - x2 ) / ( x3 - x2 )
%
%    This function returns the derivative of the factor with respect to x1:
%
%      not1d(x1,x2,x3) = 1 / ( x3 - x2 );
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X1, the evaluation point.
%
%    Input, real X2, X3, values that define the factor.
%
%    Output, real VALUE, the derivative of the basis function factor.
%
  value = 1.0 / ( x3 - x2 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef lasso
%function [Problem] = lasso(A, b, lambda, sub_mode, smooth_lambda)
% This file defines the lasso (least absolute shrinkage and selection operator) problem for L1 norm.
%
% Inputs:
%       A               dictionary matrix of size dxn.
%       b               observation vector of size dx1.
%       lambda          l1-regularized parameter.
%       sub_mode        {'prox_reg', 'subgrad_reg', 'smooth', 'const'}
%       smooth_lambda
% Output:
%       Problem     problem instance.
%
%
% The problem of interest is defined as
%
%       (1) sub_mode = 'prox_reg' for L1 regularizer solved by proximal methods
%
%           min f(w) = 1/2 * || A * w - b ||_2^2 + lambda * || w ||_1
%
%       (2) sub_mode = 'subgrad_reg' for L1 regularizer solved by subgrad methods
%
%           min f(w) = 1/2 * || A * w - b ||_2^2 + lambda * || w ||_1
%
%       (3) sub_mode = 'smooth' for smoothing
%
%           min f(w) = 1/2 * || A * w - b ||_2^2 + smooth_lambda * || w ||_1.
%
%       (43) sub_mode = 'const' for L1 constraint
%
%           min f(w) = 1/2 * || A * w - b ||_2^2,      s.t. || w ||_1 < r
%
%
% "w" is the model parameter of size d vector.
%
%
% This file is part of GDLibrary and SGDLibrary.
%
% Created by H.Kasai on Apr. 17, 2017
% Modified by H.Kasai on Mar. 25, 2018
% Modified by H.Kasai on Mar. 27, 2020
% Modified by H.Kasai on Mar. 30, 2020

    properties
        name;
        dim;
        samples;
        sub_mode;
        lambda;
        d;
        n;
        A;
        b;
        AtA;
        Atb;
        L;
        prox_flag;
        sign_flag;
        idx;
        smooth_mu;      % for smoothing
        smooth_type;    % for smoothing
    end

    methods
        function obj = lasso(A, b, lambda, sub_mode, varargin)

            obj.A = A;
            obj.b = b;

            if strcmp(sub_mode, 'prox_reg')
                obj.prox_flag = true;
            elseif strcmp(sub_mode, 'subgrad_reg')
                obj.prox_flag = false;
            elseif strcmp(sub_mode, 'const')
                obj.prox_flag = true;
            elseif strcmp(sub_mode, 'smooth')
                obj.prox_flag = false;
            else
                error('Not correctly specified for sub_mode');
            end

            obj.sub_mode = sub_mode;
            obj.lambda = lambda;

            if nargin < 6
                obj.smooth_type = 'type1';
            else
                obj.smooth_type = varargin{2};
            end

            if nargin < 5
                obj.smooth_mu = 0.1;
            else
                obj.smooth_mu = varargin{1};
            end

            obj.d = size(obj.A, 2);
            obj.n = size(obj.A, 2);

            obj.name = 'lasso';
            obj.dim = obj.d;
            obj.samples = obj.n;

            obj.AtA = obj.A'*obj.A;
            obj.Atb = obj.A'*obj.b;
            %L = max(eig(AtA));
            fprintf('Calculated Lipschitz constant (L), i.e., max(eig(AtA)), .... ')
            obj.L = eigs(obj.AtA, 1);
            fprintf('is L = %f.\n', obj.L);
        end

        function v = prox(obj, w, t)
            if strcmp(obj.sub_mode, 'prox_reg')
                v = soft_thresh(w, t * obj.lambda);
            elseif strcmp(obj.sub_mode, 'const')
                %v = sign(w) .* proj_simplex(abs(w), obj.lambda, 'ineq');
                v = proj_l1_ball(w, obj.lambda);
            else
                error('Not correctly specified for sub_mode');
            end
        end

        function f = cost(obj, w)
            reg = obj.reg(w);
            f = 1/2 * sum((obj.A * w - obj.b).^2) + obj.lambda * reg;
        end

        % calculate l1 norm
        function r = reg(obj, w)
            r = norm(w,1);
        end

        function r = residual(obj, w)
            r = - obj.A * w + obj.b;
        end

        function f = cost_batch(obj, w, indices)
            error('Not implemted yet.');
        end

        function g = full_grad(obj, w)
            %g = obj.A' * (obj.A * w - obj.obj.b);
            g = obj.AtA * w - obj.Atb;
        end

        function g = grad(obj, w, indices)
            A_partial = obj.A(:,indices);
            g = A_partial' * (A_partial * w - obj.b);
        end

        % for Subgradient descent
        function subg = full_subgrad(obj, w)
            subg = obj.lambda * sign(w);
        end

        % for Smoothing
        function smooth_g = full_smooth_grad(obj, w)
            smooth_g = zeros(obj.n, 1);

            if strcmp(obj.smooth_type, 'type1') % f_mu(w) = w^2/2mu for abs(w)<= mu, abs(w)-mu/2 for abs(w)> mu.
                for i = 1 : obj.n
                    if abs(w(i)) <= obj.smooth_mu
                        smooth_g(i) = w(i) / obj.smooth_mu;
                    elseif w(i) > obj.smooth_mu
                        smooth_g(i) =1;
                    else
                        smooth_g(i) =-1;
                    end
                end
            elseif strcmp(obj.smooth_type, 'type2') % f_mu(w) = (w^2+mu^2)^(1/2)
                for i = 1 : obj.n
                    smooth_g(i) = w(i)/sqrt(w(i)^2+obj.smooth_mu^2);
                end
            elseif strcmp(obj.smooth_type, 'type3') % f_mu(w) = mu log((exp(-w/mu) + exp(w/mu))/2)
                for i = 1 : obj.n
                    p = exp(-w(i)/obj.smooth_mu);
                    q = exp(w(i)/obj.smooth_mu);
                    smooth_g(i) = (q-p)/(p+q);
                end
                %smooth_g = smooth_g;
            end

            smooth_g = obj.lambda * smooth_g;
        end

        % for Frank-Wolfe, a.k.a. Condtional Gradient
        function [s, idx, sign_flag] = LMO(obj, grad)
            [val, idx] = max( abs(grad) );
            obj.sign_flag = -1 * sign(grad(idx));
            obj.idx = idx;

            s = zeros(size(grad));
            s(idx) = obj.sign_flag * obj.lambda;

            sign_flag = obj.sign_flag;
        end

        function h = hess(obj, w, indices)
            error('Not implemted yet.');
        end

        function h = full_hess(obj, w)
            h = obj.AtA;
        end

        function hv = hess_vec(obj, w, v, indices)
            error('Not implemted yet.');
        end

        function step = exact_line_search(obj, w, dir, idx, sign_flag)
            As = sign_flag * obj.lambda * obj.A(:,idx); % = the i-th column of the dictionary matrix A
            As_minus_Ax = As - obj.A*w;
            step = max(0, min(1, As_minus_Ax' * (obj.b-obj.A*w) / (As_minus_Ax' * As_minus_Ax)));
        end

        function w_opt = calc_solution(obj, method, options_in)

            if nargin < 1
                method = 'ista';
            end

            options.max_epoch = options_in.max_epoch;
            options.w_init = options_in.w_init;
            options.verbose = options_in.verbose;
            options.tol_optgap = 1.0e-24;
            options.tol_gnorm = 1.0e-16;
            options.step_alg = 'backtracking';

            if ~options.verbose
                fprintf('Calclation of solution started ... ');
            end

            if strcmp(method, 'sd')
                [w_opt,~] = sd(obj, options);
            elseif strcmp(method, 'cg')
                [w_opt,~] = ncg(obj, options);
            elseif strcmp(method, 'newton')
                options.sub_mode = 'INEXACT';
                options.step_alg = 'non-backtracking';
                [w_opt,~] = newton(obj, options);
            elseif strcmp(method, 'ista')
                options.step_alg = 'fix';
                options.step_init = 1/obj.L;
                options.sub_mode  = 'FISTA';
                [w_opt,~] = ag(obj, options);
            else
                options.step_alg = 'backtracking';
                options.mem_size = 5;
                [w_opt,~] = lbfgs(obj, options);
            end

            if ~options.verbose
                fprintf('done\n');
            end
        end
    end
end

author: hiroyuki-kasai

repo: SGDLibrary


text:
classdef SteerMathModelTypeEnum < matlab.mixin.SetGet
    %SteerMathModelTypeEnum Summary of this class goes here
    %   Detailed explanation goes here

    enumeration
        GenericPoly('Polynomial');
        SumOfSines('Sum of Sines');
        LinearTangent('Linear Tangent');
        FitNet('Neural Network');
    end

    properties
        name(1,:) char
    end

    methods
        function obj = SteerMathModelTypeEnum(name)
            obj.name = name;
        end
    end

    methods(Static)
        function [listBoxStr, enums] = getListBoxStr()
            m = enumeration('SteerMathModelTypeEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
            enums = m(I);
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('SteerMathModelTypeEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('SteerMathModelTypeEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function featvec = MRS4fast(im)

[featvec] = MR8fast(im);
featvec = featvec';

featvec = cat(2, featvec(:, 1:2), max(featvec(:, 3:5), [], 2), max(featvec(:, 6:8), [], 2));

author: Cloud-CV

repo: object-proposals


text:
function res = transformLine3d(line, trans)
%TRANSFORMLINE3D Transform a 3D line with a 3D affine transform
%
%   LINE2 = transformLine3d(LINE1, TRANS)
%
%   Example
%   P1 = [10 20 30];
%   P2 = [30 40 50];
%   L = createLine3d(P1, P2);
%   T = createRotationOx(P1, pi/6);
%   L2 = transformLine3d(L, T);
%   figure; hold on;
%   axis([0 100 0 100 0 100]); view(3);
%   drawPoint3d([P1;P2]);
%   drawLine3d(L, 'b');
%   drawLine3d(L2, 'm');
%
%   See also:
%   lines3d, transforms3d, transformPoint3d, transformVector3d
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2008-11-25,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2008 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

res = [...
    transformPoint3d(line(:, 1:3), trans) ...   % transform origin point
    transformVector3d(line(:,4:6), trans)];     % transform direction vect.

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [Dx,Dy] = PhysDmatrices2D(x1, y1, interp)

% function [Dr,Ds] = PhysDmatrices2D(x1, y1, interp)
% Purpose : Initialize the (x,y) differentiation matrices on the simplex
%	        evaluated at (x1,y1) at order N

Globals2D;
IDr = interp*Dr; IDs = interp*Ds;

[rx1,sx1,ry1,sy1,J1] = GeometricFactors2D(x1, y1, IDr, IDs);

n = size(interp, 1);
Dx = spdiags(rx1, 0, n, n)*IDr + spdiags(sx1, 0, n, n)*IDs;
Dy = spdiags(ry1, 0, n, n)*IDr + spdiags(sy1, 0, n, n)*IDs;
return;

repo: Aero-matlab


text:
function err = sb_test_ori(pnt,elem)
err = 1;
if(size(elem,2) == 4)
    det = sum(cross(pnt(elem(:,2),:)-pnt(elem(:,1),:),pnt(elem(:,4),:)-pnt(elem(:,1),:),2).*(pnt(elem(:,3),:)-pnt(elem(:,1),:)),2);
    if sum(det <= 0) > 0
        err = 0;
    end
elseif(size(elem,2) == 8)
    det1 = sum(cross(pnt(elem(:,6),:)-pnt(elem(:,1),:),pnt(elem(:,8),:)-pnt(elem(:,1),:),2).*(pnt(elem(:,5),:)-pnt(elem(:,1),:)),2);
    det2 = sum(cross(pnt(elem(:,3),:)-pnt(elem(:,1),:),pnt(elem(:,6),:)-pnt(elem(:,1),:),2).*(pnt(elem(:,2),:)-pnt(elem(:,1),:)),2);
    det3 = sum(cross(pnt(elem(:,8),:)-pnt(elem(:,1),:),pnt(elem(:,3),:)-pnt(elem(:,1),:),2).*(pnt(elem(:,4),:)-pnt(elem(:,1),:)),2);
    det4 = sum(cross(pnt(elem(:,8),:)-pnt(elem(:,3),:),pnt(elem(:,6),:)-pnt(elem(:,3),:),2).*(pnt(elem(:,7),:)-pnt(elem(:,3),:)),2);
    if sum(det1 <= 0) || sum(det2 <= 0) || sum(det3 <= 0) || sum(det4 <= 0)
        err = 0;
    end
else
    error('Invalid number of nodes per element!');
end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% rigid registration of frames with offsets ds
function dreg = rigidRegFrames(data, ops, ds)

% compute maximum batch size for GPU
[h, w, nFrames] = size(data);
if getOr(ops, 'useGPU', 0)
    nFramesPerBatch = getBatchSize(w*h);
else
    nFramesPerBatch = 1000; % if not on GPU, should have plenty of RAM available
end
%%
nFramesPerBatch = 8;

nBatches = ceil(nFrames/nFramesPerBatch);
startFrame = 1:nFramesPerBatch:nFrames;
endFrame = min(startFrame+nFramesPerBatch-1, nFrames);
dreg = zeros(size(data), 'single');

for iBatch = 1:nBatches
    idx = startFrame(iBatch):endFrame(iBatch);
    if ops.useGPU
        dataBatch = gpuArray(single(data(:,:,idx)));
    else
        dataBatch = single(data(:,:,idx));
    end
    [Ly, Lx, NT] = size(dataBatch);

    Ny = ifftshift([-fix(Ly/2):ceil(Ly/2)-1]);
    Nx = ifftshift([-fix(Lx/2):ceil(Lx/2)-1]);
    [Nx,Ny] = meshgrid(Nx,Ny);
    Nx = Nx / Lx;
    Ny = Ny / Ly;

    if ops.useGPU
        dsBatch = gpuArray(permute(ds(idx, :), [3, 2, 1]));
        Nx = gpuArray(single(Nx));
        Ny = gpuArray(single(Ny));
    else
        dsBatch = ds(idx, :);
    end

    if ops.useGPU % do it batch-by-batch
        dph         = 2*pi*(bsxfun(@times, dsBatch(1,1,:), Ny) + ...
            bsxfun(@times, dsBatch(:,2,:), Nx));
        fdata       = fft2(dataBatch);
        dregBatch   = gather_try(real(ifft2(fdata .* exp(1i * dph))));
    else % do it frame-by-frame
        dregBatch = zeros(size(dataBatch), 'single');
        for i = 1:NT
            dph         = 2*pi*(dsBatch(i,1)*Ny + dsBatch(i,2)*Nx);
            fdata       = fft2(single(dataBatch(:,:,i)));
            dregBatch(:,:,i) = real(ifft2(fdata .* exp(1i * dph)));
        end
    end

    dreg(:,:,idx) = dregBatch;
end

author: cortex-lab


text:
function bow2 = mydatesow_aux (epoch2)
    % beginning of day:
    epoch2_vec = mydatevec(epoch2);
    bod2_vec = epoch2_vec;  bod2_vec(:,4:6) = 0;
    bod2 = mydatenum(bod2_vec);

    % beginning of week:
    %weekday(datenum(bod2_vec))-1  % DEBUG
    bow2 = bod2 - (weekday(datenum(bod2_vec))-1)*24*3600;
end

%!test
%! % mydatesow_aux()
%! test('mydatesow')
%! test('mydatesowi')

author: Sable

repo: mcbench-benchmarks


text:
function [results] = tapas_ti_example()
%% This script simulates data from different models and checks where it comes
% from.
%
% Input
%   mi          Model to be used for inversion
%   session     Session of the experiment (either 1 or 2)
%   block       Block of the experiment, either 3, 4 or 5
%   sbj         Subject
% Output
%
% aponteeduardo@gmail.com
% copyright (C) 2015
%

DIM_THETA = sooner_gamma_ndims();

ptheta = tapas_sooner_ptheta();
ptheta.p0(13:15) = [log(1.0) log(0.5) log(0.5)];
htheta = tapas_sooner_htheta();

pars = struct();
pars.T = linspace(0.0001, 1, 16).^5;
pars.nburnin = 5000;
pars.niter = 5000;
pars.mc3it = 8;
pars.verbose = 1;
pars.samples = 1;

data = load(mat_fname());
data = data.table;
[y, u] = prepare_data_db(data);
[y, u] = select_data(y, u, 3, 1);

[pt, tfe]  = tapas_ti_estimate(y, u, ptheta, htheta, pars);
results.pt = pt;
results.fe = tfe;
results.y = y;

end % gamma_inversion

function [y, u] = select_data(y, u, block, session)
%% Selects a chuck of the data

block = block * session;

y.t = y.t(u.b == block);
y.a = y.a(u.b == block);
y.i = y.i(u.b == block);

u.tt = u.tt(u.b == block);

end % select_data

function [fname] = mat_fname()
%% Generates the name of the simulation

% Get current location
f = mfilename('fullpath');
[tdir, ~, ~] = fileparts(f);

fname = fullfile(tdir, 'data', 'tapas_sooner.mat');

end %

function [y, u] = prepare_data_db(data)
%% Generates the data ready to be feeded into the inversion method
%
% Input
%
%   data    Matrix as generated by the python script.
%
% Output
%   y       Subjects time, actions, and trial validity
%   u       Trial type and block
%

% aponteeduardo@gmail.com
% copyright (C) 2015
%

y = struct('t', [], 'a', [], 'i', []);
u = struct('tt', [], 'b', []);

nt = size(data, 1);

% Tenths of second!
y.t = data(:, 5)/100;
y.a = data(:, 6);
y.i = zeros(nt, 1);
y.i(isnan(y.t)) = 1;

u.tt = data(:, 4);
u.b =  data(:, 1) .* data(:, 2);

% Because of an inconsistency

for i = 1:nt
    if y.a(i) == 0
        y.a(i) = 1;
    else
        y.a(i) = 0;
    end
    if u.tt(i) == 1
        u.tt(i) = 0;
    else
        u.tt(i) = 1;
    end
end

end % prepare_data_db

author: translationalneuromodeling

repo: tapas


text:
function plane = medianPlane(p1, p2)
%MEDIANPLANE Create a plane in the middle of 2 points
%
%   PLANE = medianPlane(P1, P2)
%   Creates a plane in the middle of 2 points.
%   PLANE is perpendicular to line (P1 P2) and contains the midpoint of P1
%   and P2.
%   The direction of the normal of PLANE is the same as the vector from P1
%   to P2.
%
%   See also:
%   planes3d, createPlane
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 18/02/2005.
%

%   HISTORY
%   28/06/2007: add doc, and manage multiple inputs

% unify data dimension
if size(p1, 1)==1
    p1 = repmat(p1, [size(p2, 1) 1]);
elseif size(p2, 1)==1
    p2 = repmat(p2, [size(p1, 1) 1]);
elseif size(p1, 1)~=size(p2, 1)
    error('data should have same length, or one data should have length 1');
end

% middle point
p0  = (p1 + p2)/2;

% normal to plane
n   = p2-p1;

% create plane from point and normal
plane = createPlane(p0, n);


text:
function [f cumf] = sumlog2(S)
% returns log2(sum(2.^S)) in a stable way

% We compute log2(f) using our log-sum method. Precision of this
% method is approximately 2^(-1024) * (number of summands). So it is
% much less than eps = 2^(-52).
%
% Indeed, we compute (assume A>B)
%      log(A + B) = log A  + log (1 + B/A);
% B/A = 2^{log B - log A}. So we only lose precision when B/A gets below 2^{-1024}.
% But log (1 + x) <= x. Thus each time we get error less than 2^{-1024}.
%

const_log2e = 1/log(2);
log_sum = -Inf;
cumf = zeros(size(S));
idx = 1;
for cur_term = S;
	if cur_term > -Inf
		if cur_term > log_sum
			log_sum = cur_term + log1p(2^(log_sum - cur_term))*const_log2e;
		else
			log_sum = log_sum + log1p(2^(cur_term - log_sum))*const_log2e;
		end;
	end
	cumf(idx) = log_sum;
	idx = idx + 1;
end

f = log_sum;

author: yp-mit

repo: spectre

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


text:
function [mad_similarity cos_similarity r_similarity] = outliers_xval(obj)
% Similarity metrics comparing each image to the mean of others in the set using repeated k-fold cross-validation
%
% :Usage:
% ::
%
%     [mad_similarity cos_similarity r_similarity] = outliers_xval(obj)
%
% :Background:
% ::
% Assessing image outliers is crucial for detecting artifacts and other
% violations of assumptions underlying image analysis.
%
% Many methods are available, but they often use measures of variance defined relative to a sample,
% which creates problems. For example, Mahalanobis distance uses squared
% multivariate distances among images. This is useful, but detects only
% relative outliers. If there are many corrupted images, all images will
% look normal by comparison to others.  Inter-image correlations are also
% commonly used. These are also useful, but they are also sensitive to the
% overall distribution of artifacts across the set. If there are many
% corrupted images, correlations will tend to be low overall, and the bad
% images will be harder to detect. In cases where comparing individual
% images to a gold standard is desired, but there is no external standard
% image (perhaps because images look somewhat different in each sample),
% comparing images to a mean image derived from others in the set may be useful.
% In this case, however, it's desirable for the mean to be derived
% independent of the test image itself. This can be accomplished with
% cross-validation, bootstrapping, or jackknife analyses.
%
% Here, we compare each of a set of images to a mean derived from other
% images in the set using k-fold cross-validation. To average over errors
% related to the particular fold splits chosen, we average over a number of
% repeated cross-valiations.
%
% We return two error metrics, which have different desirable properties
% depending on the type of image and use case.
% - Pearson's correlation (r): correlation is insensitive to mean shift and scale
%   it is most useful when only the pattern in the image is meaningful
%   and/or images will be rescaled in analysis, and the zero-point and
%   scale are not important.
%
% - cosine similarity is sensitive to mean shift, not scale
%   This is useful for images with a meaningful zero-point that should match across images,
%   but where the scale is irrelevant (or will be removed, e.g., via
%   normalization)
%
% - Absolute agreement is most relevant when mean and scale are both relevant
%   This is the case with many images, e.g., those subjected to group statistical analysis of intensity
%   values. e.g., when testing task - control contrast images across a group,
%   We are interested in whether the group mean is 0. The zero-point is meaningful and
%   it should agree if the individual test images are replicates of the same effect.
%   The scale is also meaningful, as deviations in the values determine the variance estimate
%   of the test statistic. Any deviation is treated as error, even if it results from
%   mean shift or scale varation in the image as a whole.
%
%
% :Inputs:
%
%   **obj:**
%        an image_vector object (e.g., fmri_data) with >=5 images
%
% :Outputs:
%
%   **mad_similarity:**
%        [nimages x 1] vector of absolute agreement between each image and the mean,
%        1/(1+MAD), averaged across repeated k-fold
%
%   **cos_similarity:**
%        [nimages x 1] vector of cosine similarity between each image and the mean,
%        averaged across repeated k-fold
%
%   **r_similarity:**
%        [nimages x 1] vector of correlation between each image and the mean,
%        averaged across repeated k-fold
%
% ----------------------------------------------------------------------
% Examples:
% ----------------------------------------------------------------------
% test_images = load_image_set('emotionreg');
%
% [mad_similarity cos_similarity r_similarity ] = outliers_xval(test_images);
% figure; plot(mad_similarity)
% hold on;
% plot(cos_similarity)
% plot(r_similarity)
% legend({'Abs agreement' 'Cos sim' 'r'})
%
% % Compare with outliers based on Mahalanobis and other metrics
% figure;
% [est_outliers_uncorr, est_outliers_corr, outlier_tables] = outliers(test_images, 'notimeseries');
%
% Another, larger sample dataset
% test_images = load_image_set('kragel18_alldata');
% [mad_similarity cos_similarity r_similarity ] = outliers_xval(test_images);

% ..
%     Author and copyright information:
%
%     Copyright (C) 2023 Tor Wager
%
%     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/>.
% ..

% ----------------------------------------------------------------------
% Default parameters
% ----------------------------------------------------------------------

nfolds = 5;
nrepeats = 100;  % number of cross-validation repeats

% ----------------------------------------------------------------------
% Set up variables and functions
% ----------------------------------------------------------------------

% Analyze gray matter, where we expect meaningful variation
%obj_orig = obj;
obj = apply_mask(obj, fmri_data(which('gray_matter_mask.nii'), 'noverbose'));

nimages = size(obj.dat, 2);

if nimages < 5, error('Object must contain at least 5 images'), end

Yid = ones(nimages, 1);

[r_test, cos_sim_test, mad_test] = deal(NaN .* zeros(nimages, nrepeats));

cvpart = cvpartition(Yid, 'k', nfolds);
% [S.trIdx, S.teIdx] = xval_stratified_holdout_leave_whole_subject_out(S.Y, S.id, 'doverbose', doverbose, 'doplot', doplot, 'nfolds', nfolds);

cos_sim =  @(x, y) x' * y ./ (norm(x) * norm(y));

% ----------------------------------------------------------------------
% Run
% ----------------------------------------------------------------------

for p = 1:nrepeats

    % Draw a new k-fold cross-validation partition
    cvpart = cvpart.repartition;

    for i = 1:nfolds

        train_set = get_wh_image(obj, cvpart.training(i));
        test_set = get_wh_image(obj, cvpart.test(i));

        m = mean(train_set);

        % correlation is insensitive to mean shift and scale
        % cosine sim is sensitive to mean shift, not scale
        % for images with a meaningful zero-point that should match across images,
        % cosine sim is more meaningful.
        % if mean and scale are both relevant, as with contrast images, absolute
        % agreement may be the best metric, e.g. MAD = median absolute deviation

        r_test(cvpart.test(i), p) = corr(m.dat, test_set.dat)';

        cos_sim_test(cvpart.test(i), p) = cos_sim(m.dat, test_set.dat)';

        mad_test(cvpart.test(i), p) = median(abs(m.dat - test_set.dat))';

    end % k-fold

end % cv repeats

% ----------------------------------------------------------------------
% Calculate final similarity measures
% ----------------------------------------------------------------------

r_similarity = mean(r_test')';
cos_similarity = mean(cos_sim_test')';

mad_dissim = mean(mad_test')'; % this is actually still deviation here
mad_similarity = 1./(1+mad_dissim); % convert to similarity, scale from [0 1]
% +1 scales so that zero error (perfect agreement) will have a value of 1.

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% DOC
%
% Files
%   afemdoc              - Adaptive Finite Element Method
%   amgdoc               - ALGEBRAIC MULTIGRID METHOD
%   amgdoctest1          - AMG TEST I: DIFFERENT MESHES
%   amgdoctest2          - AMG TEST II: Different boundary conditions
%   amgdoctest3          - AMG TEST III: Robustness to time discretization
%   amgtest              - test algebraic mutligrid solver
%   amgtest2             - AMG TEST II: Different boundary conditions
%   auxstructuredoc      - Auxiliary Mesh Data Structure
%   bddoc                - Data Structure: Boundary Conditions
%   bisectdoc            - Bisection in Two Dimensions
%   coarsenAMGdoc        - Coarsening for Algebraic Multigrid
%   coarsendoc           - COARSENING in TWO DIMENSIONS
%   coloringdoc          - Coloring Vertices of a Graph
%   dof3edgedoc          - Data Structure: Lowest Order Edge Element
%   dofdoc               - Data Sturcture for Degree of Freedom (DOF)
%   dofedgedoc           - Dofedge in Two Dimensions
%   dofP2doc             - Data Structure: P2 Quadratic Element
%   fastcoding           - Vectorization
%   femcontent           - Finite Element Methods
%   femdoc               - Finite Element Methods
%   ifem                 - Display HTML documentation in the MATLAB web browser.
%   introduction         - Introduction of iFEM
%   Maxwell1doc          - Equation: Maxwell Equation Discretized by Edge Element
%   Maxwell1testdoc      - Edge Element Discretization of Maxwell Equations: Linear Element
%   Maxwell2doc          - Equation: Maxwell equation Quadratic Element in 3D
%   Maxwell2testdoc      - Edge Element Discretization of Maxwell Equations
%   Maxwelldoc           - Equation: Maxwell Equation Discretized by Edge Element
%   MaxwellNeumanBCdoc   - Local Labeling of DOFs
%   Maxwelltestdoc       - Edge Element Discretization of Maxwell Equations
%   meshamrdoc           - Adaptive Mesh Refinement/Coarsening
%   meshbasicdoc         - Basic Data Structure representing a Mesh
%   meshdoc              - Mesh
%   meshoptdoc           - Mesh Smoothing and Optimization
%   mgdoc                - MULTIGRID on BISECTION GRIDS
%   mixFEMsolverdoc      - Mixed FEM solver
%   Poisson3BDM1doc      - Equation: Poisson Equation Discretized by $BDM_1$ Element in 3D
%   Poisson3RT0doc       - Equation: Poisson Equation Discretized by $RT_0$ Element in 3D
%   PoissonBDM1doc       - Equation: Poisson Equation Discretized by $BDM_1$ Element in 2D
%   PoissonRT0doc        - Equation: Poisson Equation Discretized by $RT_0$ Element in 2D
%   sc3doc               - Simplicial Complex in Three Dimensions
%   scdoc                - Simplicial Complex in Three Dimensions
%   simplicialcomplexdoc - Simplicial Complex in 3D
%   solverdoc            - Multigrid Method
%   sparsematrixdoc      - Vectorization using Sparse Matrices
%   uniformrefine3doc    - 3-D Red Refinement


text:
clc
clear
hold on

n=100;          % number of intervals (i.e. parametric curve would be evaluted n+1 times)
k=0.5522847498; % kappa (See documentation how this value is obtained)

% First Quadrant
Px=[1 1 k 0];
Py=[0 k 1 1];

PlotBezier1(Px,Py,n);

% Second Quadrant
Px=[0 -k -1 -1];
Py=[1  1  k  0];
PlotBezier1(Px,Py,n)

% Third Quadrant
Px=[-1 -1 -k  0];
Py=[ 0 -k -1 -1];
PlotBezier1(Px,Py,n)

% Fourth Quadrant
Px=[ 0  k  1 1];
Py=[-1 -1 -k 0];
PlotBezier1(Px,Py,n)

title('Approximation of Circle using Cubic Bezier');

hold off

% % % --------------------------------
% % % Author: Dr. Murtaza Khan
% % % Email : drkhanmurtaza@gmail.com
% % % --------------------------------


text:
function i4vec_part_quick_a_test ( )

%*****************************************************************************80
%
%% I4VEC_PART_QUICK_A_TEST tests I4VEC_PART_QUICK_A.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 12;
  b = 0;
  c = n;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_PART_QUICK_A_TEST\n' );
  fprintf ( 1, '  I4VEC_PART_QUICK_A reorders an integer vector\n' );
  fprintf ( 1, '  as part of a quick sort.\n' );
  fprintf ( 1, '  Using initial random number seed = %d\n', seed );

  [ a, seed ] = i4vec_uniform_ab ( n, b, c, seed );

  i4vec_print ( n, a, '  Before rearrangement:' );

  [ a, l, r ] = i4vec_part_quick_a ( n, a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Rearranged array\n' );
  fprintf ( 1, '  Left index =  %d\n', l );
  fprintf ( 1, '  Key index =   %d\n', l+1 );
  fprintf ( 1, '  Right index = %d\n', r );

  i4vec_print ( l,     a(1:l),   '  Left half:' );
  i4vec_print ( 1,     a(l+1),   '  Key:' );
  i4vec_print ( n-l-1, a(l+2:n), '  Right half:' );

  return
end


text:
function microbeID=adaptDraftModelID(modelIDIn)
% This function adapts the name of the draft reconstruction entered into
% the pipeline into a suitable curated reconstruction ID.
%
% USAGE: microbeID = adaptKBaseModelID(modelIDIn)
%
% INPUT
% modelIDIn     name of the KBase draft reconstruction file to refine
%
% OUTPUT
% microbeID     name of the resulting refined reconstruction
%
% AUTHOR: Almut Heinken, 06/2020.

microbeID=modelIDIn;
microbeID=strrep(microbeID,'.fbamodel.sbml','');
microbeID=strrep(microbeID,'.model.sbml','');
microbeID=strrep(microbeID,'___RAST','');
microbeID=strrep(microbeID,'__.RAST','');
microbeID=strrep(microbeID,'_.RAST','');
microbeID=strrep(microbeID,'-.RAST','');
microbeID=strrep(microbeID,'..RAST','');
microbeID=strrep(microbeID,'.RAST','');
microbeID=strrep(microbeID,'_RAST','');
microbeID=strrep(microbeID,'.fna','');
microbeID=strrep(microbeID,'.fa','');
microbeID=strrep(microbeID,'_DraftModel','');
microbeID=strrep(microbeID,'.contigs','');
microbeID=strrep(microbeID,'__genome','');
microbeID=strrep(microbeID,'_genome','');
microbeID=strrep(microbeID,'_assembly','');
microbeID=strrep(microbeID,'.sbml','');
microbeID=strrep(microbeID,'.mat','');
microbeID=strrep(microbeID,' ','_');
microbeID=strrep(microbeID,'.','_');
microbeID=strrep(microbeID,'-','_');
microbeID=strrep(microbeID,'(','_');
microbeID=strrep(microbeID,')','_');
microbeID=strrep(microbeID,'\','_');
microbeID=strrep(microbeID,'/','_');
microbeID=strrep(microbeID,'+','_');
microbeID=strrep(microbeID,'©','_');
microbeID=strrep(microbeID,'____','_');
microbeID=strrep(microbeID,'___','_');
microbeID=strrep(microbeID,'__','_');
microbeID=strrep(microbeID,',','_');
end

author: opencobra

repo: cobratoolbox


text:
function comp_random_grlex_test ( )

%*****************************************************************************80
%
%% COMP_RANDOM_GRLEX_TEST tests COMP_RANDOM_GRLEX.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'COMP_RANDOM_GRLEX_TEST\n' );
  fprintf ( 1, '  A COMP is a composition of an integer N into K parts.\n' );
  fprintf ( 1, '  Each part is nonnegative.  The order matters.\n' );
  fprintf ( 1, '  COMP_RANDOM_GRLEX selects a random COMP in\n' );
  fprintf ( 1, '  graded lexicographic (grlex) order between indices RANK1 and RANK2.\n' );
  fprintf ( 1, '\n' );

  kc = 3;
  rank1 = 20;
  rank2 = 60;
  seed = 123456789;

  for test = 1 : 5
    [ xc, rank, seed ] = comp_random_grlex ( kc, rank1, rank2, seed );
    nc = sum ( xc(1:kc) );
    fprintf ( 1, '   %3d: ', rank );
    fprintf ( 1, '    %2d = ', nc );
    for j = 1 : kc - 1
      fprintf ( 1, '%2d + ', xc(j) );
    end
    fprintf ( 1, '%2d\n', xc(kc) );
  end

  return
end


text:
function hb_header_print ( title, key, totcrd, ptrcrd, indcrd, valcrd, ...
  rhscrd, mxtype, nrow, ncol, nnzero, neltvl, ptrfmt, indfmt, valfmt, ...
  rhsfmt, rhstyp, nrhs, nrhsix )

%*****************************************************************************80
%
%% HB_HEADER_PRINT prints the header of an HB file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992.
%
%  Parameters:
%
%    Input, character ( len = 72 ) TITLE, a title for the matrix.
%
%    Input, character ( len = 8 ) KEY, an identifier for the matrix.
%
%    Input, integer TOTCRD, the total number of lines of data.
%
%    Input, integer PTRCRD, the number of input lines for pointers.
%
%    Input, integer INDCRD, the number of input lines for row indices.
%
%    Input, integer VALCRD, the number of input lines for numerical values.
%
%    Input, integer RHSCRD, the number of input lines for right hand sides.
%
%    Input, character ( len = 3 ) MXTYPE, the matrix type.
%    First character is R for Real, C for complex, P for pattern only.
%    Second character is S for symmetric, U for unsymmetric, H for
%      Hermitian, Z for skew symmetric, R for rectangular.
%    Third character is A for assembled and E for unassembled
%      finite element matrices.
%
%    Input, integer NROW, the number of rows or variables.
%
%    Input, integer NCOL, the number of columns or elements.
%
%    Input, integer NNZERO.  In the case of assembled sparse matrices,
%    this is the number of nonzeroes.  In the case of unassembled finite
%    element matrices, in which the right hand side vectors are also
%    stored as unassembled finite element vectors, this is the total
%    number of entries in a single unassembled right hand side vector.
%
%    Input, integer NELTVL, the number of finite element matrix entries,
%    set to 0 in the case of assembled matrices.
%
%    Input, character ( len = 16 ) PTRFMT, the format for reading pointers.
%
%    Input, character ( len = 16 ) INDFMT, the format for reading indices.
%
%    Input, character ( len = 20 ) VALFMT, the format for reading values.
%
%    Input, character ( len = 20 ) RHSFMT, the format for reading values
%    of the right hand side.
%
%    Input, character ( len = 3 ) RHSTYP, the right hand side type.
%    First character is F for full storage or M for same as matrix.
%    Second character is G if starting "guess" vectors are supplied.
%    Third character is X if exact solution vectors are supplied.
%
%    Input, integer NRHS, the number of right hand sides.
%
%    Input, integer NRHSIX, the number of entries of storage for right
%    hand side values, in the case where RHSTYP(1:1) = 'M' and
%    MXTYPE(3:3) = 'A'.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  TOTCRD = %6d\n', totcrd );
  fprintf ( 1, '  PTRCRD = %6d\n', ptrcrd );
  fprintf ( 1, '  INDCRD = %6d\n', indcrd );
  fprintf ( 1, '  VALCRD = %6d\n', valcrd );
  fprintf ( 1, '  RHSCRD = %6d\n', rhscrd );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  KEY =    "%s".\n', key );
  fprintf ( 1, '  MXTYPE = "%s".\n', mxtype );
  fprintf ( 1, '  RHSTYP = "%s".\n', rhstyp );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  NROW   = %6d\n', nrow );
  fprintf ( 1, '  NCOL   = %6d\n', ncol );
  fprintf ( 1, '  NNZERO = %6d\n', nnzero );
  fprintf ( 1, '  NELTVL = %6d\n', neltvl );
  fprintf ( 1, '  NRHS   = %6d\n', nrhs );
  fprintf ( 1, '  NRHSIX = %6d\n', nrhsix );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  PTRFMT =    "%s".\n', ptrfmt );
  fprintf ( 1, '  INDFMT =    "%s".\n', indfmt );
  fprintf ( 1, '  VALFMT =    "%s".\n', valfmt );
  fprintf ( 1, '  RHSFMT =    "%s".\n', rhsfmt );

  return
end


text:
function [accuracy,predictlabel] = SRKDApredict(fea, gnd, model)
% SRKDApredict: Spectral Regression Kernel Discriminant Analysis Prediction
%               SRKDApredict use SRKDA as a classifier. It used the nearest
%               center rule in the SRKDA subspace for classification.
%
%       [predictlabel,accuracy,elapse] = SRKDApredict(fea, gnd, model);
%
%             Input:
%
%               fea     - data matrix. Each row is a data point.
%               gnd     - Label vector of fea.
%             model     - model trained by SRKDAtrain.m
%
%             Output:
%
%            accuracy   - classification accuracy
%         predictlabel  - predict label for fea
%
%    Examples:
%
%
% See also SRKDAtrain, KSR, KSR_caller
%
%Reference:
%
%   [1] Deng Cai, Xiaofei He, and Jiawei Han. "Speed Up Kernel Discriminant
%   Analysis", The VLDB Journal, vol. 20, no. 1, pp. 21-33, January, 2011.
%
%   [2] Deng Cai, Xiaofei He and Jiawei Han, "SRDA: An Efficient Algorithm for
%   Large Scale Discriminant Analysis" IEEE Transactions on Knowledge and
%   Data Engineering, vol. 20, no. 1, pp. 1-12, January, 2008.
%
%   [3] Deng Cai, "Spectral Regression: A Regression Framework for
%   Efficient Regularized Subspace Learning", PhD Thesis, Department of
%   Computer Science, UIUC, 2009.
%
%   version 3.0 --Jan/2012
%   version 2.0 --December/2011
%   version 1.0 --May/2006
%
%   Written by Deng Cai (dengcai AT gmail.com)
%

MAX_MATRIX_SIZE = 8000; % You can change this number based on your memory.

if ~strcmp(model.TYPE,'SRKDA')
    error('model does not match!');
end

nTrain = size(model.fea,1);
nTest = size(fea,1);
nBlock = ceil(MAX_MATRIX_SIZE*MAX_MATRIX_SIZE/nTrain);
if model.LARs
    accuracy = zeros(length(model.LassoCardi),1);
    predictlabel = zeros(nTest,length(model.LassoCardi));
    Embed_Test = cell(length(model.LassoCardi),1);
    for i=1:length(model.LassoCardi)
        Embed_Test{i} = zeros(nTest,size(model.projection{i},2));
    end
    for j = 1:ceil(nTest/nBlock)
        if j == ceil(nTest/nBlock)
            smpIdx = (j-1)*nBlock+1:nTest;
        else
            smpIdx = (j-1)*nBlock+1:j*nBlock;
        end
        KTest= constructKernel(fea(smpIdx,:),model.fea,model.options);
        if model.bSemi
            KTest = KTest*model.KtestHat;
        end
        for i=1:length(model.LassoCardi)
            if model.bSemi
                Embed_Test{i}(smpIdx,:) = KTest(:,1:model.nLabel)*model.projection{i};
            else
                Embed_Test{i}(smpIdx,:) = KTest*model.projection{i};
            end
        end
        clear KTest;
    end

    for i=1:length(model.LassoCardi)
        D = EuDist2(Embed_Test{i},model.ClassCenter{i},0);
        [dump, idx] = min(D,[],2);
        predictlabel(:,i) = model.ClassLabel(idx);
        accuracy(i) = 1 - length(find(predictlabel(:,i)-gnd))/nTest;
    end
else
    Embed_Test = zeros(nTest,size(model.projection,2));
    for i = 1:ceil(nTest/nBlock)
        if i == ceil(nTest/nBlock)
            smpIdx = (i-1)*nBlock+1:nTest;
        else
            smpIdx = (i-1)*nBlock+1:i*nBlock;
        end
        KTest= constructKernel(fea(smpIdx,:),model.fea,model.options);
        if model.bSemi
            KTest = KTest*model.KtestHat;
            Embed_Test(smpIdx,:) = KTest(:,1:model.nLabel)*model.projection;
        else
            Embed_Test(smpIdx,:) = KTest*model.projection;
        end
        clear KTest;
    end
    D = EuDist2(Embed_Test,model.ClassCenter,0);
    [dump, idx] = min(D,[],2);
    predictlabel = model.ClassLabel(idx);
    accuracy = 1 - length(find(predictlabel-gnd))/nTest;
end

author: ZJULearning

repo: MatlabFunc


text:
function hRD = plugin_rd_rx
% Reference design definition

% Call the common reference design definition function
hRD = AnalogDevices.adrv9361z7035.common.plugin_rd('ccbox_lvds','Rx');
AnalogDevices.adrv9361z7035.common.add_io(hRD, 'Rx', 'ccbox_lvds');

author: analogdevicesinc


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nickabattista[at]gmail.com
% IB2d Created: May 27th, 2015
% Institution: TCNJ
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs or non-invariant beams*)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%   .
%   .
%   .
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us to add a specific muscle model, please let Nick (nickabattista[at]gmail.co) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: updates the target point positions
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function targets = update_Target_Point_Positions(dt,current_time,targets)

%IDs = targets(:,1);                 % Stores Lag-Pt IDs in col vector
%xPts= targets(:,2);                 % Original x-Values of x-Target Pts.
%yPts= targets(:,3);                 % Original y-Values of y-Target Pts.
%kStiffs = targets(:,4);             % Stores Target Stiffnesses
%N_target = length(targets(:,1));    % Gives total number of target pts!

t1 = 0.01;                          % Period A->B (and B->A)
t2 = 0.02;                          % Period B->C (abd C->B)
period = 2*(t1+t2);                 % Time it takes to move to the right
t = rem(current_time,period);       % Time (moded out)

% Cubic Interpolation Information
p1 = 0.25;
p2 = 0.925;

% a COEFFICIENTS
a0 = 0;
a1 = 0;
a2 = 0;
a3 = 4.324324324324318;

% b COEFFICIENTS
b0 = 0.123456790123457;
b1 = -1.481481481481478;
b2 = 5.925925925925911;
b3 = -3.576910243576897;

% c COEFFICIENTS
c0 = -16.777777777777700;
c1 =  53.333333333333101;
c2 = -53.333333333333101;
c3 = 17.777777777777700;

%
% START THE INTERPOLATING BETWEEN STATES!
%
if t <= t1 % STATE A -> STATE B

    A = read_In_State('State_A.pts');
    B = read_In_State('State_B.pts');

    % Scaling time for appropriate use in interp. function so tTilde\in[0,1]
    tTilde = (t/t1);

    % Evaluate Pieceise Cubic Interpolation Poly
    if tTilde<=p1
        gFUNC = a0 + a1*tTilde + a2*tTilde^2 + a3*tTilde^3;
    elseif tTilde<=p1+p2
        gFUNC = b0 + b1*tTilde + b2*tTilde^2 + b3*tTilde^3;
    else
        gFUNC = c0 + c1*tTilde + c2*tTilde^2 + c3*tTilde^3;
    end

    targets(:,2) = A(:,1) + gFUNC*( B(:,1) - A(:,1) );
    targets(:,3) = A(:,2) + gFUNC*( B(:,2) - A(:,2) );

elseif t <= (t1+t2) % STATE B -> C

    B = read_In_State('State_B.pts');
    C = read_In_State('State_C.pts');

    % Scaling time for appropriate use in interp. function so tTilde\in[0,1]
    tTilde = (t-t1)/t2;

    % Evaluate Pieceise Cubic Interpolation Poly
    if tTilde<=p1
        gFUNC = a0 + a1*tTilde + a2*tTilde^2 + a3*tTilde^3;
    elseif tTilde<=p1+p2
        gFUNC = b0 + b1*tTilde + b2*tTilde^2 + b3*tTilde^3;
    else
        gFUNC = c0 + c1*tTilde + c2*tTilde^2 + c3*tTilde^3;
    end

    targets(:,2) = B(:,1) + gFUNC * ( C(:,1) - B(:,1) );
    targets(:,3) = B(:,2) + gFUNC * ( C(:,2) - B(:,2) );

elseif t <= (t1+2*t2) % STATE C -> B

    B = read_In_State('State_B.pts');
    C = read_In_State('State_C.pts');

    % Scaling time for appropriate use in interp. function so tTilde\in[0,1]
    tTilde = (t-t1-t2)/(t2);

    % Evaluate Pieceise Cubic Interpolation Poly
    if tTilde<=p1
        gFUNC = a0 + a1*tTilde + a2*tTilde^2 + a3*tTilde^3;
    elseif tTilde<=p1+p2
        gFUNC = b0 + b1*tTilde + b2*tTilde^2 + b3*tTilde^3;
    else
        gFUNC = c0 + c1*tTilde + c2*tTilde^2 + c3*tTilde^3;
    end

    targets(:,2) = C(:,1) + gFUNC * ( B(:,1) - C(:,1) );
    targets(:,3) = C(:,2) + gFUNC * ( B(:,2) - C(:,2) );

else % STATE B -> A

    A = read_In_State('State_A.pts');
    B = read_In_State('State_B.pts');

    % Scaling time for appropriate use in interp. function so tTilde\in[0,1]
    tTilde = (t-t1-2*t2)/(2*t2-2*t1);

    % Evaluate Pieceise Cubic Interpolation Poly
    if tTilde<=p1
        gFUNC = a0 + a1*tTilde + a2*tTilde^2 + a3*tTilde^3;
    elseif tTilde<=p1+p2
        gFUNC = b0 + b1*tTilde + b2*tTilde^2 + b3*tTilde^3;
    else
        gFUNC = c0 + c1*tTilde + c2*tTilde^2 + c3*tTilde^3;
    end

    targets(:,2) = B(:,1) + gFUNC * ( A(:,1) - B(:,1) );
    targets(:,3) = B(:,2) + gFUNC * ( A(:,2) - B(:,2) );

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: Reads in the # of vertex pts and all the vertex pts from the
%           .vertex file.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function PTS = read_In_State(struct_name)

filename = struct_name;  %Name of file to read in
fileID = fopen(filename);

% Read in the file, use 'CollectOutput' to gather all similar data together
% and 'CommentStyle' to to end and be able to skip lines in file.
C = textscan(fileID,'%f %f','CollectOutput',1);

fclose(fileID);     %Close the data file.

vertices = C{1};    %Stores all read in data in vertices (N+1,2) array

PTS = vertices(1:end,1:2);

author: nickabattista


text:
function [axons,packing] = func_axonpack_main(numelobj, d_mean, d_var, gap, iter_max)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       AxonPackin : Simulate arrangement of white matter axons
%                     author : Tom Mingasson
%             https://github.com/neuropoly/axonpacking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       CHANGE INPUTS BELOW
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% MAIN INPUTS
N = numelobj;            % number of axons i.e disks to pack
%d_mean = 3;         % theoretical mean of axon diameters in um
%d_var  = 1;         % theoretical variance of axon diameters in um
Delta  = gap;         % gap between the edge of axons in um
%iter_max = 30000;    % number of iteration i.e migrations to perform. Example: iter_max = 30000 ok if N = 1000

% SECONDARY INPUTS
threshold_high = 10;     % no diameter above 'threshold_high'
threshold_low = 0.2;     % no diameter under 'threshold_low'
iter_fvf = iter_max/10;  % to study the packing convergence the disk density i.e Fiber Volume Fraction (FVF) can be computed and displayed every 'iter_fvf' iterations

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     AxonPacking Process
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k=1:length(d_mean)

    % axons features
    axons.N{k}      = N;
    axons.d_mean{k} = d_mean(k);
    axons.d_var{k}  = d_var(k);
    axons.Delta{k}           = Delta(k);
    axons.threshold_high{k} = threshold_high;
    axons.threshold_low{k}  = threshold_low;

    % axon diameters sampling (under a gamma law or lognormal and initialization of positions 'x0' in a square area of length 'side')
    [d, x0, side] = axons_setup(axons,'gamma', k);
    axons.d{k} = d;
    axons.g_ratio{k} = compute_gratio(d);

    % packing process of the axons
    [final_positions, final_overlap, fvf_historic] = process_packing(x0, d, Delta(k), side, iter_max, iter_fvf);

    % store packing results
    % main results
    packing.initial_positions{k}    = reshape(x0,2,length(x0)/2);
    packing.final_positions{k}      = final_positions;
    % secondary results
    packing.final_overlap{k}        = final_overlap;
    packing.FVF_historic{k}         = fvf_historic;
    packing.iter_max{k}             = iter_max;

    % Statistics from the packing
    [FVF, FR, MVF, AVF] = compute_statistics(axons.d{k}, axons.Delta{k}, packing.final_positions{k}, side, axons.g_ratio{k});

    % store stats results
    stats.FVF{k}        = FVF;
    stats.FR{k}         = FR;
    stats.MVF{k}        = MVF;
    stats.AVF{k}        = AVF;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                 Save results in a folder named 'results'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% save_var = num2str(d_var);  save_var(save_var == ' ') = '';
% save_mean = num2str(d_mean); save_mean(save_mean == ' ') = '';
% save_Delta  = num2str(Delta);  save_Delta(save_Delta == ' ') = '';
% save_iter  = num2str(iter_max);
% saveName  = ['Axons', num2str(N), '_Mean', save_mean, '_Var', save_var, '_Delta', save_Delta, '_Iter',save_iter];
%
% mkdir('results')
% cd([pwd,filesep, 'results'])
%
% % save outputs
% save('axons.mat', '-struct', 'axons');
% save('packing.mat', '-struct', 'packing');
% save('stats.mat', '-struct', 'stats');
%
% % save final substrate
% saveas(figure(1000),[saveName,'.png']);
%

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text: %XOR Datafile overload

author: marianux

repo: ecg-kit

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


text:
function importIBSIPhase2DataToCERR(dataDir)
% importIBSIPhase2DataToCERR.m
% Download datasets from: https://github.com/theibsi/data_sets/
%-------------------------------------------------------------------------
% AI 10/06/2020
% AI 11/18/2022 Extnded for phase 2
% AI 12/10/2022 Extended for phase 3

%% Phase-1 : Synthetic phantoms
cerrPath = getCERRPath;
idxV = strfind(getCERRPath,filesep);
savePath = fullfile(cerrPath(1:idxV(end-1)),...
    'Unit_Testing\data_for_cerr_tests\IBSI2_synthetic_phantoms');

subDir = 'ibsi_2_digital_phantom';
synthDataDir = fullfile(dataDir,subDir,'nifti');
dirS = dir(synthDataDir);
dirS(1:2) = [];
for nDataset = 1:length(dirS)

    %Import scan
    synthDatasetName = fullfile(synthDataDir,dirS(nDataset).name,'image');
    synthDatasetDirS = dir([synthDatasetName,filesep,'*.nii.gz']);
    fileName = synthDatasetDirS.name;
    scanName = strtok(fileName,'.');
    filePath = fullfile(synthDatasetName,synthDatasetDirS.name);

    planC = nii2cerr(filePath,scanName,[],0);
    planC = updatePlanFields(planC);
    planC = quality_assure_planC(fileName,planC);

    %Create 'wholeScan' structure
    addMask3M = true(size(getScanArray(1,planC)));
    planC = maskToCERRStructure(addMask3M,0,1,...
        'wholeScan',planC);

    save_planC(planC,[],'PASSED',fullfile(savePath,...
        [scanName,'.mat']));

end

%% Phase-2 CT phantom
subDir = 'ibsi_2_ct_radiomics_phantom';
savePath = fullfile(cerrPath(1:idxV(end-1)),...
    'Unit_Testing\data_for_cerr_tests\IBSI2_CT_phantom');
scanName = subDir;
dcmMaskDir = fullfile(dataDir,subDir,'dicom','mask');
maskDirS = dir([dcmMaskDir,filesep,'*.dcm']);
dcmMaskPath = fullfile(dcmMaskDir,maskDirS.name);
outFileName = [savePath,filesep,'ibsi_2_ct_radiomics_phantom.mat'];

dcmScanDir = fullfile(dataDir,subDir,'dicom','image');
%Copy RT struct to be 'image' dir
copyfile(dcmMaskPath,dcmScanDir);

%planC = dcmdir2planC(dcmDataDir);
planC = importDICOM(dcmScanDir,savePath,true);
movefile([savePath,filesep,'image.mat'],outFileName);
planC = loadPlanC(outFileName,tempdir);
planC = updatePlanFields(planC);
planC = quality_assure_planC(outFileName,planC);

%Rename structure
indexS = planC{end};
planC{indexS.structures}(1).structureName = 'ROI';

save_planC(planC,[],'PASSED',outFileName);

%% IBSI2 phase 3
subDir = 'IBSI2_multimodal_data';
savePath = fullfile(cerrPath(1:idxV(end-1)),...
    'Unit_Testing\data_for_cerr_tests\IBSI2_multimodal_data');
allData = fullfile(dataDir,subDir);
modC = {'CT','MR_T1','PET'};

scanDirS = dir(allData);
scanDirS(1:2) = [];

for nDir = 1:length(scanDirS)

    scanName = scanDirS(nDir).name;
    outFileName = [savePath,filesep,scanName,'.mat'];

    for nMod = 1:length(modC)

        %Scan dir
        dcmScanDir = fullfile(dataDir,subDir,scanName,modC{nMod},'image');

        %Copy RT struct to be 'image' dir
        dcmMaskDir = fullfile(dataDir,subDir,scanName,modC{nMod},'mask');
        maskDirS = dir([dcmMaskDir,filesep,'*.dcm']);
        dcmMaskPath = fullfile(dcmMaskDir,maskDirS.name);
        copyfile(dcmMaskPath,dcmScanDir);

        %Import to CERR
        planD = importDICOM(dcmScanDir,savePath,true);
        movefile([savePath,filesep,'image.mat'],outFileName);
        planD = loadPlanC(outFileName,tempdir);
        planD = updatePlanFields(planD);
        planD = quality_assure_planC(outFileName,planD);
        if nMod==1
            planC = planD;
        else
            planC = planMerge(planC,planD,'all','all','all',outFileName);
        end

        %Rename structure
        indexS = planC{end};
        planC{indexS.structures}(end).structureName = [modC{nMod},'_ROI'];
    end

    save_planC(planC,[],'PASSED',outFileName);
    clear planC
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [y,fs]=readflac(filename,mode)
%READWAV  Read a .FLAC format sound file [Y,FS]=(FILENAME,MODE)
%
% Input Parameters:
%
%	FILENAME gives the name of the file (with optional .WAV extension) or alternatively
%                 can be the FIDX output from a previous call to READWAV
%	MODE		specifies the following (*=default):
%
%    Scaling: 's'    Auto scale to make data peak = +-1
%             'r'    Raw unscaled data (integer values)
%             'q'    Scaled to make 0dBm0 be unity mean square
%             'p' *	 Scaled to make +-1 equal full scale
%             'o'    Scale to bin centre rather than bin edge (e.g. 127 rather than 127.5 for 8 bit values)
%                     (can be combined with n+p,r,s modes)
%             'n'    Scale to negative peak rather than positive peak (e.g. 128.5 rather than 127.5 for 8 bit values)
%                     (can be combined with o+p,r,s modes)
%
% FLAC (Free lossless audio codec) is a compressed audio file format described here:
% http://flac.sourceforge.net/

%      Copyright (C) Mike Brookes 2008
%      Version: $Id: readflac.m 8483 2016-09-13 15:48:05Z 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<2
    mode='p';
else
    mode = [mode(:).' 'p'];
end
dirt=voicebox('dir_temp');
[fnp,fnn,fne]=fileparts(filename);
filetemp=fullfile(dirt,[fnn '.wav']);
doscom=['"' voicebox('flac') '"' ' -d -f -o "' filetemp '" "' filename '"'];
%                     fprintf(1,'Executing: %s\n',doscom);
[doss,dosr]=dos(doscom); % run the program
if doss % test for errors
    error(sprintf('Error running DOS command: %s',doscom));
end
if exist(filetemp)~=2
    error(sprintf('No output file from: %s',doscom));
end
[y,fs]=readwav(filetemp,mode);
doscom=['del /f "' filetemp '"'];
if dos(doscom) % run the program
    error(sprintf('Error running DOS command: %s',doscom));
end

author: jtkim-kaist

repo: VAD

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


text:
hb_setup();

res = rproc.read('scoresroot', ...
  fullfile(hb_path, 'matlab', 'scores', 'scores_norm_small'));

%%

categories = cellfun(@(a) a(1), res.matching.sequence, 'Uni', false);
res.matching.category = categories;
% Compute the matching scores averages
res.matching_g = varfun(@mean, res.matching, ...
  'GroupingVariables', {'descriptor', 'geom_noise', 'category'}, ...
  'InputVariables', 'ap');

%% Separate the normalisation type and split out of descriptor name

res.verification_e = rproc.postproc_norm(res.verification);
res.matching_ge = rproc.postproc_norm(res.matching_g);
res.retrieval_e = rproc.postproc_norm(res.retrieval);

%% Compute average over the norm splits
res.verification_en = varfun(@mean, res.verification_e, ...
  'GroupingVariables', {'descriptor', 'split', 'negs', 'geom_noise', 'method', 'norm_type'}, ...
  'InputVariables', 'pr_ap');

res.matching_gen = varfun(@mean, res.matching_ge, ...
  'GroupingVariables', {'descriptor', 'geom_noise', 'category', 'norm_type'}, ...
  'InputVariables', 'mean_ap');

res.retrieval_en = varfun(@mean, res.retrieval_e, ...
  'GroupingVariables', {'descriptor', 'split', 'geom_noise', 'method', 'norm_type'}, ...
  'InputVariables', 'map');

%% Pick the best normalisation per descriptor
% Preprocess the results

ver_ = res.verification_en;
ver_(~cellfun(@(a) strcmp(a, 'b'), ver_.split), :) = [];
ver_(~cellfun(@(a) strcmp(a, 'imbalanced'), ver_.method), :) = [];

ver_ = varfun(@mean, ver_, ...
  'GroupingVariables', {'descriptor', 'split', 'norm_type'}, ...
  'InputVariables', 'mean_pr_ap');
ver_.split = []; ver_.GroupCount = [];
ver_ = unstack(ver_, 'mean_mean_pr_ap', 'norm_type');
ver_.Properties.RowNames = ver_.descriptor;
ver_.descriptor = [];

match_ = res.matching_gen;
match_.GroupCount = [];
match_ = varfun(@mean, match_, ...
  'GroupingVariables', {'descriptor', 'norm_type'}, ...
  'InputVariables', 'mean_mean_ap');
match_.GroupCount = [];
match_ = unstack(match_, 'mean_mean_mean_ap', 'norm_type');
match_.Properties.RowNames = match_.descriptor;
match_.descriptor = [];

retr_ = res.retrieval_en;
retr_.GroupCount = [];
retr_(cellfun(@(a) strcmp(a, 'keepquery'), retr_.method), :) = [];
retr_ = varfun(@mean, retr_, ...
  'GroupingVariables', {'descriptor', 'norm_type'}, ...
  'InputVariables', 'mean_map');
retr_.GroupCount = [];
retr_ = unstack(retr_, 'mean_mean_map', 'norm_type');
retr_.Properties.RowNames = retr_.descriptor;
retr_.descriptor = [];

%%
ver_best = max(max(ver_{:, :}));
match_best = max(max(match_{:, :}));
retr_best = max(max(retr_{:, :}));

% Weight the scores relatively to the best score
total_score = ver_;
total_score{:,:} = (ver_{:,:} ./ ver_best + match_{:,:} ./ match_best + retr_{:,:} ./ retr_best) ./ 3;

[~, best_norm] = max(total_score{:,:}, [], 2);
best_norm_name = total_score.Properties.VariableNames(best_norm)';

best_norm = struct('descriptor', total_score.Properties.RowNames, ...
  'normstr', best_norm_name);
best_norm = struct2table(best_norm, 'asArray', true);
%%
writetable(best_norm, fullfile(hb_path, 'matlab', 'data', 'best_normalizations.csv'));

author: hpatches

repo: hpatches-benchmark


text:
filename='RVE_Square_Triangle';
ptype = 'MICRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'circleInclusion';
cost={'enforceCh_CCstar_L2','perimeter'};%enforceCh_CCstar_L2
weights=[1, 0.1];
constraint = {'volumeConstraint'};
optimizer = 'MMA'; incrementFactor = 1;%mma CON 1E-4
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.5;
Perimeter_target=3;
optimality_final =1e-4;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;

TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

%Micro
epsilon_isotropy_initial=1e-1;
epsilon_isotropy_final = 1e-3;
selectiveC_Cstar = 'Vfrac05';

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [meshStruct]=hexMeshCylinder(varargin)

% function [meshStruct]=hexMeshCylinder(cylRadius,cylLength,pointSpacing)
% ------------------------------------------------------------------------
%
%
% Kevin Mattheus Moerman
% gibbon.toolbox@gmail.com
%
% 2021/02/26 Created
% 2021/05/08 @Fireedman (Luis Antonio Aguilar) added "end" to function
%------------------------------------------------------------------------

%% Parse input

switch nargin
    case 0
        cylRadius=[];
        cylLength=[];
        pointSpacing=[];
    case 1
        cylRadius=varargin{1};
        cylLength=[];
        pointSpacing=[];
    case 2
        cylRadius=varargin{1};
        cylLength=varargin{2};
        pointSpacing=[];
    case 3
        cylRadius=varargin{1};
        cylLength=varargin{2};
        pointSpacing=varargin{3};
end

%Check for empty variables
if isempty(cylRadius)
    cylRadius=1; %Unit radius
end

if isempty(cylLength)
    cylLength=cylRadius; %Height the same as radius
end

if isempty(pointSpacing)
    pointSpacing=(cylRadius*2*pi)/10; %One tenth of circumference
end

if numel(pointSpacing)==1
    pointSpacing=pointSpacing.*ones(1,2);
end

%%
%Set number of elements for core
numElementsCore=ceil(((cylRadius*2*pi)./pointSpacing(1))/4);

if numElementsCore<1
    numElementsCore=1;
end

%Set number of nodes allong circumference
numElementsHeight=ceil(cylLength./pointSpacing(2));
if numElementsHeight<1
    numElementsHeight=1;
end

%% Create quad mesh for disc

%Raw quad mesh
[Fs,Vs]=discQuadMesh(numElementsCore,cylRadius,0.6);

%Smoothen
Eb=patchBoundary(Fs);
controlParSmooth.n=25;
controlParSmooth.Method='LAP';
controlParSmooth.RigidConstraints=unique(Eb(:));
[Vs]=patchSmooth(Fs,Vs,[],controlParSmooth);

%% Thicken to elements

%Thicken
[E,V,Fp1,Fp2]=patchThick(Fs,Vs,1,cylLength,numElementsHeight);
V(:,3)=V(:,3)-cylLength/2;

%Get element faces
[F,~]=element2patch(E,[],'hex8');

%Find boundary faces
[indFree]=freeBoundaryPatch(F);
Fb=F(indFree,:);

%Assign boundary labels (or colors)
faceBoundaryMarker=zeros(size(Fb,1),1); %Side of cylinder
faceBoundaryMarker(all(ismember(Fb,Fp1),2))=1; %Original (e.g. top)
faceBoundaryMarker(all(ismember(Fb,Fp2),2))=2; %Thickened (e.g. bottom)

%% Collect output

meshStruct.nodes=V;
meshStruct.facesBoundary=Fb;
meshStruct.boundaryMarker=faceBoundaryMarker;
meshStruct.faces=F;
meshStruct.elements=E;
meshStruct.elementMaterialID=ones(size(E,1),1);
meshStruct.faceMaterialID=ones(size(F,1),1);

end

%%
% _*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-2021 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/>.
%%
% _*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 plotPlanes(planes, figureName)
%% Plot 3d planes given plane structure
% See also: Calibration.planeStructFromPoints

%% License
% ACADEMIC OR NON-PROFIT ORGANIZATION NONCOMMERCIAL RESEARCH USE ONLY
% Copyright (c) 2018 Bingyao Huang
% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:

% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.

% If you publish results obtained using this software, please cite our paper.

% 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.

%%

nPlanes = length(planes);

% each plane has a different color
planeColors = hsv(nPlanes);

if(nargin > 1)
    if(isempty(figureName))
        error('Figure name is empty');
    end

    if(~isa(figureName, 'char'))
        error('Figure name should be a string');
    end

    % find figure by name
    figHandle = findobj( 'Type', 'Figure', 'Name', figureName );

    if(length(figHandle) > 1)
        warning(['More than 1 figures named', figureName, 'are found, use the 1st figure']);
    end

    if(length(figHandle) < 1)
        warning(['No figure named', figureName, 'is found, create one instead']);
        figHandle = figure('name', figureName);
    end

    figure(figHandle(1));
end

% plot
for i = 1:nPlanes
    corners = planes(i).corners;
    p = patch(corners(:,1), corners(:,2), corners(:,3), planeColors(i,:));
    p.FaceAlpha = 0.3;
end

end

author: BingyaoHuang

repo: single-shot-pro-cam-calib


text:
function p = pwaproj(h)

% Rewrites a0+c'*x+sum(f_i(x)) to a0+c'*x+max(aix+bi)

if prod(size(h))>1
    error('pwaproj is currently only applicable to scalars');
end

hvars = getvariables(h);
basis = getbase(h);
coefficients = basis(2:end);

Complicating = find(ismember(hvars,yalmip('extvariables')));
if length(Complicating)==length(hvars)
    LinearTerm = 0;
else
    NonComplicating = find(~ismember(hvars,yalmip('extvariables')));
    LinearTerm = recover(hvars(NonComplicating))'*coefficients(NonComplicating)';
end

if isempty(Complicating)
    p = h;
    return
end

basis = getbase(h);
coefficients = basis(2:end);

totalobjective = 0;
totalF = [];
allz = [];
allargs = [];
if ~isempty(Complicating)
    for i = 1:length(Complicating)
        extstruct = yalmip('extstruct',hvars(Complicating(i)));
        [properties{i},F{i},arguments{i}]=model(recover(hvars(Complicating(i))));
        RequiresEpi(i) = 1;
        z = sdpvar(size(arguments{i},1),size(arguments{i},2));
        allz = [allz;z];
        allargs = [allargs;arguments{i}];
        extstruct.arg{1} = z;
        t = feval(properties{i}{1}.name,z);
        f = eval(['@' properties{i}{1}.name]);
        t = f(extstruct.arg{1:end-1});
        [dummy{i},thisF]=model(t);
        totalF = totalF + thisF;
        totalobjective=coefficients(Complicating(i))*t+totalobjective;
    end
end
sdpvar t

totalF = [totalF, totalobjective < t];
P = projection([totalF,t<1234],[allz;t]);

P = full(getbase(P));
b = P(:,1);
A = P(:,2:end);
rmv = find(A(:,end)== 0);
A(rmv,:)=[];
b(rmv) = [];
rmv = find(abs(b-1234)<1e-8);
A(rmv,:)=[];
b(rmv) = [];

p = basis(1)+LinearTerm+max((A(:,1:end-1)*allargs-b)./A(:,end));

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function H = getHankelMatrix(xdat,stackmax)

H = zeros(stackmax,size(xdat,1)-stackmax);
for k=1:stackmax
    H(k,:) = xdat(k:end-stackmax-1+k,1);
end

end

author: eurika-kaiser

repo: SINDY-MPC


text:
function soln = optimTraj(problem)
% soln = optimTraj(problem)
%
% Solves a trajectory optimization problem.
%
% INPUT: "problem" -- struct with fields:
%
%   func -- struct for user-defined functions, passed as function handles
%
%       Input Notes:  square braces show size:  [a,b] = size()
%               t = [1, nTime] = time vector (grid points)
%               x = [nState, nTime] = state vector at each grid point
%               u = [nControl, nTime] = control vector at each grid point
%               t0 = scalar = initial time
%               tF = scalar = final time
%               x0 = [nState, 1] = initial state
%               xF = [nState, 1] = final state
%
%       dx = dynamics(t,x,u)
%               dx = [nState, nTime] = dx/dt = derivative of state wrt time
%
%       dObj = pathObj(t,x,u)
%               dObj = [1, nTime] = integrand from the cost function
%
%       obj = bndObj(t0,x0,tF,xF)
%               obj = scalar = objective function for boundry points
%
%       [c, ceq] = pathCst(t,x,u)
%               c = column vector of inequality constraints  ( c <= 0 )
%               ceq = column vector of equality constraints ( c == 0 )
%
%       [c, ceq] = bndCst(t0,x0,tF,xF)
%               c = column vector of inequality constraints  ( c <= 0 )
%               ceq = column vector of equality constraints ( c == 0 )
%
%       How to pass parameters to your functions:
%           - suppose that your dynamics function is pendulum.m and it
%           accepts a struct of parameters p. When you are setting up the
%           problem, define the struc p in your workspace and then use the
%           following command to pass the function:
%               problem.func.dynamics = @(t,x,u)( pendulum(t,x,u,p) );
%
%       Analytic Gradients:
%           Both the "trapezoid" and "hermiteSimpson" methods in OptimTraj
%       support analytic gradients. Type help "trapezoid" or help
%       "hermiteSimpson" for details on how to pass gradients to the
%       transcription method.
%
%
%   bounds - struct with bounds for the problem:
%
%       .initialTime.low = [scalar]
%       .initialTime.upp = [scalar]
%
%       .finalTime.low = [scalar]
%       .finalTime.upp = [scalar]
%
%       .state.low = [nState,1] = lower bound on the state
%       .state.upp = [nState,1] = lower bound on the state
%
%       .initialState.low = [nState,1]
%       .initialState.upp = [nState,1]
%
%       .finalState.low = [nState,1]
%       .finalState.upp = [nState,1]
%
%       .control.low = [nControl, 1]
%       .control.upp = [nControl, 1]
%
%
%
%   guess - struct with an initial guess at the trajectory
%
%       .time = [1, nGridGuess]
%       .state = [nState, nGridGuess]
%       .control = [nControl, nGridGuess]
%
%   options = options for the transcription algorithm (this function)
%
%       .nlpOpt = options to pass through to fmincon
%
%       .method = string to pick which method is used for transcription
%           'trapezoid'
%           'hermiteSimpson'
%           'chebyshev'
%           'rungeKutta'
%           'gpops'     ( Must have GPOPS2 installed )
%
%       .[method] = a struct to pass method-specific parameters. For
%       example, to pass the number of grid-points to the trapezoid method,
%       create a field .trapezoid.nGrid = [number of grid-points].
%
%       .verbose = integer
%           0 = no display
%           1 = default
%           2 = display warnings, overrides fmincon display setting
%           3 = debug
%
%       .defaultAccuracy = {'low','medium','high'}
%           Sets the default options for each transcription method
%
%       * if options is a struct array, the optimTraj will run the optimization
%       by running options(1) and then using the result to initialize a new
%       solve with options(2) and so on, until it runs options (end). This
%       allows for successive grid and tolerance opdates.
%
%
%
%
%
% OUTPUT: "soln"  --  struct with fields:
%
%   .grid = trajectory at the grid-points used by the transcription method
%       .time = [1, nTime]
%       .state = [nState, nTime]
%       .control = [nControl, nTime];
%
%   .interp = functions for interpolating state and control for arbitrary
%       times long the trajectory. The interpolation method will match the
%       underlying transcription method. This is particularily important
%       for high-order methods, where linear interpolation between the
%       transcription grid-points will lead to large errors. If the
%       requested time is not on the trajectory, the interpolation will
%       return NaN.
%
%       .state = @(t) = given time, return state
%           In: t = [1,n] vector of time
%           Out: x = [nState,n] state vector at each point in time
%
%       .control = @(t) = given time, return control
%           In: t = [1,n] vector of time
%           Out: u = [nControl,n] state vector at each point in time
%
%   .info = information about the optimization run
%       .nlpTime = time (seconds) spent in fmincon
%       .exitFlag = fmincon exit flag
%       .objVal = value of the objective function
%       .[all fields in the fmincon "output" struct]
%
%   .problem = the problem as it was passed to the low-level transcription,
%       including the all default values that were used
%
%
%   * If problem.options was a struct array, then soln will also be a
%   struct array, with soln(1) being the solution on the first iteration,
%   and soln(end) being the final solution.
%
%   * Both the trapezoid and hermiteSimpson method support additional
%   features. The first is that they can use analytic gradients, if
%   provided by the user. ">> help trapezoid" or ">> help hermiteSimpson"
%   for more information. These methods additionally provide error
%   estimates in two forms. The continuous collocation constraint error is
%   provided as an additional function handle (soln.interp.collCst), and
%   the absolute local error estimate for each segment is provided in
%   soln.info.error.
%

problem = inputValidation(problem);   %Check inputs
problem = getDefaultOptions(problem); % Complete options struct

% Loop over the options struct to solve the problem
nIter = length(problem.options);
soln(nIter) = struct('grid',[],'interp',[],'info',[],'problem',[]);
P = problem;  %Temp variable for passing on each iteration
for iter=1:nIter
    P.options = problem.options(iter);

    if P.options.verbose > 0    %then print out iteration count:
        disp('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~');
        disp(['Running OptimTraj, iteration ' num2str(iter)]);
    end

    if iter > 1  %Use previous soln as new guess
        P.guess = soln(iter-1).grid;
    end

    %%%% This is the key part: call the underlying transcription method:
    switch P.options.method
        case 'trapezoid'
            soln(iter) = trapezoid(P);
        case 'hermiteSimpson'
            soln(iter) = hermiteSimpson(P);
        case 'chebyshev'
            soln(iter) = chebyshev(P);
        case 'multiCheb'
            soln(iter) = multiCheb(P);
        case 'rungeKutta'
            soln(iter) = rungeKutta(P);
        case 'gpops'
            soln(iter) = gpopsWrapper(P);
        otherwise
            error('Invalid method. Type: ''help optimTraj'' for a valid list.');
    end

end
end

author: MatthewPeterKelly

repo: OptimTraj


text:
function [returns,interval] = realized_return_filter(price,time,timeType,samplingType,samplingInterval,subsamples)
% THESE COMMENTS ARE WRONG!!!
% Computes the Quantile Realized Variance of Christensen, Oomen and Podolskij (2008),
% the MinRV and MedRV estimator of Andersen, Dobrev and Shaumberg and the
% general symmetrized version suggested by Sheppard in a discussion of COP
%
% USAGE:
%   [RETURN,INTERVAL] =  realized_return_filter(PRICE,TIME,TIMETYPE,SAMPLINGTYPE,SAMPLINGINTERVAL,SUBSAMPLES)
%
% INPUTS:
%   PRICE            - m by 1 vector of high frequency prices
%   TIME             - m by 1 vector of times where TIME(i) corresponds to PRICE(i)
%   TIMETYPE         - String describing the way times are measured
%                       'wall'    24-hour clock of the form HHMMSS, e.g. 101543 or 153217
%                       'seconds' Time measured in seconds past midnight on the first day.
%                       'unit'  Unit normalized date format, e.g. .1, .234, .9
%                         Unit normalized times are more general than the other types and can be
%                         applied to data from more than one calendar day
%   SAMPLINGTYPE     - String describing the type of sampling to use when
%                        filtering PRICE
%                        'CalendarTime' - Sample in calendar time using observations separated by
%                          SAMPLINGINTERVAL seconds
%                        'CalendarUniform' - Sample in calendar time using SAMPLINGINTERVAL
%                          observations spread uniformly between TIME(1) and TIME(m)
%                        'BusinessTime' - Sample in business (tick) time using observation separated
%                          by SAMPLINGINTERVALticks
%                        'BusinessUniform' - Sample in business (tick) time using observations
%                          uniformly spaced in business time.
%                        'Fixed' - Sample at specific points in time. When using fixed,
%                          SAMPLINGINTERVAL must be a n by 1 vector of times with the same TIMETYPE
%                          as TIME (i.e. seconds if TIME is in seconds)
%   SAMPLINGINTERVAL - Scalar integer or n by 1 vector whose meaning depends on the selected SAMPLINGTYPE
%   SUBSAMPLES       - [OPTIONAL] Integer value containing the number of subsamples to use when
%                        computing RRSS. It must be the case that SUBSAMPLE<SAMPLESPERBIN. Setting
%                        SUBSAMPLES>0 will produce an estimator based on jittering the starting
%                        point over floor(linspace(0,(SUBSAMPLES)/(SUBSAMPLES+1),SUBSAMPLES+1)*SAMPLESPERBIN)+1
%                        For example, if SUBSAMPLES = 3 and SAMPLESPERBIN = 20, then the starting
%                        points for the four RR estimators will be observations 1, 6, 11 and 16.  If
%                        omitted SUBSAMPLE = 0 and RRSS=RR.
%
% OUTPUTS:
%   RETURNS          - Quantile realized variance vector (k by 1) corresponding to the values in QUANTILES
%   INTERVAL         - Subsample based version of RQ.  RQSS = RQ if SUBSAMPLES = 0
%
% COMMENTS:
%
% EXAMPLES:
%
%

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 5/1/2008

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% InputChecking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<5 || nargin>6
    error('Seven or eight inputs required.')
end

if size(price,2)>size(price,1)
    price=price';
end
if size(price,2)>1
    error('PRICE must be a m by 1 vector.')
end
if size(time,2)>size(time,1)
    time=time';
end
if any(diff(time)<0)
    error('TIME must be sorted and increasing')
end
if size(time,2)>1 || length(time)~=length(price)
    error('TIME must be a m by 1 vector.')
end
% Inserted to protect against inputing integer times
time = double(time);

timeType=lower(timeType);
if ~ismember(timeType,{'wall','seconds','unit'})
    error('TIMETYPE must be one of ''wall'', ''seconds'' or ''unit''.');
end
samplingType=lower(samplingType);
if ~ismember(samplingType,{'calendartime','calendaruniform','businesstime','businessuniform','fixed'})
    error('SAMPLINGTYPE must be one of ''CalendarTime'', ''CalendarUniform'', ''BusinessTime'', ''BusinessUniform'' or ''Fixed''.');
end

m=size(price,1);
t0=time(1);
tT=time(m);
if ismember(samplingType,{'calendartime','calendaruniform','businesstime','businessuniform'})
    % Must be a scalar integer
    if ~isscalar(samplingInterval) || floor(samplingInterval)~=samplingInterval || samplingInterval<1
        error('SAMPLINGINTERVAL must be a positive integer for the SAMPLINGTYPE selected.')
    end
else
    if size(samplingInterval,2)>size(samplingInterval,1)
        samplingInterval=samplingInterval';
    end
    if ~(any(samplingInterval>=t0) && any(samplingInterval<=tT))
        error('At least one sampling interval must be between min(TIME) and max(TIME) when using ''Fixed'' as SAMPLINGTYPE.')
    end
    if any(diff(samplingInterval)<=0)
        error('When using ''Fixed'' as SAMPLINGTYPE the vector of sampling times in SAMPLINGINTERVAL must be sorted and strictly increasing.')
    end
end

% SUBSAMPLES
if nargin==6 && ~isempty(subsamples)
    if ~isscalar(subsamples) || subsamples<0 || floor(subsamples)~=subsamples
        error('SUBSAMPLES must be a non-negative integer.')
    end
else
    subsamples=0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Convert if wall to compute the interval
if strcmp(timeType,'wall')
    time=wall2seconds(time);
    timeType = 'seconds';
end

% Initialize the return cell arrays and the interval vector
returns=cell(subsamples+1,1);
interval=zeros(subsamples+1,1);
% 1. Filter the price
logPrice = log(price);
[filteredLogPrice,filteredTimes] = realized_price_filter(logPrice ,time,timeType,samplingType,samplingInterval);
returns{1} = diff(filteredLogPrice);
interval(1) = 1;

if subsamples>0
    % Initialize the place holder for rvSS and rvAdjustment
    % Get the number of filtered times
    n = size(filteredTimes,1);
    % Compute the length of time used in computing the realized variance to
    % bias adjust the subsample RVs
    baseDifference = filteredTimes(n) - filteredTimes(1);
    switch lower(samplingType)
        case {'calendartime','calendaruniform','fixed'}
            % Uniform sampling in CT so the gap is constant
            gap=diff(filteredTimes);
            % The stepsize is constant
            step = gap /(subsamples+1);
            if ismember(samplingType,{'calendartime','calendaruniform'})
                step = [step;mean(step)];
            else
                step = [step;step(length(step))];
            end
            for i=1:subsamples
                thisSampleTime = filteredTimes;
                % Need to compute times for this subsample
                thisSampleTime = thisSampleTime(1:n)+i*step;
                % Use these times to filter the price
                filteredLogPrice = realized_price_filter(logPrice,time,'unit','fixed',thisSampleTime);
                % Compute returns
                returns{i+1} = diff(filteredLogPrice);
                % Compute the ith subsample RV
                % The adjustment depends on the amount of time used in this
                % subsample relative to the full sample RV
                interval(i+1) = (thisSampleTime(n-1)-thisSampleTime(1))/baseDifference;
            end
        case {'businesstime', 'businessuniform'}
            % Uniform sampling in business time
            if strcmp(samplingType,'businesstime')
                % If BT then samplingInterval contains the number of ticks to skip
                originalIndices = 1:samplingInterval:m;
                if ~ismember(m,originalIndices)
                    originalIndices = [originalIndices m];
                end
                gap=samplingInterval;
            else
                % If businessuniform then gap will be average
                originalIndices = floor(linspace(1,m,samplingInterval));
                if ~ismember(m,originalIndices)
                    originalIndices = [originalIndices m];
                end
                gap = m/samplingInterval;
            end
            % The step size is constant
            step = gap/(subsamples+1);
            for i=1:subsamples
                % The indices indicate which points to use since the
                % sampling is linear in ticks
                indices = floor(originalIndices + step*i);
                indices(indices>m) = m;
                returns{i+1} = diff(logPrice(indices));
                interval(i+1) = (time(max(indices))-time(min(indices)))/baseDifference;
            end
    end
end

author: bashtage

repo: mfe-toolbox


text:
function value = conex2_condition ( alpha )

%*****************************************************************************80
%
%% CONEX2_CONDITION returns the L1 condition of the CONEX1 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 February 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, the scalar defining A.
%    A common value is 100.0.
%
%    Output, real VALUE, the L1 condition.
%
  c1 = 1.0;
  c2 = abs ( 1.0 - 1.0 / alpha^2 ) + 1.0 / abs ( alpha );
  c3 = 3.0 + 1.0 / abs ( alpha );
  a_norm = max ( c1, max ( c2, c3 ) );
  c1 = 1.0;
  c2 = abs ( ( 1.0 - alpha * alpha ) / alpha ) + abs ( alpha );
  c3 = abs ( ( 1.0 + alpha * alpha ) / alpha ^ 2 ) + 2.0;
  b_norm = max ( c1, max ( c2, c3 ) );
  value = a_norm * b_norm;

  return
end

author: johannesgerer

repo: jburkardt-m