text:
classdef RegularizedPerimeterComputer < handle

    properties (Access = public)
        perimeters
        perimetersGradient
        regularizedDensity
        perimeterShapeFunctions
        epsilons
    end

    properties (Access = private)
        inputFile
        backgroundMesh
        designVariable
        scale
        nEpsilon
        epsilon
        printing
        plotting
        capturingImage
        perimeterShapeFunction
        outputFigureName
        perimeterType
        isRobinTermAdded
    end

    methods (Access = public)

        function obj = RegularizedPerimeterComputer(cParams)
            obj.init(cParams);
        end

        function compute(obj)
            obj.createEpsilonValues();
            obj.computePerimeters();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.inputFile        = cParams.inputFile;
            obj.backgroundMesh   = cParams.backgroundMesh;
            obj.scale            = cParams.scale;
            obj.designVariable   = cParams.designVariable;
            obj.outputFigureName = cParams.outputFigureName;
            obj.plotting         = cParams.plotting;
            obj.printing         = cParams.printing;
            obj.capturingImage   = cParams.capturingImage;
            obj.perimeterType    = cParams.perimeterType;
            obj.isRobinTermAdded = cParams.isRobinTermAdded;
        end

        function createEpsilonValues(obj)
            epsmin = obj.backgroundMesh.computeMeanCellSize;
            epsmax = obj.backgroundMesh.computeCharacteristicLength();
            obj.nEpsilon = min(6,ceil(log2(epsmax/epsmin)));
            obj.epsilons = epsmin*(2.^((1:obj.nEpsilon) - 1));
        end

        function computePerimeters(obj)
            obj.perimeters = zeros(obj.nEpsilon,1);
            for iepsilon = 1:obj.nEpsilon
                obj.epsilon = obj.epsilons(iepsilon);
                obj.computeRegularizedPerimeter(iepsilon);
                obj.plotDensity();
                obj.printDensity(iepsilon);
                obj.captureImage(iepsilon);
            end
        end

        function computeRegularizedPerimeter(obj,iepsilon)
            obj.createPerimeterShapeFunction();
            obj.perimeterShapeFunction.computeFunctionAndGradient();
            per    = obj.perimeterShapeFunction.value;
            dPer   = obj.perimeterShapeFunction.gradient;
            rhoReg = obj.perimeterShapeFunction.regularizedDensity;
            obj.perimeters(iepsilon) = per*obj.perimeterShapeFunction.value0;
            obj.perimetersGradient(:,iepsilon) = dPer;
            obj.regularizedDensity(:,iepsilon) = rhoReg;
            obj.perimeterShapeFunctions{iepsilon} = obj.perimeterShapeFunction;
        end

        function createPerimeterShapeFunction(obj)
            s = obj.createPerimeterParams();
            shFunc = ShFunc_Perimeter(s);
            obj.perimeterShapeFunction = shFunc;
        end

        function s = createPerimeterParams(obj)
            sC.inputFile        = obj.inputFile;
            sC.mesh             = obj.backgroundMesh;
          %  sC.designVariable.value      = obj.designVariable;
          %  sC.designVariable.nVariables = 1;
            sC.designVariable   = obj.createDesignVariable();
            sC.epsilon          = obj.epsilon;
            sC.scale            = obj.scale;
            sC.type             = obj.perimeterType;
            sC.isRobinTermAdded = obj.isRobinTermAdded;
            fCreator = PerimeterParamsCreator(sC);
            s = fCreator.perimeterParams;
        end

        function dV = createDesignVariable(obj)
          %  s.scalarProductSettings.femSettings = [];
          %  s.epsilon = [];
          %  ss.type = 'full';
          %  sLs = SettingsLevelSetCreator;
          %  s = sLs.create(ss);
            s.mesh      = obj.backgroundMesh;
            s.inputFile = obj.inputFile;
            s.scale     = obj.scale;
            s.type      = 'LevelSet';
            d = DesignVariableCreatorSettings(s);
            s = d.create();
            dV = DesignVariable.create(s);
            dV.update(obj.designVariable);
        end

        function plotDensity(obj)
            s.mesh      = obj.backgroundMesh;
            s.inputFile = obj.inputFile;
            s.scale     = obj.scale;
            s.density   = obj.perimeterShapeFunction.regularizedDensity;
            plotter = DensityPlotterForPerimeter(s);
            if obj.plotting
                plotter.plot();
            end
        end

        function printDensity(obj,iepsilon)
            s.inputFile = obj.inputFile;
            s.iter      = iepsilon;
            s.mesh      = obj.backgroundMesh;
            s.perimeter   = obj.perimeterShapeFunction;
            printer = DensityPrinterForPerimeter(s);
            printer.print();
        end

        function captureImage(obj,iepsilon)
            if obj.capturingImage
                i = iepsilon;
                f = obj.inputFile;
                outPutNameWithIter = [obj.outputFigureName,'Epsilon',num2str(i)];
                inputFileName = fullfile('Output',f,[f,num2str(i),'.flavia.res']);
                s.fileName = f;
                s.outPutImageName = outPutNameWithIter;
                s.inputFileName = inputFileName;
                imageCapturer = GiDImageCapturer(s);
                imageCapturer.capture();
            end
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function cerrFeatS = calcIBSI1PhantomRadiomics(config)
% Compute features using IBSI reference Lung CT image and configurations.
%--------------------------------------------------------------------------
% INPUTS
% config : 'C' or 'A' (to use settings corresponding to IBSI coonfig C or A).
%--------------------------------------------------------------------------
% AI 06/25/2020

%% Get config file path
switch(config)

    case 'C'
        configPath_avg = fullfile(fileparts(fileparts(getCERRPath)),...
            'Unit_Testing/settings_for_comparisons/IBSIconfigC_avg.json');
        configPath_merge = fullfile(fileparts(fileparts(getCERRPath)),...
            'Unit_Testing/settings_for_comparisons/IBSIconfigC_merge.json');

    case 'A'
        configPath_avg = fullfile(fileparts(fileparts(getCERRPath)),...
            'Unit_Testing/settings_for_comparisons/IBSIconfigA_avg.json');
        configPath_merge = fullfile(fileparts(fileparts(getCERRPath)),...
            'Unit_Testing/settings_for_comparisons/IBSIconfigA_merge.json');

end

%% Load IBSI CT phantom
fpath = fullfile(fileparts(fileparts(getCERRPath)),...
    'Unit_Testing/data_for_cerr_tests/IBSI1_CT_phantom/IBSILungCancerCTImage.mat.bz2');
planC = loadPlanC(fpath,tempdir);
planC = updatePlanFields(planC);
planC = quality_assure_planC(fpath,planC);
indexS = planC{end};

%% Compute features using CERR
%Read config files
param1S = getRadiomicsParamTemplate(configPath_avg);
param2S = getRadiomicsParamTemplate(configPath_merge);

%Get struct name
strName = param1S.structuresC;
strName = strName{1};
strC = {planC{indexS.structures}.structureName};
structNum = getMatchingIndex(strName,strC,'exact');
scanNum = getStructureAssociatedScan(structNum,planC);

%Compute features (avg)
cerrFeatS = calcGlobalRadiomicsFeatures...
    (scanNum, structNum, param1S, planC);

%Compute features (merge)
cerrFeat2S = calcGlobalRadiomicsFeatures...
    (scanNum, structNum, param2S, planC);
mergeFieldsC = intersect(fieldnames(cerrFeatS),fieldnames(cerrFeat2S));

for m = 1:length(mergeFieldsC)
    mergeFeatC = intersect( fieldnames(cerrFeatS.(mergeFieldsC{m})),...
        fieldnames(cerrFeat2S.(mergeFieldsC{m})));
    for n = 1:length(mergeFeatC)
        addFieldC = fieldnames(cerrFeat2S.(mergeFieldsC{m}).(mergeFeatC{n}));
        cerrFeatS.(mergeFieldsC{m}).(mergeFeatC{n}).(addFieldC{1}) ...
            = cerrFeat2S.(mergeFieldsC{m}).(mergeFeatC{n}).(addFieldC{1});
    end
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function drawAxis3d(varargin)
%DRAWAXIS3D Draw a coordinate system and an origin.
%
%   drawAxis3d
%	Adds 3 cylinders to the current axis, corresponding to the directions
%	of the 3 basis vectors Ox, Oy and Oz.
%	Ox vector is red, Oy vector is green, and Oz vector is blue.
%
%   drawAxis3d(L, R)
%   Specifies the length L and the radius of the cylinders representing the
%   different axes.
%
%   Example
%   drawAxis
%
%   figure;
%   drawAxis(20, 1);
%
%   See also
%   drawAxisCube
%
% ------
% Author: David Legland
% e-mail: david.legland@nantes.inra.fr
% Created: 2007-08-14,    using Matlab 7.4.0.287 (R2007a)
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

% geometrical data
origin = [0 0 0];
v1 = [1 0 0];
v2 = [0 1 0];
v3 = [0 0 1];

% default parameters
L = 1;
r = L/10;

% extract parameters from input
if ~isempty(varargin)
	L = varargin{1};
end
if length(varargin)>1
	r = varargin{2};
end

% draw 3 cylinders and a ball
hold on;
drawCylinder([origin origin+v1*L r], 16, 'facecolor', 'r', 'edgecolor', 'none');
drawCylinder([origin origin+v2*L r], 16, 'facecolor', 'g', 'edgecolor', 'none');
drawCylinder([origin origin+v3*L r], 16, 'facecolor', 'b', 'edgecolor', 'none');
drawSphere([origin 2*r], 'faceColor', 'black');

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function display_element ( prefix, element )

%*****************************************************************************80
%
%% DISPLAY_ELEMENT lists the coordinates of nodes making up an element.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 June 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string PREFIX, the common filename prefix for the node and
%    element files.
%
%    Input, integer ELEMENT, the index of the element.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'DISPLAY_ELEMENT\n' );
  fprintf ( 1, '  List coordinates of the nodes making up an element.\n' );
  fprintf ( 1, '  Element and node file prefix is "%s".\n', prefix );

  element_filename = sprintf ( '%s_elements.txt', prefix );
  node_filename = sprintf ( '%s_nodes.txt', prefix );

  element_node = load ( element_filename );
  node_xyz = load ( node_filename );

  [ element_num, element_order ] = size ( element_node );
  [ node_num, dim_num ] = size ( node_xyz );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  LEO ordering:\n' );
  fprintf ( 1, '  Node coordinates for element %d\n', element );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   I     N          X                 Y               Z\n' );
  fprintf ( 1, '\n' );
  for j = 1 : element_order
    n = element_node(element,j);
    fprintf ( 1, '  %2d  %4d  %14.6g  %14.6g  %14.6g\n', j, n, node_xyz(n,1:3) );
  end

  leo_to_gmsh = [...
     1,  2,  3,  4,  5, ...
     7, 10, 11,  8,  6, ...
    12, 18, 19, 13, 20, ...
    14,  9, 15, 16, 17 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Gmsh ordering:\n' );
  fprintf ( 1, '  Node coordinates for element %d\n', element );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   I     N          X                 Y               Z\n' );
  fprintf ( 1, '\n' );
  for j = 1 : element_order
    k = leo_to_gmsh(j);
    n = element_node(element,k);
    fprintf ( 1, '  %2d  %4d  %14.6g  %14.6g  %14.6g\n', j, n, node_xyz(n,1:3) );
  end
  return
end

author: johannesgerer

repo: jburkardt-m


text:
% ijcvMultiSegScript

imdir = '../images/stanford';
outdir = '../data';
ncv = 5;

if 1
    load '../data/stanford_imsegs2.mat';
    ijcv = load('../data/ijcv06/ijcvClassifier.mat');
    %load '../data/bsdsdata.mat';
end

cv_images = [1:numel(imsegs)];
nsegments = [5 10 15 20 25 30 35 40 45 50 60 70 80 90 100];

if ~exist('spfeatures')
    spfeatures = mcmcGetAllSuperpixelData(imdir, imsegs);
    [efeatures, adjlist] = mcmcGetAllEdgeData(spfeatures, imsegs);
    for f = 1:numel(imsegs)
        fprintf(1, '.', int32(f))
        if mod(f, 25)==0, fprintf(1, '\n'); end

        [pvSP{f}, phSP{f}, pE{f}] = mcmcInitialize(spfeatures{f}, efeatures{f}, ...
            adjlist{f}, imsegs(f), ijcv.vclassifierSP, ijcv.hclassifierSP, ijcv.eclassifier, ijcv.ecal, 'none');
        smaps{f} = generateMultipleSegmentations2(pE{f}, adjlist{f}, imsegs(f).nseg, nsegments);

        im = im2double(imread([imdir '/' imsegs(f).imname]));
        imdata = mcmcComputeImageData(im, imsegs(f));

        for k = 1:numel(nsegments)
            labdata{f, k} = mcmcGetSegmentFeatures(imsegs(f), spfeatures{f}, imdata, smaps{f}(:, k), (1:max(smaps{f}(:, k))));
            [mclab{f, k}, mcprc{f, k}, allprc{f, k}, trainw{f, k}] = segmentation2labels(imsegs(f), smaps{f}(:, k));
            unilabel{f, k} = mclab{f, k}.*(mcprc{f, k}>0.95);
            seglabel{f,k} =  1*(mcprc{f, k}>0.95) + (-1)*(mcprc{f, k}<0.95);
        end
    end

    save([outdir '/indoordata.mat'], 'spfeatures', 'efeatures', 'adjlist', ...
        'pvSP', 'phSP', 'pE', 'smaps', 'labdata', ...
        'mclab', 'mcprc', 'allprc', 'seglabel', 'unilabel', ...
        'trainw');
end
%nsegments = [3 5 6 6 7 8 9 10 11 12 13 14 15 17 18 20 23 27 32 45 50 75 100];
%nsegments = [3 4 5 6 7 7 8 8 9 9 10 11 11 12 12 13 13 14 14 15 16 16 17 19 20 21 23 25 28 31 37 51];

[vacc, hacc, vcm, hcm, pg] = testMultipleSegmentationsCV2(imsegs, ...
    labdata, labdata, smaps, ijcv.vclassifier, ijcv.hclassifier, ijcv.sclassifier, ...
    pvSP, phSP, 1);
save([outdir '/indoorResults_trainOutdoor.mat'], 'vacc', 'hacc', 'vcm', 'hcm', 'pg');
disp(num2str([vacc hacc]))
disp(num2str(vcm))
disp(num2str(hcm))

if ~exist('vclassifier')
    for k = 1:ncv
        disp(['Iteration: ' num2str(k)]);
        testind{k} = (floor((k-1)*numel(cv_images)/ncv)+1):(floor(k*numel(cv_images)/ncv));
        trainind{k} = setdiff([1:numel(cv_images)], testind{k});
        sclassifier(k) = mcmcTrainSegmentationClassifier2(...
            [labdata(trainind{k}, :) ; labdata], ...
            [seglabel(trainind{k}, :) ; seglabel], ...
            [trainw(trainind{k}, :) ; trainw]);
        [vclassifier(k), hclassifier(k)] = ...
            mcmcTrainSegmentClassifier2(...
            [labdata(trainind{k}, :) ; labdata], ...
            [unilabel(trainind{k}, :) ; unilabel], ...
            [trainw(trainind{k}, :) ; trainw], 100000);
    end
end
[vacc, hacc, vcm, hcm, pg] = testMultipleSegmentationsCV2(imsegs, ...
    labdata, labdata, smaps, vclassifier, hclassifier, sclassifier, [], [], ncv);
disp(num2str([vacc hacc]))
disp(num2str(vcm))
disp(num2str(hcm))

save([outdir '/indoorResults.mat'], 'vclassifier', 'hclassifier', 'sclassifier', 'vacc', 'hacc', 'vcm', 'hcm', 'pg');

author: Cloud-CV

repo: object-proposals


text:
function ExtractVideoCNNFeature_ucf(layer,scale,gpu_id, tag, ind)
addpath /home/lmwang/code/caffe_data_parallel/caffe/matlab
sizes =[480,640; 340,454; 240,320; 170,227; 120,160];
path_org = '/nfs/lmwang/lmwang/Data/UCF101/ucf101_org/';
folderlist = dir(path_org);
foldername = {folderlist(:).name};
foldername = setdiff(foldername,{'.','..'});

% FLOW
if  strcmp(tag,'flow')
    path2 = ['/media/RAID0/lmwang/data/UCF/ucf101_tvl1_flow_conv5_scale_',num2str(scale),'/'];
    path3 = ['/media/RAID0/lmwang/data/UCF/ucf101_tvl1_flow_conv4_scale_',num2str(scale),'/'];
    path4 = ['/media/RAID0/lmwang/data/UCF/ucf101_tvl1_flow_conv3_scale_',num2str(scale),'/'];
    path5 = ['/media/RAID0/lmwang/data/UCF/ucf101_tvl1_flow_pool2_scale_',num2str(scale),'/'];
    path_flow =  '/media/RAID0/lmwang/data/UCF/ucf101_flow_img_tvl1_gpu/';

    model_def_file = [ 'models/flow_',layer,'_scale_',num2str(scale),'_new.prototxt'];
    model_file = '10_flow_iter_90000_new.caffemodel';
    caffe.reset_all();
    caffe.set_mode_gpu();
    caffe.set_device(gpu_id);
    net = caffe.Net(model_def_file, model_file, 'test');

    for i = ind
        if ~exist([path2,foldername{i}],'dir')
            mkdir([path2,foldername{i}]);
        end
        if ~exist([path3,foldername{i}],'dir')
            mkdir([path3,foldername{i}]);
        end
        if ~exist([path4,foldername{i}],'dir')
            mkdir([path4,foldername{i}]);
        end
        if ~exist([path5,foldername{i}],'dir')
            mkdir([path5,foldername{i}]);
        end
        filelist = dir([path_org,foldername{i},'/*.avi']);
        for j = 1:length(filelist)
            filename = [path_flow,foldername{i},'/',filelist(j).name(1:end-4),'/'];
            if  ~exist([path5,foldername{i},'/',filelist(j).name(1:end-4),'.mat'],'file')
                tic;
                [feature_c5, feature_c4, feature_c3, feature_p2] = TemporalCNNFeature(filename, net, sizes(scale,1), sizes(scale,2));
                toc;
                tic;
                feature = feature_c5;
                save([path2,foldername{i},'/',filelist(j).name(1:end-4),'.mat'],'feature','-v7.3');
                feature = feature_c4;
                save([path3,foldername{i},'/',filelist(j).name(1:end-4),'.mat'],'feature','-v7.3');
                feature = feature_c3;
                save([path4,foldername{i},'/',filelist(j).name(1:end-4),'.mat'],'feature','-v7.3');
                feature = feature_p2;
                save([path5,foldername{i},'/',filelist(j).name(1:end-4),'.mat'],'feature','-v7.3');
                toc;
            end
        end
        i
    end
    caffe.reset_all();
end

author: wanglimin

repo: TDD

save_path: github-repos/MATLAB/wanglimin-TDD


text:
function checkdiff(problem, x, d, force_gradient)
% Checks the consistency of the cost function and directional derivatives.
%
% function checkdiff(problem)
% function checkdiff(problem, x)
% function checkdiff(problem, x, d)
%
% checkdiff performs a numerical test to check that the directional
% derivatives defined in the problem structure agree up to first order with
% the cost function at some point x, along some direction d. The test is
% based on a truncated Taylor series (see online Manopt documentation).
%
% Both x and d are optional and will be sampled at random if omitted.
%
% See also: checkgradient checkhessian

% If force_gradient = true (hidden parameter), then the function will call
% getGradient and infer the directional derivative, rather than call
% getDirectionalDerivative directly. This is used by checkgradient.

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%
%   March 26, 2017 (JB):
%       Detects if the approximated linear model is exact
%       and provides the user with the corresponding feedback.
%
%   April 3, 2015 (NB):
%       Works with the new StoreDB class system.
%
%   Aug. 2, 2018 (NB):
%       Using storedb.remove() to avoid unnecessary cache build-up.
%
%   Sep. 6, 2018 (NB):
%       Now checks whether M.exp() is available; uses retraction otherwise.
%
%   June 18, 2019 (NB):
%       Now issues a warning if the cost function returns complex values.

    if ~exist('force_gradient', 'var')
        force_gradient = false;
    end

    % Verify that the problem description is sufficient.
    if ~canGetCost(problem)
        error('It seems no cost was provided.');
    end
    if ~force_gradient && ~canGetDirectionalDerivative(problem)
        error('It seems no directional derivatives were provided.');
    end
    if force_gradient && ~canGetGradient(problem)
        % Would normally issue a warning, but this function should only be
        % called with force_gradient on by checkgradient, which will
        % already have issued a warning.
    end

    x_isprovided = exist('x', 'var') && ~isempty(x);
    d_isprovided = exist('d', 'var') && ~isempty(d);

    if ~x_isprovided && d_isprovided
        error('If d is provided, x must be too, since d is tangent at x.');
    end

    % If x and / or d are not specified, pick them at random.
    if ~x_isprovided
        x = problem.M.rand();
    end
    if ~d_isprovided
        d = problem.M.randvec(x);
    end

    % Compute the value f0 at f and directional derivative at x along d.
    storedb = StoreDB();
    xkey = storedb.getNewKey();
    f0 = getCost(problem, x, storedb, xkey);

    if ~force_gradient
        df0 = getDirectionalDerivative(problem, x, d, storedb, xkey);
    else
        grad = getGradient(problem, x, storedb, xkey);
        df0 = problem.M.inner(x, grad, d);
    end

    % Pick a stepping function: exponential or retraction?
    if isfield(problem.M, 'exp')
        stepper = problem.M.exp;
    else
        stepper = problem.M.retr;
        % No need to issue a warning: to check the gradient, any retraction
        % (which is first-order by definition) is appropriate.
    end

    % Compute the value of f at points on the geodesic (or approximation
    % of it) originating from x, along direction d, for stepsizes in a
    % large range given by h.
    h = logspace(-8, 0, 51);
    value = zeros(size(h));
    for k = 1 : length(h)
        y = stepper(x, d, h(k));
        ykey = storedb.getNewKey();
        value(k) = getCost(problem, y, storedb, ykey);
        storedb.remove(ykey); % no need to keep it in memory
    end

    % Compute the linear approximation of the cost function using f0 and
    % df0 at the same points.
    model = polyval([df0 f0], h);

    % Compute the approximation error
    err = abs(model - value);

    % And plot it.
    loglog(h, err);
    title(sprintf(['Directional derivative check.\nThe slope of the '...
                   'continuous line should match that of the dashed\n'...
                   '(reference) line over at least a few orders of '...
                   'magnitude for h.']));
    xlabel('h');
    ylabel('Approximation error');

    line('xdata', [1e-8 1e0], 'ydata', [1e-8 1e8], ...
         'color', 'k', 'LineStyle', '--', ...
         'YLimInclude', 'off', 'XLimInclude', 'off');

    if ~all( err < 1e-12 )
        % In a numerically reasonable neighborhood, the error should
        % decrease as the square of the stepsize, i.e., in loglog scale,
        % the error should have a slope of 2.
        isModelExact = false;
        window_len = 10;
        [range, poly] = identify_linear_piece(log10(h), log10(err), window_len);
    else
        % The 1st order model is exact: all errors are (numerically) zero
        % Fit line from all points, use log scale only in h.
        isModelExact = true;
        range = 1:numel(h);
        poly = polyfit(log10(h), err, 1);
        % Set mean error in log scale for plot.
        poly(end) = log10(poly(end));
        % Change title to something more descriptive for this special case.
        title(sprintf(...
              ['Directional derivative check.\n'...
               'It seems the linear model is exact:\n'...
               'Model error is numerically zero for all h.']));
    end
    hold all;
    loglog(h(range), 10.^polyval(poly, log10(h(range))), 'LineWidth', 3);
    hold off;

    if ~isModelExact
        fprintf('The slope should be 2. It appears to be: %g.\n', poly(1));
        fprintf(['If it is far from 2, then directional derivatives ' ...
                 'might be erroneous.\n']);
    else
        fprintf(['The linear model appears to be exact ' ...
                 '(within numerical precision),\n'...
                 'hence the slope computation is irrelevant.\n']);
    end

    if ~(isreal(value) && isreal(f0))
        fprintf(['# The cost function appears to return complex values' ...
              '.\n# Please ensure real outputs.\n']);
    end

end

author: NicolasBoumal

repo: manopt


text:
function scores = JAABADetect(expdir,varargin)
% Run classifiers trained from JAABA on experiments
% scores = JAABADetect(expdir,'classifierfiles',classifierfiles,'configfiles',configfiles)
% scores = JAABADetect(expdir,'classifierparamsfile',classifierparamsfile)
% expdir (String) -- is the experiment directory.
% classifierfiles (Cell array of Strings) -- locations of the classifiers
% for the behaviors that you want to detect.
% configfiles (Cell array of strings) -- locations of the project files
% that were used to train the classifiers.
% Example usage:
% JAABADetect('testExp','classifierfiles',{'ChaseClassifier.mat'},'configfiles',{'ChaseProject.mat'});

[blockSize,classifierfiles,configfiles,classifierparamsfile,configparams,forcecompute,DEBUG,isrelativepath,fnsrelative] = ...
  myparse(varargin,'blockSize',10000,...
  'classifierfiles',{},'configfiles',{},...
  'classifierparamsfile',0,...
  'configparams',[],...
  'forcecompute',false,...
  'debug',false,...
  'isrelativepath',false,...
  'fnsrelative',{'featureparamfilename'});

if ~isdeployed,
  SetUpJAABAPath;
end

if ischar(forcecompute),
  forcecompute = str2double(forcecompute) ~= 0;
end
if ischar(DEBUG),
  DEBUG = str2double(DEBUG);
end

if ischar(classifierparamsfile),
  configparams = ReadClassifierParamsFile(classifierparamsfile,'isrelativepath',isrelativepath,'fnsrelative',fnsrelative);
  classifierfiles = {configparams.classifierfile};
  configfiles = {configparams.configfile};
%{
%   if ~exist(classifierparamsfile,'file'),
%     error('File %s does not exist',classifierparamsfile);
%   end
%   fid = fopen(classifierparamsfile,'r');
%   classifierfiles = {};
%   configfiles = {};
%   while true,
%     l = fgetl(fid);
%     if ~ischar(l),
%       break;
%     end
%     if isempty(l),
%       continue;
%     end
%     if strcmp(l(1),'%'), continue; end
%     ws = regexp(l,',','split');
%     classifierfiles{end+1} = ws{1}; %#ok<AGROW>
%     configfiles{end+1} = ws{2}; %#ok<AGROW>
%   end
%   fclose(fid);
%}
end

nclassifiers = numel(classifierfiles);
if nclassifiers == 0,
  error('No classifiers input');
end
% read parameters
fprintf('Reading classifier parameters...\n');
if isempty(configparams),
  if nclassifiers ~= numel(configfiles),
    error('Number of classifier files and number of config files do not match');
  end
  configparams = cell(1,nclassifiers);
  for i = 1:nclassifiers,
    [~,~,ext] = fileparts(configfiles{i});
    if strcmpi(ext,'.xml');
      configparams{i} = ReadXMLConfigParams(configfiles{i});
    elseif strcmpi(ext,'.jab')
      tmp = load(configfiles{i},'-mat');
      configparams{i} = tmp.x;
    else
      configparams{i} = load(configfiles{i});
    end
  end
else
  if nclassifiers ~= numel(configparams),
    error('Number of classifier files and number of config params do not match');
  end
  if ~iscell(configparams),
    configparams = num2cell(configparams);
  end
end
classifiers = cell(1,nclassifiers);
for i = 1:nclassifiers,
  [~,~,ext] = fileparts(classifierfiles{i});
  if strcmpi(ext,'.jab'),
    tmp = load(classifierfiles{i},'-mat');
    classifiers{i} = tmp.x;
  else
    classifiers{i} = load(classifierfiles{i});
  end
end

% scores file name
scorefilenames = cell(1,nclassifiers);
for i = 1:nclassifiers,
  % obsolete, ReadXMLConfigParams does this
%   if ~isfield(configparams{i}.file,'scorefilename'),
%     if ~isfield(configparams{i},'behaviors'),
%       error('configparams %d does not have field behaviors',i);
%     elseif ~isfield(configparams{i}.behaviors,'names'),
%       error('configparams{%d}.behaviors does not have field names',i);
%     else
%       scorefilenames{i} = ['scores_',configparams{i}.behaviors.names,'.mat'];
%     end
%   else
    scorefilenames{i} = configparams{i}.file.scorefilename;
%  end
end

% find classifiers to run together
[classifiersets,featureconfigperset] = GroupClassifiers(configparams,scorefilenames);

% find classifiers that are done
docompute = true(1,nclassifiers);
if forcecompute == 0,
  for i = 1:nclassifiers,
    scorefile = fullfile(expdir,scorefilenames{i});
    if exist(scorefile,'file'),

      % check timestamp
      try
        sd = load(scorefile,'timestamp');
        if sd.timestamp == classifiers{i}.classifierTS,
          docompute(i) = false;
          behaviorname = configparams{i}.behaviors.names;
          if iscell(behaviorname),
            behaviorname = sprintf('%s_',behaviorname{:});
            behaviorname = behaviorname(1:end-1);
          end
          fprintf('Skipping classifier %s(%s) for experiment %s, scores already computed\n',...
            behaviorname,datestr(classifiers{i}.classifierTS,'yyyymmddTHHMMSS'),...
            expdir); %#ok<PFCEL>
        end
      catch ME,
        fprintf('Could not check timestamps: %s\n',getReport(ME));
        continue;
      end
    end

  end

end

scores = cell(1,nclassifiers);
if nargout >= 1,
  for i = find(~docompute),
    scorefile = fullfile(expdir,scorefilenames{i});
    try
      sd = load(scorefile,'allScores');
      scores{i} = sd.allScores.scores;
    catch ME,
      warning('Could not load in scores from %s: %s',scorefile,getReport(ME));
      docompute(i) = true;
    end
  end
end

if ~any(docompute),
  fprintf('All classifiers have been run on experiment %s. Exiting.\n',expdir);
  return;
end

for seti = 1:numel(classifiersets),

  classifierset = intersect(find(docompute),classifiersets{seti});

  if isempty(classifierset),
    continue;
  end

  nclassifiers_curr = numel(classifierset);

  fprintf('Processing %d classifiers together:\n',nclassifiers_curr);
  fprintf('%s\n',classifierfiles{classifierset});

  % get a list of all window features
  wfidx_per_classifier = cell(1,nclassifiers_curr);
  wfs_per_classifier = cell(1,nclassifiers_curr);
  wfs = {};
  fprintf('Getting list of window features...\n');
  for ii = 1:numel(classifierset),

    i = classifierset(ii);

    fprintf('Classifier %d\n',i);

    % check trans_types
    nfixtranstypes = 0;
    fns = fieldnames(classifiers{i}.windowfeaturescellparams);
    for j = 1:numel(fns),
      fn = fns{j};
      k = find(strcmp(classifiers{i}.windowfeaturescellparams.(fn),'trans_types'));
      if ~isempty(k),
        k = k+1;
        if ~iscell(classifiers{i}.windowfeaturescellparams.(fn){k}),
          classifiers{i}.windowfeaturescellparams.(fn){k} = {classifiers{i}.windowfeaturescellparams.(fn){k}}; %#ok<CCAT1>
          nfixtranstypes = nfixtranstypes+1;
        end
      end
      for k = 2:2:numel(classifiers{i}.windowfeaturescellparams.(fn)),
        if ~iscell(classifiers{i}.windowfeaturescellparams.(fn){k}),
          continue;
        end
        l = find(strcmp(classifiers{i}.windowfeaturescellparams.(fn){k},'trans_types'),1);
        if isempty(l),
          continue;
        end
        l = l+1;
        if ~iscell(classifiers{i}.windowfeaturescellparams.(fn){k}{l}),
          classifiers{i}.windowfeaturescellparams.(fn){k}{l} = {classifiers{i}.windowfeaturescellparams.(fn){k}{l}}; %#ok<CCAT1>
          nfixtranstypes = nfixtranstypes+1;
        end
      end
    end
    [~,name] = fileparts(classifierfiles{i});
    if nfixtranstypes > 0,
      warning('Fixed %d trans_types issues for classifier %s',nfixtranstypes,name);
    end

    if ~isfield(classifiers{i},'featurenames'),

      warning('featurenames field not found in classifier %s (old classifier format), computing featurenames using ComputeWindowFeatures',name);

      % all per-frame features for this classifier
      pffs_curr = fieldnames(classifiers{i}.windowfeaturesparams);

      % list of all window features for this classifier
      feature_names = {};
      for j = 1:numel(pffs_curr),
        fn = pffs_curr{j};
        [~,feature_names_curr] = ComputeWindowFeatures([0,0],...
          classifiers{i}.windowfeaturescellparams.(fn){:});
        feature_names_curr = cellfun(@(x) [{fn},x],feature_names_curr,'UniformOutput',false);
        feature_names = [feature_names,feature_names_curr]; %#ok<AGROW>
      end

      % which features are actually used
      dims = [classifiers{i}.classifier.dim];
      feature_names = feature_names(dims);

      wfs_per_classifier{ii} = feature_names;

    else
      dims = [classifiers{i}.classifier.dim];
      feature_names =  classifiers{i}.featurenames(dims);
      wfs_per_classifier{ii} = feature_names;

    end

    % put these in with the rest of the classifiers' window features
    wfidx_per_classifier{i} = nan(1,numel(dims));
    for j = 1:numel(feature_names),
      idxcurr = find(WindowFeatureNameCompare(feature_names{j},wfs),1);
      if isempty(idxcurr),
        idxcurr = numel(wfs)+1;
        wfs{idxcurr} = feature_names{j}; %#ok<AGROW>
      end
      wfidx_per_classifier{i}(j) = idxcurr;
    end

  end
  wf2pff = cellfun(@(x)x{1},wfs,'UniformOutput',false);
  [pffs,~,wf2pffidx] = unique(wf2pff);
  pfidx_per_classifier = cell(1,nclassifiers_curr);
  for ii = 1:numel(classifierset),
    pfidx_per_classifier{ii} = wf2pffidx(wfidx_per_classifier{ii});
  end

  % create parameters for ComputeWindowFeatures which will compute all these
  % features
  fprintf('Converting window feature names to ComputeWindowFeatures parameters...\n');
  windowfeaturescellparams = struct;
  for pfi = 1:numel(pffs),
    pf = pffs{pfi};
    wfidx = wf2pffidx==pfi;
    windowfeaturescellparams.(pf) = WindowFeatureName2Params(wfs(wfidx));
  end

  % make a temporary config file containing all features
  configparams_global = configparams{classifierset(1)};
  configparams_global.behaviors.names = 'JAABADetect_DummyBehavior';
  configparams_global.featureparamlist = struct;
  for i = 1:numel(pffs),
    configparams_global.featureparamlist.(pffs{i}) = struct;
  end
  configfile_global = tempname();
  % SaveXMLParams(configparams_global,configfile_global,'params');
  save(configfile_global,'-struct','configparams_global');

  % initialize with temporary configfile
  data = JLabelData(configfile_global,'openmovie',false);
  [success,msg] = data.AddExpDirNoPreload(expdir);
  if ~success,
    error(msg);
  end
  expi = 1;
  % check if all the data inecessary is there
  if ~data.filesfixable,
    fprintf('Experiment %s is missing required files that cannot be generated within this interface.',expdir);
    return;
  end

  % generate missing files if possible
  if data.filesfixable && ~data.allfilesexist,
    [success,msg] = data.GenerateMissingFiles(data.nexps,false);
    if ~success,
      fprintf('Error generating missing required files for experiment %s: %s. Removing...',expdir,msg);
      return;
    end
  end

  isfirst = true;
  for i = classifierset,
    scores{i} = cell(1,data.nflies_per_exp(expi));
  end

  for flies = 1:data.nflies_per_exp(expi),

    fprintf('Classifying flies %s...\n',mat2str(flies));

    file = data.GetPerframeFiles(expi);
    for j = 1:numel(file),
      if ~exist(file{j},'file'),
        error('Per-frame data file %s does not exist',file{j});
      end
      tmp = load(file{j});
      ndx = find(strcmp(data.allperframefns{j},pffs));
      perframedata{ndx} = tmp.data{flies(1)};
    end

    % frames for this fly
    tStart = 1;
    off = 1-data.firstframes_per_exp{expi}(flies);
    tEnd = data.endframes_per_exp{expi}(flies)+off;
    for i = classifierset,
      scores{i}{flies} = nan(1,tEnd-off);
    end

    if tEnd < 2,
      fprintf('Trajectory length too short. Not classifying\n');
      continue;
    end

    % one block of frames at a time
    for t0 = tStart:blockSize:tEnd,

      t1 = min(t0+blockSize-1,tEnd);

      % compute all window features
      X = single([]);
      feature_names = {};
      % loop through per-frame features
      for j = 1:numel(pffs),
        fn = pffs{j};

        t11 = min(t1,numel(perframedata{j}));
        [x_curr,feature_names_curr] = ...
          ComputeWindowFeatures(perframedata{j},windowfeaturescellparams.(fn){:},'t0',t0,'t1',t11);

        if t11 < t1,
          x_curr(:,end+1:end+t1-t11) = nan;
        end

        % add the window data for this per-frame feature to X
        nold = size(X,1);
        nnew = size(x_curr,2);
        if nold > nnew,
          warning('Number of examples for per-frame feature %s does not match number of examples for previous features',fn);
          x_curr(:,end+1:end+nold-nnew) = nan;
        elseif nnew > nold && ~isempty(X),
          warning('Number of examples for per-frame feature %s does not match number of examples for previous features',fn);
          X(end+1:end+nnew-nold,:) = nan;
        end
        X = [X,single(x_curr)']; %#ok<AGROW>
        % add the feature names
        feature_names = [feature_names,cellfun(@(s) [{fn},s],feature_names_curr,'UniformOutput',false)]; %#ok<AGROW>
      end

      % find indexes for each classifier into feature_names
      if isfirst,

        classifiers_indexed = cell(1,nclassifiers_curr);

        for ii = 1:nclassifiers_curr,
          i = classifierset(ii);
          wfidx_per_classifier{ii} = nan(1,numel(wfs_per_classifier{ii}));
          for j = 1:numel(wfs_per_classifier{ii}),
            idxcurr = find(WindowFeatureNameCompare(wfs_per_classifier{ii}{j},feature_names));
            if numel(idxcurr) ~= 1,
              error('Error matching wfs for classifier %d with window features computed',i);
            end
            wfidx_per_classifier{ii}(j) = idxcurr;
          end
          classifiers_indexed{ii} = classifiers{i}.classifier;
          for j = 1:numel(classifiers_indexed{ii}),
            classifiers_indexed{ii}(j).dim = j;
          end
        end

        isfirst = false;

      end

      % apply the classifiers
      for ii = 1:nclassifiers_curr,
        i = classifierset(ii);
        Xcurr = X(:,wfidx_per_classifier{ii});
        % there's some weird offseting in the stored scores files
        scores{i}{flies}(t0-off:t1-off) = myBoostClassify(Xcurr,classifiers_indexed{ii});
      end

      fprintf('%d%% done.\n',round(100*(t1-tStart)/(tEnd-tStart)));

    end

  end

  % save scores
  for ii = 1:nclassifiers_curr,
    i = classifierset(ii);
    data.SetPostprocessingParams(classifiers{i}.postprocessparams);
    allScores = struct;
    allScores.scores = scores{i};
    allScores.tStart = data.firstframes_per_exp{expi};
    allScores.tEnd = data.endframes_per_exp{expi};
    for flies = 1:data.nflies_per_exp(expi),
      [i0s i1s] = get_interval_ends(allScores.scores{flies}>0);
      allScores.t0s{flies} = i0s;
      allScores.t1s{flies} = i1s;
    end

    postprocessedscoresA = cell(1,data.nflies_per_exp(expi));
    tStartAll = data.firstframes_per_exp{expi};
    tEndAll = data.endframes_per_exp{expi};
    for flies = 1:data.nflies_per_exp(expi)
      data.windowdata.scoreNorm = classifiers{i}.scoreNorm;
      postprocessedscoresA{flies} = nan(1,tEndAll(flies));
      postprocessedscoresA{flies}(tStartAll(flies):tEndAll(flies)) = ...
        data.Postprocess(scores{i}{flies}(tStartAll(flies):tEndAll(flies)));
    end
    allScores.postprocessed = postprocessedscoresA;
    allScores.postprocessparams = data.postprocessparams;
    allScores.scoreNorm = classifiers{i}.scoreNorm; %#ok<STRNU>
    scorefilename = scorefilenames{i};
    %scorefilename = ['test_',scorefilename];
    sfn = fullfile(expdir,scorefilename);

    if exist(sfn,'file'),
      [didbak,msg] = copyfile(sfn,[sfn,'~']);
      if ~didbak,
        warning('Could not create backup of %s: %s',sfn,msg);
      end
    end
    timestamp = classifiers{i}.classifierTS; %#ok<NASGU>
    classifierfilename = classifierfiles{i}; %#ok<NASGU>
    if DEBUG == 0,
      save(sfn,'allScores','timestamp','classifierfilename');
    end
  end

end

author: kristinbranson

repo: JAABA


text:
function [datatypesout]=any2byte(datatypes)

% ANY2BYTE - Mask Initialization function for UDP Byte Packing Block
% Derived from The MathWorks, Inc. MXPCANY2BYTE
% $Revision: 1 $ $Date: 2008/03/07  $

if ~isa(datatypes,'cell')
  error('datatypes argument must be a cell array');
end

datatypesout=[];
for i=1:length(datatypes)
  datatype=datatypes{i};
  if ~isa(datatype,'char')
    error('datatypes elements must be char arrays');
  end
  if strcmp(datatype,'double')
     dtypeout=0;
  elseif strcmp(datatype,'single')
    dtypeout=1;
  elseif strcmp(datatype,'int8')
    dtypeout=2;
  elseif strcmp(datatype,'uint8')
    dtypeout=3;
  elseif strcmp(datatype,'int16')
    dtypeout=4;
  elseif strcmp(datatype,'uint16')
    dtypeout=5;
  elseif strcmp(datatype,'int32')
    dtypeout=6;
  elseif strcmp(datatype,'uint32')
    dtypeout=7;
  elseif strcmp(datatype,'boolean')
    dtypeout=8;
  else
    error('the supported data types are: double, single, int, uint8, int16, uint16, int32, uint32, boolean');
  end
  datatypesout=[datatypesout,dtypeout];
end

author: Sable

repo: mcbench-benchmarks


text:
function [prm, cst, obs] = getPrmCstObs(obj, nActCon, nIt)

% Create a copy of actual con structure (as accessing obj.con is unbelievable slow!!!)
actCon = obj.con(nActCon);

% All types of design variables
varTypes = {'prm' 'cst' 'obs'};

% For each type of design variable
for i = 1:numel(varTypes)

    % Selection of variable type (prm, cst or obs)
    varType = varTypes{i};

    % If for the actual varType values are present
    if ~isempty(actCon.(varType))

        % Names of design variable groups for selected variable type
        varNames = fields(actCon.(varType));

        % For each design variable group
        for n = 1:numel(varNames)

            % Name of selected design variable group
            varName = varNames{n};

            % Parameter, constant or observation?
            switch varType

                case 'prm'

                    % Get parameter values for selected design variable group
                    if nIt == 1, x0_or_xhat = 'x0'; else x0_or_xhat = 'xhat'; end
                    prm.(varName) = obj.prm.(x0_or_xhat)(actCon.(varType).(varName));

                case 'cst'

                    % Get constant values for selected design variable group
                    cst.(varName) = obj.cst(actCon.(varType).(varName));

                case 'obs'

                    % Get observation values for selected design variable group
                    obs.(varName) = obj.obs.b(actCon.(varType).(varName));

            end

        end

    % If for the actual varType NO values are present
    else

        switch varType

            case 'prm'

                prm = [];

            case 'cst'

                cst = [];

            case 'obs'

                obs = [];

        end

    end

end

author: pglira


text:
% aa = dir('../3D_results_plot/results/3dLandmarks_proj_resnet50/*.mat');
clear all;
close all;
grdPath = './results/AFLW-2000-3D_grdth/';
grdDir = dir('./results/AFLW-2000-3D_grdth/*.jpg');
aa = textread('./keypoints.txt', '%s');
bb = dir('./results/data/DeFA/mesh/*.mat');

pts = zeros(204, 1);
for i = 1:204
    xx = aa(i);
    xx = xx{1};
    pts(i) = str2num(xx(1:end-1));
end
pts = reshape(pts, 3, 68) + 1;
index = pts(3,:)/3;

LL = 40;
nme_list_2DASSL = zeros(2, LL);
nme_list_3DDFA = zeros(2, LL);
for ii = 1:LL
    ii
    imgName = grdDir(ii).name;
    info = load(strcat(grdPath,strrep(imgName, 'jpg', 'mat')));
    pts1 = info.pt3d_68;
    pts1 = pts1(1:3,:);
    pts1(3,:) = pts1(3,:) -min(pts1(3,:));
    img = imread(strcat(grdPath, imgName));

    pp = textread(strcat('./results/PRNet_results/kpt/', imgName(1:end-4), '_kpt.txt'));
    pts3 = pp(:,1:3)';
    pts3(3,:) = pts3(3,:) -min(pts3(3,:));

    vertex = load(strcat('./results/2DASL_results/', strrep(imgName, 'jpg', 'mat')));
    pts5 = vertex.vertex(:,index);
    pts5 = pts5(1:3,:);
    pts5(3,:) = pts5(3,:) -min(pts5(3,:));

%     imshow(img);
%     hold on
%     for i = 1:68
%         plot(pts1(1,i), pts1(2,i), 'o');
%     end

    % NME calculation
    minx = min(pts1(1,:)); miny = min(pts1(2,:));
    maxx = max(pts1(1,:)); maxy = max(pts1(2,:));
    llen = sqrt((maxx-minx)*(maxy-miny));

    dis2 = abs(pts1 - pts3);
    dis2 = sqrt(sum(dis2.^2));
    dis2 = mean(dis2)/llen;

    dis4 = abs(pts1 - pts5);
    dis4 = sqrt(sum(dis4.^2));
    dis4 = mean(dis4)/llen;

    nme_list_PRNet(1, ii) = dis2;
    nme_list_2DASSL(1, ii) = dis4;
end

% save('nme_list', 'nme_list');
dis_PRNet = nme_list_PRNet(1,:);
dis_2DASSL = nme_list_2DASSL(1,:);

[s_dis_PRNet, index] = sort(dis_PRNet);
[s_dis_2DASSL, index] = sort(dis_2DASSL);

s_dis_PRNet = s_dis_PRNet(1:LL);
s_dis_2DASSL = s_dis_2DASSL(1:LL);

x_len = 0:1:length(s_dis_PRNet);
x_len = x_len(1:length(s_dis_PRNet));
plot(s_dis_PRNet*100, x_len, 'c', 'linewidth',2);
hold on
plot(s_dis_2DASSL*100, x_len, 'g', 'linewidth',2);

axis([0 10 0 LL])
set(gca,'XLim',[0 10]);%
set(gca,'YLim',[0 LL]);%

grid on
grid minor

ax = gca;
ax.GridColor = [0 .5 .5];
ax.GridLineStyle = '--';
ax.GridAlpha = 0.5;
ax.Layer = 'top';

h = legend('PRNet', '2DASSL', 'Location','southeast')
%set(h,'Orientation','horizon', 'Fontsize',12)
set(h,'Fontsize',12)

xlabel('NME normalized by bounding box', 'fontsize', 12)
ylabel('Number of images', 'fontsize',12)

author: XgTu


text:
function A = genErdosRenyiGraph(N,pConnect)
%%GENERDOSRENYIGRAPH Generates an instance of the Erdös-Renyi model
%                    described in [1-3].
%
%INPUT:
% N: The number of nodes used for generating the graph.
% pConnect: The probability of any chosen edge being realized.
%
%OUTPUT:
% A: The adjacency matrix for the final graph.
%
%EXAMPLE: Generates an instance of the Erdös-Renyi model.
% N = 30;
% pConnect = 0.3;
% A = genErdosRenyiGraph(N,rewireProb);
% g = graph(A);
% plot(g,"MarkerSize",degree(g))
%
%REFERENCES:
%[1] P. Erdos, A. Renyi, et al., "On the evolution of random graphs,"
%    Publ. Math. Inst. Hung. Acad. Sci, vol. 5, no. 1, pp. 17-60, 1960.
%[2] P. Erdos and A. Renyi, "On random graphs i," Publ. Math. Debrecen,
%    vol. 6, pp. 290-297, 1959.
%[3] E. N. Gilbert, "Random graphs," The Annals of Mathematical
%    Statistics, vol. 30, no. 4, pp. 1141-1144, 1959.
%
%August 2022 Codie T. Lewis, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

A = zeros(N);
for i = 1:N-1
    for j = i+1:N
        A(i,j) = rand()<pConnect;
    end
end
A = A+A';
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [optN, C, N] = sshist(x,N)
% [optN, C, N] = sshist(x,N)
%
% Function `sshist' returns the optimal number of bins in a histogram
% used for density estimation.
% Optimization principle is to minimize expected L2 loss function between
% the histogram and an unknown underlying density function.
% An assumption made is merely that samples are drawn from the density
% independently each other.
%
% The optimal binwidth D* is obtained as a minimizer of the formula,
% (2K-V) / D^2,
% where K and V are mean and variance of sample counts across bins with width D.
% Optimal number of bins is given as (max(x) - min(x)) / D*.
%
% For more information, visit
% http://2000.jukuin.keio.ac.jp/shimazaki/res/histogram.html
%
% Original paper:
% Hideaki Shimazaki and Shigeru Shinomoto
% A method for selecting the bin size of a time histogram
% Neural Computation 19(6), 1503-1527, 2007
% http://dx.doi.org/10.1162/neco.2007.19.6.1503
%
% Example usage:
% optN = sshist(x); hist(x,optN);
%
% Input argument
% x:    Sample data vector.
% N (optinal):
%       A vector that specifies the number of bins to be examined.
%       The optimal number of bins is selected from the elements of N.
%       Default value is N = 2:50.
%       * Do not search binwidths smaller than a sampling resolution of data.
%
% Output argument
% optN: Optimal number of bins.
% N:    Bin numbers examined.
% C:    Cost function of N.
%
% See also SSKERNEL
%
% Copyright (c) 2009 2010, Hideaki Shimazaki All rights reserved.
% http://2000.jukuin.keio.ac.jp/shimazaki

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters Setting
x = reshape(x,1,numel(x));
x_min = min(x);
x_max = max(x);

if nargin < 2
    buf = abs(diff(sort(x)));
    dx = min(buf(logical(buf ~= 0)));
    N_MIN = 2;              % Minimum number of bins (integer)
                            % N_MIN must be more than 1 (N_MIN > 1).
    N_MAX = min(floor((x_max - x_min)/(2*dx)),50);
                            % Maximum number of bins (integer)
    N = N_MIN:N_MAX;        % # of Bins
end

SN = 30;                    % # of partitioning positions for shift average
D = (x_max - x_min) ./ N;   % Bin Size Vector

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computation of the Cost Function
Cs = zeros(length(N),SN);
for i = 1: length(N)

       shift = linspace(0,D(i),SN);
       for p = 1 : SN
               edges = linspace(x_min+shift(p)-D(i)/2,...
                        x_max+shift(p)-D(i)/2,N(i)+1);   % Bin edges

               ki = histc(x,edges);               % Count # of events in bins
               ki = ki(1:end-1);

               k = mean(ki);                      % Mean of event count
               v = sum( (ki-k).^2 )/N(i);         % Variance of event count

               Cs(i,p) = ( 2*k - v ) / D(i)^2;    % The Cost Function
       end

end
C = mean(Cs,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Optimal Bin Size Selectioin
[Cmin idx] = min(C);
optN = N(idx);                         % Optimal number of bins
%optD = D(idx);                         % *Optimal binwidth
%edges = linspace(x_min,x_max,N(idx));  % Optimal segmentation


text:
function f = bohach1_xy ( x, y )

%*****************************************************************************80
%
%% BOHACH1_XY evaluates the Bohachevsky function #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Zbigniew Michalewicz,
%    Genetic Algorithms + Data Structures = Evolution Programs,
%    Third Edition,
%    Springer Verlag, 1996,
%    ISBN: 3-540-60676-9,
%    LC: QA76.618.M53.
%
%  Parameters:
%
%    Input, real X, Y, the argument of the function.
%
%    Output, real F, the value of the function at X.
%
  f =       x * x - 0.3 * cos ( 3.0 * pi * x ) ...
    + 2.0 * y * y - 0.4 * cos ( 4.0 * pi * y ) ...
    + 0.7;

  return
end


text:
function [y,r]=smooth(x,twice)
% SMOOTH will compute the non-linear running medians of
%     the input vector.  First it smooths the data, compute
%     the residuals, smooth the residuals, and add this back
%     to the first smooth.
%
% Usage: y = smooth(x, OPTIONAL twice?)
%  where x         is the input vector of minimum length 3
%        twice?    whether to do this twice (default = 1)
%        y         is the smoothed vector
%
% [y,r] = smooth(x, OPTIONAL twice?)
%     r  is the final residuals where x = y + r
%
% See also: R3R, SR, R3RSR, H

%
% Huy Le,  Massachusetts Institute of Technology.   1997
%
% See John Wilder Tukey, "EXPLORATORY DATA ANALYSIS".
% Addison-Wesley Publishing Co. 1977. Chpt 7, 16.

format long e;
xlen = length(x); x = x(:); y = x; r = x - y; batch = 1;
if nargin<2, twice = 1; end
if nargout<2, res = 0; end
if ~all(isfinite(x)), error('Input must have finite values.'); end
if xlen<3, error('Input vector must be at least 3 elements long.'); end
% first, run the 3RSR of vector x
xtemp1 = R3RSR(x,batch);
% first hanning
xtemp2 = H(xtemp1,batch);
% second hanning
xtemp1 = H(xtemp2,batch);
if twice,
   % compute first residuals
   rtemp = x - xtemp1;
   % run 3RSR on first residuals
   rtemp1 = R3RSR(rtemp,batch);
   % hanning of first residuals
   rtemp2 = H(rtemp1,batch);
   % sum of smooths
   xtemp = xtemp1 + rtemp2;
else
   xtemp = xtemp1;
end

for ii = 1:xlen,
   y(ii) = xtemp(ii);
   r(ii) = x(ii) - y(ii);
end

return


text:
%% This file is the main file of using the SINDy-PI method to
% infer the Yeast Glycolysis Model. We will figure out what is the minimum
% data length needed for SINDy-PI to accurately discover the six-th state
% of the Yeast Glycolysis Model. This file will be used to swipe through
% different data length.
%
% Date: 2019/06/13
% Coded By: K

%% Close all, clear all, clc
close all;clear all; clc;
set(0,'defaulttextInterpreter','latex')
addpath('Functions')
addpath('Datas')
%% Define some parameters
% Define whehter you have control, if you have it, please define it
n_control=0;u=0;

% Run the ODE files and gather the simulation data.
% (We use the same data for both the iSINDy method and SINDy-PI method for better comparision)
% (Baseline data for the data length comparison of state 1,2,3,4,5,6,7)
% The data is noise clean. It is generated using 900 different initial
% conditions with 51 points of each initial condition.
load('TrainingData.mat')

% Choose whether you want to display actual ODE or not
disp_actual_ode=1;

% If the ODEs you want to display is the actual underlyting dynamics of the
% system, please set actual as 1
actual=1;

% Define how many states we have in our example
n_state=7;

% Print the actual ODE we try to discover
digits(4)
Print_ODEs(@(t,y)YeastGlycolysis_ODE(t,y),n_state,n_control,disp_actual_ode,actual);

% Create symbolic states
dz=sym('dz',[n_state,1]);

% Now we first create the parameters of the function right hand side
Highest_Poly_Order_Guess=0;
Highest_Trig_Order_Guess=0;
Highest_U_Order_Guess=0;

% Then create the right hand side library parameters
Highest_Trig_Order=0;
Highest_U_Order=0;
Highest_dPoly_Order=1;

% Set the parameter normLib=1 to normalize the librry
normLib=1;

% Set how amny iterations you want for the sparse regression
N_iter=10;

% Set whether you want to display the ODE or not
disp=0;

% Set the library size
Highest_Poly_Order_Lib=[6;6;3;3;3;6;3];

% Determine how many percent of data you want.
percent_start=1;d_percent=0.005;percent_end=12;

% Determine which states you want to discover
Which_State=6;

% Determine the left hand side guess number
if Which_State==1 || Which_State==2 || Which_State==6
    LHS_Num=2;
else
    LHS_Num=1;
end

% Get the poly-order
Highest_Poly_Order=Highest_Poly_Order_Lib(Which_State);

% Determine which LHS you want to use, the first one or the second one
LHS_Pin=1;
%%
tic
% Create the library data using all the data points
[SINDy_Data_Full,SINDy_Struct]=SINDyLib(xt,dxt(:,Which_State),Which_State,u,Highest_Poly_Order,Highest_Trig_Order,Highest_U_Order,Highest_dPoly_Order);
fprintf('\n\t Original library creation finished, using %i seconds...\n',toc)

% Create left hand side guess
[LHS_Data_Full,LHS_Sym]=GuessLib(xt,dxt(:,Which_State),Which_State,u,Highest_Poly_Order_Guess,Highest_Trig_Order_Guess,Highest_U_Order_Guess);

tic
% Create the right hand side, exclude the guess from SINDy library
[RHS_Data_Full,RHS_Struct]=ExcludeGuess(SINDy_Data_Full,SINDy_Struct,LHS_Sym{LHS_Pin});
LHS_Data_Full_Dum=LHS_Data_Full(:,LHS_Pin);
fprintf('\t Library without LHS guess creation finished, using %i seconds...\n',toc)

% Determine sparsity parameter
lambda=0.1;

% Create the new directory to save the result
FolderName=strcat('Result_DL_SINDy_Data_Length_Compare_State_',num2str(Which_State),'_LHS_Guess_',num2str(LHS_Pin),'_SwipeLength_',num2str(lambda));
[fld_status, fld_msg, fld_msgID]=mkdir(FolderName);

% Set up dummy variables for parfor
LHS_Sym_Dum=LHS_Sym{LHS_Pin};
dz_dum=dz(Which_State);
%% determine the percentage you want to use
if LHS_Pin==1
    Percent_Iter=[0.30 0.31 0.32 0.33 0.34 0.345 0.35 0.355 0.36 0.365 0.37 0.38 0.39 0.40];
else
    Percent_Iter=[0.02 0.03 0.04 0.05 0.055 0.0575 0.06 0.0625 0.065 0.0675 0.07 0.0725 0.075 0.08 0.09 0.1 0.2];
end

for per=1:size(Percent_Iter,2)
    percent=Percent_Iter(per);

    % Start!
    for Total_Run=1:10
        fprintf('\n \n Get the result for the %i time...\n',Total_Run)

        % Print the state we are testing
        fprintf('\n \t Calculating the %i expression...\n',Which_State)

        % Print the left hand side that we are testing
        fprintf('\t Testing the left hand side as %s:\n',char(LHS_Sym_Dum))

        % Define the new data length
        new_length=round((percent)*length(xt));

        % Shuffel the original data
        Sequence=randperm(size(SINDy_Data_Full,1));
        Sequence_Trimed=Sequence(1:round((percent)*length(xt)));

        RHS_Data=RHS_Data_Full(Sequence_Trimed,:);
        LHS_Data=LHS_Data_Full_Dum(Sequence_Trimed,1);

        fprintf('\n \t Data preparation finished, using %i seconds...\n',toc)

        % We fix the lambda and change the data usage.
        tic
        fprintf('\n \t Testing the percentage as %i ...\n',(percent*100))

        % Perform the sparse regression problem
        Xi=sparsifyDynamics_simplified(RHS_Data,LHS_Data,lambda,N_iter,normLib);
        fprintf('\t Uses %i seconds...\n',toc)

        % Save the calculation result of current iteration
        fprintf('\n\t Saving the result...\n')
        tic
        cc=clock;
        ResultName=strcat(FolderName,'/DL_SINDY_Data_Length_',num2str(Total_Run),'__','LHS',num2str(LHS_Pin),'_',num2str(cc(3)),'_',num2str(cc(4)),'_',num2str(cc(5)),'_',num2str(round(cc(6))),'_P3','.mat');
        save(ResultName,'Xi','lambda','percent')
        fprintf('\n\t Saving finished! Using %i seconds...\n',toc)
    end
end

author: dynamicslab

repo: SINDy-PI


text:
function [clpot, loglik] = enter_soft_evidence(engine, CPDpot, observed, pot_type, filter)
% ENTER_SOFT_EVIDENCE Add the specified soft evidence to the network (bk)
% [clpot, loglik] = enter_soft_evidence(engine, CPDpot, observed, pot_type, filter)

[ss T] = size(CPDpot);
C = length(engine.clusters);
S = length(engine.separators);
Q = length(cliques_from_engine(engine.sub_engine));
Q1 = length(cliques_from_engine(engine.sub_engine1));
clpot = cell(Q,T);
alpha = cell(C,T);

% Forwards
% The method is a generalization of the following HMM equation:
% alpha(j,t) = normalise( (sum_i alpha(i,t-1) * transmat(i,j)) * obsmat(j,t) )
% where alpha(j,t) = Pr(Q(t)=j | y(1:t))
t = 1;
[clpot(1:Q1,t), logscale(t)] = enter_soft_evidence(engine.sub_engine1, engine.clq_ass_to_node1(:), ...
					   CPDpot(:,1), find(observed(:,1)), pot_type);
for c=1:C
  k = engine.clq_ass_to_cluster1(c);
  alpha{c,t} = marginalize_pot(clpot{k,t}, engine.clusters{c});
end

%=== FH: For each separator s, divide some cluster potential by s's potential
alpha_orig = alpha(:,t);
for s=1:S
  c = engine.cluster_ass_to_separator(s);
  alpha{c,t} = divide_by_pot(alpha{c,t}, marginalize_pot(alpha_orig{c}, engine.separators{s}));
end

% For filtering, clpot{1} contains evidence on slice 1 only

%fprintf('alphas t=%d\n', t);
%for c=1:8
%  temp = pot_to_marginal(alpha{c,t});
%  temp.T
%end

% clpot{t} contains evidence from slices t-1, t for t > 1
clqs = [engine.clq_ass_to_cluster(:,1); engine.clq_ass_to_node(:,2)];
for t=2:T
  pots = [alpha(:,t-1); CPDpot(:,t)];
  [clpot(:,t), logscale(t)] = enter_soft_evidence(engine.sub_engine, clqs, pots, find(observed(:,t-1:t)),  pot_type);
  for c=1:C
    k = engine.clq_ass_to_cluster(c,2);
    cl = engine.clusters{c};
    alpha{c,t} = marginalize_pot(clpot{k,t}, cl+ss); % extract slice 2 posterior
    alpha{c,t} = set_domain_pot(alpha{c,t}, cl); % shift back to slice 1 for re-use as prior
  end
  %=== FH: For each separator s, divide some cluster potential by s's potential
  alpha_orig = alpha(:,t);
  for s=1:S
    c = engine.cluster_ass_to_separator(s);
    alpha{c,t} = divide_by_pot(alpha{c,t}, marginalize_pot(alpha_orig{c}, engine.separators{s}));
  end
end

loglik = sum(logscale);

if filter
  return;
end

% Backwards
% The method is a generalization of the following HMM equation:
% beta(i,t) = (sum_j transmat(i,j) * obsmat(j,t+1) * beta(j,t+1))
% where beta(i,t) = Pr(y(t+1:T) | Q(t)=i)
t = T;
bnet = bnet_from_engine(engine);
beta = cell(C,T);
for c=1:C
  beta{c,t} = mk_initial_pot(pot_type, engine.clusters{c} + ss, bnet.node_sizes(:), bnet.cnodes(:), ...
			     find(observed(:,t-1:t)));
end
%=== FH: For each separator s, divide some cluster potential by s's potential
beta_orig = beta(:,t);
for s=1:S
  c = engine.cluster_ass_to_separator(s);
  beta{c,t} = divide_by_pot(beta{c,t}, marginalize_pot(beta_orig{c}, engine.separators{s}+ss));
end

for t=T-1:-1:1
  clqs = [engine.clq_ass_to_cluster(:,2); engine.clq_ass_to_node(:,2)];
  pots = [beta(:,t+1); CPDpot(:,t+1)];
  temp = enter_soft_evidence(engine.sub_engine, clqs, pots, find(observed(:,t:t+1)),  pot_type);
  for c=1:C
    k = engine.clq_ass_to_cluster(c,1);
    cl = engine.clusters{c};
    beta{c,t} = marginalize_pot(temp{k}, cl); % extract slice 1
    beta{c,t} = set_domain_pot(beta{c,t}, cl + ss); % shift fwd to slice 2
  end
  %=== FH: For each separator s, divide some cluster potential by s's potential
  beta_orig = beta(:,t);
  for s=1:S
    c = engine.cluster_ass_to_separator(s);
    beta{c,t} = divide_by_pot(beta{c,t}, marginalize_pot(beta_orig{c}, engine.separators{s}+ss));
  end
end

% Combine
% The method is a generalization of the following HMM equation:
% xi(i,j,t) = normalise( alpha(i,t) * transmat(i,j) * obsmat(j,t+1) * beta(j,t+1) )
% where xi(i,j,t) = Pr(Q(t)=i, Q(t+1)=j | y(1:T))
for t=1:T-1
  clqs = [engine.clq_ass_to_cluster(:); engine.clq_ass_to_node(:,2)];
  pots = [alpha(:,t); beta(:,t+1); CPDpot(:,t+1)];
  clpot(:,t+1) = enter_soft_evidence(engine.sub_engine, clqs, pots, find(observed(:,t:t+1)),  pot_type);
end
% for smoothing, clpot{1} is undefined
for k=1:Q1
  clpot{k,1} = [];
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function E = ROBOT_DEM
% Tests routines in DEM GUI
%__________________________________________________________________________

% close and clear everything
%--------------------------------------------------------------------------
try
    cd(fullfile(spm('Dir'),'toolbox','DEM'))
catch
    cd('C:\spm\toolbox\DEM');
end
close all force
if exist('DEMO.ps','file')
    delete('DEMO.ps')
end
clc

% get scripts and functions
%--------------------------------------------------------------------------
P     = importdata('DEM_demo.m');
R     = {};
E     = {};
for i = 1:length(P)
    k = strfind(P{i},'handles,');
    if ~isempty(k)
        str          = P{i}(k + 10:end - 2);
        if exist(str,'file') == 2
            R{end + 1,1} = str;
        end
    end
end

% execute scripts and functions
%--------------------------------------------------------------------------
N     = length(R);
T     = zeros(1,N);
for i = 1:N
    try

        % run routine
        %------------------------------------------------------------------
        fprintf('\nChecking %s\n',R{i})
        t     = tic; eval(R{i}); T(i) = toc(t);

        % print graphics
        %------------------------------------------------------------------
        H     = sort(get(0,'Children'));
        for j = 1:length(H);

            figure(H(j))
            axes('position',[.05 .98 .9 .02]);
            str   = [get(gcf,'Name') ': ' R{i}];
            text(0,0.5,str,'Fontsize',10,'Fontweight','Bold')
            axis off

            spm_print('DEMO.ps',gcf);

        end
        delete(H)

    catch

        % errors
        %------------------------------------------------------------------
        E{end + 1} = lasterror;

    end

    fprintf('\n\n     --------***--------   \n\n')

end

% Show failed routines
%--------------------------------------------------------------------------
for i = 1:length(E)
    disp(E{i}.message)
    disp(E{i}.stack(end - 1))
    disp('------------------------------------------------')
end

author: spm


text:
function varargout=phon2sone(varargin)
%
% For calling details please see v_phon2sone.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_phon2sone.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: phon2sone.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_phon2sone(varargin{:});
else
    v_phon2sone(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function heapsort_complexity ( )

%*****************************************************************************80
%
%% HEAPSORT_COMPLEXITY examines the time required to perform a heap sort.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 December 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HEAPSORT_COMPLEXITY\n' );
  fprintf ( 1, '  How does the time for heapsort increase with vector length N?\n' );
%
%  Get some data for N = 1 : 200.
%  Do the loop twice to avoid startup anomalies.
%
  for i = 1 : 2

    data_size = zeros ( 200, 1 );
    data_time = zeros ( 200, 1 );

    for n = 1 : 200
      a = rand ( n, 1 );
      tic;
      a_heap = r8vec_sort_heap_a ( n, a );
      data_size(n) = n;
      data_time(n) = toc;
    end

  end

  figure ( 1 )
  plot ( data_size, data_time, 'r-' )
  grid on
  xlabel ( 'Vector length N' );
  ylabel ( 'Elapsed time T' );
  title ( 'T(N), time to heap sort a vector of length N' );

  filename = 'heapsort_1to200.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot saved as "%s".\n', filename );
%
%  Get some data for N = 1, 2, 4, ..., 2^12.
%  Do the loop twice to avoid startup anomalies.
%
  for i = 1 : 2

    data_size = zeros ( 13, 1 );
    data_time = zeros ( 13, 1 );

    n = 1;

    for nlog = 0 : 12
      a = rand ( n, 1 );
      tic;
      a_heap = r8vec_sort_heap_a ( n, a );
      data_size(nlog+1) = n;
      data_time(nlog+1) = toc;
      n = n * 2;
    end

  end

  figure ( 2 )
  loglog ( data_size, data_time, 'r-o' )
  grid on
  xlabel ( 'Vector length N' );
  ylabel ( 'Elapsed time T' );
  title ( 'T(N), time to heap sort a vector of length N' );

  filename = 'heapsort_powersoftwo.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot saved as "%s".\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HEAPSORT_COMPLEXITY:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end
function a_out = r8vec_heap_d ( n, a )

%*****************************************************************************80
%
%% R8VEC_HEAP_D reorders an R8VEC into an descending heap.
%
%  Discussion:
%
%    A descending heap is an array A with the property that, for every index J,
%    A(J) >= A(2*J) and A(J) >= A(2*J+1), (as long as the indices
%    2*J and 2*J+1 are legal).
%
%  Diagram:
%
%                  A(1)
%                /      \
%            A(2)         A(3)
%          /     \        /  \
%      A(4)       A(5)  A(6) A(7)
%      /  \       /   \
%    A(8) A(9) A(10) A(11)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    A Nijenhuis and H Wilf,
%    Combinatorial Algorithms,
%    Academic Press, 1978, second edition,
%    ISBN 0-12-519260-6.
%
%  Parameters:
%
%    Input, integer N, the size of the input array.
%
%    Input, real A(N), an unsorted array.
%
%    Output, real A_OUT(N), the array has been reordered into a heap.
%
  a_out(1:n) = a(1:n);
%
%  Only nodes N/2 down to 1 can be "parent" nodes.
%
  for i = floor ( n/2 ) : -1 : 1
%
%  Copy the value out of the parent node.
%  Position IFREE is now "open".
%
    key = a_out(i);
    ifree = i;

    while ( 1 )
%
%  Positions 2*IFREE and 2*IFREE + 1 are the descendants of position
%  IFREE.  (One or both may not exist because they exceed N.)
%
      m = 2 * ifree;
%
%  Does the first position exist?
%
      if ( n < m )
        break;
      end
%
%  Does the second position exist?
%
      if ( m + 1 <= n )
%
%  If both positions exist, take the larger of the two values,
%  and update M if necessary.
%
        if ( a_out(m) < a_out(m+1) )
          m = m + 1;
        end

      end
%
%  If the large descendant is larger than KEY, move it up,
%  and update IFREE, the location of the free position, and
%  consider the descendants of THIS position.
%
      if ( a_out(m) <= key )
        break;
      end

      a_out(ifree) = a_out(m);
      ifree = m;

    end
%
%  Once there is no more shifting to do, KEY moves into the free spot IFREE.
%
    a_out(ifree) = key;

  end

  return
end
function a_sorted = r8vec_sort_heap_a ( n, a )

%*****************************************************************************80
%
%% R8VEC_SORT_HEAP_A ascending sorts an R8VEC using heap sort.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    A Nijenhuis and H Wilf,
%    Combinatorial Algorithms,
%    Academic Press, 1978, second edition,
%    ISBN 0-12-519260-6.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the array.
%
%    Input, real A(N), the array to be sorted;
%
%    Output, real A_SORTED(N), the sorted array.
%
  if ( n < 1 )
    a_sorted = [];
    return
  end

  if ( n == 1 )
    a_sorted(1) = a(1);
    return
  end
%
%  1: Put A into descending heap form.
%
  a_sorted = r8vec_heap_d ( n, a );
%
%  2: Sort A.
%
%  The largest object in the heap is in A(1).
%  Move it to position A(N).
%
  temp = a_sorted(1);
  a_sorted(1) = a_sorted(n);
  a_sorted(n) = temp;
%
%  Consider the diminished heap of size N1.
%
  for n1 = n-1 : -1 : 2
%
%  Restore the heap structure of A(1) through A(N1).
%
    a_sorted(1:n1) = r8vec_heap_d ( n1, a_sorted );
%
%  Take the largest object from A(1) and move it to A(N1).
%
    temp = a_sorted(1);
    a_sorted(1) = a_sorted(n1);
    a_sorted(n1) = temp;

  end

  return
end


text:
function struc2xls(filename,s,varargin)
%STRUC2XLS Saves a data structure to an Excel file
% STRUC2XLS (FILE,S) writes the data structure S to an Excel workbook named
% FILE in the current directory. FILE can also be entered using path
% notation to specify a different output directory. If the file does not
% exist, it will be created in the current directory. If it does exists,
% **it will be overwritten or modified without warning!** If the output
% file exists, it cannot be open in Excel. The name of the structure can
% include indexing values to export only a subset of the structure.
%
% Optional output parameters: (...,'Sheet',SHEETNAME) specifies the name of
% the worksheet. SHEETNAME mustbe a character string. If it does not exist
% in the Excel workbook it will be added. Default is 'Sheet1'.
%
% (...,'Col',C) and (...,'Row',R) specifiy the starting column and row of
% in Excel, respectively. R must be a positive integer, and C must be a
% capital letter index from 'A' to 'IV' using the Excel's column indexing
% format. See Excel if clarification is needed. Default is column A, row 1
% (cell A1 in Excel).
%
% (...,'Orientataion','O') specifies orientation of output. O is entered as
% 'V' for vertical or 'H' for horizontal. When V is entered, field names
% appear as column headers beginning in row R beginning at column C, and
% the contents of the structure are written in the rows below them. When H
% is entered, field names appear in column C beginning at row R and the
% contents of the structure are written in the columns to the right.
% Default is V.
%
% EXAMPLES:
%   Create a sample structure 's' with two fields:
%	s=%	struct('Name',{'Jon','Jonathan','Johnny'},'Score',[85,100,75]);
%
%   Basic output to file 'demo.xls' in current directory:
%	struc2xls('demo',s)
%
%   Output beginning in cell D4 of Excel file:
%	struc2xls('demo',s,'Col','D','Row',4)
%
%   Output to worksheet 'Students':
%	struc2xls('demo',s,'Sheet','Students')
%
%   Output horizontal data (header row, data in columns):
%   struc2xls('demo',s,'Orientation','H')
%
%   Output to alternate directory
%   structxls('C:\Users\Johnny\Desktop\demo',s)
%
%
%   Modified from struct2xls.m by Francisco de Castro:
%   http://www.mathworks.com/matlabcentral/fileexchange/18530
%
%   Jeff Evans Dec-22-2008

%Defaults
sheet= 'Sheet1';
col= 'A';
fstrow= 1;
orient= 'V';

%Optional arguments
if ~isempty(varargin)
    for j= 1:2:length(varargin)
        switch varargin{j}
            case 'Sheet'
                sheet= varargin{j+1};
            case 'Col'
                col= varargin{j+1};
            case 'Row'
                fstrow= varargin{j+1};
            case 'Orientation'
                orient= varargin{j+1};
            otherwise
                error ('Unrecognized argument name');
        end
    end
end

%Output range
rangeOut=strcat(col,num2str(fstrow));

%Transform to cell
c= struct2cell(s);

%Field names
names= fieldnames(s);

%Concatenate field names and data. Create output dataset in vertical or
%horizontal orientation.
if strmatch(orient,'V')==1
    out=[names';c'];
elseif strmatch(orient,'H')==1
    out=[names, c];
end
%write
xlswrite(filename,out,sheet,rangeOut);


text:
function [gd,relres,iter]=gabfirdual(Ldual,g,a,M,varargin)
%GABFIRDUAL  Compute FIR dual window
%   Usage: gd=gabfirdual(Ldual,g,a,M);
%          gd=gabfirdual(Ldual,g,a,M, varagin);
%
%   Input parameters:
%     Ldual  : Length of dual window
%     g      : Window function
%     a      : Time shift
%     M      : Number of Channels
%     alpha1 : Weight of $l^1$-norm in the time domain
%     alpha2 : Weight of $l^1$-norm in the freq. domain
%
%   Output parameters:
%     gd     : Dual window
%
%   `gabfirdual(Ldual,g,a,M)` computes an FIR window *gd* which is an
%   approximate dual window of the Gabor system defined by *g*, *a* and
%   *M*. The FIR dual window will be supported on *Ldual* samples.
%
%   This function solve a convex optimization problem that can be written
%   as:
%
%   .. gd  = argmin_x    || alpha x||_1 +  || beta Fx||_1
%
%   ..                 + || omega (x -g_l) ||_2^2 + delta || x ||_S0
%
%   ..                 + gamma || nabla F x ||_2^2 + mu || nabla x ||_2^2
%
%   ..     such that  x is a dual windows of g
%
%   .. math:: \begin{split}  \text{gd}  = & \text{arg} \min_x    \| \alpha x \|_1 +  \| \beta \mathcal{F}x\|_1  \\ &  + \| \omega (x - g_l) \|_2^2  \\ &  \delta \| x \|_{S0}+ \mu \| \nabla x \|_2^2 +\gamma \| \nabla \mathcal{F} x \|_2^2 \\ &  \text{such that } x \text{ is a dual window of }g \end{split}
%
%   **Note**: This function require the unlocbox. You can download it at
%   `<http://unlocbox.sourceforge.net>`_
%
%   The function uses an iterative algorithm to compute the approximate
%   FIR dual. The algorithm can be controlled by the following flags:
%
%     'alpha',alpha  Weight in time. If it is a scalar, it represent the
%                  weights of the entire L1 function in time. If it is a
%                  vector, it is the associated weight assotiated to each
%                  component of the L1 norm (length: Ldual).
%                  Default value is $\alpha=0$.
%                  **Warning**: this value should not be too big in order to
%                  avoid the the L1 norm proximal operator kill the signal.
%                  No L1-time constraint: $\alpha=0$
%
%     'beta',beta  Weight in frequency. If it is a scalar, it represent the
%                  weights of the entire L1 function in frequency. If it is a
%                  vector, it is the associated weight assotiated to each
%                  component of the L1 norm in frequency. (length: Ldual).
%                  Default value is $\beta=0$.
%                  **Warning**: this value should not be too big in order to
%                  avoid the the L1 norm proximal operator kill the signal.
%                  No L1-frequency constraint: $\beta=0$
%
%     'omega',omega  Weight in time of the L2-norm. If it is a scalar, it represent the
%                  weights of the entire L2 function in time. If it is a
%                  vector, it is the associated weight assotiated to each
%                  component of the L2 norm (length: Ldual).
%                  Default value is $\omega=0$.
%                  No L2-time constraint: $\omega=0$
%
%     'glike',g_l  $g_l$ is a windows in time. The algorithm try to shape
%                  the dual window like $g_l$. Normalization of $g_l$ is done
%                  automatically. To use option omega should be different
%                  from 0. By default $g_d=0$.
%
%     'mu',mu      Weight of the smooth constraint Default value is 1.
%                  No smooth constraint: $\mu=0$
%
%     'gamma',gamma  Weight of the smooth constraint in frequency. Default value is 1.
%                    No smooth constraint: $\gamma=0$
%
%     'delta',delta  Weight of the S0-norm. Default value is 0.
%                    No S0-norm: $\delta=0$
%
%     'dual'       Look for a dual windows (default)
%
%     'painless'   Construct a starting guess using a painless-case
%                  approximation. This is the default
%
%     'zero'       Choose a starting guess of zero.
%
%     'rand'       Choose a random starting phase.
%
%     'tol',t      Stop if relative residual error is less than the
%                  specified tolerance.
%
%     'maxit',n    Do at most n iterations. default 200
%
%     'print'      Display the progress.
%
%     'debug'      Display all the progresses.
%
%     'quiet'      Don't print anything, this is the default.
%
%     'fast'       Fast algorithm, this is the default.
%
%     'slow'       Safer algorithm, you can try this if the fast algorithm
%                  is not working. Before using this, try to iterate more.
%
%     'printstep',p  If 'print' is specified, then print every p'th
%                    iteration. Default value is p=10;
%
%     'hardconstraint' Force the projection at the end (default)
%
%     'softconstaint' Do not force the projection at the end
%
%   See also: gaboptdual, gabdual, gabtight, gabfirtight, gaboptdual

% Author: Nathanael Perraudin
% Date  : 18 Feb 2014

if nargin<4
  error('%s: Too few input parameters.',upper(mfilename));
end;

% l=length(varargin);
% varargin(l+1)=cellstr('support');
% varargin(l+2)=num2cell(Ldual);

[gd,relres,iter]=gabconvexopt(g,a,M,varargin{:}, 'support',Ldual);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% EUROPEAN OPTION PRICE COMPARISON (RUN SCRIPT)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Descritpion: Script to Compare Methods For European Options Under Levy Models
%              This script compares accuracy/CPU of the following Fourier methods:
%                   1) PROJ (2015)
%                   2) Carr-Madan (2000)
%                   3) CONV (2008)
%                   4) Lewis (2001) / Lipton (2002)
%                   5) Mellin Transform (Aguilar, 2019)
%                       ... More to come
%
% Author:      Justin Kirkby
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
[folder, name, ext] = fileparts(which( mfilename('fullpath')));
cd(folder);
addpath('../../PROJ/LEVY/RN_CHF')
addpath('../../PROJ/LEVY/Helper_Functions')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Step 1) CHOOSE CONTRACT/GENERAL PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
call = 0;    %For call use 1 (else, its a put)
S_0  = 4000;  %Initial price
W    = 4000;  %Strike            %NOTE: no error handling in place for extreme values of W (increase grid if strike falls outside)
r    = .01;  %Interest rate
q    = .00;  %Dividend yield
T    = 1;    %Time (in years)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Step 2) CHOOSE MODEL PARAMETERS  (Levy Models)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model = 8;   %See Models Below (e.g. model 1 is Black Scholes), and choose specific params

params = {};
if model == 1 %BSM (Black Scholes Merton)
    params.sigmaBSM = 0.15;    %CHOOSE

elseif model == 2 %CGMY
    params.C  = 0.02;
    params.G  = 5;
    params.MM = 15;
    params.Y  = 1.2;

elseif model == 3 %NIG
    params.alpha = 8.9932;
    params.beta  = -4.5176;
    params.delta = 1.1528;

elseif model == 4 %MJD (Merton Jump Diffusion)
    params.sigma  = 0.12;
    params.lam    = 0.4;
    params.muj    = -0.12;
    params.sigmaj = 0.18;

elseif model == 5 %Kou Double Expo
    params.sigma = 0.15;
    params.lam   = 3;
    params.p_up  = 0.2;
    params.eta1  = 25;
    params.eta2  = 10;

elseif model == 6 % Heston Model
    params.v_0 = 0.0175; % initial variance
    params.theta = 0.0398;   % long term variance level
    params.kappa =1.5768;   % rate of variance mean reversion
    params.sigma_v = 0.5751;   % volatility of variance
    params.rho = -0.5711;   % correlation between Brownian motions

elseif model == 7 %KoBoL
    params.c  = 0.02;
    params.lam_p  = 15;
    params.lam_m = -5;
    params.nu  = 1.2;

elseif model == 8 % Variance Gamma
    params.sigma = 0.2;
    params.nu = 0.85;
    params.theta = -0.1;

elseif model == 9 % Bilateral Gamma
    params.alpha_p = 1.18;
    params.lam_p = 10.57;
    params.alpha_m = 1.44;
    params.lam_m = 5.57;
end

modelInput = getModelInput(model, T, r, q, params);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Carr-Madan Fourier Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath('../../Fourier/CarrMadan/')
N = 2^16;
tic
price_CM = CarrMadan_European_Price_Strikes(S_0, W, modelInput.rnCHF, N, T, r, q, call);
time_CM = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  PROJ Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath('../../PROJ/LEVY/European_Options')
logN  = 12;   %Uses N = 2^logN  gridpoint
if model == 6
    L1 = 28;
elseif model == 1
    L1 = 18;
else
    L1 = 16;
end

% ----------------------
N = 2^logN;    % grid roughly centered on [c1 - alph, c1 + alph]
alpha = getTruncationAlpha(T, L1, modelInput, model);

tic
price_PROJ = PROJ_European(3, N, alpha, r, q, T, S_0, W, call, modelInput.rnCHF, modelInput.c1*T);
time_PROJ = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  CONV Fourier Method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath('../../Fourier/CONV/')
N = 2^16;
tic
price_CONV = CONV_European_Price(S_0, W, modelInput.rnCHF, T, r, call, N, alpha);
time_CONV = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Hilbert Transform
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath('../../Fourier/HilbertTransform/')
N = 2^14; Nh = 2^5;
h = 2*pi/Nh;

tic
price_Hilb = Hilbert_European_Price(h, N, r, q, T, S_0, W, call, modelInput.rnCHF);
time_Hilb = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Lewis (2001) Fourier Method  (see also Lipton (2002))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath('../../Fourier/Lewis/')
tic
price_Lewis = Lewis_European_Price(S_0, W, modelInput.rnCHF, T, r, q, call, 500);
time_Lewis = toc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Mellin Transform (J-P Aguilar, 2020)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
addpath('../../Fourier/MellinTransform/')
has_mellin = 0;
if model == 3  % NIG
    tic
    N_terms = 50;
    price_Mellin = Mellin_NIG_European_Price( S_0, W, T, r, q, call, params.alpha, params.beta, params.delta, N_terms);
    time_Mellin = toc;
    has_mellin = 1;
end
if model == 8 %  Variance Gamma
    tic
    N_terms = 24;
    price_Mellin = Mellin_VG_European_Price( S_0, W, T, r, q, call, params.sigma, params.theta, params.nu, N_terms);
    time_Mellin = toc;
    has_mellin = 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Reference Price (Using PROJ)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nref = 2^16; L1Ref = max(16, L1+4);   % Params to obtain reference price
alphaRef = getTruncationAlpha(T, L1Ref, modelInput, model);
price_Ref = PROJ_European(3, Nref, alphaRef, r, q, T, S_0, W, call, modelInput.rnCHF, modelInput.c1*T);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% COMPARE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n---------------------------------------------\n')
fprintf('Method      |    Price      |    Err   |  CPU \n')
fprintf('---------------------------------------------\n')
fprintf('Reference   | %.8f  |          |       \n', price_Ref)
fprintf('---------------------------------------------\n')
fprintf('PROJ        | %.8f  | %.2e | %.4f \n', price_PROJ, abs(price_Ref-price_PROJ), time_PROJ)
fprintf('CONV        | %.8f  | %.2e | %.4f \n', price_CONV, abs(price_Ref-price_CONV), time_CONV)
fprintf('Carr-Madan  | %.8f  | %.2e | %.4f \n', price_CM, abs(price_Ref-price_CM), time_CM)
fprintf('Lewis(2001) | %.8f  | %.2e | %.4f \n', price_Lewis, abs(price_Ref-price_Lewis), time_Lewis)
fprintf('Hilbert     | %.8f  | %.2e | %.4f \n', price_Hilb, abs(price_Ref-price_Hilb), time_Hilb)
if has_mellin == 1
   fprintf('Mellin      | %.8f  | %.2e | %.4f \n', price_Mellin, abs(price_Ref-price_Mellin), time_Mellin)
end
fprintf('---------------------------------------------\n')


text:
% epoch() - Extract epochs time locked to specified events from continuous EEG data.
%
% Usage:
%   >> epocheddata = epoch( data, events, timelim);
%   >> [epocheddata, newtime, indices, rerefevent, rereflatencies ] = ...
%                      epoch( data, events, timelim, 'key1', value1, ... );
%
% Inputs:
%   data       - input data (chan,frames). In the case, data is a
%                3D array (chan, frames, epochs), epochs are extracted
%                only if their time windows fall within existing
%                pre-existing epochs.
%   events     - vector events (expressed in seconds or points)
%   timelim    - [init end] in second or points centered
%                on the events (i.e. [-1 2])
%
% Optional inputs:
%   'srate'      - sampling rate in Hz for events expressed in seconds
%   'valuelim'   - upper and lower limit of values that a trial should not
%                overpass. If one positive value is given, consider the
%                opposite for lower bound. Given values are also consider
%                outlier (if equal the trial is rejected). Default: none.
%   'verbose'    - ['on'|'off']. Default is 'on'.
%   'allevents'  - event vector containing the latencies of all events
%                (not only those used for epoching). The function
%                return an array 'rerefevent' which contain the latency
%                of these events in each trials (assuming the events are
%                included in the trial time window). These events can
%                be in point or second but they must be in the same format
%                as the events used for epoching.
%   'alleventrange' - for event selection, defines a time range [start end] (in
%                seconds or data points) relative to the time-locking events.
%                Default is same as 'timelim'.
%
% Outputs:
%   epocheddata - output (chan, frames, epochs)
%   indices     - indices of accepted events
%   newtime     - new time limits. See notes.
%   rerefevent  - re-referenced event cell array (size nbepochs) of array
%                 indices for each epochs (note that the number of events
%                 per trial may vary).
%   rereflatencies - re-referenced latencies event cell array (same as above
%                 but indicates event latencies in epochs instead of event
%                 indices).
%
% Note: maximum time limit will be reduced by one point with comparison to the
%       input time limits. For instance at 100 Hz, 3 seconds last 300 points,
%       but if we assign time 0 to the first point, then we must assign time
%       299 to the last point.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001
%
% See also: pop_epoch(), eeglab()

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

function [epochdat, newtime, indexes, alleventout, alllatencyout, reallim] = epoch( data, events, lim, varargin );

if nargin < 2
   help epoch;
	return;
end;
alleventout = {};

% create structure
% ----------------
if ~isempty(varargin)
   try, g = struct(varargin{:});
   catch, error('Epoch: wrong syntax in function arguments'); end;
else
    g = [];
end;

try, g.srate; 	 	     catch, g.srate = 1; end;
try, g.valuelim; 	     catch, g.valuelim =  [-Inf Inf]; end;
try, g.verbose; 	     catch, g.verbose = 'on'; end;
try, g.allevents; 	     catch, g.allevents = []; end;
try, g.alleventrange; 	 catch, g.alleventrange = lim; end;

% computing point limits
% ----------------------
reallim(1) = round(lim(1)*g.srate);   % compute offset
reallim(2) = round(lim(2)*g.srate-1); % compute offset

% epoching
% --------
fprintf('Epoching...\n');
newdatalength = reallim(2)-reallim(1)+1;

epochdat = zeros( size(data,1), newdatalength, length(events) );
g.allevents =  g.allevents(:)';
datawidth  = size(data,2)*size(data,3);
dataframes = size(data,2);
indexes = zeros(length(events),1);
alleventout = {};
alllatencyout = {};
for index = 1:length(events)
	pos0 = floor(events(index)*g.srate); % offset of time locking event
	posinit = pos0+reallim(1); % compute offset
	posend  = pos0+reallim(2); % compute offset

   if floor((posinit-1)/dataframes) == floor((posend-1)/dataframes) & posinit >= 1 & posend <= datawidth % test if within boundaries
      epochdat(:,:,index) = data(:,posinit:posend);
      if ~isinf(g.valuelim(1)) | ~isinf(g.valuelim(2))
          if (max(epochdat(:,:,index)) > g.valuelim(1)) & ...
                  (max(epochdat(:,:,index)) < g.valuelim(2))
              indexes(index) = 1;
          else
              switch g.verbose, case 'on', fprintf('Warning: event %d out of value limits\n', index); end;
          end;
      else
          indexes(index) = 1;
      end;
   else
      switch g.verbose, case 'on', fprintf('Warning: event %d out of data boundary\n', index); end;
   end;

   % rereference events
   % ------------------
   if ~isempty(g.allevents)
        posinit = pos0 + g.alleventrange(1)*g.srate; % compute offset
        posend  = pos0 + g.alleventrange(2)*g.srate; % compute offset
        eventtrial = intersect( find(g.allevents*g.srate >= posinit),  find(g.allevents*g.srate <= posend) );
        alleventout{index} = eventtrial;
        alllatencyout{index} = g.allevents(eventtrial)*g.srate-pos0;
   end;
end;
newtime(1) = reallim(1)/g.srate;
newtime(2) = reallim(2)/g.srate;

indexes = find(indexes == 1);
epochdat = epochdat(:,:,indexes);
if ~isempty(alleventout)
    alleventout = alleventout(indexes);
    alllatencyout= alllatencyout(indexes);
end;
reallim = reallim*g.srate;
return;

%% GENERATION OF NAN IN ARRAYS (old implementation)
%% ------------------------------------------------
%alleventout(index,1:length(eventtrial) ) = eventtrial+1000000;
%%then replace all zeros by Nan and subtract 1000000
%if ~isempty(alleventout)
%    alleventout( find( alleventout == 0) ) = nan;
%    alleventout = alleventout(indexes,:) - 1000000;
%   alllatencyout( find( alllatencyout == 0) ) = nan;
%   alllatencyout = alllatencyout(indexes,:) - 1000000;
%end;

function res = lat2point( lat, srate, pnts);

res = lat*srate+1 + (epoch_array-1)*pnts;

author: PatternRecognition

repo: OpenBMI


text:
function output = calloptiqsopt(interfacedata)

% Standard input interface
options = interfacedata.options;
F_struc = interfacedata.F_struc;
c       = interfacedata.c;
K       = interfacedata.K;
lb      = interfacedata.lb;
ub      = interfacedata.ub;

showprogress('Calling QSOPT',options.showprogress);

if isempty(F_struc)
    A = [];
    b = [];
else
    A =[-F_struc(K.f+1:end,2:end);-F_struc(1:K.f,2:end)];
    b =[F_struc(K.f+1:end,1);F_struc(1:K.f,1)];
end

opts = options.qsopt;
opts.nin = length(b)-K.f;
opts.display = options.verbose;

if options.savedebug
    save qsoptdebug c A b  lb ub opts
end

solvertime = tic;
[x,fval,exitflag] = qsopt(full(c), A, full(b), full(lb), full(ub),opts);
solvertime = toc(solvertime);

% No duals
D_struc = [];

switch exitflag
    case 1
        problem = 0;
    case 0
        problem = 3;
    case -1
        problem = 1;
    case -2
        problem = 11;
    otherwise
        problem = -1;
end
infostr = yalmiperror(problem,'QSOPT');

% Save all data sent to solver?
if options.savesolverinput
    solverinput.A = A;
    solverinput.b = b;
    solverinput.f = c;
    solverinput.lb = lb;
    solverinput.ub = ub;
else
    solverinput = [];
end

% Save all data from the solver?
if options.savesolveroutput
    solveroutput.x = x;
    solveroutput.fval = fval;
    solveroutput.exitflag = exitflag;
else
    solveroutput = [];
end

% Standard interface
output = createOutputStructure(x,D_struc,[],problem,infostr,solverinput,solveroutput,solvertime);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% Jason Joseph Rebello
% Carnegie Mellon University (Jan 2012 - May 2013)
% Masters in Electrical & Computer Engineering
% Logistic Regression with Regularization

% Run LogisitcRegression.m
% This program uses Logistic regression to classify handwritten digits. The
% data is comprised of a part of the MNIST dataset. 5000 examples are used.

clear all;
clc;
close all;

%% Load Training Data
t = cputime;
fprintf('Loading and Visualizing Data')

load('data.mat');

% Randomly select 100 data points to display
sel = randperm(size(X, 1));
sel = sel(1:100);

displayData(X(sel, :)); % function taken from Andrew Ng machine learning course
fprintf('...done\n');

%% Split data into training and testing
fprintf('Splitting Data into training and testing sets');
[XTrain XTest yTrain yTest] = splitData(X, y);
fprintf('...done\n\n');

%% Initialize parameters
fprintf('Initializing parameters');
m = size(X, 1); % number of examples
lambda = 0.1; % regularization parameter
numLabels = size(unique(y),1); % number of labels
fprintf('...done\n');

%% Training Logistic Regression classifier
fprintf('Training One-vs-All Logistic Regression');

theta = LRClassifier(XTrain, yTrain, numLabels, lambda);
fprintf('...done\n');

%% Predict numbers
prediction = predict(theta, XTest);

%% Calculate Accuracy over the training data
fprintf('\nTest Set Accuracy: %f\n', mean(double(prediction == yTest)) * 100);

fprintf('Program executed in %f seconds or %f minutes\n', cputime-t, (cputime-t)/60);