text:
function [eigvector, eigvalue] = LGE(W, D, options, data)
% LGE: Linear Graph Embedding
%
%       [eigvector, eigvalue] = LGE(W, D, options, data)
%
%             Input:
%               data       - data matrix. Each row vector of data is a
%                         sample vector.
%               W       - Affinity graph matrix.
%               D       - Constraint graph matrix.
%                         LGE solves the optimization problem of
%                         a* = argmax (a'data'WXa)/(a'data'DXa)
%                         Default: D = I
%
%               options - Struct value in Matlab. The fields in options
%                         that can be set:
%
%                     ReducedDim   -  The dimensionality of the reduced
%                                     subspace. If 0, all the dimensions
%                                     will be kept. Default is 30.
%
%                            Regu  -  1: regularized solution,
%                                        a* = argmax (a'X'WXa)/(a'X'DXa+ReguAlpha*I)
%                                     0: solve the sinularity problem by SVD (PCA)
%                                     Default: 0
%
%                        ReguAlpha -  The regularization parameter. Valid
%                                     when Regu==1. Default value is 0.1.
%
%                            ReguType  -  'Ridge': Tikhonov regularization
%                                         'Custom': User provided
%                                                   regularization matrix
%                                          Default: 'Ridge'
%                        regularizerR  -   (nFea x nFea) regularization
%                                          matrix which should be provided
%                                          if ReguType is 'Custom'. nFea is
%                                          the feature number of data
%                                          matrix
%
%                            PCARatio     -  The percentage of principal
%                                            component kept in the PCA
%                                            step. The percentage is
%                                            calculated based on the
%                                            eigenvalue. Default is 1
%                                            (100%, all the non-zero
%                                            eigenvalues will be kept.
%                                            If PCARatio > 1, the PCA step
%                                            will keep exactly PCARatio principle
%                                            components (does not exceed the
%                                            exact number of non-zero components).
%
%
%             Output:
%               eigvector - Each column is an embedding function, for a new
%                           sample vector (row vector) x,  y = x*eigvector
%                           will be the embedding result of x.
%               eigvalue  - The sorted eigvalue of the eigen-problem.
%               elapse    - Time spent on different steps
%
%
%
%    Examples:
%
% See also LPP, NPE, IsoProjection, LSDA.
%
%
%Reference:
%
%   1. Deng Cai, Xiaofei He, Jiawei Han, "Spectral Regression for Efficient
%   Regularized Subspace Learning", IEEE International Conference on
%   Computer Vision (ICCV), Rio de Janeiro, Brazil, Oct. 2007.
%
%   2. Deng Cai, Xiaofei He, Yuxiao Hu, Jiawei Han, and Thomas Huang,
%   "Learning a Spatially Smooth Subspace for Face Recognition", CVPR'2007
%
%   3. Deng Cai, Xiaofei He, Jiawei Han, "Spectral Regression: A Unified
%   Subspace Learning Framework for Content-Based Image Retrieval", ACM
%   Multimedia 2007, Augsburg, Germany, Sep. 2007.
%
%   4. Deng Cai, "Spectral Regression: A Regression Framework for
%   Efficient Regularized Subspace Learning", PhD Thesis, Department of
%   Computer Science, UIUC, 2009.
%
%   version 3.0 --Dec/2011
%   version 2.1 --June/2007
%   version 2.0 --May/2007
%   version 1.0 --Sep/2006
%
%   Written by Deng Cai (dengcai AT gmail.com)
%

MAX_MATRIX_SIZE = 1600; % You can change this number according your machine computational power
EIGVECTOR_RATIO = 0.1; % You can change this number according your machine computational power

if (~exist('options','var'))
   options = [];
end

ReducedDim = 30;
if isfield(options,'ReducedDim')
    ReducedDim = options.ReducedDim;
end

if ~isfield(options,'Regu') || ~options.Regu
    bPCA = 1;
    if ~isfield(options,'PCARatio')
        options.PCARatio = 1;
    end
else
    bPCA = 0;
    if ~isfield(options,'ReguType')
        options.ReguType = 'Ridge';
    end
    if ~isfield(options,'ReguAlpha')
        options.ReguAlpha = 0.1;
    end
end

bD = 1;
if ~exist('D','var') || isempty(D)
    bD = 0;
end

[nSmp,nFea] = size(data);
if size(W,1) ~= nSmp
    error('W and data mismatch!');
end
if bD && (size(D,1) ~= nSmp)
    error('D and data mismatch!');
end

bChol = 0;
if bPCA && (nSmp > nFea) && (options.PCARatio >= 1)
    if bD
        DPrime = data'*D*data;
    else
        DPrime = data'*data;
    end
    DPrime = full(DPrime);
    DPrime = max(DPrime,DPrime');
    [R,p] = chol(DPrime);
    if p == 0
        bPCA = 0;
        bChol = 1;
    end
end

%======================================
% SVD
%======================================

if bPCA
    [U, S, V] = mySVD(data);
    [U, S, V]=CutonRatio(U,S,V,options);
    eigvalue_PCA = full(diag(S));
    if bD
        data = U*S;
        eigvector_PCA = V;

        DPrime = data'*D*data;
        DPrime = max(DPrime,DPrime');
    else
        data = U;
        eigvector_PCA = V*spdiags(eigvalue_PCA.^-1,0,length(eigvalue_PCA),length(eigvalue_PCA));
    end
else
    if ~bChol
        if bD
            DPrime = data'*D*data;
        else
            DPrime = data'*data;
        end

        switch lower(options.ReguType)
            case {lower('Ridge')}
                if options.ReguAlpha > 0
                    for i=1:size(DPrime,1)
                        DPrime(i,i) = DPrime(i,i) + options.ReguAlpha;
                    end
                end
            case {lower('Tensor')}
                if options.ReguAlpha > 0
                    DPrime = DPrime + options.ReguAlpha*options.regularizerR;
                end
            case {lower('Custom')}
                if options.ReguAlpha > 0
                    DPrime = DPrime + options.ReguAlpha*options.regularizerR;
                end
            otherwise
                error('ReguType does not exist!');
        end

        DPrime = max(DPrime,DPrime');
    end
end

WPrime = data'*W*data;
WPrime = max(WPrime,WPrime');

%======================================
% Generalized Eigen
%======================================

dimMatrix = size(WPrime,2);

if ReducedDim > dimMatrix
    ReducedDim = dimMatrix;
end

if isfield(options,'bEigs')
    bEigs = options.bEigs;
else
    if (dimMatrix > MAX_MATRIX_SIZE) && (ReducedDim < dimMatrix*EIGVECTOR_RATIO)
        bEigs = 1;
    else
        bEigs = 0;
    end
end

if bEigs
    %disp('use eigs to speed up!');
    option = struct('disp',0);
    if bPCA && ~bD
        [eigvector, eigvalue] = eigs(WPrime,ReducedDim,'la',option);
    else
        if bChol
            option.cholB = 1;
            [eigvector, eigvalue] = eigs(WPrime,R,ReducedDim,'la',option);
        else
            [eigvector, eigvalue] = eigs(WPrime,DPrime,ReducedDim,'la',option);
        end
    end
    eigvalue = diag(eigvalue);
else
    if bPCA && ~bD
        [eigvector, eigvalue] = eig(WPrime);
    else
        [eigvector, eigvalue] = eig(WPrime,DPrime);
    end
    eigvalue = diag(eigvalue);

    [junk, index] = sort(-eigvalue);
    eigvalue = eigvalue(index);
    eigvector = eigvector(:,index);

    if ReducedDim < size(eigvector,2)
        eigvector = eigvector(:, 1:ReducedDim);
        eigvalue = eigvalue(1:ReducedDim);
    end
end

if bPCA
    eigvector = eigvector_PCA*eigvector;
end

for i = 1:size(eigvector,2)
    eigvector(:,i) = eigvector(:,i)./norm(eigvector(:,i));
end

function [U, S, V]=CutonRatio(U,S,V,options)
    if  ~isfield(options, 'PCARatio')
        options.PCARatio = 1;
    end

    eigvalue_PCA = full(diag(S));
    if options.PCARatio > 1
        idx = options.PCARatio;
        if idx < length(eigvalue_PCA)
            U = U(:,1:idx);
            V = V(:,1:idx);
            S = S(1:idx,1:idx);
        end
    elseif options.PCARatio < 1
        sumEig = sum(eigvalue_PCA);
        sumEig = sumEig*options.PCARatio;
        sumNow = 0;
        for idx = 1:length(eigvalue_PCA)
            sumNow = sumNow + eigvalue_PCA(idx);
            if sumNow >= sumEig
                break;
            end
        end
        U = U(:,1:idx);
        V = V(:,1:idx);
        S = S(1:idx,1:idx);
    end

author: ZJULearning

repo: MatlabFunc


text:
function result = getMotionObject()
objectname = findobj('name','Stepper Motor Simulator');
result =get(objectname,'UserData');

end

author: Sable

repo: mcbench-benchmarks


text:
function R = sep(varargin);
% SEP - An Algorithm for Converting Covariance to Spherical Error Probable
%
% Version 1.0
%
% Copyright 2004, Michael Kleder
%
% USAGE:
%
% R = sep(sx,sy,sz)
% R = sep(sx,sy,sz,prob)
%
% This function computes Spherical Error Probable radius from inputs consisting of the
% square roots of the eigenvalues of a covariance matrix (equivalently, from sigma-x,
% sigma-y, and sigma-z, of a trivariate normal distribution in a coordinate system
% where there is no cross-correlation between variables.) This means that if you have
% a covariance matrix and wish to compute the S.E.P., simply obtain the square roots
% of the eigenvalues and use these as inputs. For example, list them via "sqrt(eig(C))"
% where C is your covariance matrix.
%
% The S.E.P. is the radius of a sphere which contains a fraction of probability equal
% to the input "prob," which is asumed to be 0.5 if omitted.
%
% Note: if one of the input sigmas is significantly smaller than both others,
% calculation time may rise.
%
% By uncommenting a labeled line of code, the user can enter a diagnostic mode
% to verify the accuracy of this algorithm for whatever inputs are specified.
%
% The mathematical formulas contained herein were created by the author and
% are copyrighted. Feel free to use them provided you credit the author:
% Kleder, Michael. "An Algorithm for Converting Covariance to Spherical Error Probable"
% Mathworks Central File Exchange, 2004.

[sx,sy,sz]=deal(varargin{1:3});
target = .5;
if nargin>3
    target=varargin{4};
end
lowend=0;
highend=10*max([sx sy sz]);
R = fminbnd(@enclosederr,lowend,highend,[],sx,sy,sz,target);
% diagnostic(R,sx,sy,sz) % use this to verify results against RANDN
return
function enclerr = enclosederr(R,sx,sy,sz,target)
enclerr=dblquad(@dens,0,pi,0,2*pi,[],[],R,sx,sy,sz);
enclerr=(target-enclerr)^2;
return
function outp = dens(th,phi,R,sx,sy,sz)
sinth  = sin(th);
sinphi = sin(phi);
costh  = cos(th);
cosphi = cos(phi);
subexprA = sinth.^2.*cosphi.^2.*sy.^2.*sz.^2+sinth.^2.*...
    sinphi.^2.*sx.^2.*sz.^2+costh.^2.*sx.^2.*sy.^2;
subexprB = sx.*sy.*sz;
outp = -1./8.*sinth.*2.^(.5).*subexprB.*(2.*R.*exp(-.5.*subexprA.*...
    R.^2./subexprB.^2).*subexprA.^(.5)-pi.^(.5).*2.^(.5).*erf(.5.*...
    2.^(.5)./subexprB.*subexprA.^(.5).*R).*subexprB)./pi.^(1.5)./subexprA.^(1.5);
return
function diagnostic(R,sx,sy,sz)
disp('Running test of SEP.')
for i=10:-1:1
    fprintf('%g ',i);
    x=randn(1e6,1)*sx;;
    y=randn(1e6,1)*sy;
    z=randn(1e6,1)*sz;
    r=sqrt(sum([x y z].^2,2));
    m(i) = sum(r<R)/length(r);
end
disp(' ');
m=mean(m);
disp([num2str(m*100) ' percent of random points were located within the computed SEP.'])
return


text:
function value = base_to_i4 ( s, base )

%*****************************************************************************80
%
%% BASE_TO_I4 returns the value of an integer represented in some base.
%
%  Discussion:
%
%    BASE = 1 is allowed, in which case we allow the digits '1' and '0',
%    and we simply count the '1' digits for the result.
%
%    Negative bases between -16 and -2 are allowed.
%
%    The base -1 is allowed, and essentially does a parity check on
%    a string of 1's.
%
%  Example:
%
%        Input      Output
%    -------------  ------
%         S   BASE       I
%    ------  -----  ------
%      '101'     2       5
%    '-1000'     3     -27
%      '100'     4      16
%   '111111'     2      63
%   '111111'    -2      21
%   '111111'     1       6
%   '111111'    -1       0
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string.  The elements of S are
%    blanks, a plus or minus sign, and digits.  Normally, the digits
%    are representations of integers between 0 and |BASE-1|.  In the
%    special case of base 1 or base -1, we allow both 0 and 1 as digits.
%
%    Input, integer BASE, the base in which the representation is given.
%    Normally, 2 <= BASE <= 16.  However, there are two exceptions.
%
%    Output, integer VALUE, the integer.
%
  value = 0;
  s_length = s_len_trim ( s );

  if ( base == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BASE_TO_I4 - Serious error!\n' );
    fprintf ( 1, '  The input base is zero.\n' );
    error ( 'BASE_TO_I4 - Serious error!' )
  end

  if ( 16 < abs ( base ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BASE_TO_I4 - Serious error!\n' );
    fprintf ( 1, '  The input base is greater than 16!\n' );
    error ( 'BASE_TO_I4 - Serious error!' )
  end

  state = 0;
  isgn = 1 ;
  ichr = 1;

  while ( ichr <= s_length )

    ch = s(ichr);
%
%  Blank.
%
    if ( ch == ' ' )

      if ( state == 2 )
        break
      end
%
%  Sign, + or -.
%
    elseif ( ch == '-' )

      if ( state ~= 0 )
        break
      end

      state = 1;
      isgn = -1;

    elseif ( ch == '+' )

      if ( state ~= 0 )
        break
      end

      state = 1;

    else
%
%  Digit?
%
      digit = hex_digit_to_i4 ( ch );

      if ( abs ( base ) == 1 && ( digit == 0 || digit == 1 ) )

        value = base * value + digit;
        state = 2;

      elseif ( 0 <= digit && digit < abs ( base ) )

        value = base * value + digit;
        state = 2;

      else

        fprintf ( 1, '\n' );
        fprintf ( 1, 'BASE_TO_I4 - Serious error!\n' );
        fprintf ( 1, '  Illegal digit = "%c"\n', ch );
        fprintf ( 1, '  Conversion halted prematurely!\n' );
        error ( 'BASE_TO_I4 - Serious error!' )

      end

    end

    ichr = ichr + 1;

  end
%
%  Once we're done reading information, we expect to be in state 2.
%
  if ( state ~= 2 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'BASE_TO_I4 - Serious error!\n' );
    fprintf ( 1, '  Unable to decipher input!\n' );
    error ( 'BASE_TO_I4 - Serious error!' )
  end
%
%  Account for the sign.
%
  value = isgn * value;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [cart] = Polar2Cart(polar)
 %% Cnverts polar coordinats [alpha;r] to cartesian [x;y]

    alpha = polar(1,:);
    r = polar(2,:);
    x = zeros(1,size(polar,2));
    y = zeros(1,size(polar,2));

    for i = 1:size(polar,2)
    x(i) = r(i)*cos(alpha(i));
    y(i) = r(i)*sin(alpha(i));
    end

    cart = [x;y];
end

author: Allopart

repo: rbpf-gmapping


text:
function xyzhumanevaGenerateMovie(Z,filename)

% XYZHUMANEVAGENERATEMOVIE
%
% COPYRIGHT : Carl Henrik Ek and Neil Lawrence, 2008

% MOCAP

if(nargin<2)
  filename = 'humaneva_movie';
end

% front view
mov = avifile(strcat(filename,'_front.avi'));
h = xyzhumanevaVisualise(Z(1,:),1);
F = getframe(gca);
mov = addframe(mov,F);
for(i = 2:1:size(Z,1))
  xyzhumanevaModify(h,Z(i,:));
  F = getframe(gca);
  mov = addframe(mov,F);
end
mov = close(mov);
clf(gcf,'reset');
% isometric view
mov = avifile(strcat(filename,'_3d.avi'));
h = xyzhumanevaVisualise3d(Z(1,:),1);
F = getframe(gca);
mov = addframe(mov,F);
for(i = 2:1:size(Z,1))
  xyzhumanevaModify(h,Z(i,:));
  F = getframe(gca);
  mov = addframe(mov,F);
end
mov = close(mov);

return

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function square_hex_grid_test02 ( )

%*****************************************************************************80
%
%% SQUARE_HEX_GRID_TEST02 tests HEX_GRID_01_H.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 January 2009
%
%  Author:
%
%    John Burkardt
%
  test_num = 17;

  nodes_per_layer_test = [ ...
    1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 21, 41, 81, 101, 1001, 10001 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SQUARE_HEX_GRID_TEST02\n' );
  fprintf ( 1, '  For a hexagonal grid of points in the unit square,\n' );
  fprintf ( 1, '  given NODES_PER_LAYER, the number of grid points\n' );
  fprintf ( 1, '  along the first layer,\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  HEX_GRID_01_H computes HX and HY, the spacings\n' );
  fprintf ( 1, '  in the row and column directions.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '    NODES    LAYERS   HX          HY\n' );
  fprintf ( 1, '      PER\n' );
  fprintf ( 1, '    LAYER\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num
    nodes_per_layer = nodes_per_layer_test ( test );
    layers = hex_grid_01_layers ( nodes_per_layer );
    [ hx, hy ] = hex_grid_01_h ( nodes_per_layer );
    fprintf ( 1, '  %6d  %6d  %10f  %10f\n', nodes_per_layer, layers, hx, hy );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  LAYERS is chosen so that LAYERS-1 layers just fit\n' );
  fprintf ( 1, '  inside the unit square, but LAYERS layers do not.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  LAYERS      HY     (LAYERS-1)*HY    LAYERS*HY\n' );
  fprintf ( 1, '\n' );

  for test = 1 : test_num
    nodes_per_layer = nodes_per_layer_test ( test );
    layers = hex_grid_01_layers ( nodes_per_layer );
    [ hx, hy ] = hex_grid_01_h ( nodes_per_layer );

    temp1 =  ( layers - 1 ) * hy;
    temp2 =  ( layers ) * hy;

    fprintf ( 1, '  %6d  %10f  %10f  %10f\n', ...
      nodes_per_layer, hy, temp1, temp2 );

  end

  return
end


text:
function [mCat]=sr_makeSRC4(mCat,fTw,T)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This function removes events and creates a quiescence according to the
% chosen input variables. The lat/lon position of the quiescence is
% determined randomly.
%
% Input Variables
% mCat      Earthquake catalog in zmap-format
% fTw       Duration of PSQ
% T         Starting time of PSQ
% fR        Degree of rate change (0.75 = 75% reduction)
%
% Output Variables
% mCat      Catalog with PSQ%
%
% van Stiphout Thomas ; vanstiphout@sed.ethz.ch
% Created: 14.08.2007
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('zmap/src/thomas/seismicrates/sr_makeSRC4.m')

% rate decrease for grid node
fR=0.75;
fLonPSQ=-116.4;
fLatPSQ=34.3;
nN=300;
fLon0=fLonPSQ;
fLat0=fLatPSQ;

% get N nearest events for certain grid point
%     [caNodeIndices, vResolution_] = ex_CreateIndexCatalog(mCat, [fLon0 fLat0], 1, 0, nN, 99, 0.1, 0.1);
%     % transform cell array to matrix
%     vSel=cell2mat(caNodeIndices);

% get events within certain radius of a grid point
[caNodeIndices, vResolution] = ex_CreateIndexCatalog(mCat, [fLon0 fLat0], 1, 1, ...
    100, 10, 0.1, 0.1);
% transform cell array to matrix
caNodeIndices_=cell2mat(caNodeIndices);
vSel=caNodeIndices_;
% transform cell array to matrix
vResolution_=cell2mat(vResolution);

% control plots
% figure;plot(mCat(:,1),mCat(:,2),'.k');
% hold on;plot(mCat(vSel,1),mCat(vSel,2),'dr');

% select from N events the one in the quiescence time period
vSel3=find((mCat(vSel,3)>T-fTw) & (mCat(vSel,3)<T));
% create help matrix for events in quite period
vSel4=zeros(size(vSel3,1),1);
% determine number of events to remove
fRemove=fR;
% create vector with % events to be removed (=1)
vSel4(1:ceil(size(vSel4,1)*fRemove))=1;
vSel4=logical(vSel4(randperm(size(vSel4,1))));
% Select indices of the N events to be removed
% vSel5 are events to be removed (indices from N events)
vSel5=vSel3(vSel4,1);
% create position of events in mCatalog to be removed
vSelNaN=vSel(vSel5);
% tag events to be removed with NaN
mCat(vSelNaN,1)=NaN;
mCat=mCat(~isnan(mCat(:,1)),:);

%% calculate rate increase in box

fBvalue=1;
fInc=0.1;
fMinLat=33.6;
fMaxLat=33.75;
fMinLon=-116.75;
fMaxLon=-116.5;
mPoly=[fMaxLon fMinLat;fMaxLon fMaxLat;fMinLon fMaxLat;fMinLon fMinLat];
fMinTime=T-fTw;
fMaxTime=T;

[vSel,bnd]=inpoly(mCat(:,1:2),mPoly);
vSel=find(vSel==1);
nNSel=size(vSel,1);
mNewCat=nan(nNSel,size(mCat,2));

% Create magnitudes
[mNewCat] = syn_create_magnitudes(mNewCat, fBvalue, min(mCat(:,6)), fInc);
% Randomize
rng('shuffle');

fMaxDepth=max(mCat(:,7));
fMinDepth=min(mCat(:,7));
% Create location
mNewCat(:,1) = rand(nNSel, 1) * (fMaxLon-fMinLon) + fMinLon;
mNewCat(:,2) = rand(nNSel, 1) * (fMaxLat-fMinLat) + fMinLat;
mNewCat(:,7) = rand(nNSel, 1) * (fMaxDepth-fMinDepth) + fMinDepth;

% Randomize
rng('shuffle');

% Create focal times
mNewCat(:,3) = rand(nNSel, 1) * (fMaxTime-fMinTime) + fMinTime;
% vst: datevec does not transform mNewCat(:,3) properly. Replaced by
% decyear2mat.
mNewCat(:,3)=mNewCat(randperm(size(mNewCat,1)),3);
% [mNewCat(:,10) mNewCat(:,4) mNewCat(:,5) mNewCat(:,8) mNewCat(:,9)] = datevec(mNewCat(:,3));
[mNewCat(:,10) mNewCat(:,4) mNewCat(:,5) mNewCat(:,8) mNewCat(:,9) tmp] = decyear2mat(mNewCat(:,3));

% Remove column 10 (seconds)
% mNewCat = mNewCat(:,1:end);
mCat=[mCat; mNewCat];
[Xsort, Xi]=sort(mCat(:,3),1);
mCat=mCat(Xi,:);

author: CelsoReyes


text:
function y = cvx_subsasgn( y, varargin )
temp.type = '()';
temp.subs = varargin(1:end-1);
y = subsasgn( y, temp, varargin{end} );

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

repo: REST


text:
classdef LaunchVehicleDataValidationWarning < AbstractLaunchVehicleValidatorOutput
    %LaunchVehicleDataValidation Summary of this class goes here
    %   Detailed explanation goes here

    properties
        str  char = ''
    end

    methods
        function obj = LaunchVehicleDataValidationWarning(str)
            obj.str = str;
        end

        function writeToLabel(obj, hLbl)
            styleLabelAsWarn(hLbl,obj.str)
        end

        function [str, style] = getUiTableStringAndRowStyle(obj)
            str = string(obj.str);
            style = uistyle('BackgroundColor',[1,1,0], 'FontColor',[0,0,0], 'FontWeight','bold');
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


%% BVAR tutorial: FORECASTS
% Author:   Filippo Ferroni
% Date:     27/02/2020

clear all
close all
clc

addpath ../../cmintools/
addpath ../../bvartools/
pkg load statistics
pkg load io

%% %=========================================================================
%%% DIRECT METHODS %%%
%%=========================================================================

load('../BVAR tutorial/DataGK')
y = [logip logcpi gs1 ebp];

%% 1/ Cholesky

lags               = 12;
options.hor        = 48;
options.conf_sig   = 0.9;
dm1 = directmethods(y,lags,options);

% Define the IRF of Interest
% index of the shocks of interest (shock to gs1)
indx_sho              = [3];
% Data order
% 1. logip; 2. logcpi; 3. gs1; 4. ebp
% Change the order of the variables for the plot
% 1. gs1; 2. logcpi; 3. logip; 4. ebp
indx_var              = [3, 2, 1, 4];

% PLOT IRF
% Customize the IRF plot
% variables names for the plots
options.varnames      = {'1 year rate','CPI','IP','EBP'};
% names of the directory where the figure is saved
options.saveas_dir    = './dm_plt';
% names of the figure to save
options.saveas_strng  = 'Cholesky';
% name of the shock
options.shocksnames   = {'MP'};
% Compare with BVAR estimates
bvar_ = bvar(y,lags,options);
var_irf_sort = sort(bvar_.ir_draws,4);
% add IRF plot
options.add_irfs(:,:,:,1) = var_irf_sort(indx_var,:,indx_sho,round(bvar_.ndraws*0.95));
options.add_irfs(:,:,:,2) = var_irf_sort(indx_var,:,indx_sho,round(bvar_.ndraws*0.05));
% plot LP
plot_irfs_(dm1.ir_lp(indx_var,:,indx_sho,:),options)

options =rmfield(options,'add_irfs');

%% 2/ IV
% load the instruments
[numi,txti,rawi] = xlsread('../BVAR tutorial/factor_data.csv','factor_data');
% instrument must have the same lenght as the observed data
options.proxy  = nan(length(y),1);
% use the same instrument as GK
% instruments and data end in 2012m6
options.proxy(length(T)- length(numi)+1:end) = numi(:,4);

dm2 = directmethods(y,lags,options);

options0= options;
options0.fontsize =18;
options0.saveas_strng  = 'IV';
options0.nplots = [1 1];
options0.varnames = {'IP','CPI','1 year rate','EBP'};
% finally, the plotting command
norm = dm2.irproxy_lp(3,1,1,2)*4;
plot_irfs_(dm2.irproxy_lp(:,:,1,:)/norm,options0)

options1 = options0;
for vv = 1 :size(y,2)
    dm3 = directmethods (y(:,vv),lags,options);
    options1.saveas_strng  = ['IV_var' num2str(vv)];
    options1.varnames = options0.varnames(vv);
    % finally, the plotting command
    plot_irfs_(dm3.irproxy_lp(:,:,1,:)/norm,options1)
end

%% 3/ Bayesian LP

% run a VAR on presample data
presample = 96; % 8 years of presample
lags      = 12;
bvar_     = bvar(y(1:presample,:),lags);

% use the VAR estimates to set the priors for the LP
options.priors.name        = 'Conjugate';
% posterior mean of the VAR AR coeff and constant
options.priors.Phi.mean    = mean(bvar_.Phi_draws,3);
% average variance of the AR coeff and constant
options.priors.Phi.cov     = diag(mean(var(bvar_.Phi_draws,0,3),2));
% posterior mean of the Covariance of the VAR residuals
options.priors.Sigma.scale = mean(bvar_.Sigma_draws,3);
options.priors.Sigma.df    = size(bvar_.Phi_draws,1)-2;
options.priors.tau         = 0.5*ones(options.hor); %

options.proxy(1:presample,:) =[];

bdm = directmethods(y(presample+1:end,:),lags,options);

options.saveas_strng  = 'BLPCholesky';
% the plotting command
options.conf_sig   = 0.9;
options.fontsize   = 12;
plot_irfs_(bdm.ir_blp(indx_var,:,indx_sho,:),options)

options.saveas_strng  = 'BLPIV';

norm = median(bdm.irproxy_blp(3,1,1,:),4)*4;
% the plotting command
plot_irfs_(bdm.irproxy_blp(indx_var,:,1,:)/norm,options)

%% 4/ Optimize Shrinkage

options.priors.max_tau = 1; %
options.max_compute    = 1; % fmin search
bdm_opt                = directmethods(y(presample+1:end,:),lags,options);

% the plotting command
options.saveas_strng  = 'BLPCholeskyOpt';
plot_irfs_(bdm_opt.ir_blp(indx_var,:,indx_sho,:),options)

% the plotting command
options.saveas_strng  = 'BLPIVOpt';
plot_irfs_(bdm_opt.irproxy_blp(indx_var,:,1,:)*0.25,options)


text:
function geometry_test1746 ( )

%*****************************************************************************80
%
%% TEST1746 tests R8MAT_INVERSE_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 May 2005
%
%  Author:
%
%    John Burkardt
%
  n = 3;

  a = [ ...
    3.0, 2.0, 0.0; ...
    2.0, 2.0, 1.0; ...
    1.0, 1.0, 1.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST1746\n' );
  fprintf ( 1, '  R8MAT_INVERSE_3D inverts a 3 by 3 matrix.\n' );

  r8mat_print ( n, n, a, '  Matrix A:' );

  [ b, det ] = r8mat_inverse_3d ( a );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Determinant of A is %f\n', det );

  r8mat_print ( n, n, b, '  Inverse matrix B:' );

  return
end


text:
classdef GmsVariablesAdapter < handle

    properties (Access = private)
       adaptedDataBase
    end

    properties (Access = private)
        nnode
        nfacet
        nelem
        nodeInfo
        connecInfo
        facetsInfo
    end

    methods (Access = public)

        function obj = GmsVariablesAdapter(readDataBase)
           obj.init(readDataBase);
        end

        function aDB = getAdaptedDB(obj)
            obj.obtainCoordinates();
            obj.obtainConnectivities();
            obj.obtainSetsInfo();
            obj.obtainMasterAndSlaveRelation();
            obj.obtainCorners();
            aDB = obj.adaptedDataBase;
        end

    end

    methods (Access = private)

        function init(obj,d)
            obj.nnode      = d.nnode;
            obj.nfacet     = d.nfacet;
            obj.nelem      = d.nelem;
            obj.nodeInfo   = d.nodeInfo;
            obj.connecInfo = d.connecInfo;
            obj.facetsInfo = d.facetsInfo;
        end

        function obtainCoordinates(obj)
            coordi = obj.nodeInfo;
            obj.adaptedDataBase.coord = zeros(obj.nnode,3);
            obj.adaptedDataBase.coord = coordi(:,1:2);
        end

        function obtainSetsInfo(obj)
            connect = obj.connecInfo;
            elementInSet = connect(:,4);
            sets = unique(elementInSet);
            nsets = size(sets,1);
            isElemInThisSet = cell(nsets,1);
            for iset = 1:nsets
                setNum = sets(iset);
                isElemInThisSet{iset} = find(elementInSet == setNum);
            end
            obj.adaptedDataBase.isElemInThisSet = isElemInThisSet;
        end

        function obtainConnectivities(obj)
            connect = obj.connecInfo;
            obj.adaptedDataBase.connec = connect(:,1:3);
        end

        function obtainMasterAndSlaveRelation(obj)
            coord = obj.adaptedDataBase.coord;
            msRelator = MasterSlaveRelator(coord);
            masterSlave = msRelator.getRelation();
            obj.adaptedDataBase.masterSlave = masterSlave;
        end

        function obtainCorners(obj)
            coord = obj.adaptedDataBase.coord;
            cD = CellNodesDescriptor(coord);
            corners = cD.cornerNodes;
            obj.adaptedDataBase.corners = corners;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function save_voxel_grid_as_bin(voxel,filename)

fp = fopen(filename,'w');
fwrite(fp,uint16(size(voxel)),'uint16');
fwrite(fp,uint8(voxel),'uint8');
fclose(fp);

author: lmb-freiburg

repo: orion

save_path: github-repos/MATLAB/lmb-freiburg-orion


text:
function [x, istop, itn, normr, normAr, normA, condA, normx]...
   = lsmr(A, b, lambda, atol, btol, conlim, itnlim, localSize, show)

% LSMR   Iterative solver for least-squares problems.
%   X = LSMR(A,B) solves the system of linear equations A*X=B. If the system
%   is inconsistent, it solves the least-squares problem min ||b - Ax||_2.
%   A is a rectangular matrix of dimension m-by-n, where all cases are
%   allowed: m=n, m>n, or m<n. B is a vector of length m.
%   The matrix A may be dense or sparse (usually sparse).
%
%   X = LSMR(AFUN,B) takes a function handle AFUN instead of the matrix A.
%   AFUN(X,1) takes a vector X and returns A*X. AFUN(X,2) returns A'*X.
%   AFUN can be used in all the following syntaxes.
%
%   X = LSMR(A,B,LAMBDA) solves the regularized least-squares problem
%      min ||(B) - (   A    )X||
%          ||(0)   (LAMBDA*I) ||_2
%   where LAMBDA is a scalar.  If LAMBDA is [] or 0, the system is solved
%   without regularization.
%
%   X = LSMR(A,B,LAMBDA,ATOL,BTOL) continues iterations until a certain
%   backward error estimate is smaller than some quantity depending on
%   ATOL and BTOL.  Let RES = B - A*X be the residual vector for the
%   current approximate solution X.  If A*X = B seems to be consistent,
%   LSMR terminates when NORM(RES) <= ATOL*NORM(A)*NORM(X) + BTOL*NORM(B).
%   Otherwise, LSMR terminates when NORM(A'*RES) <= ATOL*NORM(A)*NORM(RES).
%   If both tolerances are 1.0e-6 (say), the final NORM(RES) should be
%   accurate to about 6 digits. (The final X will usually have fewer
%   correct digits, depending on cond(A) and the size of LAMBDA.)
%   If ATOL or BTOL is [], a default value of 1.0e-6 will be used.
%   Ideally, they should be estimates of the relative error in the
%   entries of A and B respectively.  For example, if the entries of A
%   have 7 correct digits, set ATOL = 1e-7. This prevents the algorithm
%   from doing unnecessary work beyond the uncertainty of the input data.
%
%   X = LSMR(A,B,LAMBDA,ATOL,BTOL,CONLIM) terminates if an estimate
%   of cond(A) exceeds CONLIM. For compatible systems Ax = b,
%   conlim could be as large as 1.0e+12 (say).  For least-squares problems,
%   conlim should be less than 1.0e+8. If CONLIM is [], the default value
%   is CONLIM = 1e+8. Maximum precision can be obtained by setting
%   ATOL = BTOL = CONLIM = 0, but the number of iterations may then be
%   excessive.
%
%   X = LSMR(A,B,LAMBDA,ATOL,BTOL,CONLIM,ITNLIM) terminates if the
%   number of iterations reaches ITNLIM.  The default is ITNLIM = min(m,n).
%   For ill-conditioned systems, a larger value of ITNLIM may be needed.
%
%   X = LSMR(A,B,LAMBDA,ATOL,BTOL,CONLIM,ITNLIM,LOCALSIZE) runs LSMR
%   with reorthogonalization on the last LOCALSIZE v_k's (v-vectors
%   generated by the Golub-Kahan bidiagonalization). LOCALSIZE = 0 or []
%   runs LSMR without reorthogonalization. LOCALSIZE = Inf specifies
%   full reorthogonalization of the v_k's.  Reorthogonalizing only u_k or
%   both u_k and v_k are not an option here, because reorthogonalizing all
%   v_k's makes the u_k's close to orthogonal. Details are given in the
%   submitted SIAM paper.
%
%   X = LSMR(A,B,LAMBDA,ATOL,BTOL,CONLIM,ITNLIM,LOCALSIZE,SHOW) prints an
%   iteration log if SHOW=true. The default value is SHOW=false.
%
%   [X,ISTOP] = LSMR(A,B,...) gives the reason for termination.
%      ISTOP  = 0 means X=0 is a solution.
%             = 1 means X is an approximate solution to A*X = B,
%                 according to ATOL and BTOL.
%             = 2 means X approximately solves the least-squares problem
%                 according to ATOL.
%             = 3 means COND(A) seems to be greater than CONLIM.
%             = 4 is the same as 1 with ATOL = BTOL = EPS.
%             = 5 is the same as 2 with ATOL = EPS.
%             = 6 is the same as 3 with CONLIM = 1/EPS.
%             = 7 means ITN reached ITNLIM before the other stopping
%                 conditions were satisfied.
%
%   [X,ISTOP,ITN] = LSMR(A,B,...) gives ITN = the number of LSMR iterations.
%
%   [X,ISTOP,ITN,NORMR] = LSMR(A,B,...) gives an estimate of the residual
%   norm: NORMR = norm(B-A*X).
%
%   [X,ISTOP,ITN,NORMR,NORMAR] = LSMR(A,B,...) gives an estimate of the
%   residual for the normal equation: NORMAR = NORM(A'*(B-A*X)).
%
%   [X,ISTOP,ITN,NORMR,NORMAR,NORMA] = LSMR(A,B,...) gives an estimate of
%   the Frobenius norm of A.
%
%   [X,ISTOP,ITN,NORMR,NORMAR,NORMA,CONDA] = LSMR(A,B,...) gives an estimate
%   of the condition number of A.
%
%   [X,ISTOP,ITN,NORMR,NORMAR,NORMA,CONDA,NORMX] = LSMR(A,B,...) gives an
%   estimate of NORM(X).
%
%   LSMR uses an iterative method requiring matrix-vector products A*v
%   and A'*u.  For further information, see
%      D. C.-L. Fong and M. A. Saunders,
%      LSMR: An iterative algorithm for sparse least-squares problems,
%      SIAM J. Sci. Comput., submitted 1 June 2010.
%      See http://www.stanford.edu/~clfong/lsmr.html.
%
% 08 Dec 2009: First release version of LSMR.
% 09 Apr 2010: Updated documentation and default parameters.
% 14 Apr 2010: Updated documentation.
% 03 Jun 2010: LSMR with local and/or full reorthogonalization of v_k.
% 10 Mar 2011: Bug fix in reorthgonalization. (suggested by David Gleich)

% David Chin-lung Fong            clfong@stanford.edu
% Institute for Computational and Mathematical Engineering
% Stanford University
%
% Michael Saunders                saunders@stanford.edu
% Systems Optimization Laboratory
% Dept of MS&E, Stanford University.
%-----------------------------------------------------------------------

  % Initialize.

  if isa(A,'numeric')
    explicitA = true;
  elseif isa(A,'function_handle')
    explicitA = false;
  else
    error('SOL:lsmr:Atype','%s','A must be numeric or a function handle');
  end

  msg = ['The exact solution is  x = 0                              '
         'Ax - b is small enough, given atol, btol                  '
         'The least-squares solution is good enough, given atol     '
         'The estimate of cond(Abar) has exceeded conlim            '
         'Ax - b is small enough for this machine                   '
         'The least-squares solution is good enough for this machine'
         'Cond(Abar) seems to be too large for this machine         '
         'The iteration limit has been reached                      '];

  hdg1 = '   itn      x(1)       norm r    norm A''r';
  hdg2 = ' compatible   LS      norm A   cond A';
  pfreq  = 20;   % print frequency (for repeating the heading)
  pcount = 0;    % print counter

  % Determine dimensions m and n, and
  % form the first vectors u and v.
  % These satisfy  beta*u = b,  alpha*v = A'u.

  u    = b;
  beta = norm(u);
  if beta > 0
    u  = u/beta;
  end

  if explicitA
    v = A'*u;
    [m n] = size(A);
  else
    v = A(u,2);
    m = size(b,1);
    n = size(v,1);
  end

  minDim = min([m n]);

  % Set default parameters.

  if nargin < 3 || isempty(lambda)   , lambda    = 0;          end
  if nargin < 4 || isempty(atol)     , atol      = 1e-6;       end
  if nargin < 5 || isempty(btol)     , btol      = 1e-6;       end
  if nargin < 6 || isempty(conlim)   , conlim    = 1e+8;       end
  if nargin < 7 || isempty(itnlim)   , itnlim    = minDim;     end
  if nargin < 8 || isempty(localSize), localSize = 0;          end
  if nargin < 9 || isempty(show)     , show      = false;      end

  if show
    fprintf('\n\nLSMR            Least-squares solution of  Ax = b')
    fprintf('\nVersion 1.11                          09 Jun 2010')
    fprintf('\nThe matrix A has %8g rows  and %8g cols', m,n)
    fprintf('\nlambda = %16.10e', lambda )
    fprintf('\natol   = %8.2e               conlim = %8.2e', atol,conlim)
    fprintf('\nbtol   = %8.2e               itnlim = %8g'  , btol,itnlim)
  end

  alpha = norm(v);
  if alpha > 0
    v = (1/alpha)*v;
  end

  % Initialization for local reorthogonalization.

  localOrtho = false;
  if localSize > 0
    localPointer    = 0;
    localOrtho      = true;
    localVQueueFull = false;

    % Preallocate storage for the relevant number of latest v_k's.

    localV = zeros(n, min([localSize minDim]));
  end

  % Initialize variables for 1st iteration.

  itn      = 0;
  zetabar  = alpha*beta;
  alphabar = alpha;
  rho      = 1;
  rhobar   = 1;
  cbar     = 1;
  sbar     = 0;

  h    = v;
  hbar = zeros(n,1);
  x    = zeros(n,1);

  % Initialize variables for estimation of ||r||.

  betadd      = beta;
  betad       = 0;
  rhodold     = 1;
  tautildeold = 0;
  thetatilde  = 0;
  zeta        = 0;
  d           = 0;

  % Initialize variables for estimation of ||A|| and cond(A).

  normA2  = alpha^2;
  maxrbar = 0;
  minrbar = 1e+100;

  % Items for use in stopping rules.
  normb  = beta;
  istop  = 0;
  ctol   = 0;         if conlim > 0, ctol = 1/conlim; end;
  normr  = beta;

  % Exit if b=0 or A'b = 0.

  normAr = alpha * beta;
  if normAr == 0, disp(msg(1,:)); return, end

  % Heading for iteration log.

  if show
    test1 = 1;
    test2 = alpha/beta;
    fprintf('\n\n%s%s'      , hdg1 , hdg2   )
    fprintf('\n%6g %12.5e'  , itn  , x(1)   )
    fprintf(' %10.3e %10.3e', normr, normAr )
    fprintf('  %8.1e %8.1e' , test1, test2  )
  end

  %------------------------------------------------------------------
  %     Main iteration loop.
  %------------------------------------------------------------------
  while itn < itnlim
    itn = itn + 1;

    % Perform the next step of the bidiagonalization to obtain the
    % next beta, u, alpha, v.  These satisfy the relations
    %      beta*u  =  A*v  - alpha*u,
    %      alpha*v  =  A'*u - beta*v.

    if explicitA
      u = A*v    - alpha*u;
    else
      u = A(v,1) - alpha*u;
    end
    beta = norm(u);

    if beta > 0
      u = (1/beta)*u;
      if localOrtho
        localVEnqueue(v);    % Store old v for local reorthogonalization of new v.
      end
      if explicitA
        v = A'*u   - beta*v;
      else
        v = A(u,2) - beta*v;
      end

      if localOrtho
        v = localVOrtho(v);  % Local-reorthogonalization of new v.
      end
      alpha  = norm(v);
      if alpha > 0,  v = (1/alpha)*v; end
    end

    % At this point, beta = beta_{k+1}, alpha = alpha_{k+1}.

    % Construct rotation Qhat_{k,2k+1}.

    alphahat = norm([alphabar lambda]);
    chat     = alphabar/alphahat;
    shat     = lambda/alphahat;

    % Use a plane rotation (Q_i) to turn B_i to R_i.

    rhoold   = rho;
    rho      = norm([alphahat beta]);
    c        = alphahat/rho;
    s        = beta/rho;
    thetanew = s*alpha;
    alphabar = c*alpha;

    % Use a plane rotation (Qbar_i) to turn R_i^T to R_i^bar.

    rhobarold = rhobar;
    zetaold   = zeta;
    thetabar  = sbar*rho;
    rhotemp   = cbar*rho;
    rhobar    = norm([cbar*rho thetanew]);
    cbar      = cbar*rho/rhobar;
    sbar      = thetanew/rhobar;
    zeta      =   cbar*zetabar;
    zetabar   = - sbar*zetabar;

    % Update h, h_hat, x.

    hbar      = h - (thetabar*rho/(rhoold*rhobarold))*hbar;
    x         = x + (zeta/(rho*rhobar))*hbar;
    h         = v - (thetanew/rho)*h;

    % Estimate of ||r||.

    % Apply rotation Qhat_{k,2k+1}.
    betaacute =   chat* betadd;
    betacheck = - shat* betadd;

    % Apply rotation Q_{k,k+1}.
    betahat   =   c*betaacute;
    betadd    = - s*betaacute;

    % Apply rotation Qtilde_{k-1}.
    % betad = betad_{k-1} here.

    thetatildeold = thetatilde;
    rhotildeold   = norm([rhodold thetabar]);
    ctildeold     = rhodold/rhotildeold;
    stildeold     = thetabar/rhotildeold;
    thetatilde    = stildeold* rhobar;
    rhodold       =   ctildeold* rhobar;
    betad         = - stildeold*betad + ctildeold*betahat;

    % betad   = betad_k here.
    % rhodold = rhod_k  here.

    tautildeold   = (zetaold - thetatildeold*tautildeold)/rhotildeold;
    taud          = (zeta - thetatilde*tautildeold)/rhodold;
    d             = d + betacheck^2;
    normr         = sqrt(d + (betad - taud)^2 + betadd^2);

    % Estimate ||A||.
    normA2        = normA2 + beta^2;
    normA         = sqrt(normA2);
    normA2        = normA2 + alpha^2;

    % Estimate cond(A).
    maxrbar       = max(maxrbar,rhobarold);
    if itn>1
      minrbar     = min(minrbar,rhobarold);
    end
    condA         = max(maxrbar,rhotemp)/min(minrbar,rhotemp);

    % Test for convergence.

    % Compute norms for convergence testing.
    normAr  = abs(zetabar);
    normx   = norm(x);

    % Now use these norms to estimate certain other quantities,
    % some of which will be small near a solution.

    test1   = normr /normb;
    test2   = normAr/(normA*normr);
    test3   =      1/condA;
    t1      =  test1/(1 + normA*normx/normb);
    rtol    = btol + atol*normA*normx/normb;

    % The following tests guard against extremely small values of
    % atol, btol or ctol.  (The user may have set any or all of
    % the parameters atol, btol, conlim  to 0.)
    % The effect is equivalent to the normAl tests using
    % atol = eps,  btol = eps,  conlim = 1/eps.

    if itn >= itnlim,   istop = 7; end
    if 1 + test3  <= 1, istop = 6; end
    if 1 + test2  <= 1, istop = 5; end
    if 1 + t1     <= 1, istop = 4; end

    % Allow for tolerances set by the user.

    if  test3 <= ctol,  istop = 3; end
    if  test2 <= atol,  istop = 2; end
    if  test1 <= rtol,  istop = 1; end

    % See if it is time to print something.

    if show
      prnt = 0;
      if n     <= 40       , prnt = 1; end
      if itn   <= 10       , prnt = 1; end
      if itn   >= itnlim-10, prnt = 1; end
      if rem(itn,10) == 0  , prnt = 1; end
      if test3 <= 1.1*ctol , prnt = 1; end
      if test2 <= 1.1*atol , prnt = 1; end
      if test1 <= 1.1*rtol , prnt = 1; end
      if istop ~=  0       , prnt = 1; end

      if prnt
	if pcount >= pfreq
	  pcount = 0;
          fprintf('\n\n%s%s'    , hdg1 , hdg2  )
	end
	pcount = pcount + 1;
        fprintf('\n%6g %12.5e'  , itn  , x(1)  )
        fprintf(' %10.3e %10.3e', normr, normAr)
        fprintf('  %8.1e %8.1e' , test1, test2 )
        fprintf(' %8.1e %8.1e'  , normA, condA )
      end
    end

    if istop > 0, break, end
  end % iteration loop

  % Print the stopping condition.

  if show
    fprintf('\n\nLSMR finished')
    fprintf('\n%s', msg(istop+1,:))
    fprintf('\nistop =%8g    normr =%8.1e'     , istop, normr )
    fprintf('    normA =%8.1e    normAr =%8.1e', normA, normAr)
    fprintf('\nitn   =%8g    condA =%8.1e'     , itn  , condA )
    fprintf('    normx =%8.1e\n', normx)
  end

% end function lsmr

%---------------------------------------------------------------------
% Nested functions.
%---------------------------------------------------------------------

  function localVEnqueue(v)

  % Store v into the circular buffer localV.

    if localPointer < localSize
      localPointer = localPointer + 1;
    else
      localPointer = 1;
      localVQueueFull = true;
    end
    localV(:,localPointer) = v;

  end % nested function localVEnqueue

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

  function vOutput = localVOrtho(v)

  % Perform local reorthogonalization of V.

    vOutput = v;
    if localVQueueFull
      localOrthoLimit = localSize;
    else
      localOrthoLimit = localPointer;
    end
    for localOrthoCount = 1:localOrthoLimit
      vtemp   = localV(:, localOrthoCount);
      vOutput = vOutput - (vOutput'*vtemp)*vtemp;
    end

  end % nested function localVOrtho

end % function lsmr

author: RoyiAvital

repo: StackExchangeCodes


text:
function [rv,rvDebiased,rvSS,rvDebiasedSS,diagnostics] = realized_variance_optimal_sampling(price,time,timeType,samplingType,samplingInterval,subsamples,options)
% Estimates realized variance using Bandi-Russell optimal sampling selection
%
% USAGE:
%   [RV,RVD,RVSS,RVSSD,DIAGNOSTICS] = realized_variance_optimal_sampling(PRICE,TIME,TIMETYPE,SAMPLINGTYPE,SAMPLINGINVERVAL)
%   [RV,RVD,RVSS,RVSSD,DIAGNOSTICS] = realized_variance_optimal_sampling(PRICE,TIME,TIMETYPE,SAMPLINGTYPE,SAMPLINGINVERVAL,SUBSAMPLES,OPTIONS)
%
% INPUTS:
%   PRICE            - m by 1 vector of high frequency prices
%   TIME             - m by 1 vector of times where TIME(i) corresponds to PRICE(i)
%   TIMETYPE         - String describing the way times are measured
%                       'wall'    24-hour clock of the form HHMMSS, e.g. 101543 or 153217
%                       'seconds' Time measured in seconds past midnight on the first day.
%                       'unit'  Unit normalized date format, e.g. .1, .234, .9
%                         Unit normalized times are more general than the other types and can be
%                         applied to data from more than one calendar day
%   SAMPLINGTYPE     - String describing the type of sampling to use when
%                        filtering PRICE
%                        'CalendarTime' - Sample in calendar time using observations separated by
%                          SAMPLINGINTERVAL seconds
%                        'CalendarUniform' - Sample in calendar time using SAMPLINGINTERVAL
%                          observations spread uniformly between TIME(1) and TIME(m)
%                        'BusinessTime' - Sample in business (tick) time using observation separated
%                          by SAMPLINGINTERVAL ticks
%                        'BusinessUniform' - Sample in business (tick) time using observations
%                          uniformly spaced in business time.
%                        'Fixed' - Sample at specific points in time. When using fixed,
%                          SAMPLINGINTERVAL must be a n by 1 vector of times with the same TIMETYPE
%                          as TIME (i.e. seconds if TIME is in seconds)
%   SAMPLINGINTERVAL  - Scalar integer or n by 1 vector whose meaning depends on SAMPLINGTYPE
%   SUBSAMPLES        - [OPTIONAL] Scalar integer indicating the number of subsample realized
%                         variance estimators to average with the original realized variance.
%                         Subsample realized variances are based on prices uniformly spaced between
%                         the times (Calendar sampling) or ticks (Business sampling).  SUBSAMPLES=1
%                         will compute a subsample realized variance using the mid-point of the
%                         price sample points, 2 will use 1/3 and 2/3, and so on. In general this
%                         number should be small so the subsample estimators will be "sparse". If
%                         the computed optimal sampling frequency is smaller than SUBSAMPLES, then
%                         SUBSAMPLES is set to the compute optimal sampling frequency.
%   OPTIONS           - [OPTIONAL] Option structure initialized by realized_options.
%                         See help realized_options for a description of fields.
%
% OUTPUTS:
%   RV                - Realized variance estimated using the Bandi-Russell estimator of the optimal
%                         sampling frequency
%   RVD               - Debiased realized variance estimated using the Bandi-Russell estimator of
%                         the optimal sampling frequency for a debiased estimator
%   RVSS              - Subsample RV using the computed optimal samplig frequency
%   RVSSD             - Debiased version of subsample RV using the computed optimal samplig frequency
%   DIAGNOSTICS       - Structure with fields
%                         OPTIMALSAMPLES          - Optimal number of samples estimated using the
%                                                   Bandi-Russell methodology
%                         OPTIMALSAMPLESDEBIASED  - Optimal number of samples estimated using the
%                                                   Bandi-Russell methodology for use with a
%                                                   debiased estimator.
%                         SAMPLES                 - Number of sampled used.  Should be equal to
%                                                   OPTIMALNUMBEROFSAMPLES unless larger than the
%                                                   number of prices or smaller than 2
%                         SAMPLESDEBIASED         - Number of sampled used in debiased estimator.
%                                                   Should be equal to
%                                                   OPTIMALNUMBEROFSAMPLESDEBIASED unless larger
%                                                   than the number of prices or smaller than 2
%                         NOISEVARIANCE           - Estimated noise variance
%                         DEBIASEDNOISEVARIANCE   - Estimated noise variance bias adjusted
%                         IQESTIMATE              - Estimated IQ
%
% COMMENTS:
%   This function estimates the optimal number of samples to use when computing realized variance
%   using the method of Bandi & Russell (2008).  The role of SAMPLINGTYPE and SAMPLINGINTERVAL are
%   to set the MAXIMUM frequency of returns.  For example, if the maximum frequency was to be 15
%   seconds, set SAMPLINGTYPE='CalendarTime' and SAMPLINGINTERVAL=15.  To use all of the data,
%   set SAMPLINGTYPE='BusinessTime' and SAMPLINGINTERVAL=1.
%
% EXAMPLE:
%  % Optimal sampling in using all data
%  rvos = realized_variance_optimal_sampling(PRICE)
%  % Optimal sampling in business time with subsampling
%  [rvos,rvosSS] = realized_variance_optimal_sampling(PRICE,TIME,'Wall','BusinessTime',1,10)
%  % Optimal sampling in business time, using every 15th trade
%  rvos = realized_variance_optimal_sampling(PRICE,TIME,'Wall','BusinessTime',15)
%  % Optimal sampling where the maximum sampling frequency is 5 seonds
%  rvos = realized_variance_optimal_sampling(PRICE,TIME,'Wall','CalendarTime',5)
%
%  See also REALIZED_VARIANCE, REALIZED_KERNEL, REALIZED_QUANTILE_VARIANCE, REALIZED_RANGE,
%  REALIZED_PRICE_FILTER, REALIZED_THRESHOLD_VARIANCE

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

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

switch nargin
    case 1
        time = linspace(wall2seconds(93000),wall2seconds(160000),length(price));
        timeType = 'seconds';
        samplingType = 'businessuniform';
        samplingInterval = 1;
        subsamples = 1;
        options = realized_options('optimal sampling');
    case 5
        subsamples = 1;
        options = realized_options('optimal sampling');
    case 6
        options = realized_options('optimal sampling');
    case 7
        if isempty(subsamples)
            subsamples = 1;
        end
        % Nothing
    otherwise
        error('5 to 7 inputs required.')
end

errorMessage = realized_variance_optimal_sampling_parameter_check(price,time,timeType,samplingType,samplingInterval,subsamples,options);
if ~isempty(errorMessage )
    error(errorMessage)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inserted to protect against inputing integer times
time = double(time);
% Filter the price
[filteredPrice,filteredTime] = realized_price_filter(price,time,timeType,samplingType,samplingInterval);

% Compute the number of prices
m = size(filteredPrice,1);

% Next need to estimate IQ using RV squared
[noiseVariance, debiasedNoiseVariance, IQEstimate] = realized_noise_estimate(filteredPrice, filteredTime, timeType, options);
% Need to double these since these are the variance of the error to returns, not to prices which is
% what is returns from realized_noise_esitmate
noiseVariance = 2*noiseVariance;
debiasedNoiseVariance = 2*debiasedNoiseVariance;

% Need to estimate the IQ,  Variance of Noise, and 4th moment of noise (Eq. (16))
optimalNumberOfSamples = round((IQEstimate/(noiseVariance)^2)^(1/3));
if optimalNumberOfSamples>m
    warning('oxfordRealized:excessiveLags','WARNING')
    SamplesUsed = m;
elseif optimalNumberOfSamples<2
    warning('oxfordRealized:tooFewLags','WARNING')
    SamplesUsed = 2;
else
    SamplesUsed = optimalNumberOfSamples;
end
% Compute the realized variance using BT since the prices have been
% filtered
actualSubsamples = min(subsamples,SamplesUsed);
[rv,rvSS] = realized_variance(filteredPrice, filteredTime, timeType,'BusinessTime',SamplesUsed,actualSubsamples);

filteredPriceNoise = realized_price_filter(filteredPrice, filteredTime, timeType, options.noiseVarianceSamplingType, options.noiseVarianceSamplingInterval);
returns = diff(log(filteredPriceNoise));
if options.useAdjustedNoiseCount
    n = sum(returns~=0);
else
    n = length(returns);
end
noise4thPower = sum(returns.^4)/n;
noiseVariance = sum(returns.^2)/n;
% FIXME: Need to change way denominator is estimated
% Eq. 18
optimalNumberOfSamplesDebiased = round((2*IQEstimate / (2*noise4thPower - 3*noiseVariance^2))^(1/2));
if optimalNumberOfSamplesDebiased >m
    warning('oxfordRealized:excessiveLags','WARNING')
    samplesDebiased = m;
elseif optimalNumberOfSamplesDebiased <2
    warning('oxfordRealized:tooFewLags','WARNING')
    samplesDebiased = 2;
else
    samplesDebiased = optimalNumberOfSamplesDebiased;
end

% TODO See if it is possible to do a debiased using 1-m*/m type change
actualSubsamples = min(subsamples,samplesDebiased);
[rvDebiased,rvDebiasedSS] = realized_variance(filteredPrice,filteredTime,timeType,samplingType,samplesDebiased,actualSubsamples);
rvDebiased = rvDebiased - optimalNumberOfSamplesDebiased * noiseVariance;
rvDebiasedSS = rvDebiasedSS - optimalNumberOfSamplesDebiased * noiseVariance;
rvDebiased = rvDebiased/(1-optimalNumberOfSamplesDebiased/m);
rvDebiasedSS = rvDebiasedSS/(1-optimalNumberOfSamplesDebiased/m);

diagnostics.m = length(filteredPrice) - 1;
diagnostics.optimalSamples = optimalNumberOfSamples;
diagnostics.samples = SamplesUsed;
diagnostics.optimalNumberOfSamplesDebiased = optimalNumberOfSamplesDebiased;
diagnostics.samplesDebiased = samplesDebiased;
diagnostics.noiseVariance = noiseVariance;
diagnostics.debiasedNoiseVariance = debiasedNoiseVariance;
diagnostics.IQEstimate = IQEstimate;

function errorMessage = realized_variance_optimal_sampling_parameter_check(price,time,timeType,samplingType,samplingInterval,subsamples,options)
% Support function for realized_variance_optimal_sampling that does input validation
%
% USAGE:
%   [ERRORMESSAGE,OPTIONS] = realized_variance_optimal_sampling_parameter_check(PRICE,TIME,TIMETYPE,SAMPLINGTYPE,SAMPLINGINVERVAL,OPTIONS)
%
% INPUTS:
%   See realized_variance_optimal_sampling
%
% OUTPUT:
%   ERRORMESSAGE - String containing a description of the error if one is detected.  Empty if no error.
%
% COMMENTS:
%   See realized_options for a description of the other fields in OPTIONS
%
%  See also REALIZED_VARIANCE_OPTIMAL_SAMPLING

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 5/1/2008
errorMessage = [];

if size(price,2)>size(price,1)
    price=price';
end
if size(price,2)>1
    errorMessage = 'PRICE must be a m by 1 vector.';
    return
end
if size(time,2)>size(time,1)
    time=time';
end
if any(diff(time)<0)
    errorMessage = 'TIME must be sorted and increasing';
    return
end
if size(time,2)>1 || length(time)~=length(price)
    errorMessage = 'TIME must be a m by 1 vector.';
    return
end

timeType=lower(timeType);
if ~ismember(timeType,{'wall','seconds','unit'})
    errorMessage = 'TIMETYPE must be one of ''wall'', ''seconds'' or ''unit''.';
    return;
end

samplingType=lower(samplingType);
if ~ismember(samplingType,{'calendartime','calendaruniform','businesstime','businessuniform','fixed'})
    errorMessage = ('SAMPLINGTYPE must be one of ''CalendarTime'', ''CalendarUniform'', ''BusinessTime'', ''BusinessUniform'' or ''Fixed''.');
    return;
end

if ~isempty(subsamples)
    if ~isscalar(subsamples) || subsamples<0 || floor(subsamples)~=subsamples
        errorMessage = 'SUBSAMPLES must be a non-negative scalar.';
        return
    end
end

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

% Options

% List of flat top kernels
flatTopKernelList = {'bartlett','twoscale','2ndorder','epanechnikov',...
    'cubic','multiscale','5thorder','6thorder','7thorder','8thorder','parzen',...
    'th1','th2','th5','th16'};

% List of non flat top kernels
nonFlatTopKernelList = {'nonflatparzen','qs','fejer','thinf','bnhls'};

% Combined kernel list
kernelList = [flatTopKernelList nonFlatTopKernelList];

% Check fields for valid values
optionsFieldNames = fieldnames(options);
% Insert any missing fiedls
defaultOptions = realized_options('Optimal Sampling');
defaultFieldNames = fieldnames(defaultOptions);
missingFields = setdiff(defaultFieldNames,optionsFieldNames);
if ~isempty(missingFields)
    for i = 1:length(missingFields)
        options.(missingFields{i}) = defaultOptions.(missingFields{i});
    end
end

for i=1:length(optionsFieldNames)
    fieldName = optionsFieldNames{i};
    fieldValue = options.(fieldName);
    if ischar(fieldValue)
        fieldValue = lower(fieldValue);
        options.(fieldName) = fieldValue;
    end

    switch fieldName
        case {'medFrequencyKernel'}
            % Member of kernelList
            if ~ismember(fieldValue,kernelList)
                errorMessage = ['OPTIONS.' fieldName ' must be one of the listed types.'];
                return
            end
        case {'medFrequencyBandwidth'}
            % Non-negative scalar
            if ~isempty(fieldValue) && ~isnonnegativescalar(fieldValue)
                errorMessage = ['OPTIONS.' fieldName ' must a non-negative scalar.'];
                return
            end
        case {'useDebiasedNoise','useAdjustedNoiseCount'}
            % Logical or scalar
            if ~islogical(fieldValue) && ~ismember(fieldValue,[0 1])
                errorMessage = 'OPTIONS.useDebiasedNoise must be a logical value.';
                return
            end
        case {'IQEstimationSamplingType','medFrequencySamplingType','noiseVarianceSamplingType'}
            % One of the sampling types
            if ~ismember(fieldValue,{'calendartime','calendaruniform','businesstime','businessuniform','fixed'})
                errorMessage = ['OPTIONS.' fieldName ' must be one of ''CalendarTime'', ''CalendarUniform'', ''BusinessTime'', ''BusinessUniform'' or ''Fixed''.'];
                return
            end
        case {'medFrequencySamplingInterval','noiseVarianceSamplingInterval','IQEstimationSamplingInterval'}
            % Non-negative scalar, less that 1 if timeType is unit
            if isempty(fieldValue) || ~isnonnegativescalar(fieldValue)
                errorMessage = ['OPTIONS.' fieldName ' must be a non-negative scalar between 0 and 1.'];
                return
            end
            if strcmp(timeType,'unit') && fieldValue>1
                errorMessage = ['OPTIONS.' fieldName ' must be less than 1 if TIMETYPE when ''unit''.'];
                return
            end
    end
end

function condition = isnonnegativescalar(x)
% Function that returns logical true if that input is a non-empty scalar >=0
condition = ~isempty(x) && isscalar(x) && x>=0;

author: bashtage

repo: mfe-toolbox


text:
function f = constructorTurbo(f, op, pref)
%   F = CONSTRUCTORTURBO(F, OP, PREF) takes the CHEBTECH F that was constructed
%   from the function handle OP using the preferences in PREF and tries to use
%   contour integrals to compute some additional higher-order coefficients to
%   high accuracy.  These coefficients are returned in the appropriate
%   positions in F.COEFFS.
%
%   This function is only meant to be called by the CHEBTECH constructor.
%
% References:
%
%   [1] Bornemann, F.  Accuracy and stability of computing high-order
%         derivatives of analytic functions by Cauchy integrals.  Found.
%         Comput. Math. 11 (2011), pp. 1-63.
%
%   [2] Wang, H. and Huybrechs, D.  Fast and accurate computation of Jacobi
%         expansion coefficients of analytic functions.  Technical Report
%         TW-645, K.U. Leuven, 2015.

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

% How many coefficients do we want to compute?
if ( ~isnan(pref.fixedLength) )
    n = pref.fixedLength;
else
    n = 2*length(f);
end

% Pick an ellipse over which to integrate.
rhoCheb = exp(abs(log(eps())) / length(f));
rho = rhoCheb^(2/3);

% Do the contour integrals.
c = chebCoeffsTurbo(op, rho, n);

% Assign the new coefficients.
if ( isreal(f) )
    f.coeffs = real(c);
elseif ( isreal(1i*f) )
    f.coeffs = imag(c);
else
    f.coeffs = c;
end

end

function c = chebCoeffsTurbo(f, rho, n)
%CHEBCOEFFSTURBO   Compute Chebyshev coefficients using contour integrals.
%   C = CHEBCOEFFSTURBO(F, RHO, N) computes the first N Chebyshev coefficients
%   of the holomorphic function represented by the function handle F using
%   Cauchy integrals over the Bernstein ellipse of parameter RHO.  F must be
%   vectorized and able to accept complex inputs.

K = 4*n;                                    % Number of quadrature nodes.
g = @(z) f((rho*z + 1./(rho*z))/2);         % Remap ellipse to unit circle.
z = exp(2*pi*1i*(0:1:(K - 1)).'/K);         % Roots of unity.

% Compute integrals with trap. rule and rescale to get coefficients.
c = bsxfun(@rdivide, fft(g(z))/K, rho.^(0:1:(K - 1)).');
c = [c(1,:) ; 2*c(2:n,:)];

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function UFlists
%UFLISTS create the web pages for each matrix list (group, name, etc.)
% Places the web pages in the matrices/ subdirectory of the current directory.
%
% Example:
%   UFlists
%
% See also UFget, UFlist

% Copyright 2006-2007, Timothy A. Davis

% create all the web pages for the lists
UFlist ('group') ;
UFlist ('name') ;
UFlist ('id') ;
UFlist ('dimension') ;
UFlist ('nnz') ;
UFlist ('symmetry') ;
UFlist ('type') ;

% do all the group pages, and the list of groups
index = UFget ;
nmat = length (index.nrows) ;

[ignore, i] = sort (index.Group) ;
g = index.Group (i) ;

% create the primary directory
[url topdir] = UFlocation ;
matrices = [topdir 'matrices'] ;
if (~exist (matrices, 'dir'))
    mkdir (matrices) ;
end

% create the groups.html file
f = fopen ([matrices filesep 'groups.html'], 'w') ;
if (f < 0)
    error ('unable to create groups.html file\n') ;
end

% add the header
fprintf (f, ...
	'<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">\n') ;
fprintf (f, '\n') ;
fprintf (f, '<meta http-equiv="content-type" content="text/html; charset=') ;
fprintf (f, 'iso-8859-1"><title>UF Sparse Matrix Collection: group list') ;
fprintf (f, '</title>\n') ;
fprintf (f, '\n') ;
fprintf (f, '
UF Sparse Matrix Collection\n') ;

fprintf (f, '
List of matrix groups in the UF Sparse Matrix Collection:\n') ;
fprintf (f, '

\n') ;
fprintf (f, 'Group	\n') ;
fprintf (f, '# matrices	\n') ;
fprintf (f, 'details	\n') ;
fprintf (f, 'description	\n') ;

% find all the groups
group = '' ;
ngroup = 0 ;
for i = 1:nmat
    if (strcmp (group, g {i}))
	continue
    end
    group = g {i} ;
    ngroup = ngroup + 1 ;
    groups {ngroup} = group ;						    %#ok
end

nmat = 0 ;
for i = 1:ngroup
    group = groups {i} ;
    UFlist ('group', group) ;
    fprintf (f, '
\n') ;

    % link to group
    fprintf (f, '%s	\n', group, group) ;

    % number of matrices
    d = dir ([topdir 'mat' filesep group filesep '*.mat']) ;
    nmat_group = size (d,1) ;
    fprintf (f, '%d	\n', nmat_group) ;
    nmat = nmat + nmat_group ;

    % link to README.txt file ("details")
    f2 = fopen ([topdir filesep 'mat' filesep group filesep 'README.txt'], 'r');
    if (f2 < 0)
	error (['no README file for group: ' group]) ;
    else
	s = fgets (f2) ;
	fprintf (f, ...
	    'details	\n', group) ;
    end

    % one-line description (first line of README.txt)
    fprintf (f, '%s	\n', s) ;
    fclose (f2);

    fprintf (f, '
\n') ;
end
fprintf (f, '
\n') ;

fprintf (f, '
Total number of matrices in UF Sparse Matrix Collection:') ;
fprintf (f, ' %d\n', nmat) ;
fprintf (f, '

Maintained by Tim Davis, last updated %s.', date) ;
fprintf (f, '
Matrix pictures by cspy, a MATLAB', ...
    url) ;
fprintf (f, ' function in the CSparse package.\n', ...
    url) ;
fprintf (f, '\n') ;
fprintf (f, '\n') ;
fclose (f) ;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef MemoryManager < handle

    methods (Access = public, Abstract)

        link(obj)
        allocateMemory(obj)
        freeSpareMemory(obj)
        transferData(obj)

    end

end


text:
function vout = Ptransp_solve(L, vin)
%Ptransp_solve  Auxi,ilary function for IR Tools
%
% vout = Ptransp_solve
%
% This computes vout = L'\vin.  If L is a function handle, then we use the
% user-supplied function, which is passed via options as 'RegMatrix'.
% Otherwise, we use the backslash operator.

% Silvia Gazzola, University of Bath
% Per Christian Hansen, Technical University of Denmark
% James G. Nagy, Emory University
% April, 2018.

% This file is part of the IR Tools package and is distributed under the
% 3-Clause BSD Licence. A separate license file should be provided as part
% of the package.

if isa(L, 'function_handle')
    transp_flag = 'transp';
    vout = L(vin, transp_flag);
else
    vout = (L')\vin;
end

repo: IRtools


text:
function u = fem1d_heat_implicit ( x_num, x, t, dt, k_fun, ...
  rhs_fun, bc_fun, element_num, element_node, quad_num, u_old )

%*****************************************************************************80
%
%% FEM1D_HEAT_IMPLICIT: finite element method, 1D heat, implicit time steps.
%
%  Discussion:
%
%    This program solves
%
%      dUdT - k * d2UdX2 = f(X,T)
%
%    over the interval [A,B] with boundary conditions
%
%      U(A,T) = UA(T),
%      U(B,T) = UB(T),
%
%    over the time interval [T0,T1] with initial conditions
%
%      U(X,T0) = U0(X)
%
%    and specified functions k(X,T) and f(X,T).
%
%    The code uses the finite element method to approximate the second derivative
%    in space, and an implicit backward Euler approximation to the first
%    derivative in time.
%
%    For a test function Vi, we write
%
%      dUdt Vi - k * d2Udx2 * Vi = f(x,t) * Vi
%
%      Int dUdt Vi - k d2Udx2 * Vi    = Int f(x,t) Vi
%      Int dUdt Vi + k dUdx   * dVidx = Int f(x,t) Vi
%
%    Take the backward Euler approximation to the derivative:
%
%      Int ( U(x,t) - U(x,t-dt) ) / dt Vi + k dUdx dVidx = Int ( f(x,t) Vi )
%
%    Now, assume the finite element projection:
%
%      U(x,t)    = sum ( uj     * Vj )
%      U(x,t-dt) = sum ( u_oldj * Vj )
%      F(x,t)    = sum ( fj     * Vj ):
%
%    Then the equation can be rewritten as:
%
%      Sum ( Int ( Vi Vj ) * ( uj - u_oldj ) / dt + k dVidx dVjdx * uj ) =
%        Sum ( Int ( Vi Vj * fj ) )
%
%    Carrying the old U to the right hand side, we have:
%
%      Sum ( Int ( Vi Vj ) * u / dt + k dVidx dVjdx * uj ) =
%        Sum ( Int ( Vi Vj * ( fj + u_oldj / dt ) ) )
%
%    or, written as a linear system for coefficients "u":
%
%      ( K + M/dT ) * u = b + c
%
%    where K is the usual stiffness matrix, M is the standard finite element
%    mass matrix, b is the usual right hand side, and c contains additions to
%    the right hand side from the time term.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer X_NUM, the number of points to use in the spatial dimension.
%
%    Input, real X(X_NUM,1), the coordinates of the nodes.
%
%    Input, real T, the current time.
%
%    Input, real DT, the size of the time step.
%
%    Input, real K_FUN(), a function to evaluate the heat conductivity.
%
%    Input, real RHS_FUN(), a function to evaluate the right hand side.
%
%    Input, real BC_FUN(), a function to set the Dirichlet conditions.
%
%    Input, integer ELEMENT_NUM, the number of elements.
%
%    Input, integer ELEMENT_NODE(2,ELEMENT_NUM), the nodes belonging to
%    each element.
%
%    Input, integer QUAD_NUM, the number of quadrature points to use.
%
%    Input, real U_OLD(X_NUM,1), the solution at time T - dt.
%
%    Output, real U(X_NUM,1), the solution at time T.
%
  [ a, b ] = assemble_fem ( x_num, x, element_num, element_node, ...
    quad_num, t, k_fun, rhs_fun );

  [ a, b ] = assemble_backward_euler ( x_num, x, element_num, ...
    element_node, quad_num, t, dt, u_old, a, b );

  [ a, b ] = assemble_bc ( x_num, x, t, bc_fun, a, b );

  u = a \ b;

  return
end


text:
function [ y, j, f ] = now_to_yjf_common ( )

%*****************************************************************************80
%
%% NOW_TO_YMDF_COMMON expresses the current date as a Common YJF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer Y, J, real F, the YJF date.
%
  c = clock ( );

  y1 = c(1);
  m1 = c(2);
  d1 = c(3);
  h1 = c(4);
  n1 = c(5);
  s1 = c(6);
  mu1 = 0.0;

  f1 = mu1;
  f1 = s1 + f1 / 1000.0;
  f1 = n1 + f1 / 60.0;
  f1 = h1 + f1 / 60.0;
  f1 = f1 / 24.0;

  [ y, j, f ] = ymdf_to_yjf_common ( y1, m1, d1, f1 );

  return
end


text:
function [X,y,isContinuous] = designMatrix(tbl,allTerms,responseVar,varargin)
% Extract a design matix from a linear model.
%  Dummy coded designs for categorical variables and zero-sum covariates for
%  contiunuous variables are returned as X
% INPUT
% tbl =  A table
% allTerms = A cell array of variable names to include in the
% design matrix ( e.g. {'ori','freq','ori:freq'} for two mains and an interaction)
% Parm/Value
% ZeroSumConstraint -  Toggle to apply the zero sum constraint to
%                   each categorical factor (This equates the marginal prior across terms in the
%                           factor; see Rouder et al) [true]
% treatAsRandom  - A char or a cell array of chars with factors that should be
%               treated as random (i.e. not fixed) effects.
%               [{}].
% forceCategorical - Logical to force each term to be treated as a
% categorical variable (useful for grouping variables like subject ID which
% may look continuous...)
% OUTPUT
% X = The complete design matrix. Cell array with one element per term.
% y = The response data
% isContinuous = Logical indicating which columns are continuous co-variates.
%
% BK - 2019.

p = inputParser;
p.addParameter('zeroSumConstraint',true,@islogical);
p.addParameter('treatAsRandom',{},@(x) ischar(x) || iscell(x));
p.addParameter('forceCategorical',false,@islogical);
p.parse(varargin{:});
treatAsRandom =p.Results.treatAsRandom;
if ischar(treatAsRandom);treatAsRandom = {treatAsRandom};end

nrAllTerms =numel(allTerms);
X = cell(1,nrAllTerms);
N = height(tbl);

%Options for modelutils.designmatrix . Its internal alg determines which
%vars are categorical quite well. But use categorical() in the data table
%to be sure (or chars/strings).
isCategorical = bf.internal.isCategorical(tbl);
if p.Results.forceCategorical
    isCategorical = true(size(isCategorical));
end
cateoricalOpts = {'model','linear','intercept',false,'DummyVarCoding','full','responseVar',responseVar,'CategoricalVars',isCategorical};
isContinuous = false(1,nrAllTerms);
for i=1:nrAllTerms
    if any(allTerms{i}==':')
        % An interaction term.
       % names = strsplit(allTerms{i},':')
        aName =extractBefore(allTerms{i},':');
        bName = extractAfter(allTerms{i},':');
        aCategorical  = bf.internal.isCategorical(tbl,aName);
        bCategorical = bf.internal.isCategorical(tbl,bName);
        bothCategorical = aCategorical &&  bCategorical;
        bothContinuous = ~aCategorical && ~bCategorical;
        if bothCategorical || p.Results.forceCategorical
            thisA = classreg.regr.modelutils.designmatrix(tbl,'PredictorVars',aName,cateoricalOpts{:});
            if ~ismember(aName,treatAsRandom) && p.Results.zeroSumConstraint
                thisA = bf.internal.zeroSumConstraint(thisA);
            end
            thisB = classreg.regr.modelutils.designmatrix(tbl,'PredictorVars',bName,cateoricalOpts{:});
            if ~ismember(bName,treatAsRandom) && p.Results.zeroSumConstraint
                thisB = bf.internal.zeroSumConstraint(thisB);
            end
            thisX = bf.internal.interaction(thisA,thisB);
        elseif bothContinuous
            thisX =  tbl.(aName).*tbl.(bName);
            thisX = thisX - mean(thisX);
            isContinuous(i) = true;
        else
            % Interaction between a categorical and a continuous covariate
            if aCategorical
                thisA = classreg.regr.modelutils.designmatrix(tbl,'PredictorVars',aName,'intercept',false,'model','linear','responseVar',tbl.(responseVar),'DummyVarCoding','full');
                if ismember(bName,allTerms)
                    % B is already included as a main effect, have to
                    % remove one level from A to avoid colinearity,
                     thisA = bf.internal.zeroSumConstraint(thisA);
%                    thisA = thisA(:,2:end); % Remove first category
                end
            else
                thisA = tbl.(aName);
            end
            if bCategorical
                thisB = classreg.regr.modelutils.designmatrix(tbl,'PredictorVars',bName,'intercept',false,'model','linear','responseVar',tbl.(responseVar),'DummyVarCoding','full');
                if ismember(aName,allTerms)
                    % A is already included as a main effect, have to
                    % remove one level from B to avoid colinearity,
                    thisB = bf.internal.zeroSumConstraint(thisB);
                    %thisB = thisB(:,2:end); % Remove first category
                end
            else
                thisB = tbl.(bName);
            end
            thisX  =  thisA.*thisB;
            thisX  = thisX-mean(thisX);
            isContinuous(i) = true;
        end
    else
        % A main term
        if  bf.internal.isCategorical(tbl,allTerms{i}) || p.Results.forceCategorical
            thisX = classreg.regr.modelutils.designmatrix(tbl,'PredictorVars',allTerms{i},cateoricalOpts{:});
            % Sum-to-zero contrasts that equates marginal priors across levels.
            if ~ismember(allTerms{i},treatAsRandom) && p.Results.zeroSumConstraint
                thisX = bf.internal.zeroSumConstraint(thisX);
            end
        else %Continuous
             thisX =  tbl.(allTerms{i});
             thisX  = thisX-mean(thisX);
            isContinuous(i)= true;
        end

    end
    X{i} =thisX; % Store for later use
end

if nargout>1
    y= tbl.(responseVar);
end
end


text:
function  s = strtrim_improve(s0)
%
% strtrim handles only a limited number of whitespace types.
% But does NOT handle 0 for example. This function turns all
% whitespaces including 0 to plain old spaces which strtrim
% does handle.
%
%
s = s0;
s(s<33) = ' ';     % Replace all non-printable chracters with spaces

% Feed the digestable new string to strtrim.
s = strtrim(s);


text:
figure
subplot(2,2,1)
hold on
plot(vTime,100 * vAggregateTFP_full,'linewidth',1.5,'linestyle','-','color',[.03,.24,.8])
plot(vTime,100 * vAggregateTFP_reduced,'linewidth',1.5,'linestyle','--','color',[.8,.24,.03])
set(gcf,'color','w')
xlim([vTime(1) vTime(end)])
grid on
title('TFP','interpreter','latex','fontsize',14)
ylabel('$\%$ deviation from s.s.','interpreter','latex')
legend('Full model','Reduced model','location','northeast')
hold off

subplot(2,2,2)
hold on
plot(vTime,100 * vAggregateOutput_full,'linewidth',1.5,'linestyle','-','color',[.03,.24,.8])
plot(vTime,100 * vAggregateOutput_reduced,'linewidth',1.5,'linestyle','--','color',[.8,.24,.03])
set(gcf,'color','w')
xlim([vTime(1) vTime(end)])
grid on
title('Output','interpreter','latex','fontsize',14)
hold off

subplot(2,2,3)
hold on
plot(vTime,100 * vAggregateConsumption_full,'linewidth',1.5,'linestyle','-','color',[.03,.24,.8])
plot(vTime,100 * vAggregateConsumption_reduced,'linewidth',1.5,'linestyle','--','color',[.8,.24,.03])
set(gcf,'color','w')
xlim([vTime(1) vTime(end)])
grid on
title('Consumption','interpreter','latex','fontsize',14)
ylabel('$\%$ deviation from s.s.','interpreter','latex')
xlabel('Quarters','interpreter','latex')
hold off

subplot(2,2,4)
hold on
plot(vTime,100 * vAggregateInvestment_full,'linewidth',1.5,'linestyle','-','color',[.03,.24,.8])
plot(vTime,100 * vAggregateInvestment_reduced,'linewidth',1.5,'linestyle','--','color',[.8,.24,.03])
set(gcf,'color','w')
xlim([vTime(1) vTime(end)])
grid on
title('Investment','interpreter','latex','fontsize',14)
xlabel('Quarters','interpreter','latex')
hold off

author: gregkaplan

repo: phact

save_path: github-repos/MATLAB/gregkaplan-phact


text:
function evalTG(pres)

iids = imgList('test');

dirname = fullfile('pb',pres,'tg');
ignore = mkdir(dirname);
for i = 1:numel(iids),
  iid = iids(i);
  fprintf(2,'Computing Pb for image %d/%d (iid=%d) using %s (%s)...\n',...
          i,numel(iids),iid,'TG',pres);
  pb = pbTG(imgRead(iid));
  imwrite(pb,fullfile(dirname,sprintf('%d.bmp',iid)),'bmp');
end

author: Cloud-CV

repo: object-proposals


text:
% KM_DEMO_KRLS_WIENER1 Wiener system identification using Kernel Recursive
% Least Square (KRLS) regression.
%
% This program implements a regression example similar to the channel
% estimation example published in
% S. Van Vaerenbergh, J. Via, and I. Santamaria. "A sliding-window
% kernel RLS algorithm and its application to nonlinear channel
% identification", 2006 IEEE International Conference on Acoustics, Speech,
% and Signal Processing (ICASSP), Toulouse, France, 2006.
%
% Author: Steven Van Vaerenbergh (steven *at* gtas.dicom.unican.es), 2010.
%
% This file is part of the Kernel Methods Toolbox for MATLAB.
% https://github.com/steven2358/kmbox

close all
clear

%% PARAMETERS

Ntrain = 1500;		% number of total train data points
Ntest = 200;		% number of test data points
Nswitch = 500;		% abrupt switch from model 1 to model 2 after N1 iterations
B1 = [1,.8668,-0.4764,0.2070]';	% model 1 linear filter
B2 = [1,-.8326,.6656,-.7153]';	% model 2 linear filter
f = @(x) tanh(x);			% Wiener system nonlinearity
SNR = 40;		% SNR in dB

pars.kernel.type = 'gauss';	% kernel type
pars.kernel.par = 2;			% kernel parameter (width in case of Gaussian kernel)
pars.M = 150;		% dictionary size

% % parameters for ALD-KRLS
% pars.algo = 'ald-krls';
% pars.thresh = 0.01;

% % parameters for SW-KRLS
% pars.algo = 'sw-krls';
% pars.c = 1E-4;

% parameters for KRLS-T
pars.algo = 'krls-t';
pars.lambda = .999;
pars.c = 1E-4;

d = 4;			% time-embedding

%% PROGRAM
tic

% generate data
fprintf('Generating Wiener system data...\n');
N = Ntrain+Ntest;
s = randn(N,1);	% Gaussian input, all data
s_mem = zeros(N,d);
for i = 1:d,
	s_mem(i:N,i) = s(1:N-i+1);	% time-embedding
end
s_train = s_mem(1:Ntrain,:);	% input train data, stored in columns
s_test = s_mem(Ntrain+1:Ntrain+Ntest,:);	% input test data, stored in columns

x1 = s_mem(1:Nswitch,:)*B1;
x2 = s_mem(Nswitch+1:Ntrain,:)*B2;
x = [x1;x2];
y = f(x);
vary = var(y);
noisevar = 10^(-SNR/10)*vary;
noise_train = sqrt(noisevar)*randn(Ntrain,1);
y_train = y + noise_train;				% noisy output train data

x_test1 = s_mem(Ntrain+1:Ntrain+Ntest,:)*B1;
x_test2 = s_mem(Ntrain+1:Ntrain+Ntest,:)*B2;
noise_test1 = sqrt(noisevar)*randn(Ntest,1);
noise_test2 = sqrt(noisevar)*randn(Ntest,1);
y_test1 = f(x_test1) + noise_test1;	% noisy output test data, model 1
y_test2 = f(x_test2) + noise_test2;	% noisy output test data, model 2

% apply KRLS
fprintf('Applying KRLS for Wiener system identification...\n');
vars = [];
MSE = zeros(Ntrain,1);
for i=1:Ntrain-1,
	if ~mod(i,Ntrain/10), fprintf('.'); end
	vars.t = i;

	% perform KRLS regression and get regression output of test signal
	switch pars.algo
		case 'ald-krls'
			vars = km_aldkrls(vars,pars,s_train(i,:),y_train(i));	% train
			y_est = km_aldkrls(vars,pars,s_test);				% evaluate
		case 'sw-krls'
			vars = km_swkrls(vars,pars,s_train(i,:),y_train(i));	% train
			y_est = km_swkrls(vars,pars,s_test);				% evaluate
		case 'krls-t'
			vars = km_krlst(vars,pars,s_train(i,:),y_train(i));	% train
			y_est = km_krlst(vars,pars,s_test);				% evaluate
		otherwise
			error('wrong algorithm')
	end

	% calcalate test error
	if i<=Nswitch, y_test = y_test1; else y_test = y_test2; 	end
	MSE(i) = mean((y_test-y_est).^2);
end
fprintf('\n');

toc
%% OUTPUT

fprintf('Mean MSE over last 500 steps = %.2d\n',mean(MSE(Ntrain-499:Ntrain)))

figure;semilogy(MSE)
title('MSE')

figure; hold on
plot(y_test,'b')
plot(y_est,'g')
legend({'Test output','Estimated output'})
title(sprintf('%s, m=%d',pars.algo,pars.M))

author: epfl-lasa


text:
function GetPlotRange(handles)
if isempty(handles.CurrentData), return; end
Current = GetCurrent(handles);
values=Current(:); values(isinf(values))=[]; values(isnan(values))=[];

if length(unique(values))>20 % it is a mask?
    values(~values)=[];
    [Min, Max] = range_outlier(values);
else
    Min=min(values);
    Max=max(values);
end

if (abs(Min - Max)<1e-3)
    Max = Max + 1;
end
if (Min > Max)
    temp = Min;
    Min = Max;
    Max = temp;
end

if (Min < 0)
    set(handles.MinSlider, 'Min',    1.5*Min);
else
    set(handles.MinSlider, 'Min',    0.5*Min);
end

if (Max < 0)
    set(handles.MaxSlider, 'Max',    0.5*Max);
else
    set(handles.MaxSlider, 'Max',    1.5*Max);
end
set(handles.MinSlider, 'Max',    Max);
set(handles.MaxSlider, 'Min',    Min);
set(handles.MinValue,  'String', Min);
set(handles.MaxValue,  'String', Max);
set(handles.MinSlider, 'Value',  Min);
set(handles.MaxSlider, 'Value',  Max);
guidata(findobj('Name','qMRLab'), handles);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function x = dipole_cdf_inv ( cdf, a, b )

%*****************************************************************************80
%
%% DIPOLE_CDF_INV inverts the Dipole CDF.
%
%  Discussion:
%
%    A simple bisection method is used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CDF, the value of the CDF.
%
%    Input, real A, B, the parameters of the PDF.
%    -1.0 <= B <= 1.0.
%
%    Output, real X, the corresponding argument of the CDF.
%
  it_max = 100;
  tol = 0.0001;
%
%  Take care of horrible input.
%
  if ( cdf <= 0.0 )
    x = -r8_huge ( )
    return
  elseif ( 1.0 <= cdf )
    x = r8_huge ( )
    return
  end
%
%  Seek X1 < X < X2.
%
  x1 = -1.0;

  while ( 1 )

    cdf1 = dipole_cdf ( x1, a, b );

    if ( cdf1 <= cdf )
      break
    end

    x1 = 2.0 * x1;

  end

  x2 = 1.0;

  while ( 1 )

    cdf2 = dipole_cdf ( x2, a, b );

    if ( cdf <= cdf2 )
      break
    end

    x2 = 2.0 * x2;

  end
%
%  Now use bisection.
%
  it = 0;

  while ( 1 )

    it = it + 1;

    x3 = 0.5 * ( x1 + x2 );
    cdf3 = dipole_cdf ( x3, a, b );

    if ( abs ( cdf3 - cdf ) < tol )
      x = x3;
      break
    end

    if ( it_max < it )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'DIPOLE_CDF_INV - Fatal error!\n' );
      fprintf ( 1, '  Iteration limit exceeded.\n' );
      error ( 'DIPOLE_CDF_INV - Fatal error!' );
    end

    if ( ( cdf3 <= cdf & cdf1 <= cdf ) | ( cdf <= cdf3 & cdf <= cdf1 ) )
      x1 = x3;
      cdf1 = cdf3;
    else
      x2 = x3;
      cdf2 = cdf3;
    end

  end

  return
end


text:
function [cellVOI] = dicomrt_loadvoi(rtstruct_filename)
% dicomrt_loadvoi(rtstruct_filename)
%
% Load Volumes Of Interests (VOIs) from dicom-rt export file.
%
% A=dicomrt_dvhcal('rtstruct_filename') load into A the VOIs
% extracted from the file 'rtstruct_filename'
%
% VOIs are stored in a cell array with the following structure:
%
%  -----------------------------
%  | [OAR 1] | [xyz contour 1] |
%  |         -------------------
%  |         | [xyz contour 2] |
%  |         -------------------
%  |         |     ...         |
%  |         -------------------
%  |         | [xyz contour n] |
%  -----------------------------
%  |   ...   |     ...         |
%  -----------------------------
%  | [OAR n] | [xyz contour 1] |
%  |         -------------------
%  |         | [xyz contour 2] |
%  |         -------------------
%  |         |     ...         |
%  |         -------------------
%  |         | [xyz contour n] |
%  -----------------------------
%
%
% See also dicomrt_loaddose, dicomrt_loadct, dicomrt_dvhcal
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Get info file

dictFlg = checkDictUse;
if dictFlg
    if isdeployed
        rtstruct = dicominfo(rtstruct_filename,'dictionary', fullfile(getCERRPath,'bin,ES - IPT4.1CompatibleDictionary.mat'));
    else
        rtstruct = dicominfo(rtstruct_filename,'dictionary', 'ES - IPT4.1CompatibleDictionary.mat');
    end
else
    rtstruct = dicominfo(rtstruct_filename);
end

%rtstruct=rtstruct_filename;
% Get number of VOIs
ROIContourSequence = fieldnames(rtstruct.ROIContourSequence);

% Define cell array
VOI=cell(size(ROIContourSequence,1),2);

% Progress bar
h = waitbar(0,'Loading progress:');
set(h,'Name','dicomrt_loadvoi: loading RTSTRUCT objects');

for i=1:size(ROIContourSequence,1) % for i=1:(number of VOIs) ...
    voilabel=getfield(rtstruct.StructureSetROISequence,char(ROIContourSequence(i)),'ROIName');
    VOI{i,1}=voilabel; % get VOI's name

    try
        ncont_temp=fieldnames(getfield(rtstruct.ROIContourSequence,char(ROIContourSequence(i)), ...
            'ContourSequence')); % get contour list per each VOI (temporary variable)
    catch
        warning(['ContourSequence not found for ROI: ', voilabel]);
        ncont_temp = [];
    end

    switch isempty(ncont_temp)
        case 0
            for j=1:size(ncont_temp,1) % for j=1:(number of contours) ...
                if j==1
                    VOI{i,2}=cell(size(ncont_temp,1),1);
                end
                try
                    NumberOfContourPoints=getfield(rtstruct.ROIContourSequence,char(ROIContourSequence(i)), ...
                        'ContourSequence', char(ncont_temp(j)),'NumberOfContourPoints');
                    ContourData=getfield(rtstruct.ROIContourSequence,char(ROIContourSequence(i)), ...
                        'ContourSequence',char(ncont_temp(j)),'ContourData');
                    x=dicomrt_mmdigit(ContourData(1:3:NumberOfContourPoints*3)*0.1,7);
                    y=dicomrt_mmdigit(ContourData(2:3:NumberOfContourPoints*3)*0.1,7);
                    z=dicomrt_mmdigit(ContourData(3:3:NumberOfContourPoints*3)*0.1,7);
                    VOI{i,2}{j,1}=cat(2,x,y,z); % this is the same as VOI{i,2}{j,1}=[x,y,z];
                end
            end
        case 1
            % set dummy values. This will be deleted later dugin the import
            NumberOfContourPoints=1;
            ContourData=[0,0,0];
            x=0;
            y=0;
            z=0;
            VOI{i,2}{1,1}=cat(2,x,y,z);
    end
    waitbar(i/size(ROIContourSequence,1),h);
    ncont_temp=[];
end
% VOI cell generation complete
% Store VOI in a cell array
cellVOI=cell(3,1);
cellVOI{1,1}=rtstruct;
cellVOI{2,1}=VOI;
cellVOI{3,1}=[];
% Close progress bar
close(h);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function y = UVtOmega(U,V,I,J,col);
% The values of U*V' in Omega, as a vector
% Omega is given by I,J, col

y = zeros(length(I), 1);
for k = 1:length(col)-1
    j = J(col(k)+1);
    Xj = U * V(j,:)';
    idx = [col(k)+1:col(k+1)];
    y(idx) = Xj(I(idx));
end

% for j = 1:length(col)-1
%     Xj = U * V(j,:)';
%     idx = [col(j)+1:col(j+1)];
%     y(idx) = Xj(I(idx));
% end

author: andrewssobral

repo: lrslibrary


text:
function [ a, b ] = p04_lim ( )

%*****************************************************************************80
%
%% P04_LIM returns the integration limits for problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 January 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A(2), B(2), the lower and upper limits of integration.
%
  a(1:2) = -1.0;
  b(1:2) =  1.0;

  return
end

author: johannesgerer

repo: jburkardt-m