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

    properties
        lvdData LvdData

        vars OptimizationVariableSet
%         objFcn(1,1) AbstractObjectiveFcn = NoOptimizationObjectiveFcn()
        objFcn(1,1) AbstractObjectiveFcn = CompositeObjectiveFcn()
        constraints(1,1) ConstraintSet =  ConstraintSet()

        %Optimization Algo Selection
        optAlgo(1,1) LvdOptimizerAlgoEnum = LvdOptimizerAlgoEnum.Fmincon;

        %Optimizers
        fminconOpt(1,1) FminconOptimizer = FminconOptimizer();
        patternSearchOpt(1,1) PatternSearchOptimizer = PatternSearchOptimizer();
        nomadOpt(1,1) NomadOptimizer = NomadOptimizer();
        ipoptOpt(1,1) IpOptOptimizer = IpOptOptimizer();
        surragateOpt(1,1) SurrogateOptimizer = SurrogateOptimizer();

        %Gradient Calc Algo Selection
        gradAlgo(1,1) LvdOptimizerGradientCalculationAlgoEnum = LvdOptimizerGradientCalculationAlgoEnum.BuiltIn;

        %Gradient Calculation Algos
        builtInGradMethod(1,1) BuiltInGradientCalculationMethod = BuiltInGradientCalculationMethod();
        customFiniteDiffsCalcMethod(1,1) CustomFiniteDiffsCalculationMethod = CustomFiniteDiffsCalculationMethod();
    end

    methods
        function obj = LvdOptimization(lvdData)
            obj.lvdData = lvdData;

            obj.vars = OptimizationVariableSet(obj.lvdData);
            obj.objFcn = CompositeObjectiveFcn(GenericObjectiveFcn.empty(1,0), ObjFcnDirectionTypeEnum.Minimize, ObjFcnCompositeMethodEnum.Sum, lvdData.optimizer, lvdData);
            obj.constraints = ConstraintSet(obj, lvdData);

            obj.optAlgo = LvdOptimizerAlgoEnum.Fmincon;
            obj.fminconOpt = FminconOptimizer();
            obj.patternSearchOpt = PatternSearchOptimizer();
            obj.nomadOpt = NomadOptimizer();
            obj.ipoptOpt = IpOptOptimizer();
            obj.surragateOpt = SurrogateOptimizer();

            obj.builtInGradMethod = BuiltInGradientCalculationMethod();
            obj.customFiniteDiffsCalcMethod = CustomFiniteDiffsCalculationMethod();
        end

        function [exitflag, message] = optimize(obj, writeOutput, callOutputFcn, hLvdMainGUI)
            obj.vars.sortVarsByEvtNum();
            optimizer = obj.getSelectedOptimizer();

            [x0All, actVars, ~] = obj.vars.getTotalScaledXVector();

            if(isempty(x0All) && isempty(actVars))
                uialert(hLvdMainGUI, 'There are no optimization variables enabled in this mission.  Optimization requires at least one variable.  Please enable at least one variable to continue with optimization.', 'Launch Vehicle Designer', 'Icon','error');

                return;
            end

            [exitflag, message] = optimizer.optimize(obj, writeOutput, callOutputFcn, hLvdMainGUI);
        end

        function [exitflag, message] = consoleOptimize(obj)
            global options_gravParamType %#ok<GVMIS>

            if(isempty(options_gravParamType))
                options_gravParamType = 'kspStockLike';
            end

            writeOutput = @(varargin) disp('');
            callOutputFcn = false;

            [exitflag, message] = obj.optimize(writeOutput, callOutputFcn, []);
        end

        function optimizer = getSelectedOptimizer(obj)
            optAlgorithm = obj.optAlgo;
            optimizer = obj.getOptimizerForEnum(optAlgorithm);
        end

        function optimizer = getOptimizerForEnum(obj, optAlgorithm)
            if(optAlgorithm == LvdOptimizerAlgoEnum.Fmincon)
                optimizer = obj.fminconOpt;
            elseif(optAlgorithm == LvdOptimizerAlgoEnum.PatternSearch)
                optimizer = obj.patternSearchOpt;
            elseif(optAlgorithm == LvdOptimizerAlgoEnum.Nomad)
                optimizer = obj.nomadOpt;
            elseif(optAlgorithm == LvdOptimizerAlgoEnum.Ipopt)
                optimizer = obj.ipoptOpt;
            elseif(optAlgorithm == LvdOptimizerAlgoEnum.Surrogate)
                optimizer = obj.surragateOpt;
            else
                error('Unknown LVD optimization algorithm!');
            end
        end

        function gradAlgo = getGradAlgoForEnum(obj, gradAlgoEnum)
            if(gradAlgoEnum == LvdOptimizerGradientCalculationAlgoEnum.BuiltIn)
                gradAlgo = obj.builtInGradMethod;
            elseif(gradAlgoEnum == LvdOptimizerGradientCalculationAlgoEnum.FiniteDifferences)
                gradAlgo = obj.customFiniteDiffsCalcMethod;
            else
                error('Unknown LVD gradient algorithm!');
            end
        end

        function tf = usesParallel(obj)
            tf = obj.getSelectedOptimizer().usesParallel();
        end

        function tf = usesStage(obj, stage)
            tf = obj.objFcn.usesStage(stage);

            tf = tf || obj.constraints.usesStage(stage);
        end

        function tf = usesEngine(obj, engine)
            tf = obj.objFcn.usesEngine(engine);

            tf = tf || obj.constraints.usesEngine(engine);
        end

        function tf = usesTank(obj, tank)
            tf = obj.objFcn.usesTank(tank);

            tf = tf || obj.constraints.usesTank(tank);
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = obj.objFcn.usesEngineToTankConn(engineToTank);

            tf = tf || obj.constraints.usesEngineToTankConn(engineToTank);
        end

        function tf = usesExtremum(obj, extremum)
            tf = obj.objFcn.usesExtremum(extremum);

            tf = tf || obj.constraints.usesExtremum(extremum);
        end

        function tf = usesGroundObj(obj, grdObj)
            tf = obj.objFcn.usesGroundObj(grdObj) || obj.constraints.usesGroundObj(grdObj);
        end

        function tf = usesCalculusCalc(obj, calculusCalc)
            tf = obj.objFcn.usesCalculusCalc(calculusCalc);

            tf = tf || obj.constraints.usesCalculusCalc(calculusCalc);
        end

        function tf = usesGeometricPoint(obj, point)
            tf = obj.constraints.usesGeometricPoint(point);
        end

        function tf = usesGeometricVector(obj, vector)
            tf = obj.constraints.usesGeometricVector(vector);
        end

        function tf = usesGeometricCoordSys(obj, coordSys)
            tf = obj.constraints.usesGeometricCoordSys(coordSys);
        end

        function tf = usesGeometricRefFrame(obj, refFrame)
            tf = obj.constraints.usesGeometricRefFrame(refFrame);
        end

        function tf = usesGeometricAngle(obj, angle)
            tf = obj.constraints.usesGeometricAngle(angle);
        end

        function tf = usesGeometricPlane(obj, plane)
            tf = obj.constraints.usesGeometricPlane(plane);
        end

        function tf = usesPlugin(obj, plugin)
            tf = obj.constraints.usesPlugin(plugin);
        end
    end

    methods(Static)
        function obj = loadobj(obj)
            if(isempty(obj.vars.lvdData))
                obj.vars.lvdData = obj.lvdData;
            end
        end
    end

    methods(Access=private)
        function evtNumToStartScriptExecAt = getEvtNumToStartScriptExecAt(obj, actVars)
            evtNumToStartScriptExecAt = obj.lvdData.script.getTotalNumOfEvents();
            for(i=1:length(actVars)) %#ok<*NO4LP>
                var = actVars(i);

                if(isVarInLaunchVehicle(var, obj.lvdData) || isVarInLaunchVehicle(var, obj.lvdData))
                    varEvtNum = 1;
                else
                    varEvtNum = getEventNumberForVar(var, obj.lvdData);

                    if(isempty(varEvtNum))
                        varEvtNum = 1;
                    end
                end

                if(varEvtNum < evtNumToStartScriptExecAt)
                    evtNumToStartScriptExecAt = varEvtNum;
                end

                if(evtNumToStartScriptExecAt == 1)
                    break; %it can't go lower than 1, so we're executing the whole thing.  No reason to keep going.
                end
            end
        end

        function evtNumToEndScriptExecAt = getEvtNumToEndScriptExecAt(obj)
            [~, activeVars, ~, ~] = obj.vars.getTotalScaledXVector();

            evtNumToEndScriptExecAt = -1;
            for(i=1:length(activeVars))
                var = activeVars(i);

                if(isVarInLaunchVehicle(var, obj.lvdData) || isVarInLaunchVehicle(var, obj.lvdData))
                    varEvtNum = 1;
                else
                    varEvtNum = getEventNumberForVar(var, obj.lvdData);

                    if(isempty(varEvtNum))
                        varEvtNum = 1;
                    end
                end

                if(varEvtNum > evtNumToEndScriptExecAt)
                    evtNumToEndScriptExecAt = varEvtNum;
                end
            end

            objFcnEvts = obj.objFcn.getObjFuncEvents();
            for(i=1:length(objFcnEvts))
                objFcnEvtNum = objFcnEvts(i).getEventNum();

                if(objFcnEvtNum > evtNumToEndScriptExecAt)
                    evtNumToEndScriptExecAt = objFcnEvtNum;
                end
            end

            constrEvts = obj.constraints.getConstrEvents();
            for(i=1:length(constrEvts))
                constrEvtNum = constrEvts(i).getEventNum();

                if(constrEvtNum > evtNumToEndScriptExecAt)
                    evtNumToEndScriptExecAt = constrEvtNum;
                end
            end
        end
    end
end


text:
function [ thetaVec,sVec ] = getEllipseLengthVector( a,b, theta0,theta1 )
%% About:
% This function is used to return the length of the ellipse as a function
% of the angle theta.

%% Inputs:
% theta0, theta1: start/end angle of the arc of the ellipse (from the parametric equation of the ellipse)
% a,b: the big and the small radious of the ellipse.

% Mohammad SAFEEA 10-Nov-2017

thetaSpan = [theta0 theta1];
s0 = 0;
opts = odeset('RelTol',1e-2,'AbsTol',1e-4);
[thetaVec,sVec] = ode45(@(theta,s) ((a*a*sin(theta)*sin(theta)+b*b*cos(theta)*cos(theta))^0.5), thetaSpan, s0,opts);

end

repo: KST-Kuka-Sunrise-Toolbox


text:
%This Matlab script can be used to reproduce Figure 7.8 in the monograph:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Empty workspace and close figures
close all;
clear;

%Maximum number of pilots
pilotMax = 40;

%Define range of average number of UEs
Kaverages = [1 10 20 40];

%Define range of effective number of UEs
Krange = 0:pilotMax;

%Prepare to store simulation results
Kdistributions = zeros(length(Krange),length(Kaverages));

%% Go through all average numbers of UEs
for n = 1:length(Kaverages)

    %Compute the probabilities of having different numbers of active UEs,
    %based on a Poission distribution
    Kdistributions(:,n) = poisspdf(Krange,Kaverages(n));

    %Enforce scheduling to make sure that there cannot be more than
    %pilotMax UEs active per cell
    Kdistributions(end,n) = 1 - sum(Kdistributions(1:end-1,n));

end

%% Plot simulation results
figure;
hold on; box on;

plot(Krange,Kdistributions(:,1),'r*-','LineWidth',1);
plot(Krange,Kdistributions(:,2),'ko-','LineWidth',1);
plot(Krange,Kdistributions(:,3),'bs-','LineWidth',1);
plot(Krange,Kdistributions(:,4),'rd--','LineWidth',1);

xlabel('Number of active UEs');
ylabel('Probability');
legend('K=1  ','K=10  ','K=20  ','K=40  ','Location','NorthWest');
ylim([0 0.6]);

author: emilbjornson

repo: massivemimobook


text:
function [ babel ] = babel( Phi )
%BABEL computes the Babel function of a dictionary
% Columns of Phi must be unit norm
% The function has been defined by Joel A Tropp

% Let us compute the Gram matrix
G = Phi' * Phi;
% Let us take element wise absolute
absG = abs(G);
% Let us sort each row in descending order
GS = sort(absG, 2,'descend');
rowSums = cumsum(GS(:, 2:end), 2);
babel = max(rowSums);

author: indigits

repo: sparse-plex


text:
function ind = meegchannels(this, modality)
% Return indices of M/EEG channels
% FORMAT ind = meegchannels(this, modality)
%
%  this      - MEEG object
%  modality  - one of EEG, LFP, MEG (excluding planar), MEGPLANAR [optional]
%
%  ind       - row vector of M/EEG channels
%
% See also eogchannels, ecgchannels, emgchannels
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel
% $Id: meegchannels.m 5045 2012-11-09 15:53:23Z vladimir $

type = chantype(this);

if nargin == 1
    ind = find(ismember(upper(type), ...
        {'EEG', 'MEG', 'MEGMAG', 'MEGGRAD', 'MEGPLANAR', 'REF', 'REFMAG', 'REFGRAD', 'LFP'}));
else
    switch modality
        case 'EEG'
            ind = find(ismember(upper(type), {'EEG'}));
        case 'LFP'
            ind = find(ismember(upper(type), {'LFP'}));
        case 'MEG'
            ind = find(ismember(upper(type), {'MEG', 'MEGMAG', 'MEGGRAD'}));
        case 'MEGMAG'
            ind = find(ismember(upper(type), {'MEGMAG'}));
        case 'MEGPLANAR'
            ind = find(ismember(upper(type), {'MEGPLANAR'}));
        case 'MEGGRAD'
            ind = find(ismember(upper(type), {'MEGGRAD'}));
        case 'MEEG'
            ind = find(ismember(upper(type), {'EEG', 'MEG', 'MEGMAG', 'MEGGRAD', 'MEGPLANAR'}));
        otherwise
            error('Unsupported modality.');
    end
end

ind = ind(:)'; % must be row to allow to use it as loop indices

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function gpsimMapEcoliResults(model, type, expNo, scale, printResults)

% GPSIMMAPECOLIRESULTS Plot the results from the MAP script.
% FORMAT
% DESC plots results from the MAP script, printing them if the
% given flag is set to true.
% ARG model : the model for which the results are being plotted.
% ARG type : gives a descriptive label for the type of the
% experiment.
% ARG expNo : number identifier for the experiment.
% ARG scale : the scaling factor applied to normalise.
% ARG printResults : function prints the plots if set to true
% (default is false).
%
% COPYRIGHT : Pei Gao, Magnus Rattray and Neil D. Lawrence, 2008
%
% SEEALSO : demEcoliMap1

% SHEFFIELDML

type(1) = upper(type(1));
switch nargin
 case 3
  scale = 1;
  printResults = false;
 case 4
  printResults = false;
end

gene = {'dinF', 'dinI', 'lexA', 'recA', 'recN', 'ruvA', 'ruvB', ...
                  'sbmC', 'sulA', 'umuC', 'umuD', 'uvrB', 'yebG', 'yjiW' };
if model.comp{1}.numGenes == 8
  targetInd = [2 3 6 7 9 11 12 14];
  gene = gene(targetInd);
end

predt = model.comp{1}.mapt;
predExprsFull = [];
for i = 1:length(model.comp)
    predExprsFull = [predExprsFull; model.comp{i}.ypred];
end
meanPredExprs = ones(size(predExprsFull, 1), 1)*mean(predExprsFull);
scalePredExprs = sqrt(var(predExprsFull - meanPredExprs));
scaleMat = ones(size(predExprsFull, 1), 1)*scalePredExprs;
% predExprsFull = meanPredExprs + (predExprsFull-meanPredExprs)./scaleMat;
startVal = 1;
endVal = length(predt);
for i = 1:length(model.comp)
    predExprs{i} = predExprsFull(startVal:endVal, :);
    startVal = startVal + length(predt);
    endVal = endVal + length(predt);
end

modelB = model.comp{1}.B.*scale;
modelS = model.comp{1}.S.*scale;
modelAlpha = model.comp{1}.alpha.*scale;

for j = 1:length(model.comp)

  scalePred = sqrt(var(model.comp{j}.f));
  truef = [0.16 0 0.19 0.75 1];
%  truef = truef/sqrt(var(truef))*scalePred;

  % Figure 2(a) histograms;
%   B = trueParams.B./scale; % From Martino paper ... but don't know the
%                            % scale
%   B = B/mean(B)*mean(modelB); % do a rough rescaling so
%                                         % that the scales match.
%   S = trueParams.S./scale;
%   S = S/mean(S)*mean(modelS);
%   D = trueParams.D; % From Martino paper, again
%                                       % we know the scale because
%                                       % p21 is fixed to 0.8.

  % Martino f from Figure 2(b), again measured with a ruler.
  if model.comp{j}.ngParam
    nonLinearity = 'multi';
  else
    nonLinearity = model.comp{j}.nonLinearity;
  end
  switch nonLinearity
   case 'linear'
    figure, lin = plot(model.comp{j}.mapt,model.comp{j}.f, '-');
    hold on;
    bh = plot(model.comp{j}.mapt, model.comp{j}.f+2*sqrt(model.comp{j}.varF),'--');
    bh =[bh plot(model.comp{j}.mapt, model.comp{j}.f-2*sqrt(model.comp{j}.varF),'--')];
    lin = [lin plot(model.comp{j}.t(1:5), truef, 'rx')];
   case 'multi'
    figure,
    predF = exp(model.comp{j}.f);
    varFplus = exp(model.comp{j}.f+2*sqrt(model.comp{j}.varf));
    varFminus = exp(model.comp{j}.f-2*sqrt(model.comp{j}.varf));
    lin = plot(model.comp{j}.mapt, predF, '-');
    hold on;
    bh = plot(model.comp{j}.mapt, varFplus, '--');
    bh = [bh plot(model.comp{j}.mapt, varFminus, '--')];
%    truef = truef/sqrt(var(truef))*sqrt(var(predF));
    title('Inferred LexA Activity', 'fontsize', 20);
%    lin = [lin plot(model.comp{j}.t(1:5), truef, 'rx')];
   otherwise
     figure,
%
%     predF = exp(model.comp{j}.f);
%     varFplus = exp(model.comp{j}.f+2*sqrt(model.comp{j}.varf));
%     varFminus = exp(model.comp{j}.f-2*sqrt(model.comp{j}.varf));
%     lin = plot(model.comp{j}.mapt, predF-min(predF), '-');
%     hold on;
%     bh = plot(model.comp{j}.mapt, varFplus-min(predF), '--');
%     bh = [bh plot(model.comp{j}.mapt, varFminus-min(predF), '--')];
%     truef = truef/sqrt(var(truef))*sqrt(var(predF));
%

    subplot(2,1,1);
    lin = plot(model.comp{j}.mapt, model.comp{j}.f, '-');
    hold on;
    bh = plot(model.comp{j}.mapt, model.comp{j}.f+2*sqrt(model.comp{j}.varf),'--');
    bh = [bh plot(model.comp{j}.mapt, model.comp{j}.f-2* ...
                                           sqrt(model.comp{j}.varf),'--')];
    lin = [lin plot(model.comp{j}.t(1:5), truef, 'rx')];
    title('f(t) Profile');

    func = str2func(model.comp{j}.nonLinearity);
    subplot(2,1,2); lin = [lin plot(model.comp{j}.mapt, func(model.comp{j}.f), '-')];
    hold on;
    bh = [bh plot(model.comp{j}.mapt, func(model.comp{j}.f+2*sqrt(model.comp{j}.varf)),'--')];
    bh = [bh plot(model.comp{j}.mapt, func(model.comp{j}.f-2* ...
                                           sqrt(model.comp{j}.varf)),'--')];
    lin = [lin plot(model.comp{j}.t(1:5), func(truef), 'rx')];
    title('g(f) Profile');
  end

  set(bh, 'lineWidth', 2);
  set(lin, 'lineWidth', 4);
  set(lin, 'markersize', 20);
  set(gca, 'fontname', 'arial', 'fontsize', 24, 'xlim', xlim)

  kernName = model.comp{j}.kern.type;
  kernName(1) = upper(kernName(1));

  switch nonLinearity
   case 'linear'
    set(gca, 'ylim', [-2 4])
   case {'exp', 'negLogLogit'}
    set(gca, 'ylim', [0 8])
%    otherwise
%     set(gca, 'ylim', [0 2])
  end

  fileName = ['dem' type num2str(expNo) kernName ...
              nonLinearity '_profile' num2str(j) '_slide'];
  if printResults
    print('-deps', ['./results/' fileName]);
  end
  pos = get(gcf, 'paperposition');
  origpos = pos;
  pos(3) = pos(3)/2;
  pos(4) = pos(4)/2;
  set(gcf, 'paperposition', pos);
  lineWidth = get(gca, 'lineWidth');
  set(gca, 'lineWidth', lineWidth*2);
  if printResults
    print('-dpng', ['./results/' fileName])
  end
  set(gca, 'lineWidth', lineWidth);
  set(gcf, 'paperposition', origpos)
  switch nonLinearity
   case 'linear'
    xStart = min(model.comp{1}.mapt);
    xEnd = max(model.comp{1}.mapt);
    axis([xStart xEnd -2 3.5]);
  end

  for index = 1:model.comp{j}.numGenes
      figure;
      lin = plot(model.comp{j}.mapt, model.comp{j}.ypred(:,index)*scale(index), '-');
      hold on;
      bh = plot(model.comp{j}.mapt, model.comp{j}.ypred(:,index)*scale(index) + ...
                2*sqrt(model.comp{j}.ypredVar(:,index)*scale(index)*scale(index)),'--');
      bh = [bh plot(model.comp{j}.mapt, model.comp{j}.ypred(:,index)*scale(index) - ...
                2*sqrt(model.comp{j}.ypredVar(:,index)*scale(index)*scale(index)) ,'--')];
      lin = [lin plot(model.comp{j}.t, model.comp{j}.y(:,index)*scale(index), 'rx')];
      titleText = [gene{index} ' mRNA'];
      title(titleText,'fontsize', 20);

      textB = ['B = ' num2str(modelB(index))];
      textD = ['D = ' num2str(model.comp{j}.D(index))];
      textS = ['S = ' num2str(modelS(index))];
      textAlpha = ['Alpha = ' num2str(modelAlpha(index))];
      textGamma = ['Gamma = ' num2str(model.comp{j}.gParam(1,index))];
      texPos = ylim;
      scalePos = (texPos(2) - texPos(1))/14;
      texh = text(35, texPos(1)+5*scalePos, textB);
      texh = [texh text(35, texPos(1)+4*scalePos,textD)];
      texh = [texh text(35, texPos(1)+3*scalePos, textS)];
      texh = [texh text(35, texPos(1)+2*scalePos, textAlpha)];
      texh = [texh text(35, texPos(1)+1*scalePos, textGamma)];
      set(bh, 'lineWidth', 3);
      set(lin, 'lineWidth', 4);
      set(lin, 'markersize', 20);
      set(texh, 'fontsize', 16);
      set(gca, 'fontname', 'arial', 'fontsize', 24, 'xlim', [min(model.comp{j}.mapt) max(model.comp{j}.mapt)])

      if printResults
          fileName = ['dem' type num2str(expNo) kernName ...
              nonLinearity '_ExprsProfile_Rep' num2str(j) '_Gene' num2str(index)];
          print('-deps', ['./results/' fileName]);
          pos = get(gcf, 'paperposition');
          origpos = pos;
          pos(3) = pos(3);
          pos(4) = pos(4);
          set(gcf, 'paperposition', pos);
          lineWidth = get(gca, 'lineWidth');
          set(gca, 'lineWidth', lineWidth*2);
          print('-dpng', ['./results/' fileName]);
          set(gca, 'lineWidth', lineWidth);
          set(gcf, 'paperposition', origpos);
      end
  end
end

order = 1:model.comp{1}.numGenes;

counter = 0;

% Plot first basal transcription rates.
figure
bar([modelB(order)]', 0.6); colormap([0 0 0; 1 1 1]);
set(gca, 'xticklabel', gene)
if printResults
  fileName = ['dem' type num2str(expNo) kernName ...
              nonLinearity '_basal'];
  print('-deps', ['./results/' fileName]);
end
pos = get(gcf, 'paperposition');
origpos = pos;
pos(3) = pos(3)/2;
pos(4) = pos(4)/2;
set(gcf, 'paperposition', pos);
lineWidth = get(gca, 'lineWidth');
set(gca, 'lineWidth', lineWidth*2);
if printResults
  print('-dpng', ['./results/' fileName]);
end
set(gcf, 'paperposition', origpos);
set(gca, 'lineWidth', lineWidth);

% Plot the sensitivities.
figure
bar([modelS(order)]', 0.6); colormap([0 0 0; 1 1 1]);
set(gca, 'xticklabel', gene);
if printResults
  fileName = ['dem' type num2str(expNo) kernName ...
              nonLinearity '_sensitivity'];
  print('-deps', ['./results/' fileName]);
end
pos = get(gcf, 'paperposition');
origpos = pos;
pos(3) = pos(3)/2;
pos(4) = pos(4)/2;
set(gcf, 'paperposition', pos);
lineWidth = get(gca, 'lineWidth');
set(gca, 'lineWidth', lineWidth*2);
if printResults
  print('-dpng', ['./results/' fileName]);
end
set(gcf, 'paperposition', origpos)
set(gca, 'lineWidth', lineWidth);

% Finally plot degradation rates.
figure
bar([model.comp{1}.D(order)]', 0.6); colormap([0 0 0; 1 1 1]);
set(gca, 'xticklabel', gene);
if printResults
  fileName = ['dem' type num2str(expNo) kernName ...
              nonLinearity '_decay'];
  print('-deps', ['./results/' fileName]);
end
pos = get(gcf, 'paperposition');
origpos = pos;
pos(3) = pos(3)/2;
pos(4) = pos(4)/2;
set(gcf, 'paperposition', pos);
lineWidth = get(gca, 'lineWidth');
set(gca, 'lineWidth', lineWidth*2);
if printResults
  print('-dpng', ['./results/' fileName]);
end
set(gcf, 'paperposition', origpos);
set(gca, 'lineWidth', lineWidth);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function y = anyCones(K)

y = any(K.q) || any(K.e) || any(K.p) || any(K.s);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [theta,axis] = dq2rot(dq)

% DQ2ROT     extracts the rotation axis and angle of a rotation dual
%            quaternion
%
%       [THETA,AXIS] = DQ2ROT(DQ) returns the rotation angle, THETA [deg],
%       and the rotation axis, AXIS, of a rotation dual quaternion DQ.
%       - DQ is a rotation dual quaternion. It is a 8-vector or an 8*N
%       array (each column is a rotation dual quaternion) where N is the
%       number of rotation dual quaternions.
%       - THETA is the rotation angle [deg]. It is comprised between 0 and
%       180deg. THETA is a scalar or an 1*N array (element i is the
%       rotation angle of dual quaternion i).
%       - AXIS is the unitary (norm equal to 1) rotation axis. It is a
%       3-vector or an 3*N array (column i is the rotation axis of
%       dual quaternion i).
%
% See also DQ2TRANS, DQ2SCREW, ROT2DQ

sdq = size(dq);
if sdq == [1 8]
    dq = dq.';
    sdq = size(dq);
end
n = sdq(2);

% wrong format
if sdq(1) ~= 8
    error('DualQuaternion:dquat2rot:wrongsize',...
        '%d rows in the DQ array. It should be 8.',sdq(1));
end

% Extraction of the rotation parameters
theta = 2*acos(dq(1,:));
axis = sym([ones(1,n) ; zeros(2,n)]);
repsin = repmat(sin(theta(n)/2),3,1);
axis(:,n) = dq(2:4,n)./repsin;

author: Sable

repo: mcbench-benchmarks


text:
% Compile libraries and download network

[home_dir, name, ext] = fileparts(mfilename('fullpath'));
% donwload network
    mkdir model
    cd model/
    if ~(exist('imagenet-vgg-m-2048.mat', 'file') == 2)
        disp('Downloading the network "imagenet-vgg-m-2048.mat" from "http://www.vlfeat.org/matconvnet/models/imagenet-vgg-m-2048.mat"...')
        urlwrite('http://www.vlfeat.org/matconvnet/models/imagenet-vgg-m-2048.mat', 'imagenet-vgg-m-2048.mat')
        disp('Done!')
    end
    if ~(exist('imagenet-vgg-verydeep-16.mat', 'file') == 2)
        disp('Downloading the network "imagenet-vgg-verydeep-16.mat" from "http://www.vlfeat.org/matconvnet/models/imagenet-vgg-verydeep-16.mat"...')
        urlwrite('http://www.vlfeat.org/matconvnet/models/imagenet-vgg-verydeep-16.mat', 'imagenet-vgg-verydeep-16.mat')
        disp('Done!')
    end
    cd(home_dir)
% compile matconvnet
    setup_path();
    vl_setupnn();
    try
        disp('Trying to compile MatConvNet with GPU support')
        vl_compilenn('enableGpu', true, ...
               'cudaRoot', '/usr/local/cuda', ...
               'cudaMethod', 'nvcc')
    catch err
        disp('Installation failed!')
        warning('ASRCF:install', 'Could not compile MatConvNet with GPU support.\nVisit http://www.vlfeat.org/matconvnet/install/ for instructions of how to compile MatConvNet.\n');
    end
    status = movefile('mex/vl_*.mex*');
    cd(home_dir)

author: Daikenan

repo: ASRCF

save_path: github-repos/MATLAB/Daikenan-ASRCF


text:
function [sumD] = mesh_fit_elec_optim(SUMD,scalp,elec)

% mesh_fit_elec_optim - optimise the fitting of electrodes to scalp vertices
%
% This function is in development.  It needs rotations and
% translations of elec vertices to fit those of the scalp
% vertices and a minimisation function for the difference
% between the vertex locations of the nearest scalp vertices
% and those of the electrodes.

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:57 $

% Licence:  GNU GPL, no implied or express warranties
% History:  09/2002, Darren.Weber_at_radiology.ucsf.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf('Still in development\n'); return

% Rotations/translations/scaling here???

[k,d] = dsearchn(scalp,elec);

sumD = sum(d) - SUMD;

return

author: PatternRecognition

repo: OpenBMI


text:
function [model, top] = rbmBB(X, numhid, varargin)
%Learn RBM with Bernoulli hidden and visible units
%This is not meant to be applied to image data
%code by Andrej Karpathy
%based on implementation of Kevin Swersky and Ruslan Salakhutdinov

%INPUTS:
%X              ... data. should be binary, or in [0,1] to be interpreted
%               ... as probabilities
%numhid         ... number of hidden layers

%additional inputs (specified as name value pairs or in struct)
%method         ... CD or SML
%eta            ... learning rate
%momentum       ... momentum for smoothness amd to prevent overfitting
%               ... NOTE: momentum is not recommended with SML
%maxepoch       ... # of epochs: each is a full pass through train data
%avglast        ... how many epochs before maxepoch to start averaging
%               ... before. Procedure suggested for faster convergence by
%               ... Kevin Swersky in his MSc thesis
%penalty        ... weight decay factor
%batchsize      ... The number of training instances per batch
%verbose        ... For printing progress
%anneal         ... Flag. If set true, the penalty is annealed linearly
%               ... through epochs to 10% of its original value

%OUTPUTS:
%model.type     ... Type of RBM (i.e. type of its visible and hidden units)
%model.W        ... The weights of the connections
%model.b        ... The biases of the hidden layer
%model.c        ... The biases of the visible layer
%model.top      ... The activity of the top layer, to be used when training
%               ... DBN's
%errors         ... The errors in reconstruction at every epoch

%Process options
%if args are just passed through in calls they become cells
if (isstruct(varargin))
    args= prepareArgs(varargin{1});
else
    args= prepareArgs(varargin);
end
[   method        ...
    eta           ...
    momentum      ...
    maxepoch      ...
    avglast       ...
    penalty       ...
    batchsize     ...
    verbose       ...
    anneal        ...
    ] = process_options(args    , ...
    'method'        ,  'CD'     , ...
    'eta'           ,  0.1      , ...
    'momentum'      ,  0.5      , ...
    'maxepoch'      ,  50       , ...
    'avglast'       ,  5        , ...
    'penalty'       , 2e-4      , ...
    'batchsize'     , 100       , ...
    'verbose'       , false     , ...
    'anneal'        , false);
avgstart = maxepoch - avglast;
oldpenalty= penalty;
[N,d]=size(X);

if (verbose)
    fprintf('Preprocessing data...\n');
end

%Create batches
numcases=N;
numdims=d;
numbatches= ceil(N/batchsize);
groups= repmat(1:numbatches, 1, batchsize);
groups= groups(1:N);
perm=randperm(N);
groups = groups(perm);
for i=1:numbatches
    batchdata{i}= X(groups==i,:);
end

%train RBM
W = 0.1*randn(numdims,numhid);
c = zeros(1,numdims);
b = zeros(1,numhid);
ph = zeros(numcases,numhid);
nh = zeros(numcases,numhid);
phstates = zeros(numcases,numhid);
nhstates = zeros(numcases,numhid);
negdata = zeros(numcases,numdims);
negdatastates = zeros(numcases,numdims);
Winc  = zeros(numdims,numhid);
binc = zeros(1,numhid);
cinc = zeros(1,numdims);
Wavg = W;
bavg = b;
cavg = c;
t = 1;
errors=zeros(1,maxepoch);

for epoch = 1:maxepoch

	errsum=0;
    if (anneal)
        %apply linear weight penalty decay
        penalty= oldpenalty - 0.9*epoch/maxepoch*oldpenalty;
    end

    for batch = 1:numbatches
		[numcases numdims]=size(batchdata{batch});
		data = batchdata{batch};

%         %go up
% 		ph = logistic(data*W + repmat(b,numcases,1));
% 		phstates = ph > rand(numcases,numhid);
%         if (isequal(method,'SML'))
%             if (epoch == 1 && batch == 1)
%                 nhstates = phstates;
%             end
%         elseif (isequal(method,'CD'))
%             nhstates = phstates;
%         end

        %%%%%%%%% START POSITIVE PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        poshidprobs = 1./(1 + exp(-data*W - repmat(b,numcases,1)));
        posprods    = data' * poshidprobs;
        poshidact   = sum(poshidprobs);
        posvisact = sum(data);

        %%%%%%%%% END OF POSITIVE PHASE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        poshidstates = poshidprobs > rand(numcases,numhid);

%         %go down
% 		negdata = logistic(nhstates*W' + repmat(c,numcases,1));
% 		negdatastates = negdata > rand(numcases,numdims);
%
%         %go up one more time
% 		nh = logistic(negdatastates*W + repmat(b,numcases,1));
% 		nhstates = nh > rand(numcases,numhid);

        %%%%%%%%% START NEGATIVE PHASE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        negdata = 1./(1 + exp(-poshidstates*W' - repmat(c,numcases,1)));
        neghidprobs = 1./(1 + exp(-negdata*W - repmat(b,numcases,1)));
        negprods  = negdata'*neghidprobs;
        neghidact = sum(neghidprobs);
        negvisact = sum(negdata);

        %%%%%%%%% END OF NEGATIVE PHASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        err= sum(sum( (data-negdata).^2 ));
        errsum = err + errsum;

        if epoch>5,
            this_momentum=1.8*momentum;
        else
            this_momentum=momentum;
        end;

%         %update weights and biases
%         dW = (data'*ph - negdatastates'*nh);
%         dc = sum(data) - sum(negdatastates);
%         db = sum(ph) - sum(nh);
% 		Winc = momentum*Winc + eta*(dW/numcases - penalty*W);
% 		binc = momentum*binc + eta*(db/numcases);
% 		cinc = momentum*cinc + eta*(dc/numcases);
% 		W = W + Winc;
% 		b = b + binc;
% 		c = c + cinc;

        %%%%%%%%% UPDATE WEIGHTS AND BIASES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Winc = this_momentum*Winc + eta*( (posprods-negprods)/numcases - penalty*W);
        binc = this_momentum*binc + (eta/numcases)*(poshidact-neghidact);
        cinc = this_momentum*cinc + (eta/numcases)*(posvisact-negvisact);

        W = W + Winc;
        c = c + cinc;
        b = b + binc;

        %%%%%%%%%%%%%%%% END OF UPDATES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%         if (epoch > avgstart)
%             %apply averaging
% 			Wavg = Wavg - (1/t)*(Wavg - W);
% 			cavg = cavg - (1/t)*(cavg - c);
% 			bavg = bavg - (1/t)*(bavg - b);
% 			t = t+1;
% 		else
			Wavg = W;
			bavg = b;
			cavg = c;
%         end

        %accumulate reconstruction error
        err= sum(sum( (data-negdata).^2 ));
		errsum = err + errsum;

    end

    errors(epoch)=errsum;
    if (verbose)
        fprintf('Ended epoch %i/%i. Reconstruction error is %f\n', ...
            epoch, maxepoch, errsum);
    end
end

model.type= 'BB';
top = logistic(X*Wavg + repmat(bavg,N,1));
model.W= Wavg;
model.b= bavg;
model.c= cavg;

author: marianux

repo: ecg-kit

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


text:
function r = multirnd(theta,k)
%MULTIRND - Random vector from multinomial distribution.
%   r = multirnd(theta,k) returns a vector randomly selected
%   from the multinomial distribution with parameter vector
%   theta, and count k (i.e. sum(r) = k).
%
%   Note: if k is unspecified, then it is assumed k=1.
%
%   Author: David Ross
%

%--------------------------------------------------------
% Check the arguments.
%--------------------------------------------------------
error(nargchk(1,2,nargin));

% make sure theta is a vector
if ndims(theta) > 2 | all(size(theta) > 1)
    error('theta must be a vector');
end

% if theta is a row vector, convert it to a column vector
if size(theta,1) == 1
    theta = theta';
end

% make sure k is a scalar?

% if the number of samples has not been provided, set
% it to one
if nargin == 1
    k = 1;
end

%--------------------------------------------------------
% Main...
%--------------------------------------------------------
n = length(theta);
theta_cdf = cumsum(theta);

r = zeros(n,1);
random_vals = rand(k,1);

for j = 1:k
    index = min(find(random_vals(j) <= theta_cdf));
    r(index) = r(index) + 1;
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function w = dualtree(x, J, Faf, af)

% Dual-tree Complex Discrete Wavelet Transform
%
% USAGE:
%    w = dualtree(x, J, Faf, af)
% INPUT:
%   x - N-point vector
%       1) N is divisible by 2^J
%       2) N >= 2^(J-1)*length(af)
%   J - number of stages
%   Faf - filters for the first stage
%   af - filters for the remaining stages
% OUTPUT:
%   w - DWT coefficients
%      w{j}{1}, j = 1..J - real part
%      w{j}{2}, j = 1..J - imaginary part
%      w{J+1}{d} - lowpass coefficients, d = 1,2
% EXAMPLE:
%    x = rand(1, 512);
%    J = 4;
%    [Faf, Fsf] = FSfarras;
%    [af, sf] = dualfilt1;
%    w = dualtree(x, J, Faf, af);
%    y = idualtree(w, J, Fsf, sf);
%    err = x - y;
%    max(abs(err))
%
% WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
% http://taco.poly.edu/WaveletSoftware/

% normalization
x = x/sqrt(2);

% Tree 1
[x1 w{1}{1}] = afb(x, Faf{1});
for j = 2:J
    [x1 w{j}{1}] = afb(x1, af{1});
end
w{J+1}{1} = x1;

% Tree 2
[x2 w{1}{2}] = afb(x, Faf{2});
for j = 2:J
    [x2 w{j}{2}] = afb(x2, af{2});
end
w{J+1}{2} = x2;

author: thomaskuestner


text:
function h=m_rectangle(long,lat,width,height,varargin)
%M_RECTANGLE Create rectangle, rounded-rectangle, or ellipse
%   M_RECTANGLE adds a default rectangle to the current map axes.
%
%   M_RECTANGLE(long,lat,width,height) creates a rectangle in 2-D coordinates.
%   The long and lat elements determine the location and the width and height
%   elements determine the size. The function plots into the current axes
%   without clearing existing content from the axes.
%   M_RECTANGLE(long,lat,width,height,flag) creates different rectangle. If
%   flag=0(default),it is a normal rectangle. If flag=1, it is a curve box
%   determined to the projection.
%
%   M_RECTANGLE(...,'Curvature',cur) adds curvature to the sides
%   of the rectangle. For different curvatures along the horizontal and
%   vertical sides, specify cur as a two-element vector of the form
%   [horizontal vertical]. For the same length of curvature along all
%   sides, specify cur as a scalar value. Specify values between 0 (no
%   curvature) and 1 (maximum curvature). Use [1 1] to create an ellipse or
%   circle.
%
%   M_RECTANGLE(...,Name,Value) specifies rectangle properties using one or
%   more Name,Value pair arguments. It is the same as function rectangle.
%
%   M_RECTANGLE(container,...) creates the rectangle in the axes, group, or
%   transform specified by container, instead of in the current axes.
%
%   R = M_RECTANGLE(...) returns the rectangle object created.
%
%   Execute GET(R), where R is a rectangle object, to see a list of
%   rectangle object properties and their current values.
%   Execute SET(R) to see a list of rectangle object properties and legal
%   property values.
%

% Shi Weiheng (tfoterye@gmail.com) 25/Mar/18
%
% Based on m_quiver written by Prof Rich Pawlowicz.
% Comments mainly come from function rectangle.
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.
%

global MAP_PROJECTION MAP_VAR_LIST

% Have to have initialized a map first

if isempty(MAP_PROJECTION)
  disp('No Map Projection initialized - call M_PROJ first!');
  return;
end

if (varargin{1}==0 || varargin{1}==1)
    flag=varargin{1};
    varargin1={varargin{2:end}};
else
    flag=0;
    varargin1=varargin;
end

[X,Y]=m_ll2xy(long,lat,'clip','off');
[X1,Y1]=m_ll2xy(long+width,lat+height,'clip','off');
[X2,Y2]=m_ll2xy(long,lat+height,'clip','off');
[X3,Y3]=m_ll2xy(long+width,lat,'clip','off');

[small_x1,small_y1]=m_ll2xy(long:width/100:long+width,lat*ones(1,101),'clip','off');
[small_x2,small_y2]=m_ll2xy(long:width/100:long+width,(lat+height)*ones(1,101),'clip','off');

if flag==0
h=rectangle('Position',[X,Y,X1-X,Y1-Y],varargin1{:});
end
if flag==1
   h=line( [X,small_x1,X1,small_x2(end:-1:1),X],[Y,small_y1,Y1,small_y2(end:-1:1),Y],varargin1{:});
end
set(h,'tag','m_rectangle');
if nargout==0
 clear h
end

author: goGPS-Project


text:
function ccvt_reflect ( )

%*****************************************************************************80
%
%% CCVT_REFLECT generates and animates a CVT dataset.
%
%  Discussion:
%
%    This program is meant to be used interactively.  It's also
%    possible to prepare a simple input file beforehand and use it
%    in batch mode.
%
%    The program requests input values from the user:
%
%    * NDIM, the spatial dimension,
%    * N, the number of points to generate,
%    * SEED, a seed to use for random number generation;
%    * INIT, initialize the points:
%      ** file, by reading data from file;
%      ** 'GRID', picking points from a grid;
%      ** 'HALTON', from a Halton sequence;
%      ** 'RAND', using MATLAB's RAND function;
%      ** 'UNIFORM', using a simple uniform RNG;
%    * IT_MAX, the maximum number of iterations;
%    * IT_FIXED, the number of iterative steps to take
%      using a fixed set of sampling points.
%    * SAMPLE, how to conduct the sampling:
%      ** 'GRID', picking points from a grid;
%      ** 'HALTON', from a Halton sequence;
%      ** 'RAND', using MATLAB's RAND function;
%      ** 'UNIFORM', using a simple uniform RNG;
%    * SAMPLE_NUM, the number of sampling points;
%    * BATCH, the number of sampling points to create at one time;
%    * OUTPUT, a file into which the data can be stored.
%
%    To indicate that no further computations are desired, it is
%    enough to input a nonsensical value, such as -1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 December 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Qiang Du, Vance Faber, and Max Gunzburger,
%    Centroidal Voronoi Tessellations: Applications and Algorithms,
%    SIAM Review, Volume 41, 1999, pages 637-676.
%
  DEBUG = 1;

  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CCVT_REFLECT\n' );
  fprintf ( 1,   MATLAB version\n' );
  fprintf ( 1, '  Generate and animate a CVT dataset.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  This program is meant to be used interactively.\n' );
  fprintf ( 1, '  It is also possible to prepare a simple input \n' );
  fprintf ( 1, '  file beforehand and use it in batch mode.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The program requests input values from the user:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  * NDIM, the spatial dimension,\n' );
  fprintf ( 1, '  * N, the number of points to generate,\n' );
  fprintf ( 1, '  * SEED, a seed to use for random number generation,\n' );
  fprintf ( 1, '  * INIT, initialize the points:\n' );
  fprintf ( 1, '    ** file, read data from a file;\n' );
  fprintf ( 1, '    ** ''GRID'', by picking points from a grid;\n' );
  fprintf ( 1, '    ** ''HALTON'', from a Halton sequence;\n' );
  fprintf ( 1, '    ** ''RAND'', using MATLAB''s RAND function;\n' );
  fprintf ( 1, '    ** ''UNIFORM'', using a simple uniform RNG;\n' );
  fprintf ( 1, '  * IT_MAX, the maximum number of iterations.\n' );
  fprintf ( 1, '  * IT_FIXED, the number of iterative steps to take\n' );
  fprintf ( 1, '    using a fixed set of sampling points.\n' );
  fprintf ( 1, '  * SAMPLE, how to conduct the sampling.\n' );
  fprintf ( 1, '    ** ''GRID'', by picking points from a grid;\n' );
  fprintf ( 1, '    ** ''HALTON'', from a Halton sequence;\n' );
  fprintf ( 1, '    ** ''RAND'', using MATLAB''s RAND function;\n' );
  fprintf ( 1, '    ** ''UNIFORM'', using a simple uniform RNG;\n' );
  fprintf ( 1, '  * SAMPLE_NUM, the number of sample points;\n' );
  fprintf ( 1, '  * BATCH, the number of sampling points to create at one time;\n' );
  fprintf ( 1, '  * OUTPUT, a file into which the data is stored.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  To indicate that no further computations are \n' );
  fprintf ( 1, '  desired, it is enough to input a nonsensical value, \n' );
  fprintf ( 1, '  such as -1.\n' );

  while ( 1 )

    fprintf ( 1, '  *\n' );
    fprintf ( 1, ' *\n' );
    fprintf ( 1, '*  Ready to generate a new dataset:\n' );
    fprintf ( 1, ' *\n' );
    fprintf ( 1, '  *\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  NDIM is the spatial dimension.\n' );
    fprintf ( 1, '  (Try ''2'' if you have no preference.)\n' );
    fprintf ( 1, '  (Any value less than 1 terminates execution.)\n' );
    ndim = [];
    ndim = input ( '  Enter NDIM:  ' );

    fprintf ( 1, '  NDIM = %12d\n', ndim );

    if ( ndim < 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_REFLECT\n' );
      fprintf ( 1, '  The input value of NDIM = %12d\n', ndim );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  N is the number of points to generate.\n' );
    fprintf ( 1, '  (Try ''25'' if you have no preference.)\n' );
    fprintf ( 1, '  (Any value less than 1 terminates execution.)\n' );
    n = [];
    n = input ( '  Enter N:  ' );

    fprintf ( 1, '  User input N = %12d\n', n );

    if ( n < 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_REFLECT\n' );
      fprintf ( 1, '  The input value of N = %12d\n', n );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  NPP is the number of boundary sample points.\n' );
    fprintf ( 1, '  (Try ''1000'' if you have no preference.)\n' );
    fprintf ( 1, '  (Any value less than 1 terminates execution.)\n' );
    npp = [];
    npp = input ( '  Enter NPP:  ' );

    fprintf ( 1, '  User input NPP = %12d\n', npp );

    if ( npp < 1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_REFLECT\n' );
      fprintf ( 1, '  The input value of NPP = %12d\n', npp );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  SEED is a seed for the random number generation.\n' );
    fprintf ( 1, '  (Try ''123456789'' if you have no preference.)\n' );
    fprintf ( 1, '  (Any value less than 0 terminates execution.)\n' );
    seed = [];
    seed = input ( '  Enter SEED:  ' );

    fprintf ( 1, '  User input SEED = %d\n', seed );

    if ( seed < 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_REFLECT\n' );
      fprintf ( 1, '  The input value of SEED = %12d\n', seed );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    seed_init = seed;

    fprintf ( 1, '\n' );
    fprintf ( 1, '  INIT is the method of initializing the data:\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  file       read data from a file;\n' );
    fprintf ( 1, '  ''GRID''     by picking points from a grid;\n' );
    fprintf ( 1, '  ''HALTON''   from a Halton sequence;\n' );
    fprintf ( 1, '  ''RAND''     using MATLAB''s RAND function;\n' );
    fprintf ( 1, '  ''UNIFORM''  using a simple uniform RNG;\n' );
    fprintf ( 1, ' \n' );
    fprintf ( 1, '  (Try ''RAND'' if you have no preference.)\n' );
    fprintf ( 1, '  (A blank value terminates execution).\n' );
    fprintf ( 1, '  (Be sure to INCLUDE QUOTES around the string!\n' );
    fprintf ( 1, ' \n' );

    init_string = [];
    init_string = input ( '  Enter INIT:  ' );

    fprintf ( 1, '  User input INIT = "%s".\n', init_string );

    if ( s_len_trim ( init_string ) <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_REFLECT\n' );
      fprintf ( 1, '  The input value of INIT \n' );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    if ( s_eqi ( init_string, 'RAND'  ) )
      init = -1;
    elseif ( s_eqi ( init_string, 'RANDOM' ) )
      init_string = 'RAND';
      init = -1;
    elseif ( s_eqi ( init_string, 'UNIFORM' ) )
      init = 0;
    elseif ( s_eqi ( init_string, 'HALTON'  ) )
      init = 1;
    elseif ( s_eqi ( init_string, 'GRID'    ) )
      init = 2;
    elseif ( 0 < s_len_trim ( init_string ) )
      init = 3;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_REFLECT\n' );
      fprintf ( 1, '  The input value of INIT \n' );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  IT_MAX is the maximum number of iterations.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  An iteration carries out the following steps:\n' );
    fprintf ( 1, '  * the Voronoi region associated with each\n' );
    fprintf ( 1, '    generator is estimated by sampling;\n' );
    fprintf ( 1, '  * the centroid of each Voronoi region is estimated.\n' );
    fprintf ( 1, '  * the generator is replaced by the centroid.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  If "enough" sampling points are used,\n' );
    fprintf ( 1, '  and "enough" iterations are taken, this process\n' );
    fprintf ( 1, '  will converge.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (Try ''50'' if you have no preference.)\n' );
    fprintf ( 1, '  (A negative value terminates execution).\n' );
    fprintf ( 1, '\n' );
    it_max = [];
    it_max = input ( '  Enter IT_MAX:  ' );

    fprintf ( 1, '  User input IT_MAX = %12d\n', it_max );

    if ( it_max < 0 )
      fprintf ( 1, ' \n' );
      fprintf ( 1, 'CCVT_REFLECT\n' );
      fprintf ( 1, '  The input value of IT_MAX = %12d\n', it_max );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  IT_FIXED is the number of consecutive iterations\n' );
    fprintf ( 1, '  to take with a fixed set of sample points.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Setting IT_FIXED to 1 means a new set of sample\n' );
    fprintf ( 1, '  points is generated on every iterative step;\n' );
    fprintf ( 1, '  Setting IT_FIXED equal to IT_MAX means a single set\n' );
    fprintf ( 1, '  of sample points is used for the entire iteration.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Any value between 1 and IT_MAX is reasonable.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (Try "%d" if you do not have a preference).\n', it_max );
    fprintf ( 1, '  (A 0 or negative value terminates execution).\n' );
    fprintf ( 1, '\n' );
    it_fixed = [];
    it_fixed = input ( '  Enter IT_FIXED:  ' );

    fprintf ( 1, '  User input IT_FIXED = %12d\n', it_fixed );

    if ( it_max < 0 )
      fprintf ( 1, ' \n' );
      fprintf ( 1, 'CCVT_REFLECT\n' );
      fprintf ( 1, '  The input value of IT_FIXED = %12d\n', it_fixed );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  SAMPLE is the method of sampling the region:\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  ''GRID''     by picking points from a grid;\n' );
    fprintf ( 1, '  ''HALTON''   from a Halton sequence;\n' );
    fprintf ( 1, '  ''RAND''     using MATLAB''s RAND function;\n' );
    fprintf ( 1, '  ''UNIFORM''  using a simple uniform RNG;\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (Try ''RAND'' if you have no preference.)\n' );
    fprintf ( 1, '  (A blank value terminates execution).\n' );
    fprintf ( 1, '  (Be sure to INCLUDE QUOTES around the string!\n' );
    fprintf ( 1, '\n' );

    sample_string = [];
    sample_string = input ( '  Enter SAMPLE:  ' );

    fprintf ( 1, '  User input SAMPLE = "%s".\n', sample_string );

    if ( s_len_trim ( sample_string ) <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_REFLECT\n' );
      fprintf ( 1, '  The input value of SAMPLE \n' );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    if ( s_eqi ( sample_string, 'RAND'  ) )
      sample = -1;
    elseif ( s_eqi ( sample_string, 'RANDOM' ) )
      sample = -1;
      sample_string = 'RAND';
    elseif ( s_eqi ( sample_string, 'UNIFORM' ) )
      sample = 0;
    elseif ( s_eqi ( sample_string, 'HALTON'  ) )
      sample = 1;
    elseif ( s_eqi ( sample_string, 'GRID'    ) )
      sample = 2;
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_REFLECT\n' );
      fprintf ( 1, '  The input value of SAMPLE \n' );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  SAMPLE_NUM is the number of sample points.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  The Voronoi regions will be explored by generating\n' );
    fprintf ( 1, '  SAMPLE_NUM points.  For each sample point, the\n' );
    fprintf ( 1, '  nearest generator is found.  Using more points\n' );
    fprintf ( 1, '  gives a better estimate of these regions.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  SAMPLE_NUM should be much larger than N, the\n' );
    fprintf ( 1, '  number of generators.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (Try ''10000'' if you have no preference.)\n' );
    fprintf ( 1, '  (A zero or negative value terminates execution.)\n' );
    fprintf ( 1, '\n' );

    sample_num = [];
    sample_num = input ( '  Enter SAMPLE_NUM:  ' );

    fprintf ( 1, '  User input SAMPLE_NUM = %12d\n', sample_num );

    if ( sample_num <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_REFLECT\n' );
      fprintf ( 1, '  The input value of SAMPLE_NUM = %12d\n', sample_num );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  BATCH is the number of sample points to create\n' );
    fprintf ( 1, '  at one time\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  BATCH should be between 1 and SAMPLE_NUM.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  It is FASTER to set BATCH to SAMPLE_NUM;\n' );
    fprintf ( 1, '  setting BATCH to 1 requires the least memory.\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (Try ''%d'' if you have no preference.)\n', ...
      min ( 1000, sample_num ) );
    fprintf ( 1, '  (A zero or negative value terminates execution.)\n' );
    fprintf ( 1, '\n' );

    batch = [];
    batch = input ( '  Enter BATCH:  ' );

    fprintf ( 1, '  User input BATCH = %12d\n', batch );

    if ( batch <= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'CCVT_REFLECT\n' );
      fprintf ( 1, '  The input value of BATCH = %12d\n', batch );
      fprintf ( 1, '  is interpreted as a request for termination.\n' );
      fprintf ( 1, '  Normal end of execution.\n' );
      break
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  MOVIE_NAME is the name of the file in which the movie is stored;\n' );
    fprintf ( 1, '\n' );
    fprintf ( 1, '  (Try ''movie.avi'' if you have no preference.)\n' );
    fprintf ( 1, '  (A blank value terminates execution).\n' );
    fprintf ( 1, '  (Be sure to INCLUDE QUOTES around the string!\n' );
    fprintf ( 1, ' \n' );

    movie_name = [];
    movie_name = input ( '  Enter MOVIE_NAME:  ' );

    fprintf ( 1, '  User input MOVIE_NAME = %s\n', movie_name );
%
%  Initialize some movie information.
%
    if ( file_exist ( movie_name ) )
      file_delete ( movie_name )
    end

    num_frames_per_second = 10;
    aviobj = avifile ( movie_name, 'fps', num_frames_per_second );
%
%  Print the header for output.
%
    if ( DEBUG )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Step    SEED          L2-Change     Energy\n' );
      fprintf ( 1, '\n' );
    end
%
%  Initialize the generators.
%
    if ( init == 3 )

      r = data_read ( sample_string, ndim, n );

    else

      seed = seed_init;
      initialize = 1;

      [ r, seed ] = cvt_sample ( ndim, n, n, init, initialize, seed );

    end
%
%  If the initialization and sampling steps use the same random number
%  scheme, then the sampling scheme does not have to be initialized.
%
    if ( init == sample )
      initialize = 0;
    else
      initialize = 1;
    end
%
%  Compute the initial energy.
%
    seed_init2 = seed_init;

    [ energy, seed ] = cvt_energy ( ndim, n, batch, sample, initialize, ...
      sample_num, seed, r );

    initialize = 0;

    it_num = 0;

    if ( DEBUG )
      fprintf ( 1, '  %4d  %12d  %14e\n', ...
        it_num, seed_init2, energy );
    end
%
%  Write out the initial points.
%
    cvt_write ( ndim, n, batch, seed_init, seed, init_string, it_max, ...
      it_fixed, it_num, energy, sample_string, sample_num, ...
      r, 'initial.txt' );

    points_eps ( 'initial.eps', ndim, n, r, 'Initial points' );

    while ( it_num < it_max )
%
%  Once every IT_FIXED steps, update the base value of SEED,
%  either to SEED_INIT before the first call to CVT_ITERATE, or to the
%  output value of SEED from the previous call to CVT_ITERATE.
%
%  Otherwise, reset the value of SEED to SEED_BASE.
%
      if ( mod ( it_num, it_fixed ) == 0 )
        seed_base = seed;
      else
        seed = seed_base;
      end

      it_num = it_num + 1;
%
%  Sample the region.
%
      seed_init2 = seed;

      r2(1:ndim,1:n) = 0.0;
      count(1:n) = 0;

      outside = 0;

      [ s, seed ] = cvt_sample ( ndim, sample_num, sample_num, sample, ...
        initialize, seed );

      initialize = 0;
      energy = 0.0;

      for js = 1 : sample_num

        distance = Inf;
        nearest = -1;

        for jr = 1 : n

          dist_sq = 0.0E+00;
          for i = 1 : ndim
            dist_sq = dist_sq + ( r(i,jr) - s(i,js) )^2;
          end

          if ( dist_sq < distance )
            distance = dist_sq;
           nearest = jr;
          end

        end

        s2(1:ndim,1) = 2.0 * r(1:ndim,nearest) - s(1:ndim,js);

        if ( ...
          0.0 <= s2(1,1) & s2(1,1) <= 1.0 & ...
          0.0 <= s2(2,1) & s2(2,1) <= 1.0 )
          r2(1:ndim,nearest) = r2(1:ndim,nearest) + s(1:ndim,js);
          count(nearest) = count(nearest) + 1;
        else
          r2(1:ndim,nearest) = r2(1:ndim,nearest) + r(1:ndim,nearest);
          count(nearest) = count(nearest) + 1;
        end

        energy = energy + sum ( ( r(1:ndim,nearest) - s(1:ndim,js) ).^2 );

      end

      energy = energy / sample_num;

      for j = 1 : n
        r2(1:ndim,j) = r2(1:ndim,j) / count(j);
      end

      r(1:ndim,1:n) = r2(1:ndim,1:n);

      handle = scatter ( r(1,1:n), r(2,1:n), [],'r', 'filled' );
      line ( [ 0.0, 1.0, 1.0, 0.0, 0.0 ], [ 0.0, 0.0, 1.0, 1.0, 0.0 ] );
      axis ( [ -0.05, 1.05, -0.05, 1.05 ] );
      axis square
%
%  Label the axes and the plot.
%
      xlabel ( 'X', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
        'FontSize', 16 );

      ylabel ( 'Y', 'FontName', 'Helvetica', 'FontWeight', 'bold', ...
        'FontSize', 16, 'Rotation', 0 );

      it_string = sprintf ( 'Step %4d', it_num );

      title ( it_string, 'FontName', 'Helvetica', 'FontWeight', ...
        'bold', 'FontSize', 16 );

      frame = getframe ( gca );
      aviobj = addframe ( aviobj, frame );

      initialize = 0;
      if ( DEBUG )
        fprintf ( 1, '  %4d  %12d  %14e\n', ...
          it_num, seed_init2, energy );
      end

    end

    aviobj = close ( aviobj );
%
%  Write the final data to a file.
%
    cvt_write ( ndim, n, batch, seed_init, seed, init_string, it_max, ...
      it_fixed, it_num, energy, sample_string, sample_num, ...
      r, 'final.txt' );

    points_eps ( 'final.eps', ndim, n, r, 'Final points' );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Press RETURN to proceed.\n' );

    pause

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Final value of SEED = %d\n', seed );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [a] = sortrow(a)
% a=ceil(rand(10,4)*3);
a=sortrows(a,1);
tmp=unique(a(:,1));
tmp_count=histc(a(:,1),tmp);
start_box=0;
end_box=0;
i=1;
startloc=1;
endloc=0;
while(i<length(tmp_count)+1)
    endloc=endloc+tmp_count(i);
    if tmp_count(i)>1
      a(startloc:endloc,:) = sortrows(a(startloc:endloc,:),2);
      start_box=[start_box startloc];
      end_box=[end_box endloc];
    end
    startloc=startloc+tmp_count(i);
    i=i+1;
end
for i=2:length(start_box)
    tmp=unique(a(start_box(i):end_box(i),2));
    tmp_count=histc(a(start_box(i):end_box(i),2),tmp);
    ii=1;
    startloc=start_box(i);
    endloc=start_box(i)-1;
    while(ii<length(tmp_count)+1)
    endloc=endloc+tmp_count(ii);
    if tmp_count(ii)>1
      a(startloc:endloc,:) = sortrows(a(startloc:endloc,:),3);
    end
    startloc=startloc+tmp_count(ii);
    ii=ii+1;
    end
end

author: Ewenwan

repo: Mathematics

save_path: github-repos/MATLAB/Ewenwan-Mathematics


text:
% Main layout-manager class
classdef GridLayout < handle

    properties
        Container

        NumRows
        NumCols
        RowHeight
        ColWidth
        HGap
        VGap
        LMargin
        RMargin
        TMargin
        BMargin

        % Sizing policies
        % 0: Absolute
        % 1: Proportional
        % 2: Automatic

        RowHeightPolicy
        RowHeightAbsolute
        RowHeightWeight

        ColWidthPolicy
        ColWidthAbsolute
        ColWidthWeight
    end

    properties (Constant)
        MinRowHeight = 5
        MinColWidth = 5
    end

    properties (Dependent)
        ActualWidth
        ActualHeight
        RequiredWidth
        RequiredHeight
    end

    properties (SetAccess = private)
        Cell  % Cell container
    end

    methods
        function Obj = GridLayout(Parent, varargin)
            Names = { ...
                'NumRows' ...
                'NumCols' ...
                'RowHeight' ...
                'ColWidth' ...
                'CellColor' ...
                'Gap' ...
                'HGap' ...
                'VGap' ...
                'Margin' ...
                'LMargin' ...
                'RMargin' ...
                'TMargin' ...
                'BMargin'};
            NamesCell = CellProps.GetCellParamNames();
            [P, PC] = ParseArgs(varargin, Names, NamesCell);

            assert(~isempty(Parent), ...
                'Parent handle cannot be empty.');

            Obj.Container = uicontainer( ...
                'Parent', Parent, ...
                'ResizeFcn', @(hsrc,ev)Update(Obj,true), ...
                'BackgroundColor', P.CellColor);
            Position = getpixelposition(Obj.Container);
            setappdata(Obj.Container,'Position',Position);

            % NumRows
            DefaultNumRows = 1;
            if ~isempty(P.RowHeight) && iscell(P.RowHeight)
                DefaultNumRows = length(P.RowHeight);
            end
            Obj.NumRows = GetArg(P,'NumRows',DefaultNumRows);
            % NumCols
            DefaultNumCols = 1;
            if ~isempty(P.ColWidth) && iscell(P.ColWidth)
                DefaultNumCols = length(P.ColWidth);
            end
            Obj.NumCols = GetArg(P,'NumCols',DefaultNumCols);
            % RowHeight
            Value = GetArg(P,'RowHeight','*');
            if ~iscell(Value)
                Obj.RowHeight = cell(1,Obj.NumRows);
                [Obj.RowHeight{1:end}] = deal(Value);
            else
                Obj.RowHeight = P.RowHeight;
            end
            % ColWidth
            Value = GetArg(P,'ColWidth','*');
            if ~iscell(Value)
                Obj.ColWidth = cell(1,Obj.NumCols);
                [Obj.ColWidth{1:end}] = deal(Value);
            else
                Obj.ColWidth = P.ColWidth;
            end
            % HGap
            Obj.HGap = GetArg(P,'HGap',5);
            % VGap
            Obj.VGap = GetArg(P,'VGap',5);
            % Gap (overrides HGap, VGap)
            if ~isempty(P.Gap)
                Obj.HGap = P.Gap;
                Obj.VGap = P.Gap;
            end
            % Margin
            if isempty(P.Margin)
                DefaultLMargin = Obj.HGap;
                DefaultRMargin = Obj.HGap;
                DefaultTMargin = Obj.VGap;
                DefaultBMargin = Obj.VGap;
            else
                DefaultLMargin = P.Margin;
                DefaultRMargin = P.Margin;
                DefaultTMargin = P.Margin;
                DefaultBMargin = P.Margin;
            end
            Obj.LMargin = GetArg(P,'LMargin',DefaultLMargin);
            Obj.RMargin = GetArg(P,'RMargin',DefaultRMargin);
            Obj.TMargin = GetArg(P,'TMargin',DefaultTMargin);
            Obj.BMargin = GetArg(P,'BMargin',DefaultBMargin);

            % Create array of cell containers
            Obj.Cell = zeros(Obj.NumRows,Obj.NumCols);
            for RIdx = 1:Obj.NumRows
                for CIdx = 1:Obj.NumCols
                    ThisCell = uicontainer( ...
                        'Parent', Obj.Container, ...
                        'Units', 'pixels');
                    if ~isempty(P.CellColor)
                        set(ThisCell, 'BackgroundColor', P.CellColor);
                    end
                    Props = CellProps(PC);
                    setappdata(ThisCell,'Props',Props);
                    setappdata(ThisCell,'RSpan',[RIdx RIdx]);
                    setappdata(ThisCell,'CSpan',[CIdx CIdx]);
                    Obj.Cell(RIdx,CIdx) = ThisCell;
                end
            end

            % RowHeightPolicy
            Obj.RowHeightPolicy   = zeros(1,Obj.NumRows);
            Obj.RowHeightAbsolute = zeros(1,Obj.NumRows);
            Obj.RowHeightWeight   = zeros(1,Obj.NumRows);
            for RIdx = 1:Obj.NumRows
                RHeight = Obj.RowHeight{RIdx};
                if isnumeric(RHeight)
                    % Absolute
                    Obj.RowHeightPolicy(RIdx) = 0;
                    Obj.RowHeightAbsolute(RIdx) = max(RHeight,Obj.MinRowHeight);
                elseif ischar(RHeight) && RHeight(end) == '*'
                    % Proportional
                    Obj.RowHeightPolicy(RIdx) = 1;
                    Weight = RHeight(1:end-1);
                    if isempty(Weight)
                        Weight = '1';
                    end
                    Obj.RowHeightWeight(RIdx) = str2double(Weight);
                elseif ischar(RHeight) && isequal(RHeight,'a')
                    % Automatic
                    Obj.RowHeightPolicy(RIdx) = 2;
                else
                    error('Illegal RowHeight definition.');
                end
            end

            % ColWidthPolicy
            Obj.ColWidthPolicy   = zeros(1,Obj.NumCols);
            Obj.ColWidthAbsolute = zeros(1,Obj.NumCols);
            Obj.ColWidthWeight   = zeros(1,Obj.NumCols);
            for CIdx = 1:Obj.NumCols
                CWidth = Obj.ColWidth{CIdx};
                if isnumeric(CWidth)
                    % Absolute
                    Obj.ColWidthPolicy(CIdx) = 0;
                    Obj.ColWidthAbsolute(CIdx) = max(CWidth,Obj.MinColWidth);
                elseif ischar(CWidth) && CWidth(end) == '*'
                    % Proportional
                    Obj.ColWidthPolicy(CIdx) = 1;
                    Weight = CWidth(1:end-1);
                    if isempty(Weight)
                        Weight = '1';
                    end
                    Obj.ColWidthWeight(CIdx) = str2double(Weight);
                elseif ischar(CWidth) && isequal(CWidth,'a')
                    % Automatic
                    Obj.ColWidthPolicy(CIdx) = 2;
                else
                    error('Illegal ColWidth definition.');
                end
            end

            if ~isempty(Obj.RequiredHeight)
                HeightDiff = Obj.RequiredHeight-Obj.ActualHeight;
                Position = getpixelposition(Parent);
                Position(4) = Position(4) + HeightDiff;
                setpixelposition(Parent, Position);
            end
        end

        function MergeCells(Obj, RSpan, CSpan)
            if isempty(RSpan)
                RSpan = [1 Obj.NumRows];
            end
            if isscalar(RSpan)
                RSpan = [RSpan RSpan];
            end
            if isempty(CSpan)
                CSpan = [1 Obj.NumCols];
            end
            if isscalar(CSpan)
                CSpan = [CSpan CSpan];
            end
            for RIdx = RSpan(1):RSpan(2)
                for CIdx = CSpan(1):CSpan(2)
                    ThisCell = Obj.Cell(RIdx,CIdx);
                    if RIdx == RSpan(1) && CIdx == CSpan(1)
                        setappdata(ThisCell,'RSpan',RSpan);
                        setappdata(ThisCell,'CSpan',CSpan);
                    else
                        delete(ThisCell);
                    end
                end
            end
        end

        function RemoveCells(Obj, RSpan, CSpan)
            if isempty(RSpan)
                RSpan = [1 Obj.NumRows];
            end
            if isscalar(RSpan)
                RSpan = [RSpan RSpan];
            end
            if isempty(CSpan)
                CSpan = [1 Obj.NumCols];
            end
            if isscalar(CSpan)
                CSpan = [CSpan CSpan];
            end
            for RIdx = RSpan(1):RSpan(2)
                for CIdx = CSpan(1):CSpan(2)
                    ThisCell = Obj.Cell(RIdx,CIdx);
                    if ishghandle(ThisCell)
                        delete(ThisCell);
                    end
                end
            end
        end

        function FormatCells(Obj, RSpan, CSpan, varargin)
            if isscalar(RSpan)
                RSpan = [RSpan RSpan];
            end
            if isscalar(CSpan)
                CSpan = [CSpan CSpan];
            end
            for RIdx = RSpan(1):RSpan(2)
                for CIdx = CSpan(1):CSpan(2)
                    ThisCell = Obj.Cell(RIdx,CIdx);
                    if ~ishghandle(ThisCell)
                        continue; % Cell is removed
                    end
                    GridLayout.FormatCell(ThisCell, varargin{:});
                end
            end
        end

        function Update(Obj, IsCallback)
            if nargin < 2
                IsCallback = false;
            end
            % Layout size
            NewPosition = getpixelposition(Obj.Container);
            OldPosition = getappdata(Obj.Container,'Position');
            if isequal(NewPosition,OldPosition) && IsCallback
                return;
            end
            setappdata(Obj.Container,'Position',NewPosition);
            LayoutWidth  = NewPosition(3);
            LayoutHeight = NewPosition(4);

            % Get children cells; as a vector!
            Cells = get(Obj.Container,'Children');
            % Get their row and column spans
            RSpan = zeros(length(Cells),2);
            CSpan = zeros(length(Cells),2);
            for i = 1:length(Cells)
                RSpan(i,:) = getappdata(Cells(i),'RSpan');
                CSpan(i,:) = getappdata(Cells(i),'CSpan');
            end

            % Determine column widths
            CWidthSum = LayoutWidth-(Obj.NumCols-1)*Obj.VGap-Obj.LMargin-Obj.RMargin;
            CWidth = zeros(1,Obj.NumCols);
            % 1. Absolute
            CMask = Obj.ColWidthPolicy == 0;
            CWidth(CMask) = Obj.ColWidthAbsolute(CMask);
            % 2. Automatic
            for CIdx = find(Obj.ColWidthPolicy == 2)
                MaxWidth = 0;
                Mask = CSpan(:,1) == CIdx & CSpan(:,2) == CIdx;
                assert(sum(Mask) > 0, ...
                    'Any column with automatic width must have at least one non-merged cell.');
                for MyCell = Cells(Mask)'
                    Props = getappdata(MyCell,'Props');
                    Width = Props.LMargin+Props.RMargin;
                    Child = get(MyCell,'Children');
                    if ~isempty(Child)
                        ChildSize = get(Child,'Position');
                        Width = Width + ChildSize(3);
                    end
                    MaxWidth = max(MaxWidth,Width);
                end
                CWidth(CIdx) = MaxWidth;
            end
            % 3. Proportional
            CWidthProportionalMask = Obj.ColWidthPolicy == 1;
            CWidthSumProportional = CWidthSum - sum(CWidth);
            CWidthSumProportional = max(CWidthSumProportional,0);
            CWidthWeight = Obj.ColWidthWeight(CWidthProportionalMask);
            CWidth(CWidthProportionalMask) = CWidthSumProportional*CWidthWeight/sum(CWidthWeight);
            CWidth = max(CWidth, Obj.MinColWidth);

            % Determine row heights
            RHeightSum = LayoutHeight-(Obj.NumRows-1)*Obj.HGap-Obj.TMargin-Obj.BMargin;
            RHeight = zeros(1,Obj.NumRows);
            % 1. Absolute
            RMask = Obj.RowHeightPolicy == 0;
            RHeight(RMask) = Obj.RowHeightAbsolute(RMask);
            % 2. Automatic
            for RIdx = find(Obj.RowHeightPolicy == 2)
                MaxHeight = 0;
                Mask = RSpan(:,1) == RIdx & RSpan(:,2) == RIdx;
                assert(sum(Mask) > 0, ...
                    'Any row with automatic height must have at least one non-merged cell.');
                for MyCell = Cells(Mask)'
                    Props = getappdata(MyCell,'Props');
                    Height = Props.TMargin+Props.BMargin;
                    Child = get(MyCell,'Children');
                    if ~isempty(Child)
                        ChildSize = get(Child,'Position');
                        Height = Height + ChildSize(4);
                    end
                    MaxHeight = max(MaxHeight,Height);
                end
                RHeight(RIdx) = MaxHeight;
            end
            % 3. Proportional
            RMask = Obj.RowHeightPolicy == 1;
            RHeightSumProportional = RHeightSum - sum(RHeight);
            RHeightSumProportional = max(RHeightSumProportional,0);
            Weight = Obj.RowHeightWeight(RMask);
            RHeight(RMask) = RHeightSumProportional*Weight/sum(Weight);
            RHeight = max(RHeight, Obj.MinRowHeight);

            % Horizontal and vertical offsets for each cell
            CellOffsetH = cumsum([Obj.LMargin CWidth(1:end-1)+Obj.VGap]);
            CellOffsetV = LayoutHeight-cumsum([RHeight(1)+Obj.BMargin RHeight(2:end)+Obj.HGap]);

            % Set cell positions
            for i = 1:length(Cells)
                RIdx  = RSpan(i,1):RSpan(i,2);
                CIdx  = CSpan(i,1):CSpan(i,2);

                set(Cells(i), ...
                    'Position', [ ...
                        CellOffsetH(CIdx(1)) ...
                        CellOffsetV(RIdx(end)) ...
                        sum(CWidth(CIdx))+Obj.VGap*(length(CIdx)-1) ...
                        sum(RHeight(RIdx))+Obj.HGap*(length(RIdx)-1)]);

                GridLayout.ResizeCell(Cells(i));
            end
        end

        function Value = get.ActualWidth(Obj)
            Position = getpixelposition(Obj.Container);
            Value = Position(3);
        end
        function Value = get.ActualHeight(Obj)
            Position = getpixelposition(Obj.Container);
            Value = Position(4);
        end

        function Value = get.RequiredWidth(Obj)
            if any(Obj.ColWidthPolicy == 1)
                Value = [];
                return;
            end
            Mask = Obj.ColWidthPolicy == 0;
            Value = sum(Obj.ColWidthAbsolute(Mask)) + Obj.HGap*(Obj.NumCols-1) + Obj.LMargin + Obj.RMargin;
        end
        function Value = get.RequiredHeight(Obj)
            if any(Obj.RowHeightPolicy == 1)
                Value = [];
                return;
            end
            Mask = Obj.RowHeightPolicy == 0;
            Value = sum(Obj.RowHeightAbsolute(Mask)) + Obj.VGap*(Obj.NumRows-1) + Obj.TMargin + Obj.BMargin;
        end
    end

    methods (Static)
        function FormatCell(Cell, varargin)
            assert(isscalar(Cell) && ishghandle(Cell), ...
                'Cell must be a hghandle.');
            Props = getappdata(Cell,'Props');
            for i = 1:2:length(varargin)
                Props.(varargin{i}) = varargin{i+1};
            end
            setappdata(Cell,'Props',Props);
        end
    end

    methods (Static, Access = private)
        function ResizeCell(Src)
            Child = get(Src,'Children');
            if isempty(Child)
                return;
            end

            % Special case for axes with legend
            % The legend has the same parent as the axes!
            if ~isscalar(Child)
                assert(numel(Child) == 2 && ...
                    ishghandle(Child(1),'axes') && strcmp(get(Child(1),'Tag'),'legend') && ...
                    ishghandle(Child(2),'axes'), ...
                    'If Cell has two children, they must be an ''axes'' and its legend.');
                Child = Child(2);
            end

            % Special treatment for axes objects
            IsChildAxes = ishghandle(Child,'axes');
            if IsChildAxes
                OldAxesUnits = get(Child,'Units');
                set(Child,'Units','Pixels');
                UseOuterPosition = strcmp(get(Child,'ActivePositionProperty'),'outerposition');
                if UseOuterPosition
                    ChildPosition = get(Child,'OuterPosition');
                else
                    ChildPosition = get(Child,'Position');
                end
            else
                ChildPosition = get(Child,'Position');
            end

            % Child offset and size
            ChildOffset = ChildPosition(1:2);
            ChildSize = ChildPosition(3:4);
            % Cell size
            CellPosition = getpixelposition(Src);
            CellSize = CellPosition(3:4);
            % Cell propoerties
            Props = getappdata(Src,'Props');
            % Horizontal alignment
            HAlign  = Props.HAlign;
            LMargin = Props.LMargin;
            RMargin = Props.RMargin;
            if strcmp(HAlign,'Left')
                ChildOffset(1) = LMargin;
            elseif strcmp(HAlign,'Right')
                ChildOffset(1) = CellSize(1)-ChildSize(1)-RMargin;
            elseif strcmp(HAlign,'Center')
                ChildOffset(1) = (CellSize(1)-ChildSize(1))/2;
            else % Stretch
                ChildOffset(1) = LMargin;
                ChildSize(1) = CellSize(1)-LMargin-RMargin;
            end
            % Vertical alignment
            VAlign  = Props.VAlign;
            TMargin = Props.TMargin;
            BMargin = Props.BMargin;
            if strcmp(VAlign,'Bottom')
                ChildOffset(2) = BMargin;
            elseif strcmp(VAlign,'Top')
                ChildOffset(2) = CellSize(2)-ChildSize(2)-TMargin;
            elseif strcmp(VAlign,'Center')
                ChildOffset(2) = (CellSize(2)-ChildSize(2))/2;
            else % Stretch
                ChildOffset(2) = BMargin;
                ChildSize(2) = CellSize(2)-TMargin-BMargin;
            end
            % Prevent non-positive sizes
            ChildSize = max(ChildSize,.1);

            % Update child position
            ChildPosition = [ChildOffset ChildSize];
            if IsChildAxes
                if UseOuterPosition
                    set(Child,'OuterPosition',ChildPosition);
                else
                    set(Child,'Position',ChildPosition);
                end
                set(Child,'Units',OldAxesUnits);
            else
                set(Child,'Position',ChildPosition);
            end
        end

    end
end

% EOF

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
function [Fopt,OPT_stats_out]=obj_DEMO_FEBio_iFEA_uniaxial_transiso_02(Pn,objectiveStruct)

%% Unnormalize and constrain parameters

P=Pn.*objectiveStruct.parNormFactors; %Scale back, undo normalization
P_in=P; %Proposed P

%Constraining parameters
for q=1:1:numel(P);
    [P(q)]=parLimNat(objectiveStruct.Pb_struct.xx_c(q),objectiveStruct.Pb_struct.xxlim(q,:),P(q));
end

%% SETTING MATERIAL PARAMETERS

mat_struct.par_names={{'solid','Ogden unconstrained','c1'},...
     {'solid','Ogden unconstrained','m1'},...
    {'solid','Ogden unconstrained','c2'},...
     {'solid','Ogden unconstrained','m2'},...
    {'solid','Ogden unconstrained','cp'},...
    {'solid','ellipsoidal fiber distribution','ksi'},...
    {'solid','ellipsoidal fiber distribution','beta'},...
    };

%Acces material parameters
mat_struct.id=1;

c1=P(1);
m1=P(2);
ksi=[P(3) P(3) P(3)*P(4)];
beta=P(5)*ones(1,3);
cp=(2*c1+mean(ksi))*objectiveStruct.k_factor;

mat_struct.par_values={c1 m1 c1 -m1 cp ksi beta};

disp('SETTING MATERIAL PARAMETERS...');
disp(['Proposed (norm.): ',sprintf(repmat('%6.16e ',[1,numel(Pn)]),Pn)]);
disp(['Proposed        : ',sprintf(repmat('%6.16e ',[1,numel(P_in)]),P_in)]);
disp(['Set (constr.)   : ',sprintf(repmat('%6.16e ',[1,numel(P)]),P)]);

%Assign material parameters
docNode=set_mat_par_FEBIO(objectiveStruct.FEB_struct.run_filename,objectiveStruct.FEB_struct.run_filename,{mat_struct});

disp('Done')

%%

for q=1:1:2 %Direction cases

    disp('SETTING FIBRE DIRECTIONS...');

    switch q
        case 1
            alphaFib=0;
        case 2
            alphaFib=0.5*pi;
    end

    [R,~]=euler2DCM([0,alphaFib,0]);
    v_fib=(R*[0 0 1]')';
    V_fib=v_fib(ones(numel(objectiveStruct.FEB_struct.Geometry.ElementData.MatAxis.ElementIndices),1),:);

    %Adding fibre direction, construct local orthonormal basis vectors
    [a,d]=vectorOrthogonalPair(V_fib);
    VF_E=zeros(size(V_fib,1),size(V_fib,2),2);
    VF_E(:,:,1)=a; %a1 ~ e1 ~ X or first direction
    VF_E(:,:,2)=d; %a2 ~ e2 ~ Y or second direction
    %Vf_E %a3 ~ e3 ~ Z, third direction, or fibre direction
    objectiveStruct.FEB_struct.Geometry.ElementData.MatAxis.Basis=VF_E;
    docNode=addMatAxisFibreElementData_FEB(docNode,objectiveStruct.FEB_struct);
    write_XML_no_extra_lines(objectiveStruct.FEB_struct.run_filename,docNode)% Saving XML file
    disp('Done')

    %% START FEBio NOW

    [runFlag]=runMonitorFEBio(objectiveStruct.FEBioRunStruct);

    stretch_exp=objectiveStruct.stretch_exp;
    stress_cauchy_exp=objectiveStruct.stress_cauchy_exp(:,q);

    if runFlag==1

        %Importing stress
        [~,S_mat,~]=importFEBio_logfile(objectiveStruct.FEB_struct.run_output_names{2}); %Element Cauchy stresses
        S_mat=S_mat(:,2:end,:); %Final stress, no element labels
        S_mat=squeeze(mean(S_mat,1)); %Mean across elements
        stress_cauchy_sim=[0; S_mat(:)]; %Cauchy stress
        stress_cauchy_sim=stress_cauchy_sim.*1e3; %Scale to kPa

        %Import principal strains
        [~, E_mat,~]=importFEBio_logfile(objectiveStruct.FEB_struct.run_output_names{1}); %Nodal displacements
        E_mat=E_mat(:,2:end,:);

        %Derive x stretch
        EX_mat=E_mat(:,1,:);
        EX_mat=squeeze(mean(EX_mat,1)); %Mean across elements
        stretch_sim_X=sqrt(2*EX_mat+1);
        stretch_sim_X=[1; stretch_sim_X(:)];
        stretch_sim_X_end=(stretch_sim_X(end));

        %Derive y stretch
        EY_mat=E_mat(:,2,:);
        EY_mat=squeeze(mean(EY_mat,1)); %Mean across elements
        stretch_sim_Y=sqrt(2*EY_mat+1);
        stretch_sim_Y=[1; stretch_sim_Y(:)];
        stretch_sim_Y_end=(stretch_sim_Y(end));

        %Derive z stretch
        EZ_mat=E_mat(:,3,:);
        EZ_mat=squeeze(mean(EZ_mat,1)); %Mean across elements
        stretch_sim_Z=sqrt(2*EZ_mat+1);
        stretch_sim_Z=[1; stretch_sim_Z(:)];
        stretch_sim_Z_end=(stretch_sim_Z(end));

        %Interpolate experiment onto simulated points
        stress_cauchy_sim_exp = interp1(stretch_sim_Z,stress_cauchy_sim,stretch_exp,'pchip');

        %Derive Fopt
        stressDev=stress_cauchy_exp-stress_cauchy_sim_exp;
        Fopt_stress=mean(abs(stressDev)./max(abs(stress_cauchy_exp)));

        switch q
            case 1
                stretchDev=[1./(sqrt(0.7)) 1./(sqrt(0.7))]-[stretch_sim_X_end stretch_sim_Y_end];
            case 2
                stretchDev=[exp(log(stretch_sim_Z_end)*-0.36) exp(log(stretch_sim_Z_end)*-0.65)]-[stretch_sim_X_end stretch_sim_Y_end];
        end

        Fopt_stretch=mean(abs(stretchDev));

        [R_sq]=R_squared(stress_cauchy_sim_exp,stress_cauchy_exp);

    else %Output NaN
        stress_cauchy_sim=NaN(size(stretch_exp));
        stretch_sim_Z=NaN(size(stretch_exp));
        stressDev=NaN(size(stretch_exp));
        stretchDev=NaN(1,2);
        Fopt_stress=Nan;
        Fopt_stretch=NaN;
    end

    OPT_stats_out.stress_cauchy_sim{q}=stress_cauchy_sim;
    OPT_stats_out.stretch_sim{q}=stretch_sim_Z;
    OPT_stats_out.stretch_sim_end(q,:)=[stretch_sim_X_end stretch_sim_Y_end stretch_sim_Z_end];
    OPT_stats_out.stressDev{q}=stressDev;
    OPT_stats_out.stretchDev{q}=stretchDev;
    OPT_stats_out.Fopt(:,q)=[Fopt_stress Fopt_stretch];
    OPT_stats_out.R_sq(q)=R_squared(stress_cauchy_exp,stress_cauchy_sim_exp);

end

switch objectiveStruct.method
    case 1
        Fopt=OPT_stats_out.Fopt(1,1)+OPT_stats_out.Fopt(1,2)+OPT_stats_out.Fopt(2,2); %Scalar objective function
    case 2
        Fopt=[OPT_stats_out.Fopt(1,1) OPT_stats_out.Fopt(1,2) OPT_stats_out.Fopt(2,2)]; %Objective function vectors
end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% ACQUIRE_LECROY_SCOPE_DATA M-Code for communicating with an instrument.
%
%   This is the machine generated representation of an instrument control
%   session using a device object. The instrument control session comprises
%   all the steps you are likely to take when communicating with your
%   instrument. These steps are:
%
%       1. Create a device object
%       2. Connect to the instrument
%       3. Configure properties
%       4. Invoke functions
%       5. Disconnect from the instrument
%
%   To run the instrument control session, type the name of the M-file,
%   acquire_LeCroy_scope_data, at the MATLAB command prompt.
%
%   The M-file, ACQUIRE_LECROY_SCOPE_DATA must be on your MATLAB PATH. For additional information
%   on setting your MATLAB PATH, type 'help addpath' at the MATLAB command
%   prompt.
%
%   Example:
%       acquire_LeCroy_scope_data;
%
%   See also ICDEVICE.

% Copyright 2009 - 2010 The MathWorks, Inc.
%   Creation time: 08-Jan-2009 14:37:45

% Create a TCPIP object.
interfaceObj = instrfind('Type', 'tcpip', 'RemoteHost', '127.0.0.1', 'RemotePort', 1861, 'Tag', '');

% Create the TCPIP object if it does not exist
% otherwise use the object that was found.
if isempty(interfaceObj)
    interfaceObj = tcpip('172.31.146.86', 1861);
else
    fclose(interfaceObj);
    interfaceObj = interfaceObj(1);
end

% Create a device object.
deviceObj = icdevice('lecroy_basic_driver.mdd', interfaceObj);

% Connect device object to hardware.
connect(deviceObj);

% Execute device object function(s).
groupObj = get(deviceObj, 'Waveform');
groupObj = groupObj(1);
[Y, T, XUNIT, YUNIT, HEADER] = invoke(groupObj, 'readwaveform', 'channel1');

% Disconnect device object from hardware.
disconnect(deviceObj);

% Delete objects.
delete([deviceObj interfaceObj]);


text:
function tid = timerinit(par1,par2)
%TIMERINIT Initialize a new timer.
%   TID = TIMERINIT() initializes a new timer for counting elapsed time,
%   and returns its id.
%
%   TID = TIMERINIT('TIMERNAME') sets the timer name to the specified
%   string for display purposes.
%
%   TID = TIMERINIT(ITERNUM) initializes a new ETA timer for a process with
%   ITERNUM iterations. An ETA timer can be used for both counting elapsed
%   time and estimating remaining time.
%
%   TID = TIMERINIT('TIMERNAME',ITERNUM) sets the ETA timer name to the
%   specified string for display purposes.
%
%   Example:
%
%     tid = timerinit(100);
%     for i = 1:100
%       pause(0.07);
%       timereta(tid,i,1);
%     end
%     timereta(tid,i);
%
%   See also TIMERETA, TIMERCLEAR.

%  Ron Rubinstein
%  Computer Science Department
%  Technion, Haifa 32000 Israel
%  ronrubin@cs
%
%  June 2008

global utiltbx_timer_start_times    % start times
global utiltbx_time_lastdisp        % last display times
global utiltbx_timer_iternums       % iteration numbers
global utiltbx_timer_lastiter       % last queried iteration numbers
global utiltbx_timer_name           % timer names
global utiltbx_timer_callfun        % timer calling functions

% parse function arguments %

if (nargin==0)

  iternum = -1;
  timername = '';

elseif (nargin==1)

  if (ischar(par1))
    iternum = -1;
    timername = par1;

  elseif (isnumeric(par1) && numel(par1)==1 && par1>0)
    iternum = par1;
    timername = '';

  else
    error('Invalid number of iterations');
  end

elseif (nargin==2)

  if (ischar(par1) && isnumeric(par2))
    if (numel(par2)==1 && par2>0)
      timername = par1;
      iternum = par2;
    else
      error('Invalid number of iterations');
    end
  else
    error('Invalid function syntax');
  end

else
  error('Too many function parameters');
end

% register the timer %

if (isempty(utiltbx_timer_start_times))
  utiltbx_timer_start_times = clock;
  utiltbx_time_lastdisp = utiltbx_timer_start_times;
  utiltbx_timer_iternums = double(iternum);
  utiltbx_timer_lastiter = 0;
  utiltbx_timer_name = { timername };
  utiltbx_timer_callfun = {};
  tid = 1;
else
  utiltbx_timer_start_times(end+1,:) = clock;
  utiltbx_time_lastdisp(end+1,:) = utiltbx_timer_start_times(end,:);
  utiltbx_timer_iternums(end+1) = double(iternum);
  utiltbx_timer_lastiter(end+1) = 0;
  utiltbx_timer_name{end+1} = timername;
  tid = size(utiltbx_timer_start_times,1);
end

% detect timer calling function %

st = dbstack;
if (length(dbstack) >= 2)
  utiltbx_timer_callfun{end+1} = st(2).name;
else
  utiltbx_timer_callfun{end+1} = '';
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB