text:
% This script estimates the prior of a hedge fund return and processes extreme views on CVaR
% according to Entropy Pooling
% This script complements the article
%	"Fully Flexible Extreme Views"
%	by A. Meucci, D. Ardia, S. Keel
%	available at www.ssrn.com
% The most recent version of this code is available at
% MATLAB Central - File Exchange

% IMPORTANT - This script is about the methodology, not the input data, which has been modified

clear; clc; close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fetch data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load pseudodata;

xi = 100 * cell2mat(data(:, 2));
n = length(xi);

% bandwith
bw = kernelbw(xi);

% weights
lambda = log(2) / (n / 2);
wi = exp(-lambda * (n - (n:-1:1)'));
wi = flipud(wi) / sum(wi);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Prior market model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% kernel density
market = [];
market.mu  = mean(xi);
market.pdf = @(x) kernelpdf(x, xi, bw, wi);
market.cdf = @(x) kernelcdf(x, xi, bw, wi);
market.inv = @(x) kernelinv(x, xi, bw, wi);
market.VaR95 = market.inv(0.05);
market.CVaR95 = quadgk(@(x) (x .* market.pdf(x)), -100, market.VaR95) / 0.05;

% numerical (Gauss-Hermite grid) prior
load ghq1000.mat % load mesh of GH zeros
tmp = (ghqx-min(ghqx))/(max(ghqx)-min(ghqx)); % rescale GH zeros so they belong to [0,1]
epsilon = 1e-10;
Lower = market.inv(epsilon);
Upper = market.inv(1-epsilon);
X  = Lower + tmp * (Upper-Lower); % rescale mesh

p = integrateSubIntervals(X, market.cdf);
p = normalizeProb(p);
J = length(X);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Entropy posterior from extreme view on mean and CVaR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
view = [];
view.mu   = mean(xi) - 1.0;
view.CVaR95 = market.CVaR95 - 1.0;

% Netwton Raphson
s = [];
idx = cumsum(p) <= 0.05 ;
s(1) = sum(idx);
[p_, KLdiv] = optimizeEntropy(p, idx', 0.05, [ones(1, J); X'; (idx .* X)'], [1; view.mu; 0.05 * view.CVaR95]);
%[p_, KLdiv] = optimizeEntropy(p, [idx'; (idx .* X)'], [0.05; 0.05 * view.CVaR95], [ones(1, J); X'], [1; view.mu]);

doStop = 0;
i = 1;
while ~doStop
   i = i + 1;

   idx = [ones(1, s(i-1)), zeros(1, J - s(i-1))]';
   [dummy, KLdiv_s]  = optimizeEntropy(p, idx', 0.05, [ones(1, J); X'; (idx .* X)'], [1; view.mu; 0.05 * view.CVaR95]);
   %[dummy, KLdiv_s]  = optimizeEntropy(p, [idx'; (idx .* X)'], [0.05; 0.05 * view.CVaR95], [ones(1, J); X'], [1; view.mu]);

   idx = [ones(1, s(i-1) + 1), zeros(1, J - s(i-1) - 1)]';
   [dummy, KLdiv_s1]  = optimizeEntropy(p, idx', 0.05, [ones(1, J); X'; (idx .* X)'], [1; view.mu; 0.05 * view.CVaR95]);
   %[dummy, KLdiv_s1] = optimizeEntropy(p, [idx'; (idx .* X)'], [0.05; 0.05 * view.CVaR95], [ones(1, J); X'], [1; view.mu]);

   idx = [ones(1, s(i-1) + 2), zeros(1, J - s(i-1) - 2)]';
   [dummy, KLdiv_s2]  = optimizeEntropy(p, idx', 0.05, [ones(1, J); X'; (idx .* X)'], [1; view.mu; 0.05 * view.CVaR95]);
   %[dummy, KLdiv_s2] = optimizeEntropy(p, [idx'; (idx .* X)'], [0.05; 0.05 * view.CVaR95], [ones(1, J); X'], [1; view.mu]);

   % first difference
   DE  = KLdiv_s1 - KLdiv_s;
   % second difference
   D2E = KLdiv_s2 - 2 * KLdiv_s1 + KLdiv_s;
   % optimal s
   s = [s; round( s(i-1) - (DE / D2E) )];

   tmp = [];
   idx  = [ones(1, s(i)), zeros(1, J - s(i))]';
   [tmp.p_, tmp.KLdiv]  = optimizeEntropy(p, idx', 0.05, [ones(1, J); X'; (idx .* X)'], [1; view.mu; 0.05 * view.CVaR95]);
   %[tmp.p_, tmp.KLdiv] = optimizeEntropy(p, [idx'; (idx .* X)'], [0.05; 0.05 * view.CVaR95], [ones(1, J); X'], [1; view.mu]);
   p_ = [p_, tmp.p_];
   KLdiv = [KLdiv; tmp.KLdiv]; %#ok<*AGROW>

   % if change in KLdiv less than one percent, stop
   if abs((KLdiv(i) - KLdiv(i-1)) / KLdiv(i-1)) < 0.01
      doStop = 1;
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Display results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;
plot(X, p, 'b', 'LineWidth', 1.5);
hold on;
plot(X, p_(:, end), 'r', 'LineWidth', 1.5);
hold on; plot(market.mu, 0.0, 'o', 'Color', 'Blue', 'MarkerSize', 8, 'MarkerFaceColor', 'Blue');
hold on; plot(market.CVaR95, 0.0, '^', 'Color', 'Blue', 'MarkerSize', 8, 'MarkerFaceColor', 'Blue');
hold on; plot(view.mu, 0.0, 'o', 'Color', 'Red', 'MarkerSize', 8, 'MarkerFaceColor', 'Red');
hold on; plot(view.CVaR95, 0.0, '^', 'Color', 'Red', 'MarkerSize', 8, 'MarkerFaceColor', 'Red');
grid on;
xlabel('Returns [%]');
h = legend('Prior', 'Posterior');
set(h, 'FontSize', 8);
set(gca, 'YTickLabel', '');


text:
function [BB vx] = spm_get_bbox(V, thr, premul)
% Compute volume's bounding box, for full field of view or object bounds
% FORMAT [BB vx] = spm_get_bbox(V, thr)
% V   - mapped image volume(s) (from spm_vol) or filename (empty for GUI)
% thr - threshold, such that BB contains voxels with intensities > thr
%       or strings 'nz', 'nn', fv', for non-zero, non-NaN, or field of view
%       where 'fv' (the default) uses only the image's header information.
%
% BB  - a [2 x 3] array of the min and max X, Y, and Z coordinates {mm},
%       i.e. BB = [minX minY minZ; maxX maxY maxZ].
% vx  - a [1 x 3] vector of voxel dimensions {mm}.
%__________________________________________________________________________
% Copyright (C) 2011 Wellcome Trust Centre for Neuroimaging

% Ged Ridgway
% $Id: spm_get_bbox.m 4205 2011-02-21 15:39:08Z guillaume $

% Undocumented expert options:
% V           - can be a 4D @nifti object (but not 5D), image-based BBs
%               will be computed using "all" along the 4th dimension.
% thr = 'old' - reproduce spm_write_sn/bbvox_from_V (and elsewhere)
% premul      - a matrix that premultiplies V.mat, as used in spm_orthviews

%-Get an SPM volume structure
%--------------------------------------------------------------------------
if nargin < 1 || isempty(V)
    [V, sts] = spm_select(1, 'image', 'Select Image');
    if ~sts, error('Must select an image'), end
end
if ischar(V), V = spm_vol(V); end

%-Get volume structure from @nifti object if given
%--------------------------------------------------------------------------
if isa(V, 'nifti')
    V = spm_vol(V.dat.fname); % (potentially a struct array of volumes)
end

%-Compute voxel dimensions (for compatibility with bbvox_from_V)
%--------------------------------------------------------------------------
P = spm_imatrix(V(1).mat);
vx = P(7:9);
% the above agrees with sqrt(sum(V.mat(1:3,1:3).^2)) for simple rotations,
% and seems more appropriate if there are reflections and/or skews.
% Note that spm_imatrix(diag([-1 1 1 1])) is [-1 1 1] as expected.

%-Compute bounding box
%--------------------------------------------------------------------------
if nargin < 2 || isempty(thr) || strcmpi(thr, 'fv')
    % overall field-of-view bounding box from header information
    d = V(1).dim;
    corners = [
        1    1    1    1
        1    1    d(3) 1
        1    d(2) 1    1
        1    d(2) d(3) 1
        d(1) 1    1    1
        d(1) 1    d(3) 1
        d(1) d(2) 1    1
        d(1) d(2) d(3) 1
        ]';
    XYZ = V(1).mat(1:3, :) * corners;
elseif strcmpi(thr, 'old')
    % code from spm_write_sn/bbvox_from_V (and other places)
    % NB: main difference is that vx(1)<0 gives descending BB(:,1),
    % shouldn't be used if V.mat contains rotations or skews.
    o  = V(1).mat\[0 0 0 1]';
    o  = o(1:3)';
    BB = [-vx.*(o-1) ; vx.*(V(1).dim(1:3)-o)];
    if exist('premul', 'var')
        warning('spm_get_bbox:old_and_premul', 'old method ignores premul')
    end
else
    % image-based bounding box using voxel intensities
    [img XYZ] = spm_read_vols(V);
    if ischar(thr)
        switch lower(thr)
            case 'nn'  % non-NaN, though include +/- Inf in computation
                img = ~isnan(img);
            case 'nz'  % special case of non-zero (rather than > 0)
                img = ~isnan(img) & img ~= 0;
            otherwise
                error('Unknown threshold type %s', thr)
        end
    else
        % treat thr as numeric threshold
        img = img > thr;
    end
    if ndims(img) == 4
        img = all(img, 4);
    end
    if nnz(img) == 0
        warning('spm_get_bbox:nothing', ...
            'Threshold leaves no voxels, returning full field of view');
    else
        XYZ = XYZ(:, img); % keep only coords that satisfy condition
    end
end

if ~exist('BB', 'var') % exists already if 'old' case chosen above
    if exist('premul', 'var')
        XYZ = premul(1:3, :) * [XYZ; ones(1, size(XYZ, 2))];
    end
    BB = [
        min(XYZ, [], 2)'
        max(XYZ, [], 2)'
        ];
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%  Figure 3.3      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
% script to generate Fig. 3.3
%% fig3_03.m Example 3.4
% Frequency response
clf;
k=1;
num=1;                      % form numerator
den=[1 k];                  % form denominator
sys=tf(num,den);            % form system
% frequency range
w=logspace(-2,2);
[mag,phase]=bode(sys,w);    % compute frequency response
% plot frequency response
subplot(2,1,1); loglog(w,mag(:));xlabel('\omega (rad/sec)');
ylabel('M');title('Figure 3.3: Magnitude, phase');grid;
subplot(2,1,2); semilogx(w,phase(:));xlabel('\omega (rad/sec)');
ylabel('\phi (deg)')
grid;
% Bode grid
bodegrid

author: Sable

repo: mcbench-benchmarks


text:
function value = i1mach ( i )

%*****************************************************************************80
%
%% I1MACH returns integer machine constants.
%
%  Discussion:
%
%    Input/output unit numbers.
%
%      I1MACH(1) = the standard input unit.
%      I1MACH(2) = the standard output unit.
%      I1MACH(3) = the standard punch unit.
%      I1MACH(4) = the standard error message unit.
%
%    Words.
%
%      I1MACH(5) = the number of bits per integer storage unit.
%      I1MACH(6) = the number of characters per integer storage unit.
%
%    Integers.
%
%    Assume integers are represented in the S digit base A form:
%
%      Sign * (X(S-1)*A^(S-1) + ... + X(1)*A + X(0))
%
%    where 0 <= X(1:S-1) < A.
%
%      I1MACH(7) = A, the base.
%      I1MACH(8) = S, the number of base A digits.
%      I1MACH(9) = A^S-1, the largest integer.
%
%    Floating point numbers
%
%    Assume floating point numbers are represented in the T digit
%    base B form:
%
%      Sign * (B^E) * ((X(1)/B) + ... + (X(T)/B^T) )
%
%    where 0 <= X(I) < B for I=1 to T, 0 < X(1) and EMIN <= E <= EMAX.
%
%      I1MACH(10) = B, the base.
%
%    Single precision
%
%      I1MACH(11) = T, the number of base B digits.
%      I1MACH(12) = EMIN, the smallest exponent E.
%      I1MACH(13) = EMAX, the largest exponent E.
%
%    Double precision
%
%      I1MACH(14) = T, the number of base B digits.
%      I1MACH(15) = EMIN, the smallest exponent E.
%      I1MACH(16) = EMAX, the largest exponent E.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 April 2007
%
%  Author:
%
%    Original FORTRAN77 version by Phyllis Fox, Andrew Hall, Norman Schryer
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Phyllis Fox, Andrew Hall, Norman Schryer,
%    Algorithm 528,
%    Framework for a Portable Library,
%    ACM Transactions on Mathematical Software,
%    Volume 4, Number 2, June 1978, page 176-188.
%
%  Parameters:
%
%    Input, integer I, chooses the parameter to be returned.
%    1 <= I <= 16.
%
%    Output, integer VALUE, the value of the chosen parameter.
%
  if ( i < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I1MACH - Fatal error!\n' );
    fprintf ( 1, '  The input argument I is out of bounds.\n' );
    fprintf ( 1, '  Legal values satisfy 1 <= I <= 16.\n' );
    fprintf ( 1, '  I = %d\n', i );
    value = 0;
    error ( 'I1MACH - Fatal error!' );
  elseif ( i == 1 )
    value = 5;
  elseif ( i == 2 )
    value = 6;
  elseif ( i == 3 )
    value = 7;
  elseif ( i == 4 )
    value = 6;
  elseif ( i == 5 )
    value = 32;
  elseif ( i == 6 )
    value = 4;
  elseif ( i == 7 )
    value = 2;
  elseif ( i == 8 )
    value = 31;
  elseif ( i == 9 )
    value = 2147483647;
  elseif ( i == 10 )
    value = 2;
  elseif ( i == 11 )
    value = 24;
  elseif ( i == 12 )
    value = -125;
  elseif ( i == 13 )
    value = 128;
  elseif ( i == 14 )
    value = 53;
  elseif ( i == 15 )
    value = -1021;
  elseif ( i == 16 )
    value = 1024;
  elseif ( 16 < i )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I1MACH - Fatal error!\n' );
    fprintf ( 1, '  The input argument I is out of bounds.\n' );
    fprintf ( 1, '  Legal values satisfy 1 <= I <= 16.\n' );
    fprintf ( 1, '  I = %d\n', i );
    value = 0;
    error ( 'I1MACH - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [u,v,un,vn]=CalculateVelocityFromVortex2D(s,xc,yc,x,y)

%--------------------------------------------------------------------------
%CalculateVelocityFromVortex2D
%Version 1.00
%Created by Stepen
%Created 8 January 2011
%--------------------------------------------------------------------------
%CalculateVelocityFromVortex2D calculates the flow velocity influenced by
%several given vortex flow at several given locations in 2D domain.
%CalculateVelocityFromVortex2D assumes that the vortexs that influence the
%flow are point vortexes.
%--------------------------------------------------------------------------
%Syntax:
%[u,v,un,vn]=CalculateVelocityFromvortex2D(s,xc,yc,x,y)
%Input argument:
%- s (m x 1 num) specifies the strength of all point vortexes.
%- xc (m x 1 num) specifies the x axis location of all point vortexes.
%- yc (m x 1 num) specifies the y axis location of all point vortexes.
%- x (p x q num) specifies the x axis location of the corresponding
%  velocity vector.
%- y (p x q num) specifies the y axis location of the corresponding
%  velocity vector.
%Output argument:
%- u (m x n) specifies the horizontal velocity component at the
%  corresponding location.
%- u (m x n) specifies the vertical velocity component at the corresponding
%  location.
%- un (m x n num) specifies the normalized horizontal velocity component of
%  the corresponding location. un output was generated for plotting
%  purpose. Use this output instead of u for velocity quiver plotting.
%- vn (m x n num) specifies the normalized vertical velocity component of
%  the corresponding location. vn output was generated for plotting
%  purpose. Use this output instead of v for velocity quiver plotting.
%--------------------------------------------------------------------------

%CodeStart-----------------------------------------------------------------
%Checking input point vortexes
    s_size=size(s);
    xc_size=size(xc);
    yc_size=size(yc);
    if (s_size(2)~=1)||(xc_size(2)~=1)||(yc_size(2)~=1)
        error('Input point vortexs must be a i x 1 array!')
    end
    if (numel(s)~=numel(xc))||(numel(s)~=numel(yc))
        error('Input s, xc, and yc do not have the same size!')
    end
    vortexcount=numel(s);
%Checking input location
    x_size=size(x);
    y_size=size(y);
    if (x_size(1)~=y_size(1))||(x_size(2)~=y_size(2))
        error('Input x and y do not have the same size!')
    end
    rowcount=x_size(1);
    colcount=x_size(2);
%Preallocating array for speed
    u=zeros(rowcount,colcount);
    v=zeros(rowcount,colcount);
    un=zeros(rowcount,colcount);
    vn=zeros(rowcount,colcount);
%Calculating influenced velocity
    for i=1:1:rowcount
        for j=1:1:colcount
            for k=1:1:vortexcount
                if (x(i,j)==xc(k))&&(y(i,j)==yc(k))
                    u(i,j)=u(i,j);
                    v(i,j)=v(i,j);
                else
                    r=norm([(x(i,j)-xc(k)),(y(i,j)-yc(k))]);
                    vtheta=-s(k)/(2*pi()*r);
                    u(i,j)=u(i,j)+(vtheta*((y(i,j)-yc(k))/r));
                    v(i,j)=v(i,j)+(vtheta*((xc(k)-x(i,j))/r));
                end
            end
        end
    end
%Calculating normalized vector
    for i=1:1:rowcount
        for j=1:1:colcount
            V=norm([u(i,j),v(i,j)]);
            un(i,j)=u(i,j)/V;
            vn(i,j)=v(i,j)/V;
        end
    end
%CodeEnd-------------------------------------------------------------------

end


text:
%kkinverse 'Inverse Filtering (Restoration) in Fourier Frequency Domain'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kinverse.pane file
%
% Parameters:
% InputFile: i1 'Source data object (FFT)', required: 'source data object (FFT)'
% InputFile: i2 ' PSF data object (FFT)', required: 'PSF data object (FFT)'
% OutputFile: o 'Filtered data object (FFT)', required: 'filtered data object (FFT)'
%
% Example: o = kkinverse({i1, i2}, {'i1','';'i2','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kinverse - Inverse Filtering (Restoration) in Fourier Frequency Domain
%
%  DESCRIPTION
% .I kinverse
% performs restoration of data in the spatial domain by inverse filtering in
% the frequency domain.  The idea is to deconvolve the blurring function
% from the blurred data.  A few problems exist in performing this operation,
% as described in great detail in any good image processing book.
%
% The major difficulty lies in how small the denominator is allowed to be
% before you call it garbage.  The value chosen is supplied as the
% -t thresh value.  Ideally, this will be just above the noise floor in the
% image, but experimentation is the best guide here.
% .I kinverse
% places zeros in the output data for those data points where the PSF data object
% magnitude is below the threshold. In this respect, kinverse is actually an
% implementation of a "pseudoinverse\fR filter.
%
% Note that the input data objects will be treated simply as a field of complex
% data, presumably obtained by application of an FFT, although the data could
% be generated by many other means. The output data will be either KCOMPLEX if
% the input data is any type other than KDCOMPLEX. Otherwise, the output will be
% of type KDCOMPLEX. Data objects with maps attached are not supported.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
% R.C. Gonzalez, P. Wintz, "Digital Image Processing", Addison-Wesley 1987,
% ISBN 0-201-11026-1, pp. 221-224.
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkinverse(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,..] = kkinverse(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'i2', '__input';'o', '__output'};
maxval={0,0,0};
minval={0,0,0};
istoggle=[0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; 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 'kinverse"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function select3dtool(arg)
%SELECT3DTOOL A simple tool for interactively obtaining 3-D coordinates
%
% SELECT3DTOOL(FIG) Specify figure handle
%
% Example:
%   surf(peaks);
%   select3dtool;
%   % click on surface

if nargin<1
   arg = gcf;
end

if ~ishandle(arg)
   feval(arg);
   return;
end

%% initialize gui %%
fig = arg;
figure(fig);

uistate = uiclearmode(fig);
[tool, htext] = createUI;
hmarker1 = line('marker','o','markersize',10,'markerfacecolor','k','erasemode','xor','visible','off');
hmarker2 = line('marker','o','markersize',10,'markerfacecolor','r','erasemode','xor','visible','off');

state.uistate = uistate;
state.text = htext;
state.tool = tool;
state.fig = fig;
state.marker1 = hmarker1;
state.marker2 = hmarker2;
setappdata(fig,'select3dtool',state);
setappdata(state.tool,'select3dhost',fig);

set(fig,'windowbuttondownfcn','select3dtool(''click'')');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function off

state = getappdata(gcbf,'select3dtool');

if ~isempty(state)
    delete(state.tool);
end

fig = getappdata(gcbf,'select3dhost');

if ~isempty(fig) && ishandle(fig)
    state = getappdata(fig,'select3dtool');
    uirestore(state.uistate);
    delete(state.marker1);
    delete(state.marker2);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function click

[p, v, vi] = select3d;
state = getappdata(gcbf,'select3dtool');

if ~ishandle(state.text)
    state.text = createUI;
end

if ~ishandle(state.marker1)
    state.marker1 = [];
end

if ~ishandle(state.marker2)
    state.marker2 = [];
end

setappdata(state.fig,'select3dtool',state);

if isempty(v)
    v = [nan nan nan];
    vi = nan;
    set(state.marker2,'visible','off');
else
    set(state.marker2,'visible','on','xdata',v(1),'ydata',v(2),'zdata',v(3));
end

if isempty(p)
    p = [nan nan nan];
    set(state.marker1,'visible','off');
else
    set(state.marker1,'visible','on','xdata',p(1),'ydata',p(2),'zdata',p(3));
end

% Update tool and markers
set(state.text,'string',createString(p(1),p(2),p(3),v(1),v(2),v(3),vi));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [fig, h] = createUI

pos = [200 200 200 200];

% Create selection tool %
fig = figure('handlevisibility','off','menubar','none','resize','off',...
    'numbertitle','off','name','Select 3-D Tool','position',pos,'deletefcn','select3dtool(''off'')');

h = uicontrol('style','text','parent',fig,'string',createString(0,0,0,0,0,0,0),...
    'units','norm','position',[0 0 1 1],'horizontalalignment','left');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [str] = createString(px,py,pz,vx,vy,vz,vi)

str = sprintf('  Position:\n  X  %f\n  Y:  %f\n  Z:  %f  \n\n  Vertex:\n  X:  %f\n  Y:  %f\n  Z:  %f  \n\n  Vertex Index:\n  %d',px,py,pz,vx,vy,vz,vi);


text:
function y=renyi_entro(DATA,q)
  [M,N]=size(DATA);
       y=zeros(1,N);
       for n=1:N
           y(1,n)=log(sum(DATA(:,n).^q))/(1-q);
       end


text:
function varargout = svd(f)
%SVD   Singular value decomposition of a DISKFUN.
%   SVD(F) returns the singular values of F. The number of singular values
%   returned is equal to the rank of the DISKFUN.
%
%   S = SVD(F) returns S, a vector of singular values in non-increasing
%   order.
%
%   [U, S, V] = SVD(F) returns the SVD of F. V is a quasimatrix of
%   orthogonal CHEBFUN objects, U is a quasimatrix of CHEBFUN objects that
%   are orthogonal with respect to the r weight on [0,1] (derived from the
%   measure on the disk) and S is a diagonal matrix with the singular
%   values on the diagonal.
%
%   The length and rank of a DISKFUN are slightly different quantities.
%   LENGTH(F) is the number of pivots used by the constructor, and
%   RANK(F) is the number of significant singular values of F. The relation
%   RANK(F) <= LENGTH(F) should always hold.

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

if ( isempty(f) )
    varargout = { [ ] };
    return
end

% Get CDR decomposition of f:
[C, D, R] = cdr(f);
C = restrict(C, [0, 1]);

% Do QR in both variables, one with the weighted inner-product and one with
% the standard L2 inner-product.
[QwC, RwC] = diskQR(C);
[QwR, RwR] = qr(R);

% Use the QR factorizations of the columns and rows to make up the SVD of
% the DISKFUN object.  Since
%
%        C * D * R' = QwC * ( RwC * D * RwR' ) * QwR'
%
% we compute the SVD of ( RwC * D * RwR' ).
[U, S, V] = svd( RwC * D * RwR.' );
U = QwC * U;
V = QwR * V;

% Output just like the svd of a matrix.
if ( nargout > 1 )
    varargout = { U, S, V };
else
    varargout = { diag(S) };
end

end

function [Q, R] = diskQR(A)
% Fast version of the abstractQR code, specifically for the disk.

n = length(A) + 1;           % Can probably get away with a smaller n.
[r, w] = legpts(n, [0,1]);    % Legpts

% Do a weighted QR, and then unweight the QR:
WR = spdiags(sqrt(w.'.*r), 0, n, n);
invWR = spdiags(1./sqrt(w.'.*r), 0, n, n);

% Discrete QR with inner product <u,v> = sum(r*conj(u).*v):
[discreteQ, discreteR] = qr(WR*A(r, :), 0);

s = sign(diag(discreteR));    % }
s(~s) = 1;                    % } Enforce diag(R) >= 0
S = diag(s);                  % }

% Undo the weighting:
discreteQ = invWR*discreteQ*S;
% Correct the signs in R:
discreteR = S*discreteR;

% Go back to continuous land:
Q = chebfun(legvals2chebvals(discreteQ), [0, 1]);
R = discreteR;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%fill_prmm Compute the null-space and fill PRMM.
%
%   Parameters:

function [P,X, u1,u2, lambda, info] = fill_prmm(M, Idepths, central, ...
                                                opt, info)

[NULLSPACE, result] = create_nullspace(M, Idepths, central, ...
                                       opt.create_nullspace);

info.create_nullspace             = opt.create_nullspace;
info.sequence{end}.tried          = result.tried;
info.sequence{end}.tried_perc     = result.tried/comb(size(M,2), 4) *100;
info.sequence{end}.used           = result.used;
info.sequence{end}.used_perc      = result.used/result.tried *100;
info.sequence{end}.failed         = result.failed;
info.sequence{end}.size_nullspace = size(NULLSPACE);

if opt.verbose
  disp(sprintf('Tried/used: %d/%d (%.1e %%/ %2.1f %%)', result.tried, ...
               result.used, info.sequence{end}.tried_perc, ...
               info.sequence{end}.used_perc));
  disp(sprintf('%d x %d'' is size of the nullspace', size(NULLSPACE,2), ...
               size(NULLSPACE,1))); end

[m,n] = size(M); m = m/3;

if size(NULLSPACE,2) == 0
  P = []; X = []; u1 = 1:m; u2 = 1:n; lambda=[];
else
  r = 4;
  [L, S] = nullspace2L(NULLSPACE, r, opt);
  clear NULLSPACE;

  if isempty(opt.create_nullspace), threshold = .01;
  else, threshold = opt.create_nullspace.threshold; end
  if svd_suff_data(S, r, threshold)
    if opt.verbose
      dS = diag(S); %disp(diag(S));
      fprintf(1,'Smallest 2r singular values:%s.\n', sprintf(' %f', ...
                                                  dS(end-2*r:end))); end
    [Mdepths, lambda] = L2depths(L, M, Idepths, opt); info.Mdepths = Mdepths;
    [P,X, u1b,u2] = approximate(Mdepths, r, L, opt);
    u1 = union(ceil(u1b/3),[]); killb = setdiff(k2i(u1),u1b);
    if ~isempty(killb), r1b = setdiff(1:3*m,u1b); kill = killb;
      for ib = killb(1:end-1), lower = find(killb > ib);
        if kill(lower(1)-1) < kill(lower(1)) - 1,
          kill(lower) = kill(lower) -1; end; end
      P = P(setdiff(1:length(r1b),kill),:);end
    lambda = lambda(setdiff(1:m,u1),setdiff(1:n,u2)); % to fit P*X
  else P=[]; X=[]; u1=1:m; u2=1:n; lambda=[]; end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [L,S] = nullspace2L(NULLSPACE, r, opt)
% Compute the basis of MM (L) from the null-space.

if opt.verbose, fprintf(1,'Computing the basis...'); tic; end
if size(NULLSPACE,2) < 10*size(NULLSPACE,1) % orig:[A,S,U] = svd(NULLSPACE',0);
  [U,S,V]             = svd(NULLSPACE);
else
  [U,SS]              = eig(NULLSPACE*NULLSPACE');
  dSS                 = diag(SS); l = length(dSS);
  [sortdSS(l:-1:1),I] = sort(dSS); I(l:-1:1) = I;  U = U(:,I);
  sortdSS             = max(sortdSS,zeros(1,l)); sortdSS=sqrt(sortdSS);
  for i=1:l, S(i,i) = sortdSS(i); end
end

lenU = length(U);
L = U(:,lenU+1-r:lenU);
if opt.verbose, disp(['(' num2str(toc) ' sec)']); end

function r = comb(n,k)

% returns combination number n over k

r=1;
for i=1:k
  r=r*(n-i+1)/i;
end

function y = svd_suff_data(S,r, threshold)

% S is the singular value part of the svd of the nullspaces of the column
% r-tuples.  We'll want to be able to take the r least significant columns
% of U.  This is right because the columns of U should span the whole space
% that M's columns might span.  That is, M is FxP.  The columns of U should
% span the F-dimensional Euclidean space, since U is FxF.  However, we want
% to make sure that the F-r-1'th singular value of S isn't tiny.  If it is,
% our answer is totally unreliable, because the nullspaces of the column
% r-tuples don't have sufficient rank.  If this happens, it means that the
% intersection of the column cross-product spaces is bigger than r-dimensional,
% and randomly choosing an r-dimensional subspace of that isn't likely to
% give the right answer.

Snumcols = size(S,2);
Snumrows = size(S,1);
if (Snumrows == 0 | Snumcols + r < Snumrows | Snumrows <= r)
  y = 0;
else
  y = S(Snumrows-r,Snumrows-r) > threshold;
end

author: strawlab

repo: MultiCamSelfCal


text:
function yred=rednoise(N,g,a);
% REDNOISE: A fast rednoise generator using filter.
%
% example: rednoise(N,g);
%
% Description: generates a rednoise series
% with zero process mean.
% note: statistical mean&variance will be different.
%
% Inputs: n - desired length of time series
%         g - lag-1 autocorrelation
%         a - noise innovation variance parameter (optional, default=1)
%
% Example:
%   plot(rednoise(1000,.95))
%
% Aslak Grinsted 2006-2014

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

if nargin<3
    a=1;
end

if g==0
    yred=randn(N,1)*a;
    return
end
tau=ceil(-2/log(abs(g))); %2 x de-correlation time

yred=filter([1 0],[1;-g],randn(tau+N,1)*a);
yred=yred(tau+1:end);

author: grinsted

repo: wavelet-coherence


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

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

% Last Modified by GUIDE v2.5 29-Jun-2010 21:17:19

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

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

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

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

% Update handles structure
guidata(hObject, handles);
img=imread('ascii.png');    %Leggo la tabella ascii da disco
axes(handles.axes2)         %Fisso l'axes per il plot
imshow(img)                 %Mostro la tabella ascii
% UIWAIT makes ascii wait for user response (see UIRESUME)
% uiwait(handles.figure1);

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

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

% --- Executes when figure1 is resized.
function figure1_ResizeFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


text:
% OPTIMI toolbox
% Version 0.13		Monday 21 Nov 2005 at 08:49
% Copyright (c) 2005 Neil D. Lawrence
%
% CMPNDTIEPARAMETERS Tie parameters together.
% EXPTRANSFORM Constrains a parameter to be positive through exponentiation.
% NEGLOGLOGIT Constrains a parameter to be positive.
% OPTIMISEPARAMS Optimise parameters.
% OPTOPTIONS Give optimisation options for NETLAB.
% SIGMOIDTRANSFORM Constrains a parameter to be between 0 and 1.

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function test_bug3174

% WALLTIME 00:20:00
% MEM 3gb
% DEPENDENCY ft_volumerealign

load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug3174.mat'));

% co-register the CT to the MRI
cfg             = [];
cfg.method      = 'spm';
cfg.spmversion  = 'spm12';
cfg.viewresult  = 'yes'; % view realignment result
cfg.coordsys    = 'tal';
ct2              = ft_volumerealign(cfg, cttal, mri);


text:
classdef LIRCMOP1 < PROBLEM
% <multi> <real> <large/none> <constrained>
% Constrained benchmark MOP with large infeasible regions

%------------------------------- Reference --------------------------------
% Z. Fan, W. Li, X. Cai, H. Huang, Y. Fang, Y. You, J. Mo, C. Wei, and E.
% Goodman, An improved epsilon constraint-handling method in MOEA/D for
% CMOPs with large infeasible regions, Soft Computing, 2019, 23:
% 12491-12510.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Wenji Li

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values and constraint violations
        function Population = Evaluation(obj,varargin)
            X = varargin{1};
            X = max(min(X,repmat(obj.upper,size(X,1),1)),repmat(obj.lower,size(X,1),1));
            x_odd       = X(:,3:2:end);
            x_even      = X(:,2:2:end);
            g_1         = sum((x_odd - sin(0.5 * pi * X(:,1))).^2,2);
            g_2         = sum((x_even - cos(0.5 * pi * X(:,1))).^2,2);
            PopObj(:,1) = X(:,1) + g_1;
            PopObj(:,2) = 1 - X(:,1) .^ 2 + g_2;
            PopCon(:,1) = (0.5 - g_1).*(0.51 - g_1);
            PopCon(:,2) = (0.5 - g_2).*(0.51 - g_2);
            Population  = SOLUTION(X,PopObj,PopCon,varargin{2:end});
            obj.FE      = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - R(:,1).^2;
            R      = R + 0.5;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
% Digital Video Stabilization and Rolling Shutter Correction using Gyroscopes
% Copyright (C) 2011 Alexandre Karpenko
%
% 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
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

function new_theta = integrate_and_resample(w, old_time, new_time)

dT = old_time(2:end) - old_time(1:end-1);
theta = cumsum( w(2:end) .* dT );

new_theta = interp1(old_time(2:end), theta(:,1), new_time);
new_theta(isnan(new_theta)) = 0;

author: alex-golts

repo: Video-Stabilization


text:
function uI = faceinterpolate(u,node,elem,elemType)
%% FACEINTERPOLATE interpolate to face elements (RT0 or BDM1).
%
% uI = faceinterpolate(u,node,elem,elemType) interpolates a given function u
% into the lowesr order RT0 or BDM1 finite element spaces. The coefficient
% is given by the line integral int_e u*n ds. The input elemType can be 'RT0'
% or 'BDM1'.
%
% Example: RT0
%
%      maxIt = 5;
%      node = [1,0; 1,1; 0,1; -1,1; -1,0; -1,-1; 0,-1; 0,0]; % nodes
%      elem = [1,2,8; 3,8,2; 8,3,5; 4,5,3; 7,8,6; 5,6,8];    % elements
%      bdFlag = setboundary(node,elem,'Dirichlet');
%      pde = mixBCdata;
%      err = zeros(maxIt,2);
%      h = zeros(maxIt,1);
%      for k =1:maxIt
%        [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
%        [u,sigma,eqn] = PoissonRT0(node,elem,bdFlag,pde);
%        sigmaI = faceinterpolate(pde.Du,node,elem,'RT0');
%        err(k,1) = getL2errorRT0(node,elem,pde.Du,sigmaI);
%        err(k,2) = sqrt((sigma-sigmaI)'*eqn.M*(sigma-sigmaI));
%        h(k) = 1./(sqrt(size(node,1))-1);
%      end
%      figure;
%      showrateh2(h,err(:,1),2,'r-+','|| \sigma - \sigma_I ||',...
%                 h,err(:,2),2,'b-+','|| \sigma_h - \sigma_I ||');
%
% Example: BDM1
%
%      maxIt = 5;
%      node = [1,0; 1,1; 0,1; -1,1; -1,0; -1,-1; 0,-1; 0,0]; % nodes
%      elem = [1,2,8; 3,8,2; 8,3,5; 4,5,3; 7,8,6; 5,6,8];    % elements
%      bdFlag = setboundary(node,elem,'Dirichlet');
%      pde = mixBCdata;
%      err = zeros(maxIt,2);
%      h = zeros(maxIt,1);
%      for k =1:maxIt
%        [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
%        [u,sigma,eqn] = PoissonBDM1(node,elem,bdFlag,pde);
%        sigmaI = faceinterpolate(pde.Du,node,elem,'BDM1');
%        err(k,1) = getL2errorBDM1(node,elem,pde.Du,sigmaI);
%        err(k,2) = sqrt((sigma-sigmaI)'*eqn.M*(sigma-sigmaI));
%        h(k) = 1./(sqrt(size(node,1))-1);
%      end
%      figure;
%      showrateh2(h,err(:,1),2,'r-+','|| \sigma - \sigma_I ||',...
%                 h,err(:,2),2,'b-+','|| \sigma_h - \sigma_I ||');
%
%
% Example: RT1
%
%      maxIt = 5;
%      node = [1,0; 1,1; 0,1; -1,1; -1,0; -1,-1; 0,-1; 0,0]; % nodes
%      elem = [1,2,8; 3,8,2; 8,3,5; 4,5,3; 7,8,6; 5,6,8];    % elements
%      bdFlag = setboundary(node,elem,'Dirichlet');
%      pde = mixBCdata;
%      err = zeros(maxIt,2);
%      h = zeros(maxIt,1);
%      for k = 1:maxIt
%        [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
%        sigmaI = faceinterpolate(pde.Du,node,elem,'RT1');
%        err(k,1) = getL2errorRT1(node,elem,pde.Du,sigmaI);
%        h(k) = 1./(sqrt(size(node,1))-1);
%      end
%      figure;
%      showrateh2(h,err(:,1),2,'r-+','|| \sigma - \sigma_I ||',...
%                 h,err(:,2),2,'b-+','|| \sigma_h - \sigma_I ||');
%
%
% See also edgeinterpolate, edgeinterpolate1, edgeinterpolate2
%
% Created by Ming Wang at Mar 29, 2011, M-lint modified at May 14, 2011.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('elemType','var'), elemType = 'RT0'; end
%% edge
if size(elem,2) == 2 % the input elem is edge
    edge = elem;
else
    [elem2edge,edge] = dofedge(elem);
end
NE = size(edge,1);
NT = size(elem,1);
edgeVec = node(edge(:,2),:) - node(edge(:,1),:);
nVec = zeros(NE,2);
nVec(:,1) = edgeVec(:,2);
nVec(:,2) = -edgeVec(:,1);
clear edgeVec

%% dof for RT0
[lambda,weight] = quadpts1(4);
nQuad = size(lambda,1);
uI = zeros(NE,1);
for i = 1:nQuad
    pxy = lambda(i,1)*node(edge(:,1),:)+lambda(i,2)*node(edge(:,2),:);
    flux = u(pxy);
    uI = uI + weight(i)*dot(flux,nVec,2);
end

%% dof for BDM1
if strcmp(elemType,'BDM1') || strcmp(elemType,'RT1')
   uI(NE+1:2*NE) = zeros(NE,1);
   for i = 1:nQuad
        pxy = lambda(i,1)*node(edge(:,1),:)+lambda(i,2)*node(edge(:,2),:);
        uI(NE+1:2*NE) = uI(NE+1:2*NE)+ ...
                     weight(i)*3*(lambda(i,1)-lambda(i,2))*dot(u(pxy),nVec,2);
   end
end

%% dof for RT1
if strcmp(elemType,'RT1')
    mid = (node(edge(:,1),:) + node(edge(:,2),:))/2;
    umid = u(mid);
    uI(2*NE+2*NT) = 0;
    % Face dof coefficients
    uquadpts = umid(elem2edge(:,1),:)+umid(elem2edge(:,2),:)+umid(elem2edge(:,3),:);
    lf = [4*dot(nVec(elem2edge(:,2),:),uquadpts,2) ...
          4*dot(nVec(elem2edge(:,3),:),uquadpts,2)];
    elem2edgeDofValue = [uI(elem2edge(:,1:3)) uI(elem2edge(:,1:3)+NE)];
    localMatrix = [4 8 4 -4  0 4; ...
                   8 4 -4 0 -4 4]';
    lf = lf - elem2edgeDofValue*localMatrix;
    uI((2*NE+1):end) = [2*lf(:,1) - lf(:,2); ...
                        2*lf(:,2) - lf(:,1)]/3;
end

repo: ifem


text:
function [rDisambig,sourceInfo]=disambiguateClust1D(rMeas,ampMeas,deltaR,rMax,clustDist,twoWayDisambig,threshold)
%%DISAMBIGUATECLUST1D This function uses a clustering algorithm to
%                   disambiguate measurements in one dimension. It can
%                   handle measurements from multiple simultaneous targets.
%                   Information on which detections went into the
%                   disambiguated measurements is also returned. The
%                   snapshots all have a different number of bins; the
%                   lengths in deltaR should be relatively prime (coprime)
%                   integers (possibly all times a scaling constant). When
%                   the "true" value of a detection is outside of the range
%                   0<=rMeas<deltaR(i) for the ith interval, it is aliased
%                   back in. This function tries to dealias the
%                   measurements and determine how many detections are
%                   present. This function only considers aliasing in a
%                   positive (+deltaR) direction and is thus best suited
%                   for dealiasing ranges.
%
%INPUTS: rMeas A numIntX1 or 1XnumInt cell matrix with the detections in
%              each interval. Each cell contains a numMeasX1 or 1XnumMeas
%              vector of the detections in that interval. Detections should
%      ampMeas A numIntX1 or 1XnumInt cell matrix with the positive real
%              amplitudes associated with the detections in each interval.
%              These are used to perform a weighted average to determine
%              the true target location.
%       deltaR The numIntX1 or 1XnumInt vector of the unambiguous range in
%              each interval. These values must be relatively prime
%              integers (times a common possibly non-integer scale factor)
%              for the algorithm to work.
%         rMax The maximum unambiguous range to consider for the
%              completely unaliased measurements (all dealiased values
%              are less than this. If this value is omitted or an empty
%              matrix is passed, then the default of prod(deltaR) is used,
%              which is the maximum from the Chinese remainder theorem.
%              However, as demonstrated below, the presence of multiple
%              targets can cause extra distant "ghost" detections to appear
%              due to the interference of the targets. Thus, it is often
%              good to design the maximum unambiguous value to be very big
%              and to limit the maximum number of bins actually considered.
%              Put another way, use more coprime intervals.
%    clustDist This is the distance to determine that two aliased-out
%              values are in the same cluster. This depends on the scaling
%              of the intervals. The default if this parameter is omitted
%              or an empty matrix is passed is 0.5, which is good if deltaR
%              are all integers.
% twoWayDisambig If this is false, then the integers over which aliasing
%              are performed are only positive. Otherwise, they can be
%              positive and negative, which can be more convenient for
%              Doppler disambiguation.
%    threshold An optional threshold for declaring a detection. At least
%              this many common intervals must cluster together. The
%              default if this parameter is omitted or an empty matrix is
%              passed is numInt.
%
%OUTPUTS: rDisambig The values of the detections, a numDetectX1 vector.
%              These are weighted averages of the values in rMeas that went
%              into the detections. If there are no detections, then this
%              will be an empty matrix.
%        sourceInfo Information allowing one to reconstruct which values in
%              rMeas went into forming each detection. This is a
%              numDetectX1 cell array. Each cell ion the array has a
%              numMeasX2 matrix providing information on the measurements
%              in that cluster. The first column gives the number of the
%              interval that produced the measurement (index of rMeas) and
%              the second column provides information on which measurement
%              in that interval led to the detection.
%
%The general idea of such clustering is described in [1] and [2]. However,
%it is done differently here. The authors in [2] described a
%disambiguiation algorithm that is unusable if one of the intervals does
%not produce detections and the basic algorithm in [1] is only focussed on
%a single target. Here, we alias out the measurement across all ambiguity
%intervals and then determine which aliased out measurements gate. Two
%values gate if they are from different intervals and are within clustDist
%of each other. A DisjointSet data structure is used to take the gating
%ifnromation and cluster the aliased out measurements. Clusters of a
%sufficient size are counted as detections. Note that this clustering is
%not globally optimal as it is still possible for two detections from a
%single PRI to end up in one cluster, which is only particularly likely
%with closely spaced targets. An optimal clustering algorithm would be a
%multiframe assignment problem.
%
%Note that because linear averaging is performed, points very near the
%edges will not average to the correct values. However, since measurement
%noises are typically small and unambiguous intervals large, this is seldom
%a problem.
%
%EXAMPLE 1:
%This is similar to the toy example given in [2], except we add some noise
%to the measurements In this instance, the algorithm will produce results
%that are close to the true unaliased bins.
% deltaR=[7;8;11];
% maxR=24;
%
% trueBins=[6;13];
%
% rMeas=cell(3,1);
% for curInt=1:3
%     rMeas{curInt}=unique(wrapRange(trueBins,0,deltaR(curInt)));
%     %Add noise.
%     rMeas{curInt}=rMeas{curInt}+0.1*randn(size(rMeas{curInt}));
% end
%
%[rDisambig,sourceInfo]=disambiguateClust1D(rMeas,[],deltaR,maxR)
%
%EXAMPLE 2:
%This example is better suited for Doppler disambiguation Here, positive
%and negaitve values are present, which could correspond to positive and
%negative range rates. The disambiguation is thus two-sided and the maximum
%number of bins used is limited so that irrelevant extra-fast targets are
%not detected.
% deltaR=[7;8;11;13];
% maxR=100;
% %-50 to 49
% trueBins=[6;-13];
%
% rMeas=cell(3,1);
% for curInt=1:4
%     rMeas{curInt}=unique(wrapRange(trueBins,0,deltaR(curInt)));
%     %Add noise.
%     rMeas{curInt}=rMeas{curInt}+0.1*randn(size(rMeas{curInt}));
% end
% disambiguateClust1D(rMeas,[],deltaR,maxR,[],true)
%
%REFERENCES:
%[1] G. Trunk and S. Brockett, "Range and velocity ambiguity resolution,"
%    in Record of the IEEE National Radar Conference, Lynnfield, MA, 20-22
%    Apr. 1993, pp. 146-149.
%[2] P. Stinco, M. Greco, F. Gini, A. Farina, and L. Timmoneri, "Analysis
%    and comparison of two disambiguity algorithms: The modified CA and
%    CRT," in Proceedings of the International Radar Conference -
%    Surveillance for a Safer World, Bordeaux, France, 12-16 Oct. 2009.
%
%December 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<4||isempty(rMax))
   rMax=prod(deltaR);
end

numPRI=length(rMeas);
if(nargin<5||isempty(clustDist))
   clustDist=0.5;
end

if(nargin<6||isempty(twoWayDisambig))
   twoWayDisambig=false;
end

if(nargin<7||isempty(threshold))
   threshold=numPRI;
end

%Determine the maximum number of aliased out measurements in both
%directions.
numAliased=0;
for curPRI=1:numPRI
    aliasedMax=ceil(rMax/deltaR(curPRI));
    numAliased=numAliased+length(rMeas{curPRI})*aliasedMax;
end

%We will now alias out all of the measurements
measList=zeros(numAliased,1);
%souceData(i,1) holds the PRI that produced the aliased-out measurement.
%sourceData(i,2) holds the index of the measurement in the PRI.
sourceData=zeros(numAliased,2);
numAdded=0;
for curPRI=1:numPRI
    curMeas=rMeas{curPRI};
    numMeas=length(curMeas);

    aliasedMax=ceil(rMax/deltaR(curPRI));
    if(twoWayDisambig)
        valRange=(-ceil((aliasedMax)/2)):fix((aliasedMax)/2);
    else
        valRange=(0:(aliasedMax-1));
    end

    aliasedVals=bsxfun(@plus,curMeas(:),deltaR(curPRI)*valRange);
    origIdx=repmat(1:numMeas,1,aliasedMax+1);
    aliasedVals=aliasedVals(:);
    origIdx=origIdx(:);
    if(twoWayDisambig)
        %Get rid of values that aliased to values that are too low.
        %Get rid of values that aliased too low or too high
        sel=(aliasedVals>=-rMax/2)&(aliasedVals<=rMax/2);
        aliasedVals=aliasedVals(sel);
    else
        %Get rid of values that aliased to values that are too high.
        sel=aliasedVals<=rMax;
        aliasedVals=aliasedVals(sel);
    end

    origIdx=origIdx(sel);

    num2Add=length(aliasedVals);
    idx2Add=(numAdded+1):(numAdded+num2Add);

    measList(idx2Add)=aliasedVals;

    sourceData(idx2Add,1)=curPRI;
    sourceData(idx2Add,2)=origIdx;

    numAdded=numAdded+num2Add;
end

%Shrink to fit.
measList=measList(1:numAdded);
sourceData=sourceData(1:numAdded,:);

%Sort the values.
[measList,idx]=sort(measList,'ascend');
sourceData=sourceData(idx,:);

%Create a DisjointSet and determine which measurements gate with each other.
%The DisjointSet will cluster them. To gate, they must be sufficiently
%close and also not within the same interval. The loops go through all
%pairs of aliased-out measurements. We take advantage of having sorted the
%measurements so that once a distance comparison fails, no measurements are
%farther ranges are considered.
gateSet=DisjointSet(numAdded);
for idx1=1:(numAdded-1)
    for idx2=(idx1+1):numAdded
        %If they originate in different PRIs.
        if(sourceData(idx1,1)~=sourceData(idx2,1))
            if(abs(measList(idx1)-measList(idx2))<clustDist)
                %The two values gate together.
                gateSet.unionFromList([idx1;idx2]);
            else
                break;
            end
        end
    end
end

%Form clusters.
theClust=gateSet.createClusterSet();

%Count how many clusteres are large enough to keep.
numClusters=theClust.numClusters();

numPastThresh=sum(theClust.clustSize>=threshold);
rDisambig=zeros(numPastThresh,1);
sourceInfo=cell(numPastThresh,1);

curDetectClust=0;
for curClust=1:numClusters
    numMeasCur=theClust.clustSize(curClust);
    if(numMeasCur>=threshold)
        sourceDataCur=sourceData(theClust(curClust,:),:);

        rAmpCur=zeros(numMeasCur,1);
        if(~isempty(ampMeas))
            for curMeas=1:numMeasCur
                PRI=sourceDataCur(curMeas,1);
                idx=sourceDataCur(curMeas,2);
                idxR=theClust(curClust,curMeas);

                rAmpCur(curMeas,1)=measList(idxR);
                rAmpCur(curMeas,2)=ampMeas{PRI}(idx);
            end
        else
             for curMeas=1:numMeasCur
                idxR=theClust(curClust,curMeas);

                rAmpCur(curMeas,1)=measList(idxR);
                rAmpCur(curMeas,2)=1;
            end
        end

        curDetectClust=curDetectClust+1;

        %A simple weighted average.
        rDisambig(curDetectClust)=sum(rAmpCur(:,1).*rAmpCur(:,2))/sum(rAmpCur(:,2));
        sourceInfo{curDetectClust}=sourceDataCur;
    end
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
 function out = fatrix2_subsref_dot(st, args)
%function out = fatrix2_subsref_dot(st, args)
% handle subscript references like ob.ref or ob.ref(arg[s])
% This will called from ../subsref with struct(ob)
% and it will have arg(1).type = '.'
% Based on @strum/subsref
% Copyright 2006-1-19, Jeff Fessler, University of Michigan

%if args(1).type ~= '.'
%	error 'only ob.data and ob.meth(args) done'
%end

name = args(1).subs;
%printm('here, name=%s', name)

switch name
case 'imask_array'
	out = fatrix2_imask_array(st);
	return

case 'omask_array'
	out = fatrix2_omask_array(st);
	return

end

% ob.name
if isfield(st, name)
	if length(args) == 1
		out = st.(name);
	else
		out = subsref(st.(name), args(2:end));
	end

% ob.name or ob.name(arg[s]), where name is in meth field
elseif isfield(st.meth, name)
	fun = st.meth.(name); % function handle for method
	if length(args) == 1
		out = fun(st.arg);
	elseif length(args) == 2
		out = fun(st.arg, args(2).subs{:});
	else
		printm('unknown subs')
		keyboard
	end

% ob.name, where name is in data field 'arg' (last resort)
elseif isfield(st.arg, name)
	out = st.arg.(name);
	if length(args) > 1

		% ob.name{:} does not work properly because there seems
		% to be no way to return an arbitrary number of arguments
		% (in a list) to the caller.
		% This is unfortunate because
		% a.b = {1,2} and a.b{:} works for ordinary structures.
% todo: see this solution involving numel()
% http://www.mathworks.com/support/solutions/data/1-1ABOD.html?solution=1-1ABOD
		if iscell(out) && ...
		(isequal(args(2).subs, {':'}) || length(args(2).subs{1}) > 1)
			warn 'instead of strum.data{:} etc., use two steps:'
			warn 'tmp = strum.data; tmp{:} or tmp{1:end} etc.'
			fail 'strum.data{:} etc. unsupported thanks to matlab'
		end

		% fix: strum1.strum2.plot() does not work here
		out = subsref(out, args(2:end));
	end

else
	disp(st)
	fail(['unknown field name: ' name])
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [u v] = coarse_to_fine(I1, I2, settings, show_flow, h)

[M N] = size(I1);

% computes the maximum number of pyramid levels; the coarsest image should
% have a width or height around 10
pyramid_levels = min(...
    ceil(log(10/M)/log(settings.pyramid_factor)), ...
    ceil(log(10/N)/log(settings.pyramid_factor)));

pyrI1 = cell(pyramid_levels, 1);
pyrI2 = cell(pyramid_levels, 1);

pyrI1{1} = I1;
pyrI2{1} = I2;

% build the pyramids
for i = 2:pyramid_levels
  pyrI1{i} = imresize(I1, (settings.pyramid_factor)^(i-1), settings.resampling_method);
  pyrI2{i} = imresize(I2, (settings.pyramid_factor)^(i-1), settings.resampling_method);
end

% start coarse to fine processing
for level = pyramid_levels:-1:1;

  [M N] = size(pyrI1{level});
  if level == pyramid_levels

    % initialization
    u = zeros(M, N);
    v = zeros(M, N);

    pu = zeros(M, N, 2);
    pv = zeros(M, N, 2);

  else
    % previous dimensions
    [Mp Np] = size(pyrI1{level+1});

    % upsample the flow to next level
    u = imresize(u, [M N], settings.resampling_method) * N/Np;
    v = imresize(v, [M N], settings.resampling_method) * M/Mp;

    pu_tmp = pu;
    pv_tmp = pv;

    pu = zeros(M, N, 2);
    pv = zeros(M, N, 2);

    for i=1:2
      pu(:,:,i) = imresize(pu_tmp(:,:,i), [M N], settings.resampling_method);
      pv(:,:,i) = imresize(pv_tmp(:,:,i), [M N], settings.resampling_method);
    end
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % solve the optical flow on the current level
  [u, v, pu, pv] = solve_flow_on_level(u, v, pu, pv, pyrI1{level}, pyrI2{level}, settings, level, show_flow, h);
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end

author: thomas-koehler

repo: SupER


text:
function llh = GFM_llh( theta, tree )

% Evaluate the likelihood of data under the estimated GFM model
% This is *exactly* the same procedure for likelihood computation as in
% GFM_EM
%
% Syntax:
%   llh = GFM_llh( theta, tree )
%
%
% Input:
%   theta : Parameter cell for GFM model.
%
%   tree  : Wavelet tree from DWT2_TO_CELL
%
%
% Output:
%   llh   : Vector with log-likelihood for each direction and all trees

if size(theta, 1) ~= numel(tree)-1
    error('Wavelet transform must have the same size as theta');
end

L = length(tree) - 1;
M = numel( theta{1} );

% Naming convention:
% beta1 is beta_u
% beta2 is beta_{rho(u),u}

[no_rows, no_cols] = cellfun( @(x) size(x{1}), tree(2:end) );
[beta1, lh, cond_joint_probs] = deal( cell(L,1) );

for l = 1:L
    beta1{l} = ones(no_rows(l), no_cols(l), M);
    lh{l} = zeros(no_rows(l), no_cols(l));
    cond_joint_probs{l} = ones(no_rows(l), no_cols(l), M, M);
end

[beta2, prod_beta2, state_probs] = deal( beta1 );
[beta1_rep] = deal( cond_joint_probs );
llh = zeros(3, 1);

for d = 1:3
    for l = L:-1:1
        w{l} = repmat(tree{l+1}{d}, [1 1 M]);
        likelihood = normpdf( w{l}, 0, ...
            repmat(reshape(theta{l,3,d}, [1 1 M]), [no_rows(l) no_cols(l) 1]) );

        % Make sure numerically large coefficients doesn't evaluate to all
        % zeros
        idx = find( sum(likelihood, 3) == 0 );
        if ~isempty( idx )
            likelihood(idx) = 1/M;
        end

        state_probs{l} = repmat( reshape( theta{l,1,d}, [1 1 M] ), [no_rows(l) no_cols(l) 1] );

        % For the non-leaf nodes beta1 needs an additional factor
        if l < L
            % Compute beta2
            trans_probs{l+1} = repmat(reshape( theta{l+1,2,d}, [1 1 M M] ), [no_rows(l+1) no_cols(l+1) 1]);
            beta1_rep{l+1} = repmat( beta1{l+1}, [1 1 1 M] );
            state_probs_rep{l+1} = repmat( state_probs{l+1}, [1 1 1 M] );

            beta2_full = trans_probs{l+1} .* beta1_rep{l+1} ./ state_probs_rep{l+1};
            beta2{l+1} = squeeze( sum(beta2_full, 3) );

            % Product of beta2's with common parent
            prod_beta2{l} = ...
                beta2{l+1}(1:2:end, 1:2:end, :) .* ...
                beta2{l+1}(2:2:end, 1:2:end, :) .* ...
                beta2{l+1}(1:2:end, 2:2:end, :) .* ...
                beta2{l+1}(2:2:end, 2:2:end, :);

            % Unnormalized beta1
            tmp = prod_beta2{l} .* likelihood .* state_probs{l};
        else
            % Unnormalized beta1
            tmp = likelihood .* state_probs{l};
        end

        % Likelihood and beta1
        lh{l} = sum( tmp, 3 );
        beta1{l} = tmp ./ repmat( lh{l}, [1 1 M] );
    end

    % The total likelihood for this direction
    llh(d) = sum( cellfun(@(x) sum(log( max(x(:), eps) )), lh) );
end

end


text:
function [w,u] = biharmonicP2(node,elem,pde,bdFlag,option)
%   [w,u] = biharmonicP1(node,elem,pde,bdFlag) produces the mixed cubic finite element
%   approximation of the biharmonic equation, where w = laplace u
%   See also biharmonicP1, biharmonicP2, biharmonicP3.
%   Created by Jie Zhou.
%   Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if nargin<5, option = []; end
tic;
%% Construct Data Structure
[elem2dof,edge,bdDof] = dofP2(elem);
N = size(node,1);  NT = size(elem,1);  Ndof = N+size(edge,1);

%% Compute geometric quantities and gradient of local basis
[Dlambda,area] = gradbasis(node,elem);

%% Assemble stiffness matrix
% Since Dphi_i*Dphi_j is quadratic, numerical quadrature rule is used here
if ~isfield(option,'quadorder')
    option.quadorder = 4;   % default order
end
[lambda, weight] = quadpts(option.quadorder);
nQuad = size(lambda,1);
ii = zeros(21*NT,1); jj = zeros(21*NT,1); sA = zeros(21*NT,1);sB = zeros(21*NT,1);
index = 0;
for i = 1:6
    for j = i:6
        Bij = 0;
        Aij = 0;
        for p = 1:nQuad
                 Bij = Bij + weight(p)*dot(Dphi(p,i),Dphi(p,j),2);
                 Aij = Aij + weight(p)*dot(phi(p,i),phi(p,j),2);
        end
        Bij = Bij.*area;
        Aij = Aij.*area;

        ii(index+1:index+NT) = double(elem2dof(:,i));
        jj(index+1:index+NT) = double(elem2dof(:,j));
        sB(index+1:index+NT) = Bij;
        sA(index+1:index+NT) = Aij;
        index = index + NT;
    end
end
clear Aij Bij
diagIdx = (ii == jj);   upperIdx = ~diagIdx;
B = sparse(ii(diagIdx),jj(diagIdx),sB(diagIdx),Ndof,Ndof);
A = sparse(ii(diagIdx),jj(diagIdx),sA(diagIdx),Ndof,Ndof);
% A = spdiags(accumarray(ii(diagIdx),sA(diagIdx),[Ndof 1]),0,Ndof,Ndof);
BU = sparse(ii(upperIdx),jj(upperIdx),sB(upperIdx),Ndof,Ndof);
AU = sparse(ii(upperIdx),jj(upperIdx),sA(upperIdx),Ndof,Ndof);
B = B + BU + BU';
A = A + AU + AU';

%% boundary condition
%%
    fixedDof = [];
    isFixedDof = false(Ndof,1);
    if ~isempty(bdFlag)
        elem2edge = elem2dof(:,4:6)-N;
        isDirichlet(elem2edge(bdFlag(:)==1)) = true;
        isFixedDof(edge(isDirichlet,:)) = true;
        isFixedDof(N + find(isDirichlet')) = true;
        fixedDof = find(isFixedDof);
        freeDof = find(~isFixedDof);
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction Dphi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function s = Dphi(p,i) % gradient of basis phi
    switch i
        case 1
            s = (4*lambda(p,1)-1).*Dlambda(:,:,1);
        case 2
            s = (4*lambda(p,2)-1).*Dlambda(:,:,2);
        case 3
            s = (4*lambda(p,3)-1).*Dlambda(:,:,3);
        case 4
            s = 4*(lambda(p,2)*Dlambda(:,:,3)+lambda(p,3)*Dlambda(:,:,2));
        case 5
            s = 4*(lambda(p,3)*Dlambda(:,:,1)+lambda(p,1)*Dlambda(:,:,3));
        case 6
            s = 4*(lambda(p,1)*Dlambda(:,:,2)+lambda(p,2)*Dlambda(:,:,1));
    end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction Dphi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function s = phi(p,i) % gradient of basis phi
    switch i
        case 1
            s = (2*lambda(p,1)-1).*lambda(p,1);
        case 2
            s = (2*lambda(p,2)-1).*lambda(p,2);
        case 3
            s = (2*lambda(p,3)-1).*lambda(p,3);
        case 4
            s = 4*lambda(p,3).*lambda(p,2);
        case 5
            s = 4*lambda(p,1).*lambda(p,3);
        case 6
            s = 4*lambda(p,1).*lambda(p,2);
    end
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Assemble right hand side by high order quadrature rule
% To reduce the effect of the error introduced by the numerical quadrature,
% the load term is computed using the 3rd order qudrature rule.
b = zeros(Ndof,1);
u=zeros(Ndof,1);
w=zeros(Ndof,1);

if ~isfield(option,'fquadorder')
    option.fquadorder = 3;   % default order
end
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if ~isempty(pde.f)
    % quadrature points in the barycentric coordinate
    [lambda,w] = quadpts(option.fquadorder);
    nQuad = size(lambda,1);
%     phi(:,1) = lambda(:,1).*(2*lambda(:,1)-1);
%     phi(:,2) = lambda(:,2).*(2*lambda(:,2)-1);
%     phi(:,3) = lambda(:,3).*(2*lambda(:,3)-1);
%     phi(:,4) = 4*lambda(:,2).*lambda(:,3);
%     phi(:,5) = 4*lambda(:,3).*lambda(:,1);
%     phi(:,6) = 4*lambda(:,1).*lambda(:,2);
    bt = zeros(NT,6);
    for p = 1:nQuad
        % quadrature points in the x-y coordinate
        pxy = lambda(p,1)*node(elem(:,1),:) ...
            + lambda(p,2)*node(elem(:,2),:) ...
            + lambda(p,3)*node(elem(:,3),:);
        if isfield(pde,'f') && isnumeric(pde.f)
            fp = pde.f;        % piecewise constant
        else
            fp = pde.f(pxy);   % function handle
        end
        for j = 1:6
            bt(:,j) = bt(:,j) + w(p)*phi(p,j)*fp;
        end
    end
    bt = bt.*repmat(area,1,6);
    b = accumarray(elem2dof(:),bt(:),[Ndof 1]);
end

 function [b1,u] = getbdP2(b)
    %% Boundary conditions for Poisson equation: P2 quadratic FEM.
    %
    % The set up of boundary condition consists of two parts:
    %

    %
    %  Modify the right hand side b. The Neumann boundary integral is added
    %  to b. For Dirichlet boundary ndoes, b(fixedDof) is the evaluation of
    % pde.g_D.
    %

    u = zeros(Ndof,1);

    %% Part 1: Find boundary edges and modify the load b
    % Neumann boundary condition
        if ~isfield(option,'gNquadorder')
            option.gNquadorder = 6;
        end
                idxN = (bdFlag(:) == 1);      % all Neumann edges in bdFlag
        Neumannidx = elem2edge(idxN ); % index of Neumann and Robin edges
        % since boundary integral is also needed for Robin edges
        Neumann   = edge(Neumannidx,:);
        b1 = zeros(Ndof,1);
        [lambdagN,weightgN] = quadpts1(option.gNquadorder);
        nQuadgN = size(lambdagN,1);
        % quadratic bases (1---3---2)
        bdphi = zeros(nQuadgN,3);
        bdphi(:,1) = (2*lambdagN(:,1)-1).*lambdagN(:,1);
        bdphi(:,2) = (2*lambdagN(:,2)-1).*lambdagN(:,2);
        bdphi(:,3) = 4*lambdagN(:,1).*lambdagN(:,2);
        % length of edge
        el = sqrt(sum((node(Neumann(:,1),:) - node(Neumann(:,2),:)).^2,2));
        ge = zeros(size(Neumann,1),3);
        for pp = 1:nQuadgN
            ppxy = lambdagN(pp,1)*node(Neumann(:,1),:) ...
                 + lambdagN(pp,2)*node(Neumann(:,2),:);
            gNu = pde.g_N(ppxy);
            ge(:,1) = ge(:,1) + weightgN(pp)*gNu*bdphi(pp,1);
            ge(:,2) = ge(:,2) + weightgN(pp)*gNu*bdphi(pp,2);
            ge(:,3) = ge(:,3) + weightgN(pp)*gNu*bdphi(pp,3); % interior bubble
        end
        % update RHS
        ge = ge.*repmat(el,1,3);
        b1(1:N) = accumarray(Neumann(:), [ge(:,1); ge(:,2)],[N,1]);
        b1(N+Neumannidx) = b1(N+Neumannidx) + ge(:,3);

    %% Part 2: Find Dirichlet boundary edges and compute the boundary value
    % Dirichlet boundary conditions

        isDirichlet(elem2edge(bdFlag(:)==1)) = true;
        % interpolation
                idx = (fixedDof > N);                         % index of edge nodes
        u(fixedDof(~idx)) = pde.g_D(node(fixedDof(~idx),:)); % bd value at vertex dofs
        bdEdgeIdx = fixedDof(idx) - N;
        bdEdgeMid = (node(edge(bdEdgeIdx,1),:) + node(edge(bdEdgeIdx,2),:))/2;
        u(fixedDof(idx)) = pde.g_D(bdEdgeMid);
        b1 = b1 - B*u;

    end % end of getbdP2

[b1,u] = getbdP2(b);
B(:,fixedDof)=[];
Nu=size(freeDof,1);

b(fixedDof)=[];
     bigA = [A, B; ...
             B', sparse(Nu,Nu)];
     bigF = [b1; -b];

% Solver
tic;

%bigU=PFGMRES(bigA, bigF,sparse(Ndof+Nu,1), Ndof, 10, 1e-6, [],[]);

   bigU=bigA\bigF;
   w=bigU(1:Ndof);
   u(freeDof)=bigU(Ndof+1:Ndof+Nu);
toc
end                 % end of function PoissonP2

repo: Aero-matlab


text:
function parse_data_int(filename, max_memo_GB)
% max_memo_GB is an idea of the number of GB allocated for the data to be
% stored in RAM, so it is used to compute the number of segments in which
% the data should be split for processing
% Based in read_intan_data.m version 1.1, June 26, 2010
% (c) 2010, Intan Technologies, LLC
% For more information, see http://www.intantech.com
% For updates and latest version, see http://www.intantech.com/software.html

[unused, fname, ext] = fileparts(filename);
ext = lower(ext(2:end));

if strcmp(ext,'.int')
    error('Incorrect extension.');
end

with_memory=true;
try
	memory;
catch
	with_memory=false;
end
if with_memory
	[userview,systemview] = memory;
	memo_avaible = floor(systemview.PhysicalMemory.Available*0.80);
	if exist('max_memo_GB','var') && ~isempty(max_memo_GB)
        max_memo = max_memo_GB*(1024)^3;
		if max_memo > memo_avaible
			error('max_memo_GB > 80% of Physical Memory Available')
		end
	else
		max_memo = memo_avaible;
	end
else
	max_memo = max_memo_GB*(1024)^3;
end

fid = fopen(filename, 'r');

% Read first three header bytes encoding file version
for i=1:3
    header(i) = fread(fid, 1, 'uint8');
end

if (header(1) ~= 128)
    error('Improper data file format.');
end

if (header(2) ~= 1 || header(3) ~= 1)
    warning('Data file version may not be compatible with this m-file.');
end

% Now see which amplifier channels are saved in this file.

for i=1:64
    amp_on(i) = fread(fid, 1, 'uint8');
end

num_amps = sum(amp_on);
channels = find(amp_on);
s = dir(filename);
filesize = s.bytes;

lts = (filesize - 67)/(num_amps*4 + 1);

sr = 25000;
save('intan_meta_data','sr','lts','channels')

%-----------------------------------

fprintf(1, '\nData file contains %0.2f seconds of data from %d amplifier channels.\n', lts/sr, num_amps);
fprintf(1, 'Channels: %s \n',num2str(channels));

%--------------------------------------

% Go back to the beginning of the file...
frewind(fid);
% ...skip the header this time...
fread(fid, 3+64, 'uint8');

samples_per_channel = ceil(max_memo/(num_amps*4+1));
num_segments = ceil(lts/samples_per_channel);

outfile_handles = cell(1,num_amps);

for i = 1:num_amps
    outfile_handles{i} = fopen([fname '_' num2str(channels(i)) '.intch'],'w','l');
end

for j=1:num_segments
    ini = (j-1)*samples_per_channel;
    fin = min(j*samples_per_channel,lts);
    data2 = fread(fid,(fin-ini)*(num_amps*4+1),'uint8=>uint8');
    data2((num_amps*4)+1:num_amps*4+1:end) = [];
    data2 = typecast(data2,'single');
    for ind = 1:num_amps
        fwrite(outfile_handles{ind},data2(ind:num_amps:end),'single');
    end
    fprintf('Segment %d out of %d processed.\n',j,num_segments);
end

fclose('all');

author: csn-le


text:
function rec=iEWT1D(ewt,mfb)

% ======================================================================
% function rec=iEWT1D(ewt,mfb)
%
% Perform the Inverse Empirical Wavelet Transform of ewt accordingly to
% the filter bank mfb
%
% Inputs:
%   -ewt: cell containing the EWT components
%   -mfb: filter bank used during the EWT
%
% Outputs:
%   -rec: reconstructed signal
%
% Author: Jerome Gilles
% Institution: UCLA - Department of Mathematics
% Year: 2012
% Version: 1.0
% ======================================================================

l=round(length(ewt{1})/2);

%We perform the adjoint operator to get the reconstruction
for k=1:length(ewt)
    ewt{k}=[ewt{k}(l-1:-1:1);ewt{k};ewt{k}(end:-1:end-l+1)];
    if k==1
        rec=zeros(length(ewt{1}),1);
    end
    rec=rec+real(ifft(fft(ewt{k}).*mfb{k}));
end
rec=rec(l:end-l);


text:
%%*****************************************************************
%% NTrhsfun: compute the right-hand side vector of the
%%           Schur complement equation for the NT direction.
%%*****************************************************************
%% SDPT3: version 4.0
%% Copyright (c) 1997 by
%% Kim-Chuan Toh, Michael J. Todd, Reha H. Tutuncu
%% Last Modified: 16 Sep 2004
%%*****************************************************************

    function [rhs,EinvRc,hRd] = NTrhsfun(blk,At,par,X,Z,rp,Rd,sigmu,hRd,dX,dZ);

    spdensity = par.spdensity;
    m = length(rp);
    if (nargin > 8)
       corrector = 1;
    else
       corrector = 0;
       hRd = zeros(m,1);
    end
    hEinvRc = zeros(m,1);
    EinvRc  = cell(size(blk,1),1);
    rhsfree = [];
%%
    for p = 1:size(blk,1)
        pblk = blk(p,:);
        n = sum(pblk{2});  numblk = length(pblk{2});
        if strcmp(pblk{1},'l')
	   if iscell(sigmu)
              EinvRc{p} = sigmu{p}./Z{p} -X{p};
	   else
              EinvRc{p} = sigmu./Z{p} -X{p};
           end
           Rq = sparse(n,1);
           if (corrector) & (norm(par.parbarrier{p})==0)
              Rq = dX{p}.*dZ{p}./Z{p};
           else
              tmp  = par.dd{p}.*Rd{p};
              tmp2 = mexMatvec(At{p},tmp,1);
              hRd = hRd + tmp2;
           end
	   EinvRc{p} = EinvRc{p} - Rq;
           tmp2 = mexMatvec(At{p},EinvRc{p},1);
           hEinvRc = hEinvRc + tmp2;
	elseif strcmp(pblk{1},'q')
	   if iscell(sigmu)
              EinvRc{p} = qops(pblk,-sigmu{p}./(par.gamz{p}.*par.gamz{p}),Z{p},4) -X{p};
	   else
              EinvRc{p} = qops(pblk,-sigmu./(par.gamz{p}.*par.gamz{p}),Z{p},4) -X{p};
           end
           Rq = sparse(n,1);
   	   if (corrector) & (norm(par.parbarrier{p})==0)
              w = sqrt(par.gamz{p}./par.gamx{p});
              hdx = qops(pblk,w,par.ff{p},5,dX{p});
              hdz = qops(pblk,w,par.ff{p},6,dZ{p});
              hdxdz = Arrow(pblk,hdx,hdz);
              vv = qops(pblk,w,par.ff{p},5,X{p});
              Vihdxdz = Arrow(pblk,vv,hdxdz,1);
              Rq = qops(pblk,w,par.ff{p},6,Vihdxdz);
           else
              tmp  = par.dd{p}.*Rd{p} + qops(pblk,qops(pblk,Rd{p},par.ee{p},1),par.ee{p},3);
              tmp2 = mexMatvec(At{p},tmp,1);
              hRd = hRd + tmp2;
           end
	   EinvRc{p} = EinvRc{p} - Rq;
           tmp2 = mexMatvec(At{p},EinvRc{p},1);
           hEinvRc = hEinvRc + tmp2;
        elseif strcmp(pblk{1},'s')
           n2 = pblk{2}.*(pblk{2}+1)/2;
	   if iscell(sigmu)
     	      %%ss = [0,cumsum(pblk{2})];
              %%sigmuvec = zeros(n,1);
              %%for k = 1:length(pblk{2});
              %%   sigmuvec(ss(k)+1:ss(k+1)) = sigmu{p}(k)*ones(pblk{2}(k),1);
              %%end
              sigmuvec = mexexpand(pblk{2},sigmu{p});
              tmp = spdiags(sigmuvec./par.sv{p} -par.sv{p},0,n,n);
           else
              tmp = spdiags(sigmu./par.sv{p} -par.sv{p},0,n,n);
           end
           EinvRc{p} = Prod3(pblk,par.G{p}',tmp,par.G{p},1);
           Rq = sparse(n,n);
           if (corrector) & (norm(par.parbarrier{p})==0)
              hdZ = Prod3(pblk,par.G{p},dZ{p},par.G{p}',1);
              hdX = spdiags(qops(pblk,par.parbarrier{p}',1./par.sv{p},3)-par.sv{p},0,n,n)-hdZ;
              tmp = Prod2(pblk,hdX,hdZ,0);
              tmp = 0.5*(tmp+tmp');
              if (numblk == 1)
                 d = par.sv{p};
                 e = ones(pblk{2},1);
                 Rq = 2*tmp./(d*e'+e*d');
                 if (nnz(Rq) <= spdensity*n2); Rq = sparse(Rq); end
              else
                 Rq = sparse(n,n);
                 ss = [0, cumsum(pblk{2})];
                 for i = 1:numblk
                    pos = [ss(i)+1 : ss(i+1)];
                    d = par.sv{p}(pos); e = ones(length(pos),1);
                    Rq(pos,pos) = 2*tmp(pos,pos)./(d*e' + e*d');
                 end
              end
              Rq = Prod3(pblk,par.G{p}',Rq,par.G{p},1);
           else
              tmp = Prod3(pblk,par.W{p},Rd{p},par.W{p},1,par.nzlistAy{p});
              tmp2 = AXfun(pblk,At(p,:),par.permA(p,:),{tmp});
              hRd = hRd + tmp2;
           end
           EinvRc{p} = EinvRc{p} - Rq;
           tmp2 = AXfun(pblk,At(p,:),par.permA(p,:),EinvRc(p));
           hEinvRc = hEinvRc + tmp2;
        elseif strcmp(pblk{1},'u')
           rhsfree = [rhsfree; Rd{p}];
        end
    end
%%
    rhs = rp + hRd - hEinvRc;
    rhs = full([rhs; rhsfree]);
%%*******************************************************************

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [MF, VF, KF, MC, VC] = garchfor2(data, resids, ht, kt, parameters, model,ar, ma, p, q, max_forecast)
%{
-----------------------------------------------------------------------
 PURPOSE:
 Mean, Volatility and Kurtosis Forecasting
-----------------------------------------------------------------------
 USAGE:
 [MF, VF, KF, MC, VC] = garchfor2(data, resids, ht, kt, parameters, model,ar, ma, p, q, max_forecast)

 INPUTS:
 data:          a vector of series
 resids:        a vector of residuals
 ht:            a vector of conditional variances
 kt:            a vector of conditional kurtosis
 parameters:    a vector of parameters
 model:         'GARCH', 'GJR', 'AGARCH', 'NAGARCH'
 ar:            positive scalar integer representing the order of AR
 am:            positive scalar integer representing the order of MA
 p:             positive scalar integer representing the order of ARCH
 q:             positive scalar integer representing the order of GARCH
 max_forecasts: maximum number of forecasts (i.e. 1-trading months 22 days)

 OUTPUTS:
 MF:            a vector of mean forecasts
 VF:            a vector of volatility forecasts
 KF:            a vector of kurtosis forecasts
 MC:            a vector of cumulative mean forecasts
 VC:            a vector of cumulative volatility forecasts
-----------------------------------------------------------------------
 Author:
 Alexandros Gabrielsen, a.gabrielsen@city.ac.uk
 Date:     08/2011
-----------------------------------------------------------------------
%}

if nargin == 0
    error('Data, Residuals, Variance, Kurtosis, GARCH Model, AR, MA, ARCH, GARCH, Maximum Number of Forecasts')
end

if size(data,2) > 1
   error('Data vector should be a column vector')
end

if (length(ar) > 1) | (length(ma) > 1) | ar < 0 | ma < 0
    error('AR and MA should be positive scalars')
end

if (length(p) > 1) | (length(q) > 1) | p < 0 | q < 0
    error('P and Q should be positive scalars')
end

z=ar+ma;
MF=[];
VF=[];
KF=[];

% Verifying that the vector of parameters is a column vector
[r,c]=size(parameters);
if r<c
    parameters=parameters';
end

% Forecasting the Mean
MF = parameters(1:1+z)'*[1; data(end-(1:ar)); resids(end-(0:ma-1))]; % 1-period ahead forecast
 for i = 2:max_forecast
     MF(i,1) = sum([parameters(1); ones(1,ar)*parameters(2:2+ar)*MF(i-1,1); ones(1,ma)*parameters(3+ar:2+ar+ma)*resids(end-(0:ma-1-i))]);
 end
 clear i

 % Forecasting Volatility and Kurtosis
 % Please note that for a given set of parameters it is possible to forecast
 % a smaller than 3 kurtosis, and therefore we floor the forecast at 3.
if isequal(strcat(model), 'GARCH')
    VF(1,1) = parameters(2+z:2+z+p+q)'*[1; resids(end-(0:(p-1))).^2;ht(end-(0:(q-1)))]; % 1-period ahead forecast
    KF(1,1) = parameters(3+z+p+q:3+z+2*p+2*q)'*[1; (resids(end-(0:(p-1))).^4)./(ht(end-(0:(p-1))).^2); kt(end-(0:(q-1)))]; % 1-period ahead forecast
    if KF(1,1) < 3
        KF(1,1) = 3
    end
    for i = 2:max_forecast
       VF(i,1) = parameters(2+z) + ones(1,p+q)*parameters(3+z:2+z+p+q)*VF(i-1,1);
       KF(i,1) = parameters(3+z+p+q) + ones(1,p+q)*parameters(4+z+p+q:3+z+2*p+2*q)*KF(i-1,1);
    end
clear i
elseif  isequal(strcat(model), 'GJR')
    VF(1,1) = parameters(2+z:2+z+2*p+q)'*[1; resids(end-(0:(p-1))).^2;  ht(end-(0:(q-1))); (resids(end-(0:(p-1)))<0).*resids(end-(0:(p-1))).^2];
    KF(1,1) = parameters(3+z+p+q:3+z+2*p+2*q)'*[1; (resids(end-(0:(p-1))).^4)./(ht(end-(0:(p-1))).^2); kt(end-(0:(q-1)))]; % 1-period ahead forecast
    if KF(1,1) < 3
        KF(1,1) = 3
    end
    for i = 2:max_forecast
         VF(i,1) = parameters(2+z) + ones(1,2*p+q)*parameters(3+z:2+z+2*p+q)*VF(i-1,1);
         KF(i,1) = parameters(3+z+p+q) + ones(1,p+q)*parameters(4+z+p+q:3+z+2*p+2*q)*KF(i-1,1);
    end
clear i
elseif isequal(strcat(model), 'AGARCH')
    asym=parameters(3+z+p+q:2+z+2*p+q);
    VF(1,1) = parameters(2+z:2+z+p+q)'*[1; (resids(end-(0:(p-1))) - asym).^2; ht(end-(0:(q-1)))]; % 1-period ahead forecast
    KF(1,1) = parameters(4+z+p+q:4+z+2*p+2*q)'*[1; (resids(end-(0:(p-1))).^4)./(ht(end-(0:(p-1))).^2); kt(end-(0:(q-1)))]; % 1-period ahead forecast
    if KF(1,1) < 3
        KF(1,1) = 3
    end
    for i = 2:max_forecast
       VF(i,1) = parameters(2+z:2+z+p+q)'*[1; (sqrt(VF(i-1,1)) - asym).^2; VF(i-1,1)];
       KF(i,1) = parameters(4+z+p+q) + ones(1,p+q)*parameters(5+z+p+q:4+z+2*p+2*q)*KF(i-1,1);
    end
clear i
elseif isequal(strcat(model), 'NAGARCH')
    asym=parameters(3+z+p+q:2+z+2*p+q);
    VF(1,1) = parameters(2+z:2+z+p+q)'*[1; (resids(end-(0:(p-1)))*sqrt(ht(end-(0:(p-1)))) - asym).^2; ht(end-(0:(p-1)))]; % 1-period ahead forecast
    KF(1,1) = parameters(4+z+p+q:4+z+2*p+2*q)'*[1; (resids(end-(0:(p-1))).^4)./(ht(end-(0:(p-1))).^2); kt(end-(0:(q-1)))]; % 1-period ahead forecast
    if KF(1,1) < 3
        KF(1,1) = 3
    end
    for i = 2:max_forecast
       VF(i,1) = parameters(2+z:2+z+p+q)'*[1; (sqrt(VF(i-1,1)) - asym).^2; VF(i-1,1)]
       KF(i,1) = parameters(4+z+p+q) + ones(1,p+q)*parameters(5+z+p+q:4+z+2*p+2*q)*KF(i-1,1);
    end
    clear i
else error('Invalid GARCH Model');
end

% Estimating cumulative measures
MC = cumsum(MF);
VC = sqrt(cumsum(VF));

end


text:
function sigma_n = sigma ( n )

%*****************************************************************************80
%
%% SIGMA returns the value of SIGMA(N), the divisor sum.
%
%  Definition:
%
%    SIGMA(N) is the sum of the distinct divisors of N, including 1 and N.
%
%  First values:
%
%     N  SIGMA(N)
%
%     1    1
%     2    3
%     3    4
%     4    7
%     5    6
%     6   12
%     7    8
%     8   15
%     9   13
%    10   18
%    11   12
%    12   28
%    13   14
%    14   24
%    15   24
%    16   31
%    17   18
%    18   39
%    19   20
%    20   42
%
%  Formula:
%
%    SIGMA(U*V) = SIGMA(U) * SIGMA(V) if U and V are relatively prime.
%
%    SIGMA(P^K) = ( P^(K+1) - 1 ) / ( P - 1 ) if P is prime.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 July 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the value to be analyzed.
%
%    Output, integer SIGMA_N, the value of SIGMA(N).  If N is less than
%    or equal to 0, SIGMA_N will be returned as 0.  If there is not
%    enough room for factoring N, SIGMA_N is returned as -1.
%
  maxfactor = 20;

  if ( n <= 0 )
    sigma_n = 0;
    return
  end

  if ( n == 1 )
    sigma_n = 1;
    return
  end
%
%  Factor N.
%
  [ nfactor, factor, power, nleft ] = i4_factor ( n );

  if ( nleft ~= 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SIGMA - Fatal error!\n' );
    fprintf ( 1, '  Not enough factorization space.\n' );
    sigma_n = -1;
    error ( 'SIGMA - Fatal error!' );
  end

  sigma_n = 1;
  for i = 1 : nfactor
    sigma_n = ( sigma_n * ( factor(i)^( power(i) + 1 ) - 1 ) ) ...
      / ( factor(i) - 1 );
  end

  return
end


text:
classdef GradientCirclePerimeterExperiment < handle

    properties (Access = private)
        radius
        nameCase
        inputFile
        scale
        imesh

        backgroundMesh
        boundaryMesh
        unfittedMesh

        domainLength
        levelSet

        regularizedPerimeter

        xNodesToPlot
        gradientInNodesToPlot

        outputFolder

    end

    methods (Access = public)

        function obj = GradientCirclePerimeterExperiment()
            obj.init();

            for imesh = 1:numel(obj.inputFile)
                obj.imesh = imesh;
                obj.createBackgroundAndBoundaryMesh();
                obj.computeDomainLength();
                obj.createLevelSet();
                obj.computeRegularizedPerimeters();
                obj.computeGradientSurfaces();
            end

        end

    end

    methods (Access = private)

        function init(obj)
            obj.radius = 0.2499999;
            obj.nameCase = 'GradientCirclePerimeterExperiment';
            obj.inputFile = {'SquareMacroTriangle'; ...
                'SquareMacroTriangleFine';...
                'SquareMacroTriangleFineFine'};
            obj.scale     = 'MACRO';
            obj.outputFolder = '/home/alex/Dropbox/Perimeter/';
        end

        function createBackgroundAndBoundaryMesh(obj)
            s.inputFile = obj.inputFile{obj.imesh};
            s.isBackgroundMeshRectangularBox = true;
            mCreator = BackgroundAndBoundaryMeshCreatorFromInputFile(s);
            obj.backgroundMesh = mCreator.backgroundMesh;
            obj.boundaryMesh   = mCreator.boundaryMesh;
        end

        function computeDomainLength(obj)
            x = obj.backgroundMesh.coord;
            d = max(x(:,1)) - min(x(:,1));
            obj.domainLength = d;
        end

        function createLevelSet(obj)
            s.type = 'circleInclusion';
            halfSide = obj.domainLength/2;
            s.fracRadius = obj.radius/halfSide;
            s.coord      = obj.backgroundMesh.coord;
            s.ndim       = obj.backgroundMesh.ndim;
            lsCreator = LevelSetCreator.create(s);
            obj.levelSet = lsCreator.getValue();
        end

        function computeGradientSurfaces(obj)
            s.outPutFolder = obj.outputFolder;
            s.mesh         = obj.backgroundMesh;
            s.rPerimeter   = obj.regularizedPerimeter;
            s.iMesh        = obj.imesh;
            s.domainLength = obj.domainLength;
            gComputer = GradientSurfPerimeterComputer(s);
            gComputer.compute();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan