text:
function view_aValue(lab1, re3) % autogenerated function wrapper
    % This routine is the equivalent to the view_bva.m but for calc_avalgrid.m
    % view_maxz plots the maxz LTA values calculated
    % with maxzlta.m or other similar values as a color map
    % needs re3, gx, gy, stri
    %
    % define size of the plot etc.
    %
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    if ~exist('Prmap','var') || isempty(Prmap)
        Prmap = nan(size(re3));
    end

    if isempty(name)
        name = '  ';
    end
    think
    report_this_filefun(mfilename('fullpath'));
    co = 'w';

    % Find out if figure already exists
    %
    bmap=findobj('Type','Figure','-and','Name','a-value-map');

    % Set up the Seismicity Map window Enviroment
    %
    if isempty(bmap)
        bmap = figure_w_normalized_uicontrolunits( ...
            'Name','a-value-map',...
            'NumberTitle','off', ...
            'NextPlot','new', ...
            'backingstore','on',...
            'Visible','off', ...
            'Position',[ (fipo(3:4) - [600 400]) ZmapGlobal.Data.map_len]);

        lab1 = 'a-value:';
        create_my_menu();

        re4 = re3;

        colormap(jet)
        tresh = nan;
        minpe = nan;
        Mmin = nan;

    end   % This is the end of the figure setup

    % Now lets plot the color-map of the z-value
    %
    figure(bmap);
    delete(findobj(bfig,'Type','axes'));
    % delete(sizmap);
    reset(gca)
    cla
    hold off
    watchon;
    set(gca,'visible','off','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','normal',...
        'FontWeight','normal','LineWidth',1.,...
        'Box','on','SortMethod','childorder')

    rect = [0.18,  0.10, 0.7, 0.75];
    rect1 = rect;

    % find max and min of data for automatic scaling
    %
    ZG.maxc = max(max(re3));
    ZG.maxc = fix(ZG.maxc)+1;
    ZG.minc = min(min(re3));
    ZG.minc = fix(ZG.minc)-1;

    % set values gretaer tresh = nan
    %
    % re4 = re3;
    % l = r > tresh;
    % re4(l) = nan(1,length(find(l)));
    % l = Prmap < minpe;
    % re4(l) = nan(1,length(find(l)));
    % l = MaxCMap <  Mmin;
    % re4(l) = nan(1,length(find(l)));

    % plot image
    %
    orient landscape
    %set(gcf,'PaperPosition', [0.5 1 9.0 4.0])

    axes('position',rect)
    hold on
    pco1 = pcolor(gx,gy,re3);

    axis([ min(gx) max(gx) min(gy) max(gy)])
    set(gca,'dataaspect',[1 cosd(mean(ZG.primeCatalog.Latitude,'omitnan')) 1]);
    hold on
    if sha == 'fl'
        shading flat
    else
        shading interp
    end
    % make the scaling for the recurrence time map reasonable
    if lab1(1) =='T'
        l = isnan(re3);
        re = re3;
        re(l) = [];
        caxis([min(re) 5*min(re)]);
    end
    if fre == 1
        caxis([fix1 fix2])
    end

    title([name ';  '   num2str(t0b) ' to ' num2str(teb) ],'FontSize',ZmapGlobal.Data.fontsz.s,...
        'Color','r','FontWeight','normal')

    xlabel('Longitude [deg]','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s)
    ylabel('Latitude [deg]','FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s)

    % plot overlay
    %
    hold on
    update(mainmap())
    ploeq = plot(ZG.primeCatalog.Longitude,ZG.primeCatalog.Latitude,'k.');
    set(ploeq,'Tag','eq_plot','MarkerSize',ZG.ms6,'Marker',ty,'Color',co,'Visible',vi)

    set(gca,'visible','on','FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','normal',...
        'FontWeight','normal','LineWidth',1.,...
        'Box','on','TickDir','out')

    h1 = gca;
    hzma = gca;

    % Create a colorbar
    %
    h5 = colorbar('horiz');
    set(h5,'Pos',[0.35 0.07 0.4 0.02],...
        'FontWeight','normal','FontSize',ZmapGlobal.Data.fontsz.s,'TickDir','out')

    rect = [0.00,  0.0, 1 1];
    axes('position',rect)
    axis('off')
    %  Text Object Creation
    txt1 = text(...
        'Color',[ 0 0 0 ],...
        'EraseMode','normal',...
        'Units','normalized',...
        'Position',[ 0.2 0.06 0 ],...
        'HorizontalAlignment','right',...
        'Rotation',[ 0 ],...
        'FontSize',ZmapGlobal.Data.fontsz.s,....
        'FontWeight','normal',...
        'String',lab1);

    % Make the figure visible
    %
    set(gca,'FontSize',ZmapGlobal.Data.fontsz.s,'FontWeight','normal',...
        'FontWeight','normal','LineWidth',1.,...
        'Box','on','TickDir','out')
    figure(bmap);
    %sizmap = signatur('ZMAP','',[0.01 0.04]);
    %set(sizmap,'Color','k')
    axes(h1)
    set(gcf,'color','w');
    watchoff(bmap)
    %whitebg(gcf,[ 0 0 0 ])
    done

    %% ui functions
    function create_my_menu()
        add_menu_divider();

        add_symbol_menu('eq_plot');

        options = uimenu('Label',' Select ');
        uimenu(options,'Label','Refresh ', 'callback',@callbackfun_001)
        uimenu(options,'Label','Select EQ in Circle',...
            'callback',@callbackfun_002)
        uimenu(options,'Label','Select EQ in Circle - Constant R',...
            'callback',@callbackfun_003)
        uimenu(options,'Label','Select EQ in Circle - Overlay existing plot',...
            'callback',@callbackfun_004)

        uimenu(options,'Label','Select EQ in Polygon -new ',...
            'callback',@callbackfun_005)
        uimenu(options,'Label','Select EQ in Polygon - hold ',...
            'callback',@callbackfun_006)

        % Menu 'Maps'
        op1 = uimenu('Label',' Maps ');
        % A-Value map calculated by the MaxLikelihoodA...
        uimenu(op1,'Label','a-value map ',...
            'callback',@callbackfun_007)
        % B-Value map (fixed b-value by input from calc_avalgrid.m (NEVER CALLED)
        uimenu(op1,'Label','b-value map ',...
            'callback',@callbackfun_008)
        % Magnitude of completeness calculated by MaxCurvature
        uimenu(op1,'Label','Magnitude of completness map ',...
            'callback',@callbackfun_009)
        % Resolution estimation by mapping the needed radius to cover ni
        % earthquakes
        uimenu(op1,'Label','Resolution map',...
            'callback',@callbackfun_010)
        % Earthquake density map
        uimenu(op1,'Label','Earthquake density map',...
            'callback',@callbackfun_011)
        % Mu-value of the normal CDF
        uimenu(op1,'Label','Mu-value of the normal CDF',...
            'callback',@callbackfun_012)
        %  Sigma-value of the normal CDF
        uimenu(op1,'Label','Sigma-value of the normal CDF',...
            'callback',@callbackfun_013)

        % Menu 'Display'
        add_display_menu(4);
    end

    %% callback functions

    function callbackfun_001(mysrc,myevt)
        % automatically created callback function from text
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        view_aValue(lab1, re3);
    end

    function callbackfun_002(mysrc,myevt)
        % automatically created callback function from text
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        h1 = gca;
        met = 'ni';
        ZG=ZmapGlobal.Data;
        ZG.hold_state=false;
        cirbva;
        watchoff(bmap);
    end

    function callbackfun_003(mysrc,myevt)
        % automatically created callback function from text
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        h1 = gca;
        met = 'ra';
        ZG=ZmapGlobal.Data;
        ZG.hold_state=false;
        cirbva;
        watchoff(bmap);
    end

    function callbackfun_004(mysrc,myevt)
        % automatically created callback function from text
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        h1 = gca;
        ZG.hold_state=true;
        cirbva;
        watchoff(bmap);
    end

    function callbackfun_005(mysrc,myevt)
        % automatically created callback function from text
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        cufi = gcf;
        ZG.hold_state=false;
        selectp;
    end

    function callbackfun_006(mysrc,myevt)
        % automatically created callback function from text
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        cufi = gcf;
        ZG.hold_state=true;
        view_aValue(lab1, re3);
    end

    function callbackfun_007(mysrc,myevt)
        % automatically created callback function from text
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1 ='a-value';
        re3 = aValueMap;
        view_aValue(lab1, re3);
    end

    function callbackfun_008(mysrc,myevt)
        % automatically created callback function from text
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1='b-value';
        re3 = bValueMap;
        view_aValue(lab1, re3);
    end

    function callbackfun_009(mysrc,myevt)
        % automatically created callback function from text
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1 = 'Mcomp';
        re3 = MaxCMap;
        view_aValue(lab1, re3);
    end

    function callbackfun_010(mysrc,myevt)
        % automatically created callback function from text
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1='Radius in [km]';
        re3 = reso;
        view_aValue(lab1, re3);
    end

    function callbackfun_011(mysrc,myevt)
        % automatically created callback function from text
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1='log(EQ per km^2)';
        re3 = log10(ni./(reso.^2*pi));
        view_aValue(lab1, re3);
    end

    function callbackfun_012(mysrc,myevt)
        % automatically created callback function from text
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1='Mu-Value';
        re3 = MuMap;
        view_aValue(lab1, re3);
    end

    function callbackfun_013(mysrc,myevt)
        % automatically created callback function from text
        callback_tracker(mysrc,myevt,mfilename('fullpath'));
        lab1='Sigma-Value';
        re3 = SigmaMap;
        view_aValue(lab1, re3);
    end
end

author: CelsoReyes


text:
%GENREGSIN Generate sinusoidal regression data
%
%     X = GENDATSIN(N,SIGMA)
%
% INPUT
%    N      Number of objects to generate
%    SIGMA  Standard deviation of the noise
%
% OUTPUT
%    X      Regression dataset
%
% DESCRIPTION
% Generate an artificial regression dataset [X,Y] with:
%
%    y = sin(4x) + noise.
%
% where noise is Gaussian distributed with standard deviation sigma.
%
%  X = GENDATSIN(100)
% generates 100 x,y pairs with data and default noise (sigma = 0.1).
%
%  x = (0:0.01:1)';
%  y = genregsin(x,0);
% generates the true function along the x-axis, with zero noise.
%
% SEE ALSO (PRTools Guide)
% GENDATR, GENDATSINC

% Copyright: D.M.J. Tax, D.M.J.Tax@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands
function x = gendatsin(nrx,noise)
if nargin<2
  noise = 0.1;
end

if (length(nrx)>1)
  x = nrx;
  nrx = size(x);
  x = sin(4*x) + noise*randn(nrx);
else
  x = rand(nrx,1);
  y = sin(4*x) + noise*randn(nrx,1);
end
x = gendatr(x,y);

return

author: marianux

repo: ecg-kit

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


text:
function CstDefineFrequencydomainSolver(mws,startFreq,endFreq,samples)

ChangeSolverType = invoke(mws,'ChangeSolverType','HF Frequency Domain');

%'@ define frequency domain solver parameters
Mesh = invoke(mws,'Mesh');
invoke(Mesh,'SetCreator','High Frequency');
FDSolver = invoke(mws,'FDSolver');
invoke(FDSolver,'Reset');
invoke(FDSolver,'SetMethod','Tetrahedral','General purpose');
invoke(FDSolver,'OrderTet','Second');
invoke(FDSolver,'OrderSrf','First');
invoke(FDSolver,'Stimulation','All','All');
invoke(FDSolver,'ResetExcitationList');
invoke(FDSolver,'AutoNormImpedance','True');
invoke(FDSolver,'NormingImpedance','50');
invoke(FDSolver,'ModesOnly','False');
invoke(FDSolver,'ConsiderPortLossesTet','True');
invoke(FDSolver,'SetShieldAllPorts','False');
invoke(FDSolver,'AccuracyHex','1e-6');
invoke(FDSolver,'AccuracyTet','1e-4');
invoke(FDSolver,'AccuracySrf','1e-3');
invoke(FDSolver,'LimitIterations','False');
invoke(FDSolver,'MaxIterations','0');
invoke(FDSolver,'SetCalculateExcitationsInParallel','True','True','');
invoke(FDSolver,'StoreAllResults','False');
invoke(FDSolver,'StoreResultsInCache','False');
invoke(FDSolver,'UseHelmholtzEquation','True');
invoke(FDSolver,'LowFrequencyStabilization','True');
invoke(FDSolver,'Type','Auto');
invoke(FDSolver,'MeshAdaptionHex','False');
invoke(FDSolver,'MeshAdaptionTet','True');
invoke(FDSolver,'AcceleratedRestart','True');
invoke(FDSolver,'FreqDistAdaptMode','Distributed');
invoke(FDSolver,'NewIterativeSolver','True');
invoke(FDSolver,'TDCompatibleMaterials','False');
invoke(FDSolver,'ExtrudeOpenBC','False');
invoke(FDSolver,'SetOpenBCTypeHex','Default');
invoke(FDSolver,'SetOpenBCTypeTet','Default');
invoke(FDSolver,'AddMonitorSamples','True');
invoke(FDSolver,'CalcStatBField','False');
invoke(FDSolver,'CalcPowerLoss','True');
invoke(FDSolver,'CalcPowerLossPerComponent','False');
invoke(FDSolver,'StoreSolutionCoefficients','True');
invoke(FDSolver,'UseDoublePrecision','False');
invoke(FDSolver,'UseDoublePrecision_ML','True');
invoke(FDSolver,'MixedOrderSrf','False');
invoke(FDSolver,'MixedOrderTet','False');
invoke(FDSolver,'PreconditionerAccuracyIntEq','0.15');
invoke(FDSolver,'MLFMMAccuracy','Default');
invoke(FDSolver,'MinMLFMMBoxSize','0.3');
invoke(FDSolver,'UseCFIEForCPECIntEq','true');

invoke(FDSolver,'UseFastRCSSweepIntEq','false');
invoke(FDSolver,'UseSensitivityAnalysis','False');
% invoke(FDSolver,'SweepErrorThreshold','True','0.01');
invoke(FDSolver,'RemoveAllStopCriteria','Hex');
invoke(FDSolver,'AddStopCriterion','All S-Parameters','0.01','2','Hex','True');
invoke(FDSolver,'AddStopCriterion','Reflection S-Parameters','0.01','2','Hex','False');
invoke(FDSolver,'AddStopCriterion','Transmission S-Parameters','0.01','2','Hex','False');
invoke(FDSolver,'RemoveAllStopCriteria','Tet');

invoke(FDSolver,'AddStopCriterion','All S-Parameters','0.01','2','Tet','True');
invoke(FDSolver,'AddStopCriterion','Reflection S-Parameters','0.01','2','Tet','False');
invoke(FDSolver,'AddStopCriterion','Transmission S-Parameters','0.01','2','Tet','False');
% invoke(FDSolver,'All Probes','0.05','2','Tet','True');
invoke(FDSolver,'RemoveAllStopCriteria','Srf');

invoke(FDSolver,'AddStopCriterion','All S-Parameters','0.01','2','Srf','True');
invoke(FDSolver,'AddStopCriterion','Reflection S-Parameters','0.01','2','Srf','False');
invoke(FDSolver,'AddStopCriterion','Transmission S-Parameters','0.01','2','Srf','False');
invoke(FDSolver,'SweepMinimumSamples','3');
invoke(FDSolver,'SetNumberOfResultDataSamples','1001');

invoke(FDSolver,'SetResultDataSamplingMode','Automatic');
invoke(FDSolver,'SweepWeightEvanescent','1.0');
invoke(FDSolver,'AccuracyROM','1e-4');
invoke(FDSolver,'AddSampleInterval','','','1','Automatic','True');
invoke(FDSolver,'AddSampleInterval',num2str(startFreq),num2str(endFreq),num2str(samples),'Automatic','False');
% invoke(FDSolver,'AddSampleInterval','','','','Automatic','False');
invoke(FDSolver,'MPIParallelization','False');
invoke(FDSolver,'UseDistributedComputing','False');
invoke(FDSolver,'NetworkComputingStrategy','RunRemote');
invoke(FDSolver,'NetworkComputingJobCount','3');
invoke(FDSolver,'LimitCPUs','True');
invoke(FDSolver,'MaxCPUs','96');
invoke(FDSolver,'MaximumNumberOfCPUDevices','2');

IESolver = invoke(mws,'IESolver');
invoke(IESolver,'Reset');
invoke(IESolver,'UseFastFrequencySweep','True');
invoke(IESolver,'UseIEGroundPlane','False');
% invoke(IESolver,'SetRealGroundMaterialName','');
invoke(IESolver,'CalcFarFieldInRealGround','False');
invoke(IESolver,'RealGroundModelType','Auto');
invoke(IESolver,'PreconditionerType','Auto');
invoke(IESolver,'ExtendThinWireModelByWireNubs','False');

IESolver = invoke(mws,'IESolver');
invoke(IESolver,'SetFMMFFCalcStopLevel','0');
invoke(IESolver,'SetFMMFFCalcNumInterpPoints','6');
invoke(IESolver,'UseFMMFarfieldCalc','True');
invoke(IESolver,'SetCFIEAlpha','0.500000');
invoke(IESolver,'LowFrequencyStabilization','False');
invoke(IESolver,'LowFrequencyStabilizationML','True');
invoke(IESolver,'Multilayer','False');
invoke(IESolver,'SetiMoMACC_I','0.0001');

invoke(IESolver,'SetiMoMACC_M','0.0001');
invoke(IESolver,'DeembedExternalPorts','True');
invoke(IESolver,'SetOpenBC_XY','True');
invoke(IESolver,'OldRCSSweepDefintion','False');
invoke(IESolver,'SetAccuracySetting','Custom');
invoke(IESolver,'CalculateSParaforFieldsources','True');

invoke(IESolver,'ModeTrackingCMA','True');
invoke(IESolver,'NumberOfModesCMA','3');
invoke(IESolver,'StartFrequencyCMA','-1.0');
invoke(IESolver,'SetAccuracySettingCMA','Default');
invoke(IESolver,'FrequencySamplesCMA','0');
invoke(IESolver,'SetMemSettingCMA','Auto');

FDSolver = invoke(mws,'FDSolver');
invoke(FDSolver,'Start');

end

repo: CST-MATLAB-API


text:
function fx = p47_fun ( n, x )

%*****************************************************************************80
%
%% P47_FUN evaluates the integrand for problem 47.
%
%  Discussion:
%
%    The function is singular at the left endpoint.
%
%  Interval:
%
%    0 <= x <= 1
%
%  Integrand:
%
%    sqrt ( x ) * ln ( x )
%
%  Exact Integral:
%
%    -4/9 = -0.4444...
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Piessens, Elise de Doncker-Kapenga,
%    Christian Ueberhuber, David Kahaner,
%    QUADPACK: A Subroutine Package for Automatic Integration,
%    Springer, 1983, page 101.
%
%  Parameters:
%
%    Input, integer N, the number of evaluation points.
%
%    Input, real X(N), the evaluation points.
%
%    Output, real FX(N), the integrand values.
%
  fx(1:n) = sqrt ( x(1:n) ) .* log ( x(1:n) );

  i = find ( x == 0.0 );
  fx(i) = 0.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function Z=SoftThreshold_GS( A, thres )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Soft-thresholding for group lasso
%
% Reference:
% Yuan, Ming, and Yi Lin.
% "Model selection and estimation in regression with grouped variables."
% Journal of the Royal Statistical Society: Series B
% (Statistical Methodology) 68.1 (2006): 49-67.
%
% Provider:
% Hongteng Xu @ Georgia Tech
% June 12, 2017
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Z = zeros(size(A));
for u = 1:size(A, 3)
    for v = 1:size(A, 1)
        tmp = 1 - thres/norm(A(v,:,u));
        if tmp>0
            Z(v,:,u) = tmp*A(v,:,u);
        end
    end
end

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
function [N, eta, sigma_g] = scd_noise_fit_histo(data,varargin)
% [N, eta, sigma_g] = scd_noise_fit_histo(data)
% [N, eta, sigma_g] = scd_noise_fit_histo(__,Name,Value)

data=data(:);
p=inputParser;
addRequired(p,'data',@isnumeric)
addOptional(p,'nbins',max(15,sqrt(length(data))));
[~, xmax] = range_outlier(data);
addOptional(p,'maxval',xmax);
addOptional(p,'maxy',0);
addOptional(p,'color',[0 0 1],@isnumeric);
addOptional(p,'fig',1,@isnumeric);
addOptional(p,'plotfit',0,@isnumeric);
addOptional(p,'distrib','Non-central Chi',@(x) any(validatestring(x,{'Non-central Chi','Rician'})));

if logical(exist('OCTAVE_VERSION', 'builtin'))
    in.fig = varargin{2};
    in.distrib = varargin{4};
    in.nbins = max(15,sqrt(length(data)));
    [~, in.maxval] = range_outlier(data);
    in.maxy = 0;
    in.color = [0 0 1];
    in.plotfit = 0;
else
    parse(p,data,varargin{:});
    in=p.Results;
end

%% PREPARE DATA
data(data==0)=[];
data(data>in.maxval)=[];

if isinteger(data)
    xout=double(unique(data(:)));
    n=double(histc(data(:),xout));
else
    [n,xout]=hist(data,in.nbins);
end
xout(n==0)=[]; n(n==0)=[];  n(xout==0)=0;
n=n/trapz(xout,n);

xout = xout(:); n = n(:);
maxnoise=cumtrapz(xout,n);
maxnoise=xout(find(maxnoise>0.999,1,'first'));
data(data>maxnoise)=[];

if isinteger(data)
    xout=double(unique(data(:)));
    n=double(histc(data(:),xout));
else
    [n,xout]=hist(data,in.nbins);
end

xout(n==0)=[]; n(n==0)=[];  n(xout==0)=0;
n=n/trapz(xout,n);
n(xout<0)=[];
xout(xout<0)=[];

xout = xout(:); n = n(:);
xi=0:max(xout)/in.nbins:max(xout);
yi = interp1(xout,n,xi);
yi(isnan(yi))=0;
yi=yi/trapz(xi,yi);
xi=xout; yi=n;

%% Fit: 'non-central chi'.
%eta = sum(xi.*yi.*mean(diff(xi)));
%    [N      eta     sigma]
var0=[1           eps            std(double(data))];
lb = [1           0                  0.1        ];
ub= [10    double(max(data)) std(double(data)*2)];
switch in.distrib
    case 'Non-central Chi'
        fx = [0 0 0];
    case 'Rician'
        fx = [1 1 0];
end

options = optimset('Algorithm','trust-region-reflective','TolFun',1e-8,'Display','final','MaxIter',10,'Display','off');
disp('     N        eta      sigma_g')
[varfit] = lsqnonlin(@(x1) noncentralchi_error(addfixparameters(var0,x1,fx),xi,yi,in.plotfit), var0(~fx), lb(~fx), ub(~fx),options);
varfit = addfixparameters(var0,varfit,fx);
disp(varfit)
N=varfit(1); eta=varfit(2); sigma_g=varfit(3);
fval = noncentralchi(xi, N,sigma_g,eta);

if in.fig
    figure(73)
    set(73,'Name','Noise histogram','NumberTitle','off')
    plot(xi,yi,'+','Color',in.color,'MarkerSize',10); hold on, plot(xi,fval,'Color',[1 0 0],'Linewidth',2);
    xlim([0 in.maxval]);
    if ~in.maxy, maxy=2*max(yi); else maxy=in.maxy; end
    ylim([0 maxy])
    title(['N = ' num2str(N) ', eta = ' num2str(eta) ', sigma = ' num2str(sigma_g)])
end

end

function Error = noncentralchi_error(var,x,y,plotfit)
N=var(1); eta=var(2); sigma_g=var(3);
f= noncentralchi(x, N,sigma_g,eta);
f(isnan(f))=0;
f=f/trapz(x,f);
if rand<plotfit
    disp(num2str(var))
    figure(4)
    hold off, plot(x,y); hold on, plot(x,f,'r');
    xlim([0 max(x)]);
    ylim([0 2*max(y)])
    pause(0.1)
end
Error = abs(f-y);
end

function p= noncentralchi(x, N,sigma_g,eta)
p = abs(x).^N./(sigma_g^2*eta^(N-1)).*exp(-(x.^2+eta^2)/(2*sigma_g^2)).*besseli(N-1,x*eta/sigma_g^2);
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function sphere_voronoi_test ( )

%*****************************************************************************80
%
%% SPHERE_VORONOI_TEST tests the SPHERE_VORONOI library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 April 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_VORONOI_TEST\n' );
  fprintf ( 1, '  MATLAB version:\n' );
  fprintf ( 1, '  Test the SPHERE_VORONOI library.\n' );

  sphere_voronoi_test01 ( );
  sphere_voronoi_test02 ( );
  sphere_voronoi_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_VORONOI_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
% Rank aware thresholding needs much more measurements than MUSIC to work
% properly.

close all;
clear all;
clc;
rng('default');
png_export = true;
pdf_export = false;
% Create the directory for storing images
[status_code,message,message_id] = mkdir('bin');

mf = spx.graphics.Figures();

% Signal space
N = 256;
% Sparsity level
K = 6;
% Number of measurements
M = 32;
% Number of signals
S = 4;
% Construct the signal generator.
gen  = spx.data.synthetic.SparseSignalGenerator(N, K, S);
% Generate bi-uniform signals
X = gen.biUniform(1, 2);
% Sensing matrix
Phi = spx.dict.simple.gaussian_dict(M, N);
% Measurement vectors
Y = Phi.apply(X);
%  Rank Aware Thresholding solver instance
solver = spx.pursuit.joint.RankAwareThresholding(Phi, K);
% Solve the sparse recovery problem
result = solver.solve(Y);
% Solution vector
Z = result.Z;
% Comparison
cs = spx.commons.SparseSignalsComparison(X, Z, K);
cs.summarize();

for s=1:S
    mf.new_figure(sprintf('MMV signal: %d', s));
    subplot(411);
    stem(X(:, s), '.');
    title('Sparse vector');
    subplot(412);
    stem(Z(:, s), '.');
    title('Recovered sparse vector');
    subplot(413);
    stem(abs(X(:, s) - Z(:, s)), '.');
    title('Recovery error');
    subplot(414);
    stem(Y(:, s), '.');
    title('Measurement vector');
end

author: indigits

repo: sparse-plex


text:
function mcolon_install
% function mcolon_install
% installation of the package mcolon
%
% Author: Bruno Luong <brunoluong@yahoo.com>
% Date: 03-Jan-2011: fix bug with pathname having space

[pname] = fileparts(mfilename('fullpath'));

if ispc()
    cpcmd = 'COPY';
else
    cpcmd = 'cp';
end

CMD = [cpcmd ' "' ...
       pname filesep() 'mcolon_MatlabBased.m" "' ...
       pname filesep() 'mcolon.m"'];
fprintf('%s\n', CMD);
system(CMD);

CMD = [cpcmd ' "' ...
       pname filesep() 'castarrays_MatlabBased.m" "' ...
       pname filesep() 'castarrays.m"'];
fprintf('%s\n', CMD);
system(CMD);

if strfind(computer(), '64')
    mex -v -O -largeArrayDims mcolonmex.c
    mex -v -O -largeArrayDims getclassidmex.c
else
    mex -v -O mcolonmex.c
    mex -v -O -largeArrayDims getclassidmex.c
end

CMD = [cpcmd ' "' ...
       pname filesep() 'mcolon_MexBased.m" "' ...
       pname filesep() 'mcolon.m"'];
fprintf('%s\n', CMD);
system(CMD);

CMD = [cpcmd ' "' ...
       pname filesep() 'castarrays_MexBased.m" "' ...
       pname filesep() 'castarrays.m"'];
fprintf('%s\n', CMD);
system(CMD);

fprintf('\n');
fprintf([repmat('-',1,100) '\n']);
fprintf(['Last step: Plase add manually "%s" \n' ...
         'in your Matlab path using (menu File -> Set Path)\n'], pname);

end

author: Sable

repo: mcbench-benchmarks


text:
function runprop = pbsNonSharedParallelDecodeFcn( runprop )
%PBSNONSHAREDPARALLELDECODEFCN Prepares a worker to run a MATLAB task.
% This function is referenced by pbsNonSharedParallelSubmitFcn.
% THIS FUNCTION MUST BE ON THE PATH OF THE MATLAB WORKERS;
% typically, this is accomplished by putting this function into
% $MATLABROOT/toolbox/local on the workers, or changing
% to the directory where this function is, before starting those
% MATLAB workers.

% Copyright 2006 The MathWorks, Inc.

% Read environment variables into local variables. The names of
% the environment variables were determined by the submit function
storageConstructor = getenv('MDCE_STORAGE_CONSTRUCTOR');
storageLocation = getenv('MDCE_STORAGE_LOCATION');
jobLocation = getenv('MDCE_JOB_LOCATION');

% For a parallel job, the task location is derived from labindex
taskLocation = [jobLocation filesep 'Task' num2str( labindex ) ];

% Set runprop properties from the local variables:
set(runprop, ...
        'StorageConstructor', storageConstructor, ...
        'StorageLocation', storageLocation, ...
        'JobLocation', jobLocation, ...
        'TaskLocation', taskLocation);


text:
function [ x, seed ] = hypersphere_uniform_sample ( m, n, r, c, seed )

%*****************************************************************************80
%
%% HYPERSPHERE_SURFACE_UNIFORM samples hypersphere surface.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Russell Cheng,
%    Random Variate Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley, 1998, pages 168.
%
%    George Marsaglia,
%    Choosing a point from the surface of a sphere,
%    Annals of Mathematical Statistics,
%    Volume 43, Number 2, April 1972, pages 645-646.
%
%    Reuven Rubinstein,
%    Monte Carlo Optimization, Simulation, and Sensitivity
%    of Queueing Networks,
%    Wiley, 1986, page 234.
%
%  Parameters:
%
%    Input, integer M, the dimension of the space.
%
%    Input, integer N, the number of points.
%
%    Input, real R, the radius of the sphere.
%
%    Input, real C(M,1), the center of the sphere.
%
%    Input/output, integer SEED, a seed for the random number generator.
%
%    Output, real X(M,N), the points.
%
  [ x, seed ] = r8mat_normal_01 ( m, n, seed );

  v(1,1:n) = sqrt ( sum ( x.^2, 1 ) );
  vv = repmat ( v, m, 1 );
  x = x ./ vv;
%
%  Scale by the sphere radius.
%
  x = r * x;
%
%  Shift to the sphere center.
%
  x = x + repmat ( c, 1, n );

  return
end


text:
function I=gausslaguerr(f2,a,b,n)

%I=gausslaguerr(f,a,b,n)
%Aproximates integral using Gauss-Laguerre method

%Laguerre polynomial
p=polelague(n);
%Roots
x=roots(p(n+1,:));

G=feval(f2,x);		%Function evaluation on the nodes

%Coeficients
for i=1:n
   C(i)=((factorial(n).^2).*x(i))./(polyval(p(n+1,:),(x(i))).^2);
end

I=dot(C,G);


text:
function [x] = getjet(i)
file    = matfile('jet_data/jetLES.mat');
x       = squeeze(file.p(i,:,:));
end

author: SpectralPOD


text:
%
% Default configuration for Redundant Robotic Manipulator nDOF
%
% by Olzhas Adiyatov
% 6/10/2013

conf = struct;

conf.delta_goal_point = 3;                                          % Radius of goal point
conf.num_link = 6;                                                  % Number of links in the manipulator model
conf.max_ang = [pi; 3*pi/4*ones(conf.num_link-1,1)];                % Maximum joint angle for each joint
conf.min_ang = [0 ;-3*pi/4*ones(conf.num_link-1,1)];                % Minimum joint angle for each joint
conf.delta_ang_max = 10 * pi / 180;                                 % Maximum joint angle change at each iteration
conf.delta_ang_neighbor = 30 * pi / 180;                            % Maximum joint angle change at each neighbor search
conf.disp_interval = 100;                                           % Interval for progress monitoring
conf.init_conf = [pi/2; zeros(conf.num_link-1,1)];                  % Initial configuration of the manipulator model
conf.step_div = 4;                                                  % Number of transition states between nodes, used for collision detection

conf.len_link    = 14*ones(conf.num_link, 1);                       % Array storing the length of the links
conf.height_link = 12 *ones(conf.num_link, 1);                      % Array storing the height of the links
conf.width_link  = 4 *ones(conf.num_link, 1);                       % Array storing the width of the links

author: olzhas


text:
 function warn(varargin)
%function warn(varargin)
%
% print concise warning message

[name line] = caller_name(1);
printf(['Warn: %s %d: ' varargin{1}], name, line, varargin{2:end})

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function y = tdist_unnorm(x, p, type)
%TDIST_UNNORM   Unnormalized t-distribution robust function.
%   TDIST_UNNORM(X, P, TYPE) evaluates the unnormalize t-distribution
%   robust function with parameters P at point(s) X.  P(1) corresponds to
%   the exponent alpha; P(2) corresponds to the scaling sigma.
%   TYPE selects the evaluation type:
%    - 0: function value
%    - 1: first derivative
%    - 2: second derivative
%
%   This is a private member function of the class 'robust_function'.
%
%   Author:  Stefan Roth, Department of Computer Science, TU Darmstadt
%   Contact: sroth@cs.tu-darmstadt.de
%   $Date:  $
%   $Revision: $

% Copyright 2004-2007, Brown University, Providence, RI. USA
% Copyright 2007-2010 TU Darmstadt, Darmstadt, Germany.
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

  r = p(1);
  s = p(2);

  switch (type)
   case 0
    y = r * log(1 + 0.5 * (x / s).^2);
   case 1
    y = r * x ./ (s^2 * (1 + 0.5 * (x / s).^2));
   case 2
    y = r ./ (s^2 * (1 + 0.5 * (x / s).^2));
  end

author: kristinbranson

repo: JAABA


text:
classdef NodalFieldPrinter < FieldRepresenterPrinter

    properties (Access = private)
        nodes
    end

    methods (Access = public)

        function obj = NodalFieldPrinter(cParams)
            obj.init(cParams);
            obj.nodes(:,1) = 1:size(obj.fieldValues,1);
        end

        function printResultsLineHeader(obj)
            iD = obj.fileID;
            fN = obj.fieldName;
            sC = obj.simulationStr;
            it = obj.iter;
            fT = obj.fieldType;
            rL = obj.fieldPosition;
            fprintf(iD,'\nResult "%s" "%s" %.0f %s %s\n',fN,sC,it,fT,rL);
        end

        function printFieldLines(obj)
            iD = obj.fileID;
            pformat = ['%6.0f ',repmat('%12.5d ',1,obj.nComp),' \n'];
            printV  = [obj.nodes, obj.fieldValues]';
            fprintf(iD,pformat,printV);
        end

    end

     methods (Access = private)

        function init(obj,d)
            fieldsNames = fieldnames(d);
            for ifield = 1:length(fieldsNames)
                ifieldName = fieldsNames{ifield};
                fieldValue = d.(ifieldName);
                obj.(fieldsNames{ifield}) = fieldValue;
            end
        end

     end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [vert,conn,tria,tnum] = refine2(varargin)
%REFINE2 (Frontal)-Delaunay-refinement for two-dimensional,
%polygonal geometries.
%   [VERT,EDGE,TRIA,TNUM] = REFINE2(NODE,EDGE) returns a co-
%   nstrained Delaunay triangulation of the polygonal region
%   {NODE,EDGE}. NODE is an N-by-2 array of polygonal verti-
%   ces and EDGE is an E-by-2 array of edge indexing. Each
%   row in EDGE represents an edge of the polygon, such that
%   NODE(EDGE(JJ,1),:) and NODE(EDGE(JJ,2),:) are the coord-
%   inates of the endpoints of the JJ-TH edge. If the argum-
%   ent EDGE is omitted it assumed that the vertices in NODE
%   are connected in ascending order.
%
%   [...] = REFINE2(NODE,EDGE,PART) computes a triangulation
%   for a multiply-connected geometry. PART is a cell-array
%   of polygonal "parts", where each element PART{KK} is an
%   array of edge indices defining a given polygonal region.
%   EDGE(PART{KK}, :) is the set of edges in the KK-TH part.
%
%   VERT is a V-by-2 array of XY coordinates in the triangu-
%   lation, EDGE is an array of constrained edges, TRIA is a
%   T-by-3 array of triangles, and TNUM is a T-by-1 array of
%   part indices. Each row of TRIA and EDGE define an eleme-
%   nt. VERT(TRIA(II,1),:), VERT(TRIA(II,2),:) and VERT(TRIA
%   (II,3),:) are the coordinates of the II-TH triangle. The
%   edges in EDGE are defined in a similar manner. NUM is an
%   array of part indexing, such that TNUM(II) is the index
%   of the part in which the II-TH triangle resides.
%
%   [...] = REFINE2(..., OPTS) passes an additional options
%   structure OPTS, containing various user-defined paramet-
%   ers, including:
%
% - OPTS.KIND = {'DELFRONT'}, 'DELAUNAY' -- the type of ref-
%   inement employed. The 'DELFRONT' algorithm is typically
%   slower, but produces higher quality output.
%
% - OPTS.RHO2 = {1.025} -- the maximum allowable radius-edge
%   ratio. Refinement proceeds until all interior triangles
%   satisfy the radius-edge threshold. Smaller radius-edge
%   ratios lead to improved triangle shape, with RHO2=1 req-
%   uiring that all angles exceed 30 degrees. Setting RHO2<1
%   may lead to non-convergence.
%
% - OPTS.REF1 = {'REFINE'}, 'PRESERVE' -- refinement 'flag'
%   for 1-dimensional faces (i.e. edges). The 'PRESERVE' op-
%   tion results in minimal refinement, attempting to retain
%   the initial edges without further subdivision. Edges are
%   split only to satisfy basic geomertical conformance.
%
% - OPTS.REF2 = {'REFINE'}, 'PRESERVE' -- refinement 'flag'
%   for 2-dimensional faces (i.e. trias). The 'PRESERVE' op-
%   tion results in minimal refinement, attempting to retain
%   the initial trias without further subdivision. Trias are
%   split only to satisfy basic geomertical conformance.
%
% - OPTS.SIZ1 = {1.333} -- the normalised rel.-length th-
%   reshold for edge-elements. Each exterior edge is refined
%   until LL/HH<SIZ1, where LL is the edge-length, HH is the
%   edge-centred mesh-size value.
%
% - OPTS.SIZ2 = {1.300} -- the normalised rel.-length th-
%   reshold for tria-elements. Each interior tria is refined
%   until RE/HH<SIZ2, where RE is an effective tria length,
%   based on the circumradius, HH is the tria-centred mesh-
%   size value.
%
% - OPTS.DISP = { +10 } -- refinement verbosity. Set to INF
%   for quiet execution.
%
%   [...] = REFINE2(..., HFUN,ARGS) also passes an optional
%   mesh-size function argument. Setting HFUN = HMAX, where
%   HMAX is a scalar value, imposes a constant size constra-
%   int over the full domain. HFUN can also be defined as a
%   general function handle [HH] = HFUN(PP), where PP is an
%   N-by-2 array of XY coordinates and HH is the associated
%   vector of mesh-size values. User-defined HFUN must be
%   fully vectorised. Additional arguments {A1,A2,...AN} for
%   HFUN can be passed as trailing parameters to REFINE2. In
%   such cases, HFUN must adopt a signature [HH] = HFUN(PP,
%   A1,A2,...,AN). HFUN must return positive values.
%
%   See also SMOOTH2, TRIDIV2, TRICOST, TRIDEMO

%   This routine implements a "multi-refinement" variant of
%   Delaunay-refinement type mesh-generation. Both standard
%   Delaunay-refinement and Frontal-Delaunay type algorithms
%   are available. The Frontal-Delaunay approach is a simpl-
%   ified version of the JIGSAW algorithm, described in:
%
% * D. Engwirda, (2014): "Locally-optimal Delaunay-refineme-
%   nt and optimisation-based mesh generation", Ph.D. Thesis
%   School of Mathematics and Statistics, Univ. of Sydney.
%   http://hdl.handle.net/2123/13148
%
% * D. Engwirda & D. Ivers, (2016): "Off-centre Steiner poi-
%   nts for Delaunay-refinement on curved surfaces", Comput-
%   er-Aided Design, (72), 157--171.
%   http://dx.doi.org/10.1016/j.cad.2015.10.007

%   This work is an extension of the "off-centre" type tech-
%   niques introduced in:
%
% * H. Erten & A. Ungor, (2009): "Quality triangulation with
%   locally optimal Steiner points", SIAM Journal on Scient-
%   ific Comp. 31(3), 2103--2130.
%   http://doi.org/10.1137/080716748
%
% * S. Rebay, (1993): "Efficient Unstructured Mesh Generati-
%   on by Means of Delaunay Triangulation and Bowyer-Watson
%   Algorithm, J. Comp. Physics 106(1), 125--138.
%   http://dx.doi.org/10.1006/jcph.1993.1097

%   Generally speaking, the Delaunay-refinement method impl-
%   emented here is a variantion of the "classical" algorit-
%   hm introduced in:
%
% * J. Ruppert, (1995): "A Delaunay refinement algorithm for
%   quality 2-dimensional mesh generation." Journal of Algo-
%   rithms 18(3), 548--585.
%   http://dx.doi.org/10.1006/jagm.1995.1021
%
%   See also: S. Cheng, T. Dey & J. Shewchuk, (2012): "Dela-
%   unay mesh generation", CRC Press, for comprehensive cov-
%   erage of Delaunay-based meshing techniques.

%   A much more advanced, and fully three-dimensional imple-
%   mentation is available in the JIGSAW library. For addit-
%   ional information, see:
%   https://github.com/dengwirda/jigsaw-matlab

%-----------------------------------------------------------
%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 13/02/2020
%-----------------------------------------------------------

    node = []; PSLG = []; part = {}; opts = [] ;
    hfun = []; harg = {};

%---------------------------------------------- extract args
    if (nargin>=+1), node = varargin{1}; end
    if (nargin>=+2), PSLG = varargin{2}; end
    if (nargin>=+3), part = varargin{3}; end
    if (nargin>=+4), opts = varargin{4}; end
    if (nargin>=+5), hfun = varargin{5}; end
    if (nargin>=+6), harg = varargin(6:end); end

   [opts] = makeopt(opts) ;

%---------------------------------------------- default EDGE
    nnod = size(node,1) ;

    if (isempty(PSLG))
        PSLG = [(1:nnod-1)',(2:nnod)'; nnod,1] ;
    end

%---------------------------------------------- default PART
    ncon = size(PSLG,1) ;

    if (isempty(part)), part{1} = (1:ncon)'; end

%---------------------------------------------- basic checks
    if (~isnumeric(node) || ~isnumeric(PSLG) || ...
        ~iscell   (part) || ~isstruct (opts) )
        error('refine2:incorrectInputClass' , ...
            'Incorrect input class.') ;
    end

%---------------------------------------------- basic checks
    if (ndims(node) ~= +2 || ndims(PSLG) ~= +2)
        error('refine2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end
    if (size(node,2) < +2 || size(PSLG,2) < +2)
        error('refine2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

%---------------------------------------------- basic checks
    if (min([PSLG(:)])<+1 || max([PSLG(:)])>nnod)
        error('refine2:invalidInputs', ...
            'Invalid EDGE input array.') ;
    end

    pmin = cellfun(@min,part);
    pmax = cellfun(@max,part);

    if (min([pmin(:)])<+1 || max([pmax(:)])>ncon)
        error('refine2:invalidInputs', ...
            'Invalid PART input array.') ;
    end

%-------------------------------- prune any non-unique topo.
   [ivec,ivec,jvec] = ...
        unique(sort(PSLG,+2),'rows') ;

    PSLG = PSLG(ivec,:) ;

    for ppos = +1:length(part)

        if ( ~isnumeric(part{ppos}) )
            error (  ...
            'refine2:incorrectInputClass', ...
                'Incorrect input class. ') ;
        end

        part{ppos} = ...
            unique(jvec(part{ppos})) ;

    end

%-------------------------------- check part "manifold-ness"
    for ppos = +1:length(part)

        eloc = PSLG(part{ppos},:) ;
        nadj = ...
            accumarray(eloc(:),1) ;

        if (any(mod(nadj,2) ~= 0) )
        error('refine2:nonmanifoldInputs', ...
            'Non-manifold PART detected.') ;
        end

    end

%---------------------------------------------- output title
    if (~isinf(opts.disp))
        fprintf(1,'\n') ;
        fprintf(1,' Refine triangulation...\n') ;
        fprintf(1,'\n') ;
        fprintf(1,[...
' -------------------------------------------------------\n', ...
'      |ITER.|          |CDT1(X)|          |CDT2(X)|     \n', ...
' -------------------------------------------------------\n', ...
             ] ) ;
    end

%-------------------------------- PASS 0: inflate box bounds
    vert = node; tria = []; tnum = []; iter = 0 ;
    conn = PSLG;

    vmin = min(vert,[],1);      % inflate bbox for stability
    vmax = max(vert,[],1);

    vdel = vmax - 1.*vmin;
    vmin = vmin - .5*vdel;
    vmax = vmax + .5*vdel;

    vbox = [
        vmin(1), vmin(2)
        vmax(1), vmin(2)
        vmax(1), vmax(2)
        vmin(1), vmax(2)
           ] ;
    vert = [vert ; vbox] ;

%-------------------------------- PASS 0: shield sharp feat.
   [vert,conn,tria,tnum,iter] = ...
        cdtbal0(vert,conn,tria,tnum, ...
            node,PSLG,part,opts,hfun,harg,iter);

%-------------------------------- PASS 1: refine 1-simplexes
   [vert,conn,tria,tnum,iter] = ...
        cdtref1(vert,conn,tria,tnum, ...
            node,PSLG,part,opts,hfun,harg,iter);

%-------------------------------- PASS 2: refine 2-simplexes
   [vert,conn,tria,tnum,iter] = ...
        cdtref2(vert,conn,tria,tnum, ...
            node,PSLG,part,opts,hfun,harg,iter);

    if (~isinf(opts.disp)), fprintf(1,'\n'); end

%-------------------------------- trim extra adjacency info.
    tria = tria( :,1:3) ;

%-------------------------------- trim vert. - deflate bbox.
    keep = false(size(vert,1),1);
    keep(tria(:)) = true;
    keep(conn(:)) = true;

    redo = zeros(size(vert,1),1);
    redo(keep) = ...
        (+1:length(find(keep)))';

    conn = ...
        reshape(redo(conn),[],2);
    tria = ...
        reshape(redo(tria),[],3);

    vert = vert(keep,:) ;

end

function [vert,conn,tria,tnum,iter] = ...
            cdtbal0(vert,conn,tria,tnum, ...
                node,PSLG,part,opts,hfun,harg,iter)
%CDTBAL0 constrained Delaunay-refinement for "sharp" 0-dim.
%features at PSLG vertices.
%   [...] = CDTBAL0(...) refines the set of 1-simplex eleme-
%   nts incident to "sharp" features in the PSLG. Specifica-
%   lly, edges that subtend "small" angles are split about a
%   set of new "collar" vertices, equi-distributed about the
%   centre of "sharp" features. Collar size is computed as a
%   min. of the incident edge-len. and local mesh-size cons-
%   traints.

    if (iter <= opts.iter)

    %------------------------------------- build current CDT
       [vert,conn, ...
        tria,tnum] = deltri2(vert,conn, ...
                             node,PSLG, ...
                             part, ...
                             opts.dtri) ;

    %------------------------------------- build current adj
       [edge,tria] = tricon2(tria,conn) ;

       [feat,ftri] = isfeat2(vert, ...
                             edge,tria) ;

        apex = false(size(vert,1), 1) ;
        apex(tria(ftri)) =  true ;

    %------------------------------------- eval. length-fun.
        if (~isempty(hfun))
            if (isnumeric(hfun))
            vlen = hfun * ...
              ones(size(vert,1),1) ;
            else
            vlen = feval( ...
                hfun,vert,harg{:}) ;
            vlen = vlen(:) ;
            end
        else
            vlen = +inf * ...
              ones(size(vert,1),1) ;
        end

    %------------------------------------- form edge vectors
        evec = vert(conn(:,2),:) ...
             - vert(conn(:,1),:) ;
        elen = sqrt(sum(evec.^2,2));
        evec = evec./[elen,elen] ;

    %------------------------------------- min. adj. lengths
        for epos = +1 : size(conn,1)

            ivrt = conn(epos,1) ;
            jvrt = conn(epos,2) ;

            vlen(ivrt) = min( ...
            vlen(ivrt), .67*elen(epos)) ;
            vlen(jvrt) = min( ...
            vlen(jvrt), .67*elen(epos)) ;

        end

    %------------------------------------- mark feature edge
        iref = apex(conn(:,1)) ...      %- refine at vert. 1
            & ~apex(conn(:,2)) ;
        jref = apex(conn(:,2)) ...      %- refine at vert. 2
            & ~apex(conn(:,1)) ;
        dref = apex(conn(:,1)) ...      %- refine at both!
            &  apex(conn(:,2)) ;

        keep =~apex(conn(:,1)) ...      %- refine at neither
            & ~apex(conn(:,2)) ;

    %------------------------------------- protecting collar
        ilen = vlen(conn(iref,1)) ;
        inew = vert(conn(iref,1),:) ...
        + [ilen,ilen].*evec(iref,:) ;

        jlen = vlen(conn(jref,2)) ;
        jnew = vert(conn(jref,2),:) ...
        - [jlen,jlen].*evec(jref,:) ;

        Ilen = vlen(conn(dref,1)) ;
        Inew = vert(conn(dref,1),:) ...
        + [Ilen,Ilen].*evec(dref,:) ;

        Jlen = vlen(conn(dref,2)) ;
        Jnew = vert(conn(dref,2),:) ...
        - [Jlen,Jlen].*evec(dref,:) ;

        vnew = [inew; jnew; Inew; Jnew] ;

    %------------------------------------- add new vert/edge
        iset = (1:size(inew,1))' ...
                + size(vert,1) ;

        jset = (1:size(jnew,1))' ...
                + size(inew,1) + ...
                + size(vert,1) ;

        Iset = (1:size(Inew,1))' ...
                + size(inew,1) + ...
                + size(jnew,1) + ...
                + size(vert,1) ;

        Jset = (1:size(Jnew,1))' ...
                + size(inew,1) + ...
                + size(jnew,1) + ...
                + size(Inew,1) + ...
                + size(vert,1) ;

        vert = [vert ; vnew] ;

        cnew = [conn(iref,1), iset ;
                conn(iref,2), iset ;
                conn(jref,2), jset ;
                conn(jref,1), jset ;
                conn(dref,1), Iset ;
                conn(dref,2), Jset ;
                Iset, Jset] ;
        conn = [conn(keep,:); cnew ] ;

    end

end

function [vert,conn,tria,tnum,iter] = ...
            cdtref1(vert,conn,tria,tnum, ...
                node,PSLG,part,opts,hfun,harg,iter)
%CDTREF1 constrained Delaunay-refinement for 1-simplex elem-
%nts embedded in R^2.
%   [...] = CDTREF1(...) refines the set of 1-simplex eleme-
%   nts embedded in the triangulation until all constraints
%   are satisfied. Specifically, edges are refined until all
%   local mesh-spacing and encroachment conditions are met.
%   Refinement proceeds according to either a Delaunay-refi-
%   nement or Frontal-Delaunay type approach, depending on
%   user-settings. In either case, new steiner vertices are
%   introduced to split "bad" edges - those that violate the
%   set of prescribed constraints. In the "-DR" type process
%   edges are split about their circumballs (midpoints). In
%   the "-FD" approach, new vertices are positioned such th-
%   at mesh-spacing constraints are satisfied in a "locally-
%   optimal" fashion.

    tcpu.full = +0. ;
    tcpu.ball = +0. ;
    tcpu.hfun = +0. ;
    tcpu.encr = +0. ;
    tcpu.offc = +0. ;

    vidx = (1:size(vert,1))';     %- "new" vert list to test

    tnow =  tic ;

    ntol = +1.55;

    while (strcmpi(opts.ref1,'refine'))

        iter = iter + 1 ;

        if (iter>=opts.iter),break; end

    %------------------------------------- calc. circumballs
        ttic = tic ;

        bal1 = cdtbal1(vert,conn) ;

        tcpu.ball = ...
            tcpu.ball + toc(ttic) ;

    %------------------------------------- eval. length-fun.
        ttic = tic ;

        if (~isempty(hfun))
            if (isnumeric(hfun))
            fun0 = hfun * ...
              ones(size(vert,1),1);
            fun1 = hfun ;
            else
            fun0(vidx) = ...
                feval(hfun, ...
            vert(vidx,:), harg{:});
            fun0 = fun0(:) ;
            fun1 = fun0(conn(:,1))...
                 + fun0(conn(:,2));
            fun1 = fun1 / +2. ;
            end
        else
            fun0 = +inf * ...
              ones(size(vert,1),1);
            fun1 = +inf ;
        end

        siz1 = ...
         +4. * bal1(:,3)./(fun1.*fun1) ;

        tcpu.hfun = ...
            tcpu.hfun + toc(ttic) ;

    %------------------------------------- test encroachment
        ttic = tic ;

        bal1(:,3) = ...
            (1.-eps^.75) * bal1(:,3) ;

       [vp,vi] = ...
           findball(bal1,vert(:,1:2));

    %------------------------------------- near=>[vert,edge]
        next = +0;
        ebad = false(size(conn,1),1) ;
        near = zeros(size(conn,1),1) ;
        for ii = +1 : size(vp,1)
            for ip = vp(ii,1):vp(ii,2)
                jj = vi(ip);
                if (ii ~= conn(jj,1) ...
                &&  ii ~= conn(jj,2) )
                next = next + 1;
                near(next,1) = ii;
                near(next,2) = jj;
                end
            end
        end

        near = near(1:next-0,:);

        if (~isempty(near))
    %-- mark edge "encroached" if there is a vert within its
    %-- dia.-ball that is not joined to either of its vert's
    %-- via an existing edge...
            ivrt = conn(near(:,2),1);
            jvrt = conn(near(:,2),2);

            pair = [near(:,1), ivrt];
            ivec = setset2(pair,conn) ;

            pair = [near(:,1), jvrt];
            jvec = setset2(pair,conn) ;

            okay = ~ivec & ~jvec ;

            ebad(near(okay,2))=true ;

        end

        tcpu.encr = ...
            tcpu.encr + toc(ttic);

    %------------------------------------- refinement queues
        ref1 = false(size(conn,1),1);
        ref1(ebad)           = true ;   %- edge encroachment
        ref1(siz1>opts.siz1* ...        %- bad equiv. length
                  opts.siz1) = true ;

        num1 = find(ref1)  ;

    %------------------------------------- dump-out progess!
        if (mod(iter,opts.disp)==0)
            numc = size(conn,1) ;
            numt = size(tria,1) ;
            fprintf(+1, ...
            '%11i %18i %18i\n', ...
            [iter,numc,numt]) ;
        end

    %------------------------------------- nothing to refine
        if (isempty(num1)), break; end

    %------------------------------------- refine "bad" tria
        switch (lower(opts.kind))
        case 'delaunay'
    %------------------------------------- do circ-ball pt's
        new1 = bal1(ref1, 1:2) ;

        vidx = (1:size(new1,1))' ...
                + size(vert,1) ;

        cnew = [conn( ref1,1), vidx
                conn( ref1,2), vidx];
        conn = [conn(~ref1,:); cnew];

    %------------------------------------- update vertex set
        vert = [vert; new1(:,1:2)];

        case 'delfront'
    %-- symmetric off-centre scheme:- refine edges from both
    %-- ends simultaneously, placing new vertices to satisfy
    %-- the worst of mesh-spacing and local voronoi constra-
    %-- ints.

        ttic = tic ;

        evec = vert(conn(ref1,2),:) ...
             - vert(conn(ref1,1),:) ;
        elen = sqrt(sum(evec.^2,2)) ;
        evec = evec ./ [elen, elen] ;

    %------------------------------------- "voro"-type dist.
        vlen = sqrt(bal1(ref1,3));

    %------------------------------------- "size"-type dist.
        ihfn = fun0(conn(ref1,1));
        jhfn = fun0(conn(ref1,2));

    %------------------------------------- bind "safe" dist.
        ilen = min(vlen,ihfn) ;
        jlen = min(vlen,jhfn) ;

    %------------------------------------- locate offcentres
        inew = vert(conn(ref1,1),:) ...
             + [ilen,ilen].*evec ;
        jnew = vert(conn(ref1,2),:) ...
             - [jlen,jlen].*evec ;

    %------------------------------------- iter. "size"-type
        for ioff = +1 : +3
    %------------------------------------- eval. length-fun.
        if (~isempty(hfun))
            if (isnumeric(hfun))
            iprj = hfun * ...
              ones(size(inew,1),1);
            jprj = hfun * ...
              ones(size(jnew,1),1);
            else
            iprj = feval( ...
                hfun,inew,harg{:});
            jprj = feval( ...
                hfun,jnew,harg{:});
            iprj = iprj(:);
            jprj = jprj(:);
            end
        else
            iprj = +inf * ...
              ones(size(inew,1),1);
            jprj = +inf * ...
              ones(size(jnew,1),1);
        end

        iprj = 0.5*ihfn + 0.5*iprj;
        jprj = 0.5*jhfn + 0.5*jprj;

    %------------------------------------- bind "safe" dist.
        ilen = min(vlen,iprj) ;
        jlen = min(vlen,jprj) ;

    %------------------------------------- locate offcentres
        inew = vert(conn(ref1,1),:) ...
             + [ilen,ilen].*evec ;
        jnew = vert(conn(ref1,2),:) ...
             - [jlen,jlen].*evec ;

        end

    %------------------------------------- merge i,j if near
        near = ...
            ilen+jlen>=vlen*ntol ;

        znew = inew(near,:) * .5 ...
             + jnew(near,:) * .5 ;

        inew = inew(~near,1:2) ;
        jnew = jnew(~near,1:2) ;

    %------------------------------------- split constraints
        zset = (1:size(znew,1))' ...
                + size(vert,1) ;

        iset = (1:size(inew,1))' ...
                + size(znew,1) + ...
                + size(vert,1) ;

        jset = (1:size(jnew,1))' ...
                + size(znew,1) + ...
                + size(inew,1) + ...
                + size(vert,1) ;

        set1 = num1( near);
        set2 = num1(~near);

        cnew = [conn( set1,1), zset
                conn( set1,2), zset
                conn( set2,1), iset
                conn( set2,2), jset
                iset, jset ] ;
        conn = [conn(~ref1,:); cnew];

    %------------------------------------- update vertex set
        vert = [vert; znew(:,1:2)];
        vert = [vert; inew(:,1:2)];
        vert = [vert; jnew(:,1:2)];

        vidx = [zset; iset; jset] ;

        tcpu.offc = ...
            tcpu.offc + toc(ttic) ;

        end % switch(lower(opts.kind))

    end

    tcpu.full = ...
        tcpu.full + toc(tnow) ;

    if (~isinf(opts.disp) )
    %------------------------------------- print final stats
        numc = size(conn,1) ;
        numt = size(tria,1) ;
        fprintf(+1, ...
        '%11i %18i %18i\n', ...
        [iter,numc,numt]) ;
    end

    if (opts.dbug)
    %------------------------------------- print debug timer
        fprintf(1,'\n') ;
        fprintf(1,' 1-simplex REF. timer...\n');
        fprintf(1,'\n') ;
        fprintf(1, ...
        ' FULL: %f \n', tcpu.full);
        fprintf(1, ...
        ' BALL: %f \n', tcpu.ball);
        fprintf(1, ...
        ' HFUN: %f \n', tcpu.hfun);
        fprintf(1, ...
        ' ENCR: %f \n', tcpu.encr);
        fprintf(1, ...
        ' OFFC: %f \n', tcpu.offc);
        fprintf(1,'\n') ;
    end

end

function [vert,conn,tria,tnum,iter] = ...
            cdtref2(vert,conn,tria,tnum, ...
                node,PSLG,part,opts,hfun,harg,iter)
%CDTREF2 constrained Delaunay-refinement for 2-simplex elem-
%nts embedded in R^2.
%   [...] = CDTREF2(...) refines the set of 2-simplex eleme-
%   nts embedded in the triangulation until all constraints
%   are satisfied. Specifically, triangles are refined until
%   all local mesh-spacing and element-shape conditions are
%   met. Refinement proceeds according to either a Delaunay-
%   refinement or Frontal-Delaunay type approach, depending
%   on user-settings. In either case, new steiner points are
%   introduced to split "bad" triangles - those that violate
%   the set of prescribed constraints. In the "-DR" type pr-
%   ocess triangles are split about their circumballs. In
%   the "-FD" approach, new vertices are positioned such th-
%   at mesh-spacing and element-shape constraints are satis-
%   fied in a "locally-optimal" fashion.

    tcpu.full = +0. ;
    tcpu.dtri = +0. ;
    tcpu.tcon = +0. ;
    tcpu.ball = +0. ;
    tcpu.hfun = +0. ;
    tcpu.offc = +0. ;
    tcpu.filt = +0. ;

    vidx = (1:size(vert,1))';     %- "new" vert list to test

    tnow =  tic ;

    near = +.775;

    while (strcmpi(opts.ref2,'refine'))

        iter = iter + 1 ;

    %------------------------------------- build current CDT
        ttic = tic ;

        nold = size(vert,1) ;

       [vert,conn, ...
        tria,tnum]= deltri2(vert,conn, ...
                            node,PSLG, ...
                            part, ....
                            opts.dtri) ;

        nnew = size(vert,1) ;

        vidx = ...
       [vidx; (nold:nnew)'] ;

        tcpu.dtri = ...
            tcpu.dtri + toc(ttic) ;

    %------------------------------------- build current adj
        ttic = tic ;

       [edge,tria]= tricon2(tria,conn) ;

        tcpu.tcon = ...
            tcpu.tcon + toc(ttic) ;

        if (iter>=opts.iter),break; end

    %------------------------------------- calc. circumballs
        ttic = tic ;

        bal1 = cdtbal1(vert,conn) ;
        bal2 = cdtbal2(vert, ...
                       edge,tria) ;
        len2 = minlen2(vert,tria) ;

        rho2 = bal2(:,+3) ./ len2 ;

    %------------------------------------- refinement scores
        scr2 = rho2 .* bal2(:,+3) ;

        tcpu.ball = ...
            tcpu.ball + toc(ttic) ;

    %------------------------------------- eval. length-fun.
        ttic = tic ;

        if (~isempty(hfun))
            if (isnumeric(hfun))
            fun0 = hfun * ...
              ones(size(vert,1),1);
            fun2 = hfun ;
            else
            fun0(vidx) = ...
                feval(hfun, ...
            vert(vidx,:), harg{:});
            fun0 = fun0(:) ;
            fun2 = fun0(tria(:,1))...
                 + fun0(tria(:,2))...
                 + fun0(tria(:,3));
            fun2 = fun2 / +3. ;
            end
        else
            fun0 = +inf * ...
              ones(size(vert,1),1);
            fun2 = +inf ;
        end

        siz2 = ...
         +3. * bal2(:,3)./(fun2.*fun2) ;

        tcpu.hfun = ...
            tcpu.hfun + toc(ttic) ;

    %------------------------------------- refinement queues
        ref1 = false(size(conn,1),1);
        ref2 = false(size(tria,1),1);

        stri = isfeat2(vert,edge,tria) ;

        ref2(rho2>opts.rho2* ...        %- bad rad-edge len.
                  opts.rho2) = true ;
        ref2(stri) = false ;
        ref2(siz2>opts.siz2* ...        %- bad equiv. length
                  opts.siz2) = true ;

        num2 = find(ref2);

    %------------------------------------- dump-out progess!
        if (mod(iter,opts.disp)==0)
            numc = size(conn,1) ;
            numt = size(tria,1) ;
            fprintf(+1, ...
            '%11i %18i %18i\n', ...
            [iter,numc,numt]) ;
        end

    %------------------------------------- nothing to refine
        if (isempty(num2)), break; end

       [scr2,idx2] = sort( ...
            scr2(num2),'descend');
        num2 = num2(idx2);

    %------------------------------------- refine "bad" tria
        switch (lower(opts.kind))
        case 'delaunay'
    %------------------------------------- do circ-ball pt's
        new2 = zeros(length(num2),3);
        new2(:,1:2) = bal2(num2,1:2);

        rmin = ...                      %- min. insert radii
            len2(num2)*(1.-eps^.75)^2 ;

        new2(:,  3) = max( ...
            bal2(num2,3)*near^2,rmin) ;

        case 'delfront'
    %-- off-centre scheme -- refine triangles by positioning
    %-- new vertices along a local segment of the voronoi
    %-- diagram, bounded by assoc. circmballs. New points
    %-- are placed to satisfy the worst of local mesh-length
    %-- and element-shape constraints.

        ttic = tic ;

    %------------------------------------- find frontal edge
       [lmin,emin] = ...
            minlen2(vert,tria(num2,:)) ;

        ftri = false(length(num2),1) ;
        epos = zeros(length(num2),1) ;
        tadj = zeros(length(num2),1) ;

        for ii = +1 : length(epos)
            epos(ii) = tria( ...
                num2(ii),emin(ii)+3) ;
        end

    %------------------------------------- find frontal tria
        for enum = +1 : +3

            eidx = tria(num2,enum+3) ;

            ftri = ...
            ftri | edge(eidx,5) > +0 ;

            ione = ...
                num2 ~= edge(eidx,3) ;
            itwo = ~ione ;

            tadj(ione) = ...
                edge(eidx(ione),3);
            tadj(itwo) = ...
                edge(eidx(itwo),4);

            okay = tadj > +0 ;
            tidx = tadj(okay);

            ftri(okay) = ...
            ftri(okay) | ~ref2(tidx) ;

        end

        if (~any(ftri))                 %- can this happen!?
        ftri = true(length(num2),+1) ;
        end

    %------------------------------------- locate offcentres
        emid = vert(edge(epos,+1),:) ...
             + vert(edge(epos,+2),:) ;
        emid = emid * +0.50 ;

        elen = sqrt(lmin(:));

    %------------------------------------- "voro"-type dist.
        vvec = bal2(num2,1:2)-emid ;
        vlen = sqrt(sum(vvec.^2,2));
        vvec = vvec ./ [vlen,vlen] ;

        hmid = fun0(edge(epos,+1),:) ...
             + fun0(edge(epos,+2),:) ;
        hmid = hmid * +0.50 ;

    %------------------------------------- "ball"-type dist.
        rtri = elen * opts.off2 ;
        rfac = elen * +0.50 ;
        dsqr = rtri.^2 - rfac.^2;
        doff = rtri + ...
            sqrt(max(+0.,dsqr)) ;

    %------------------------------------- "size"-type dist.
        dsiz = +sqrt(3.)/2. * hmid ;

    %------------------------------------- bind "safe" dist.
       [dist,ioff] = ...
          min([dsiz,doff,vlen],[],2) ;

    %------------------------------------- locate offcentres
        off2 = ...
        emid + [dist,dist] .* vvec ;

    %------------------------------------- iter. "size"-type
        for isub = +1 : +3
    %------------------------------------- eval. length-fun.
        if (~isempty(hfun))
            if (isnumeric(hfun))
            hprj = hfun * ...
              ones(size(off2,1),1) ;
            else
            hprj = feval( ...
                hfun,off2,harg{:}) ;
            hprj = hprj(:) ;
            end
        else
            hprj = +inf * ...
              ones(size(off2,1),1) ;
        end

    %------------------------------------- "size"-type dist.
        hprj = .33*hmid + .67*hprj ;

        dsiz = +sqrt(3.)/2. * hprj ;

        dsiz(dsiz<elen*.50) = +inf ;    %- edge-ball limiter
        dsiz(dsiz>vlen*.95) = +inf ;    %- circ-ball limiter

    %------------------------------------- bind "safe" dist.
       [dist,ioff] = ...
          min([dsiz,doff,vlen],[],2) ;

    %------------------------------------- locate offcentres
        off2 = ...
        emid + [dist,dist] .* vvec ;

        end

        orad = ...
        sqrt((elen*.5).^2 + dist.^2) ;

    %------------------------------------- do offcentre pt's
        new2 = ...
        zeros(length(find(ftri)),+3) ;
        new2(:,1:2) = off2(ftri,1:2) ;

        rmin = ...                      %- min. insert radii
            lmin(ftri)*(1.-eps^.75)^2 ;

        new2(:,  3) = max( ...
            (orad(ftri)*near).^2,rmin);

        tcpu.offc = ...
            tcpu.offc + toc (ttic) ;

        end % switch(lower(opts.kind))

    %------------------------------------- inter.-ball dist.
        ttic = tic ;

    %------------------------------------- proximity filters
       [vp,vi] = ...
          findball(new2,new2(:,1:2)) ;

        keep = true (size(new2,1),1) ;
        for ii = size(vp,1):-1:+1
            for ip = vp(ii,1) ...
                   : vp(ii,2)
                jj = vi(ip);
                if (keep(jj) && ...
                    keep(ii) && ...
                    jj < ii )

                keep(ii) = false ;
                break;

                end
            end
        end

        new2 = new2(keep,:);

    %------------------------------------- test encroachment
        bal1(:,3) = ...
            (1.-eps^.75) * bal1(:,3);

       [vp,vi] = ...
          findball(bal1,new2(:,1:2));

        keep = true (size(new2,1),1);
        for ii = +1:+1:size(vp,1)
            for ip = vp(ii,1) ...
                   : vp(ii,2)
                jj = vi(ip);
                ref1(jj) =  true ;
                keep(ii) = false ;
            end
        end

    %------------------------------------- leave sharp edges
        ebnd = false(size(edge,1),1);
        ebnd(tria(stri,4:6)) = true ;

        enot = ...
        setset2(conn,edge(ebnd,1:2));

        ref1(enot) = false ;

    %------------------------------------- preserve boundary
        if (strcmp(lower(opts.ref1),...
            'preserve'))
        ref1(:)    = false ;
        end

    %------------------------------------- refinement points
        new2 = new2(keep,:);
        new1 = bal1(ref1,:);

        tcpu.filt = ...
            tcpu.filt + toc(ttic) ;

    %------------------------------------- split constraints
        idx1 = ...
       (1:size(new1))'+size(vert,1) ;

        idx2 = ...
       (1:size(new2))'+size(new1,1) ...
                      +size(vert,1) ;

        cnew = [conn( ref1,1), idx1
                conn( ref1,2), idx1];
        conn = [conn(~ref1,:); cnew];

        vidx = [idx1; idx2];

    %------------------------------------- update vertex set
        nold = size(vert,1);
        vert = [vert; new1(:,1:2)];
        vert = [vert; new2(:,1:2)];
        nnew = size(vert,1);

        if (nnew == nold), break; end   %- we *must* be done

    end

    tcpu.full = ...
        tcpu.full + toc(tnow) ;

    if (~isinf(opts.disp) )
    %------------------------------------- print final stats
        numc = size(conn,1) ;
        numt = size(tria,1) ;
        fprintf(+1, ...
        '%11i %18i %18i\n', ...
        [iter,numc,numt]) ;
    end

    if (opts.dbug)
    %------------------------------------- print debug timer
        fprintf(1,'\n') ;
        fprintf(1,' 2-simplex REF. timer...\n');
        fprintf(1,'\n') ;
        fprintf(1, ...
        ' FULL: %f \n', tcpu.full);
        fprintf(1, ...
        ' DTRI: %f \n', tcpu.dtri);
        fprintf(1, ...
        ' TCON: %f \n', tcpu.tcon);
        fprintf(1, ...
        ' BALL: %f \n', tcpu.ball);
        fprintf(1, ...
        ' HFUN: %f \n', tcpu.hfun);
        fprintf(1, ...
        ' OFFC: %f \n', tcpu.offc);
        fprintf(1, ...
        ' FILT: %f \n', tcpu.filt);
        fprintf(1,'\n') ;
    end

end

function [opts] = makeopt(opts)
%MAKEOPT setup the options structure for REFINE2.

    if (~isfield(opts,'dtri'))
        opts.dtri = 'constrained';
    else
    if (~strcmpi(opts.dtri, 'conforming') && ...
        ~strcmpi(opts.dtri,'constrained') )
        error( ...
    'refine2:invalidOption','Invalid constraint DTRI.');
    end
    end

    if (~isfield(opts,'kind'))
        opts.kind = 'delfront';
    else
    if (~strcmpi(opts.kind, 'delfront') && ...
        ~strcmpi(opts.kind, 'delaunay') )
        error( ...
    'refine2:invalidOption','Invalid refinement KIND.');
    end
    end

    if (~isfield(opts,'ref1'))
        opts.ref1 = 'refine';
    else
    if (~strcmpi(opts.ref1,   'refine') && ...
        ~strcmpi(opts.ref1, 'preserve') )
        error( ...
    'refine2:invalidOption','Invalid refinement REF1.');
    end
    end

    if (~isfield(opts,'ref2'))
        opts.ref2 = 'refine';
    else
    if (~strcmpi(opts.ref2,   'refine') && ...
        ~strcmpi(opts.ref2, 'preserve') )
        error( ...
    'refine2:invalidOption','Invalid refinement REF2.');
    end
    end

    if (~isfield(opts,'iter'))
        opts.iter = +inf;
    else
    if (~isnumeric(opts.iter))
        error('refine2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.iter)~= +1)
        error('refine2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.iter <= +0)
        error('refine2:invalidOptionValues', ...
            'Invalid OPT.ITER selection.') ;
    end
    end

    if (~isfield(opts,'disp'))
        opts.disp = +10 ;
    else
    if (~isnumeric(opts.disp))
        error('refine2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.disp)~= +1)
        error('refine2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.disp <= +0)
        error('refine2:invalidOptionValues', ...
            'Invalid OPT.DISP selection.') ;
    end
    end

    if (~isfield(opts,'rho2'))
        opts.rho2 = 1.025;
    else
    if (~isnumeric(opts.rho2))
        error('refine2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.rho2)~= +1)
        error('refine2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.rho2 < +1.)
        error('refine2:invalidOptionValues', ...
            'Invalid OPT.RHO2 selection.') ;
    end
    end

    if (~isfield(opts,'off2'))
        opts.off2 = 0.933;
    else
    if (~isnumeric(opts.off2))
        error('refine2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.off2)~= +1)
        error('refine2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.off2 < +.7)
        error('refine2:invalidOptionValues', ...
            'Invalid OPT.OFF2 selection.') ;
    end
    end

    if (~isfield(opts,'siz1'))
        opts.siz1 = 1.333;
    else
    if (~isnumeric(opts.siz1))
        error('refine2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.siz1)~= +1)
        error('refine2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.siz1 <= 0.)
        error('refine2:invalidOptionValues', ...
            'Invalid OPT.SIZ1 selection.') ;
    end
    end

    if (~isfield(opts,'siz2'))
        opts.siz2 = 1.300;
    else
    if (~isnumeric(opts.siz2))
        error('refine2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.siz2)~= +1)
        error('refine2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.siz2 <= 0.)
        error('refine2:invalidOptionValues', ...
            'Invalid OPT.SIZ2 selection.') ;
    end
    end

    if (~isfield(opts,'dbug'))
        opts.dbug = false;
    else
    if (~islogical(opts.dbug))
        error('refine2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.dbug)~= +1)
        error('refine2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    end

end

author: dengwirda


text:
function [widths, maxDepth, nodePositions] = treeGetWidths(tree)

% TREEGETWIDTHS give width of each level of tree.
% FORMAT
% DESC gives the width of a tree at each level of the hierarchy and
% the maximum depth of a tree as well as the node stored at a give
% depth and breadth.
% ARG tree : the tree for which the dimensions are required.
% RETURN widths : stores the width at each depth level.
% RETURN maxDepth : the maximum depth of the tree.
% RETURN nodePositions : stores the nodeIndex present at depth i
% and breadth j.
%
% COPYRIGHT : Andrew J. Moore, 2006
%
% SEEALSO : treeFindParents, treeFindChildren

% NDLUTIL

maxDepth = 0;
%widths stores the width at each depth level. Since the max depth isn't yet
%known, allocate enough memory to widths for the worst case scenario where
%each node has only 1 child and the depth is equal to the number of nodes.
widths = zeros(length(tree), 1);
%nodePositions
nodePositions = zeros(length(tree), length(tree));
rootInd = treeFindRoots(tree);
indAlreadySeen = [];
for i = 1:length(rootInd)
  traverseTree(rootInd(i), 1);
end
widths = widths(1:maxDepth, :); %trim off excess rows

  function traverseTree(nodeIndex, depthLevel)
    if ~any(indAlreadySeen == nodeIndex)
      widths(depthLevel) = widths(depthLevel) + 1;
    end
    indAlreadySeen = [indAlreadySeen nodeIndex];

    nodePositions(depthLevel, widths(depthLevel)) = nodeIndex;
    if length(tree(nodeIndex).children) > 0
      for i=1:length(tree(nodeIndex).children)
        traverseTree(tree(nodeIndex).children(i), depthLevel + 1);
      end
    else
      if depthLevel > maxDepth
        maxDepth = depthLevel;
      end
    end
  end

end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [cpulse, verbose] = tapas_physio_findpeaks_template_xcorr(...
    c, pulseCleanedTemplate, cpulseSecondGuess, averageHeartRateInSamples, ...
    verbose, varargin)
% Finds peaks of a time series via pre-determined template via maxima of
% matlab cross correlations (xcorr) via going backward from search starting
% point in time series, and afterwards forward again
%
%   [cpulse, verbose] = tapas_physio_findpeaks_template_correlation(...
%       c, pulseCleanedTemplate, cpulseSecondGuess, ...
%           averageHeartRateInSamples, verbose)
%
% IN
%   varargin    property name/value pairs for additional options
%
%
% OUT
%
% EXAMPLE
%   tapas_physio_findpeaks_template_correlation
%
%   See also

% Author: Steffen Bollmann, cleanup, xcorr: Lars Kasper
% Created: 2014-08-05
% Copyright (C) 2014 TNU, Institute for Biomedical Engineering,
%                         University of Zurich and ETH Zurich.
%
% This file is part of the TAPAS PhysIO Toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Determine starting peak for the search:
%   search for a representative R-peak a range of peaks

nSamples = size(c,1);

debug = verbose.level >= 4;

idxStartPeakSearch = [0 20];

halfTemplateWidthInSamples = floor(numel(pulseCleanedTemplate)/2);

[~,zTransformedTemplate] = tapas_physio_corrcoef12(pulseCleanedTemplate,...
    pulseCleanedTemplate);
isZTransformed = [0 1];

% start and end point of search for representative start cycle
centreSampleStart = round(2*halfTemplateWidthInSamples+1);

if idxStartPeakSearch(1) > 0
    centreSampleStart = centreSampleStart + ...
        cpulseSecondGuess(idxStartPeakSearch(1));
end
centreSampleEnd = cpulseSecondGuess(idxStartPeakSearch(2));

iSignalStart = centreSampleStart - halfTemplateWidthInSamples;
iSignalEnd = centreSampleEnd + halfTemplateWidthInSamples;
signalPart = c(iSignalStart:iSignalEnd);

similarityToTemplate = xcorr(zTransformedTemplate, flipud(signalPart));

% not needed, since only zero-filled
% similarityToTemplate = similarityToTemplate(1:centreSampleEnd);

[C_bestMatch, I_bestMatch] = max(similarityToTemplate);

%% now compute backwards to the beginning:
% go average heartbeat by heartbeat back and look (with
% decreasing weighting for higher distance) for highest
% correlation with template heartbeat

n = I_bestMatch;
bestPosition = n; % to capture case where 1st R-peak is best

peakNumber = 1;

similarityToTemplate = zeros(nSamples,1);

searchStepsTotal    = round(0.5*averageHeartRateInSamples);
searchPositionArray = -searchStepsTotal:searchStepsTotal;
nSamplesSignalPart  = 2*searchStepsTotal+1;
locationWeight      = ones(nSamplesSignalPart,1);

while n > 1+searchStepsTotal+halfTemplateWidthInSamples

    % Nested function, needs c, zTransformedTemplate, n, halfTemplateWidthInSamples,
    % searchStepsTotal, locationWeight
    similarityToTemplate(n+searchPositionArray) = ...
        get_similarity_to_template();

    %find biggest correlation-peak from the last search
    indexSearchStart    = n-searchStepsTotal;
    indexSearchEnd      = n+searchStepsTotal;

    indexSearchRange    = indexSearchStart:indexSearchEnd;
    searchRangeValues   = similarityToTemplate(indexSearchRange);
    [C_bestMatch,I_bestMatch] = max(searchRangeValues);
    bestPosition = indexSearchRange(I_bestMatch);

    cpulse(peakNumber) = bestPosition;
    peakNumber = peakNumber+1;

    n=bestPosition-averageHeartRateInSamples;
end % END: going backwards to beginning of time course

%% Now go forward through the whole time series
n           = bestPosition; % 1st R-peak
peakNumber  = 1;
clear cpulse;

% Now correlate template with PPU signal at the positions
% where we would expect a peak based on the average heartrate and
% search in the neighborhood for the best peak, but weight the peaks
% deviating from the initial starting point by a gaussian
searchStepsTotal = round(0.5*averageHeartRateInSamples);

% for weighted searching of max correlation
locationWeight = tapas_physio_gausswin(nSamplesSignalPart);

n = max(n, searchStepsTotal + halfTemplateWidthInSamples + 1);

% zero-pad c at end to allow for detection of last peak by
% template-matching up to the last sample of c
c = [c; zeros(searchStepsTotal + halfTemplateWidthInSamples + 1, 1)];

while n < nSamples % -searchStepsTotal - halfTemplateWidthInSamples

     similarityToTemplate(n+searchPositionArray) = ...
         get_similarity_to_template();

    %find biggest correlation-peak from the last search
    indexSearchStart    = n - searchStepsTotal;
    indexSearchEnd      = n + searchStepsTotal;

    indexSearchRange    = indexSearchStart:indexSearchEnd;
    searchRangeValues   = similarityToTemplate(indexSearchRange);
    [C_bestMatch,I_bestMatch] = max(searchRangeValues);
    bestPosition        = indexSearchRange(I_bestMatch);

    cpulse(peakNumber) = bestPosition;
    peakNumber = peakNumber+1;

    %only take the last 20 cpulses to compute the current HeartRate
    foundCpulses = size(cpulse,2);

    if  foundCpulses < 3
        currentHeartRateInSamples=averageHeartRateInSamples;
    end

    if (foundCpulses < 21) && (foundCpulses >= 3)
        currentHeartRateInSamples = round(mean(diff(cpulse)));
    end

    if foundCpulses >= 21
        currentCpulses = cpulse (foundCpulses-20:foundCpulses);
        currentHeartRateInSamples = round(mean(diff(currentCpulses)));
    end

    %check currentHeartRate
    checkSmaller    = currentHeartRateInSamples > 0.5*averageHeartRateInSamples;
    checkLarger     = currentHeartRateInSamples < 1.5*averageHeartRateInSamples;

    %jumpToNextPeakSearchArea
    if (checkSmaller && checkLarger)
        n = bestPosition + currentHeartRateInSamples;
    else
        n = bestPosition + averageHeartRateInSamples;
    end
end

%% Nested function,
% computes point-wise similarity (cross-correlation) of time course snippet
% to given peak template (z-transformed)
% nested function to improve performance
%
% needs c, zTransformedTemplate, n, halfTemplateWidthInSamples,
% searchStepsTotal, searchPositionArray, locationWeight
    function similarityToTemplateTmp = get_similarity_to_template()
        % (c, zTransformedTemplate, n, halfTemplateWidthInSamples, searchStepsTotal, locationWeight);
        iSignalStart    = n - halfTemplateWidthInSamples - searchStepsTotal;
        iSignalEnd      = n + halfTemplateWidthInSamples + searchStepsTotal;

        signalPart = c(iSignalStart:iSignalEnd);
        similarityToTemplateTmp = ...
            xcorr(zTransformedTemplate, flipud(signalPart));
        similarityToTemplateTmp = similarityToTemplateTmp(1:nSamplesSignalPart);
        % crop beginning and end
        %similarityToTemplateTmp(1:halfTemplateWidthInSamples) = [];
        %similarityToTemplateTmp(end-halfTemplateWidthInSamples+1:end) = [];

        % reweight correlations with distance from expected heart beat
        amplitudeWeight = abs(c((n+1)-searchPositionArray));
        similarityToTemplateTmp =  locationWeight.*amplitudeWeight .* ...
            similarityToTemplateTmp;

    end

end

author: translationalneuromodeling

repo: tapas


text:
function [ W, W1,S, index_mat,IDX_vote,S_tan,Data_PCA ] = createAffMatrix( feats,params,affinity_type )
%% W1 affinity matrix using Euclidean distances and RBF.

switch affinity_type
    case 'local_PCA'
       % [W1,index_mat, S] = RBF_affnity(feats,params);
        [W1,index_mat, S,IDX_vote,S_tan] = RBF_affnity_knn_fast(feats,params);
        %% W2 affinity using local tangent space with local PCA
        [Data_PCA,W2,theta] = locPCA_affnity_fast(feats, IDX_vote, params);
        %index_mat=IDX_vote;
        %% Connecting all the points
        Wc = W1.*W2;
        W = Wc.*Wc';

        %W_large=W_large.*W2;
        %W_large=W_large.*W_large;

    case 'Hillinger_affinity'
        %params.affinity_type='cosine';
      %  [W1, index_mat, S] = affinity.RBF_affnity(feats,params);
        [W1,index_mat, S,IDX_vote,S_tan] = RBF_affnity_knn(feats,params);
        params.affinity_type= 'Hillinger_affinity';

        [W2, W_l]= Hillinger_affinity( feats, IDX_vote,params);
        Wc = W1.*W2;
        W = Wc.*Wc';
    case 'sqeuclidean'
       % [W1,index_mat, S] = affinity.RBF_affnity(feats,params);
      %  [W1,index_mat, S,IDX_vote,S_tan] = affinity.RBF_affnity_knn(feats,params);
        [W1,index_mat, S,IDX_vote,S_tan] = RBF_affnity_knn_fast(feats,params);
        W=W1;
        case 'outliers'
        [W1,index_mat, S,IDX_vote,S_tan] = RBF_affnity_knn_fast(feats,params);
         W=W1;

    case 'cosine'
       % [W1,index_mat, S] = affinity.RBF_affnity(feats,params);
      % [W1,index_mat, S,IDX_vote,S_tan] = affinity.RBF_affnity_knn(feats,params);
       [W1,index_mat, S,IDX_vote,S_tan] = RBF_affnity_knn_fast(feats,params);

        W = W1;
    case 'Tensor_Voting'
        %[W1,index_mat, S] = affinity.RBF_affnity(feats,params);
        [W1,index_mat, S,IDX_vote,S_tan] = RBF_affnity_knn(feats,params);
        %% W2 affinity using the Tensor Voting Graph
        [Normal_Space,W2,theta] = TVG_affinity_knn(feats, IDX_vote, params,S_tan);
        %% Connecting all the points
        Wc = W1.*W2;
        %W=(Wc + Wc')./2;
        W = Wc.*Wc';
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [parameters, LL, ht, VCVrobust, VCV, scores, diagnostics] = agarch(epsilon, p, q, model_type, error_type, startingvals, options)
% AGARCH(P,Q) and NAGARCH(P,Q) with different error distributions:
% Normal, Students-T, Generalized Error Distribution, Skewed T
%
% USAGE:
%   [PARAMETERS] = agarch(EPSILON,P,Q)
%   [PARAMETERS,LL,HT,VCVROBUST,VCV,SCORES,DIAGNOSTICS]
%                                = agarch(EPSILON,P,Q,MODEL_TYPE,ERROR_TYPE,STARTINGVALS,OPTIONS)
%
% INPUTS:
%   EPSILON      - A column of mean zero data
%   P            - Positive, scalar integer representing the number of symmetric innovations
%   Q            - Non-negative, scalar integer representing the number of lags of conditional
%                    variance (0 for ARCH-type model)
%   MODEL_TYPE   - [OPTIONAL] The type of variance process, either
%                    'AGARCH'  - Asymmetric GARCH, Engle (1990) [DEFAULT]
%                    'NAGARCH' - Nonlinear Asymmetric GARCH, Engle & Ng (1993)
%   ERROR_TYPE   - [OPTIONAL] The error distribution used, valid types are:
%                    'NORMAL'    - Gaussian Innovations [DEFAULT]
%                    'STUDENTST' - T distributed errors
%                    'GED'       - Generalized Error Distribution
%                    'SKEWT'     - Skewed T distribution
%   STARTINGVALS - [OPTIONAL] A (2+p+q), plus 1 for STUDENTST OR GED (nu),  plus 2 for SKEWT
%                    (nu,lambda), vector of starting values.
%                  [omega alpha(1) ... alpha(p) gamma beta(1) ... beta(q) [nu lambda]]'.
%   OPTIONS      - [OPTIONAL] A user provided options structure. Default options are below.
%
% OUTPUTS:
%   PARAMETERS   - A 2+p+q column vector of parameters with
%                  [omega alpha(1) ... alpha(p) gamma beta(1) ... beta(q) [nu lambda]]'.
%   LL           - The log likelihood at the optimum
%   HT           - The estimated conditional variances
%   VCVROBUST    - Robust parameter covariance matrix
%   VCV          - Non-robust standard errors (inverse Hessian)
%   SCORES       - Matrix of scores (# of params by t)
%   DIAGNOSTICS  - Structure of optimization output information.  Useful to check for convergence
%                     problems
% COMMENTS:
%   The following (generally wrong) constraints are used:
%    (1) omega > 0
%    (2) alpha(i) >= 0 for i = 1,2,...,p
%    (3) beta(i)  >= 0 for i = 1,2,...,q
%    (4) -q(.01,EPSILON)<gamma<q(.99,EPSILON) for AGARCH
%    (5) sum(alpha(i) + beta(k)) < 1 for i = 1,2,...p and k=1,2,...,q for
%    AGARCH and sum(alpha(i)*(1+gamma^2) + beta(k)) < 1 for NAGARCH
%    (6) nu>2 of Students T and nu>1 for GED
%    (7) -.99<lambda<.99 for Skewed T
%
%    The conditional variance, h(t), of a AGARCH(P,Q) process is given by:
%
%     h(t)  = omega
%             + alpha(1)*(r_{t-1}-gamma)^2 + ... + alpha(p)*(r_{t-p}-gamma)^2
%             + beta(1)*h(t-1) +...+ beta(q)*h(t-q)
%
%    The conditional variance, h(t), of a NAGARCH(P,Q) process is given by:
%
%     h(t)  = omega
%             + alpha(1)*(r_{t-1}-gamma*sqrt(h(t-1)))^2 + ... + alpha(p)*(r_{t-p}-gamma*sqrt(h(t-p)))^2
%             + beta(1)*h(t-1) +...+ beta(q)*h(t-q)
%
%   Default Options
%     options  =  optimset('fminunc');
%     options  =  optimset(options , 'TolFun'      , 1e-005);
%     options  =  optimset(options , 'TolX'        , 1e-005);
%     options  =  optimset(options , 'Display'     , 'iter');
%     options  =  optimset(options , 'Diagnostics' , 'on');
%     options  =  optimset(options , 'LargeScale'  , 'off');
%     options  =  optimset(options , 'MaxFunEvals' , '200*numberOfVariables');
%
%  See also AGARCH_LIKELIHOOD, AGARCH_CORE, AGARCH_PARAMETER_CHECK, AGARCH_TRANSFORM, AGARCH_ITRANSFORM
%
%  You should use the MEX files (or compile if not using Win64 Matlab) as they provide speed ups of
%  approx 10 times relative to the m file

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 7/12/2009

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch nargin
    case 3
        [p,q,model_type,error_type,startingvals,options]=agarch_parameter_check(epsilon, p, q);
    case 4
        [p,q,model_type,error_type,startingvals,options]=agarch_parameter_check(epsilon, p, q, model_type);
    case 5
        [p,q,model_type,error_type,startingvals,options]=agarch_parameter_check(epsilon, p, q, model_type, error_type);
    case 6
        [p,q,model_type,error_type,startingvals,options]=agarch_parameter_check(epsilon, p, q, model_type, error_type, startingvals);
    case 7
        [p,q,model_type,error_type,startingvals,options]=agarch_parameter_check(epsilon, p, q, model_type, error_type, startingvals, options);
    otherwise
        error('Number of inputs must be between 3 and 7');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Initial setup
m  =  max([p q]);

% Augment the epsilon with local back casts to avoid costly memory allocations
back_cast_length = max(floor(length(epsilon)^(1/2)),1);
back_cast_weights = .05*(.9.^(0:back_cast_length ));
back_cast_weights = back_cast_weights/sum(back_cast_weights);
back_cast = back_cast_weights*((epsilon(1:back_cast_length+1)).^2);
if back_cast==0
     back_cast=cov(epsilon);
end
epsilon_augmented=[sqrt(back_cast)*ones(m,1);epsilon];
%Compute the length of the augmented epsilon
T = size(epsilon_augmented,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Starting values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This flag is for the robustness check below, if the user supplies starting
%values, there will be no robustness check
[startingvals,nu,lambda]=agarch_starting_values(startingvals,epsilon,p,q,model_type,error_type);
%Finally, initialize the starting values
startingvals = [startingvals; nu; lambda];
% Compute transform bounds
transform_bounds = quantile(epsilon,[.01 .99]);
%Transform the starting vals
[garch_params_transformed,nu_transformed,lambda_transformed]=agarch_transform(startingvals,p,q,model_type,error_type,transform_bounds);
%Re-append nu, lambda
startingvals_transformed = [garch_params_transformed; nu_transformed; lambda_transformed];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Starting values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Estimate the parameters. Note the 1 in the last argument to indicate it
% is a constrained optimization

%LL0 is used to make sure the log likelihood improves
LL0=agarch_likelihood(startingvals_transformed,epsilon_augmented,p,q,model_type,error_type,transform_bounds ,back_cast,T,1);
%Parameter estimation
[parameters,LL,exitflag,output]=fminunc('agarch_likelihood',startingvals_transformed,options,epsilon_augmented,p,q,model_type,error_type,transform_bounds,back_cast,T,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Estimation Robustness
% This portion of the code is to make sure that the optimization converged
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is the case where the optimization did not converge, but improved on
% the initial log likelihood
if  exitflag<=0 && LL<LL0
    % Try more iterations, only do more iterations if the final likelihood is
    % actually better than the initial

    % Increase the max iterations and max fun evals
    % Also switch to steepest descent
    if ischar(options.MaxFunEvals)
        options.MaxIter=2*100*length(parameters);
    else
        options.MaxIter=2*options.MaxIter;
    end
    if ischar(options.MaxFunEvals)
        options.MaxFunEvals=4*100*length(parameters);
    else
        options.MaxFunEvals=2*options.MaxFunEvals;
    end
    options.HessUpdate='steepdesc';
    % Estimate the parameters.
    [parameters,LL,exitflag,output]=fminunc('agarch_likelihood',parameters,options,epsilon_augmented,p,q,model_type,error_type,transform_bounds,back_cast,T,1);
end

% Transform the parameters from the real line to the restricted space
[parameters,nu,lambda]=agarch_itransform(parameters,p,q,model_type,error_type,transform_bounds);
parameters=[parameters;nu;lambda];
% Compute the log likelihood if needed
if nargout>1
    [LL, likelihoods, ht]=agarch_likelihood(parameters,epsilon_augmented,p,q,model_type,error_type,transform_bounds,back_cast,T);
    LL=-LL;
end

%Compute standard errors using RobustVCV if needed.
if nargout>3
    nw=0; %No newey west on scores
    [VCVrobust,A,B,scores,hess]=robustvcv('agarch_likelihood',parameters,nw,epsilon_augmented,p,q,model_type,error_type,transform_bounds,back_cast,T);
    VCV=hess^(-1)/(T-m);
end

%Report diagnostics in case requested
diagnostics.EXITFLAG=exitflag;
diagnostics.ITERATIONS=output.iterations;
diagnostics.FUNCCOUNT=output.funcCount;
diagnostics.MESSAGE=output.message;

author: bashtage

repo: mfe-toolbox


text:
classdef MaF15 < PROBLEM
% <multi/many> <real> <large/none>
% Inverted LSMOP8

%------------------------------- Reference --------------------------------
% R. Cheng, M. Li, Y. Tian, X. Zhang, S. Yang, Y. Jin, and X. Yao, A
% benchmark test suite for evolutionary many-objective optimization,
% Complex & Intelligent Systems, 2017, 3(1): 67-81.
%------------------------------- 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".
%--------------------------------------------------------------------------

    properties(Access = private)
        sublen;	% Number of variables in each subcomponent
        len;    % Cumulative sum of lengths of variable groups
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            % Parameter setting
            if isempty(obj.M); obj.M = 3; end
            if isempty(obj.D); obj.D = 20*obj.M; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = [ones(1,obj.M-1),10.*ones(1,obj.D-obj.M+1)];
            obj.encoding = ones(1,obj.D);
            % Calculate the number of variables in each subcomponent
            nk = 2;
            c  = 3.8*0.1*(1-0.1);
            for i = 1 : obj.M-1
                c = [c,3.8.*c(end).*(1-c(end))];
            end
            obj.sublen = floor(c./sum(c).*(obj.D-obj.M+1)/nk);
            obj.len    = [0,cumsum(obj.sublen*nk)];
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            [N,D] = size(PopDec);
            M     = obj.M;
            nk    = 2;
            PopDec(:,M:D) = (1+repmat(cos((M:D)./D*pi/2),N,1)).*PopDec(:,M:D) - repmat(PopDec(:,1)*10,1,D-M+1);
            G = zeros(N,M);
            for i = 1 : 2 : M
                for j = 1 : nk
                    G(:,i) = G(:,i) + Griewank(PopDec(:,obj.len(i)+M-1+(j-1)*obj.sublen(i)+1:obj.len(i)+M-1+j*obj.sublen(i)));
                end
            end
            for i = 2 : 2 : M
                for j = 1 : nk
                    G(:,i) = G(:,i) + Sphere(PopDec(:,obj.len(i)+M-1+(j-1)*obj.sublen(i)+1:obj.len(i)+M-1+j*obj.sublen(i)));
                end
            end
            G      = G./repmat(obj.sublen,N,1)./nk;
            PopObj = (1+G+[G(:,2:end),zeros(N,1)]).*(1-fliplr(cumprod([ones(N,1),cos(PopDec(:,1:M-1)*pi/2)],2)).*[ones(N,1),sin(PopDec(:,M-1:-1:1)*pi/2)]);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R = UniformPoint(N,obj.M);
            R = 1 - R./repmat(sqrt(sum(R.^2,2)),1,obj.M);
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            if obj.M == 2
                R = obj.GetOptimum(100);
            elseif obj.M == 3
                a = linspace(0,pi/2,10)';
                R = {1-sin(a)*cos(a'),1-sin(a)*sin(a'),1-cos(a)*ones(size(a'))};
            else
                R = [];
            end
        end
    end
end

function f = Griewank(x)
    f = sum(x.^2,2)./4000 - prod(cos(x./repmat(sqrt(1:size(x,2)),size(x,1),1)),2) + 1;
end

function f = Sphere(x)
    f = sum(x.^2,2);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO