text:
function [F,V,regionInd]=multiRegionTriMesh2D(regionSpec,pointSpacing,resampleCurveOpt,plotOn)

% function [F,V,regionInd]=multiRegionTriMesh2D(regionSpec,pointSpacing,plotOn)
% ------------------------------------------------------------------------
% This function creates a 2D triangulation for each of the regions
% specified in the variable regionSpec. The mesh aims to obtain a point
% spacing as defined by the input pointSpacing. The triangulation is based
% on a 2D constrained Delaunary triangulation. The constraints are formed
% by the boundary curves inside the regionSpec cell. Large areas, with
% respect to the pointSpacing, will contain a near homogeneous and
% aproximately equilateral triangulations. Other regions (e.g. at
% boundaries and thin/complex shapes) will contain other triangulation
% types.
% The function output contains the triangular faces in F, the vertices in V
% and the per triangle region indices in regionInd. By setting plotOn to 0
% or 1 plotting can be switched on or off.
%
% More on the specification of regions:
% Regions are defined as cell entries in the input variable regionSpec for
% instnace region 1 is found in regionSpec{1}. Each region entry is itself
% also a cell array containing all the boundary curves, e.g. for a two
% curve region 1 we would have something like regionSpec{1}={V1,V2} where
% V1 and V2 are the boundary curves. Multiple curves may be given here. The
% first curve should form the outer boundary of the entire region, the
% curves that follow should define holes inside this boundary and the
% space inside them is therefore not meshed. The boundary vertices for
% regions that share boundaries are merged and will share these boundary
% vertices.
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 2013/14/08
%------------------------------------------------------------------------

%% PLOT SETTINGS
if plotOn==1
    fontSize=20;
    hf1=cFigure;
    title('Smoothened triangulated mesh','FontSize',fontSize);
    xlabel('X','FontSize',fontSize);ylabel('Y','FontSize',fontSize);zlabel('Z','FontSize',fontSize);
    hold on;
end

%% MESHING REGIONS

%The total vertex, face and color (=region number) matrices
V=[]; F=[]; regionInd=[];
for qRegion=1:1:numel(regionSpec)

    %Define region cell
    regionCell=regionSpec{qRegion};

    %Meshing region
    [Fs,Vs]=regionTriMesh2D(regionCell,pointSpacing,resampleCurveOpt,0);

    %Joining regions
    F=[F;Fs+size(V,1)]; %Add new faces and fix vertex indices
    V=[V;Vs]; %Add points
    regionInd=[regionInd; qRegion*ones(size(Fs,1),1)]; %Create region index for faces
end

%% PLOTTING
if plotOn==1
    figure(hf1);
    gpatch(F,V,regionInd);
    colormap(gjet(numel(regionSpec))); icolorbar;
    axisGeom(gca,fontSize);
    drawnow;
end

%% REMOVING DOUBLE POINTS
%Removing double points (region curve points may appear multiple times)
[F,V]=mergeVertices(F,V);

%%
% _*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:
function setOptsDoubleValueInObject(handles, options, optPropName, handlesFieldName)
    strValue = strtrim(get(handles.(handlesFieldName),'String'));

    if(~isempty(strValue)) %allow empty, implies default
        value = str2double(strValue);
    else
        value = NaN;
    end
    options.(optPropName) = value;

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function lars_init()

global RESOLUTION_OF_LARS;
global REGULARIZATION_FACTOR;
RESOLUTION_OF_LARS = 1e-12;
REGULARIZATION_FACTOR = 1e-11;
% RESOLUTION_OF_LARS = 1e-10;
% REGULARIZATION_FACTOR = 1e-9;
end

author: Sable

repo: mcbench-benchmarks


text:
function computeAllPrediction_LGG(pathModels,pathResults,training,param,maxOrder,nBoot,imbalance,seed)
% -------------------------------------------------------------------------
% function computeAllPrediction_LGG(pathModels,pathResults,training,param,maxOrder,nBoot,imbalance,seed)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes prediction performance estimation for a given
% feature set type, and for all model orders of all experiments with
% different degrees of freedom. See ref. [1] for more details.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% 1. pathModels: Full path to the Models folder of the corresponding experiment.
%                --> Ex: '/myProject/WORKSPACE/LOGISTIC_REGRESSION/MODELS'
% 2. pathResults: Full path to the Results folder of the corresponding experiment.
%                --> Ex: '/myProject/WORKSPACE/LOGISTIC_REGRESSION/RESULTS'
% 3. training: Structure defining all parameters for the given experiment
%              to perform. See organizeRadiomicsExp_LGG.m for more details.
%              --> IMPORTANT VARIABLE DEFINING ALL PARAMETERS IN TRAINING
% 4. param: Cell of two strings, defining 1) The feature set type/name; and
%           2) the outcome to model in the 'training' structure.
% 5. maxOrder: Integer specifying the maximal model order to construct.
%              --> Ex: 10
% 6. nBoot: Number of bootstrap samples to use.
%           --> Ex: 100
% 7. imbalance: String specifying the type of imbalance-adjustement strategy
%               employed. Either 'IABR' for imbalance-adjusted bootstrap
%               resampling (see ref.[1]), or 'IALR' for imbalance-adjusted
%               logistic regression (formal reference to come).
%               --> Ex: 'IALR'
% 7. seed: Random generator seed for reproducibility of experiments.
% -------------------------------------------------------------------------
% OUTPUTS: Prediction performance results are saved in a folder named 'RESULTS' in the
% corresponding folder of each experiment (e.g. 'Experiment1', 'Experiment2', etc.)
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2017
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015-2017  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath = pwd;
warning off

cd(pathModels)
fSetName = param{1};
outcomeName = param{2};
outcome = training.(outcomeName).outcome;
models = load(['MODELS_',fSetName,'_',outcomeName]); models = struct2cell(models); models = models{1};

tic
fprintf(['\n --> COMPUTING PREDICTION PERFORMANCE (MODEL ORDERS OF 1 to % u) FOR  "',fSetName,'" FEATURE SET AND "',outcomeName,'" OUTCOME ... '],maxOrder)
for j = 1:maxOrder
    orderName = ['Order',num2str(j)];
    data = models.(orderName).Data;
    [orderResults] = predictionPerformanceEstimation(data,outcome,nBoot,imbalance,seed);
    results.(orderName) = orderResults;
    results.(orderName).Data = models.(orderName).Data;
    results.(orderName).Name = models.(orderName).Name;
end
cd(pathResults), save(['RESULTS_',fSetName,'_',outcomeName],'results')
fprintf('DONE!\n'), toc

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
% Test file read
% File:  c:\M-files\shortcut_updates\readtxt_TwinT.m
% 2/15/07
%
clc;clear;
ns=1e-9;
[time,v1,v3]=textread('c:\Spiceapps\datfiles\xformer_tran.txt','%f %f %f');
%
h=plot(time/ns,v1,'k',time/ns,v3,'r');
set(h,'LineWidth',2);
axis([0 200 -5 15]);
grid on
xlabel('Time (ns)');
ylabel('Volts');
title('Transformer Pulse Response');
legend('V1','V3',0);
figure(1)
%


text:
function eAltitude = altitudeefficiency(h,M,throttle,assumptions)
% Change in gas turbine core thermal efficiency as function of altitude.
%
%   relativeEfficiency = altitudeefficiency(h,M,throttle,assumptions)
%
%   h is in meters.
%
%   See also CALCULATEPSFC.

%% Interpolation method
%{
h0 = [0
    3474.72
    6797.04
    9083.04
    11277.6
    14630.4
    17068.8];

e0 = [0.846
    0.93
    0.978
    0.995
    1
    0.992
    0.966];

eAltitude = interp1(h0,e0,h,'linear');
%}

%% Quadratic curve fit from interpolation data
if nargin < 4
    assumptions.jnk = nan;
end
if ~isfield(assumptions,'efficiencyAtSeaLevel')
    assumptions.efficiencyAtSeaLevel = .846;
end
if ~isfield(assumptions,'hMaxEfficiency')
    assumptions.hMaxEfficiency = 11277.6; % meters (37000 ft)
end

k = (1-assumptions.efficiencyAtSeaLevel)./assumptions.hMaxEfficiency.^2;
eAltitude = 1-k.*(h-assumptions.hMaxEfficiency).^2;

end


text:
% add relevant paths
hodi = fileparts(which('zmap'));
path_list = [
    {hodi};
    fullfile(hodi, {'src';
    'help';
    'zmapwww';
    'importfilters';
    'resrc';
    });
    split(genpath(fullfile(hodi,'dem'))); %location for DEM data files
    fullfile(hodi, 'resrc', {'features';'sample';'focalmech'});
    fullfile(hodi, 'src',{'afterrate';
    'cgr_utils';
    'as_class';
    'declus';
    'pvals';
    'synthetic';
    'utils';
    'distance2curve';
    fullfile('system_dependent',computer)});
    fullfile(hodi,'src','cgr_utils',{'grids';'gui';'selections'});
    {fullfile(hodi, 'src', 'danijel')};
    fullfile(hodi, 'src', 'danijel', {'calc';
    'ex';
    'focal';
    'gui';
    'plot';
    'probfore'});
    {fullfile(hodi, 'src', 'jochen')};
    fullfile(hodi, 'src', 'jochen', {'auxfun';
    'ex';
    'plot';
    'seisvar';
    'stressinv'});
    fullfile(hodi, 'src', 'jochen', 'seisvar', 'calc');
    {fullfile(hodi, 'src', 'thomas')};
    fullfile(hodi, 'src', 'thomas', {'decluster';
    'etas';
    'slabanalysis';
    'seismicrates';
    });
    fullfile(hodi, 'src','thessa');
    fullfile(hodi, 'src', 'thomas', 'decluster', 'reasen');
    ];
addpath(path_list{:});

author: CelsoReyes


text:
function uninstallGibbon

% function uninstallGibbon
%-------------------------------------------------------------------------
% Change log:
% 2022/03/16 Created based on installGibbon
%-------------------------------------------------------------------------

%% Add GIBBON library path so functions are known to use here
gibbonPath=fileparts(mfilename('fullpath')); %Get the GIBBON path
addpath(fullfile(gibbonPath,'lib')); %Add gibbon lib path so gibbon functions used here are known

%% Settings

W=60; %Scrollbar width
squareSizeMax=1200;

%% Open figure

%Force groot units to be pixels
graphicalRoot=groot;
grootUnits=graphicalRoot.Units;
if ~strcmp(grootUnits,'pixels')
    graphicalRoot.Units='pixels';
end

screenSizeGroot = get(groot,'ScreenSize');

figStruct.Name='uninstalling GIBBON';

hf = cFigure(figStruct);
hf.NumberTitle='off';

squareSize=round([min(screenSizeGroot(3:4)) min(screenSizeGroot(3:4))]); % width, height
squareSize(squareSize>squareSizeMax)=squareSizeMax;

hf.Units='pixels';
windowSize=[2*W (screenSizeGroot(4)-squareSize(2))-2*W squareSize(1) squareSize(2)]; % left bottom width height

hf.Position=windowSize; % left bottom width height

%Remove tool and menu bars
ht = findobj(allchild(hf),'flat','Type','uitoolbar');
delete(ht);
ht = findobj(allchild(hf),'flat','Type','uimenu');
delete(ht);

set(hf,'ResizeFcn',{@resizeAll,{hf}});

hf.UserData.W=W;
hf.UserData.gibbonPath=gibbonPath;

% Reset groot units if a change was needed
if ~strcmp(grootUnits,'pixels')
    graphicalRoot.Units=grootUnits;
end

%% Initilize text fields

textPositionInstall=[W hf.Position(4)-W hf.Position(3)-W*2 round(W/1.5)];

top_title='uninstalling GIBBON';
hTextTitle = uicontrol(hf,'Style','text','String',top_title,...
    'Position',textPositionInstall,...
    'BackgroundColor',[1 1 1],'HorizontalAlignment','Left','FontSize',18,'FontWeight','normal');

hf.UserData.uihandles.hTextTitle=hTextTitle;

top_statement='Removing gibbon paths. Please wait...';
textPositionStatement=textPositionInstall;
textPositionStatement(2)=textPositionInstall(2)-W;

hTextStatement = uicontrol(hf,'Style','text','String',top_statement,...
    'Position',textPositionStatement,...
    'BackgroundColor',[1 1 1],'HorizontalAlignment','Left','FontSize',12,'FontWeight','bold');

hf.UserData.uihandles.hTextStatement=hTextStatement;

%% Add image

imagePath=fullfile(gibbonPath,'docs','img','GIBBON_overview.jpg');
M=importdata(imagePath);

figure(hf);
image(M);
axis tight; axis equal; axis off;
hAxis=gca;

hAxis.Units='pixels';
% [left bottom width height]
axisHeight=hf.Position(4)-7*W;
axisWidth=hf.Position(3)-2*W;
hAxis.Position=[round((hf.Position(3)-axisWidth)/2) round((hf.Position(4)-axisHeight))-6*W axisWidth axisHeight];

hf.UserData.uihandles.hAxis=hAxis;

%%
drawnow;

%% Adding temp folder if it does not exist

tempPath=fullfile(gibbonPath,'data','temp');

if exist(tempPath,'file')
    rmdir(tempPath,'s');
end

%% Adding paths

[pathNames]=getSubPaths(gibbonPath);

%Remove "hidden" folders with . in path name
logicKeep=~gcontains(pathNames,'.');
pathNames=pathNames(logicKeep);

for q=1:1:numel(pathNames)
     pathNameNow=pathNames{q};
     rmpath(pathNameNow); %Remove path
     hTextStatement.String=[top_statement,' ',num2str(round(100*q/numel(pathNames))),'% complete'];
     drawnow;
end
rmpath(fullfile(gibbonPath,'lib'));

hTextStatement.String='Done removing toolbox paths';
drawnow;
pause(0.5);

%% Unzip compressed data
hf.UserData.uihandles.hTextStatement.String='Removing unzipped data';
drawnow;

zipFolder=fullfile(hf.UserData.gibbonPath,'data','DICOM/0001_human_calf');

files = dir(fullfile(zipFolder,'*.dcm'));
files={files(1:end).name};
files=sort(files(:));
if ~isempty(files)
    for q=1:numel(files)
        delete(fullfile(zipFolder,files{q}));
    end
end

hf.UserData.uihandles.hTextStatement.String='Done removing unzipped data'; drawnow;
pause(0.5);

%% Removing documentation integrating

% hf.UserData.uihandles.hTextStatement.String='Removing documentation'; drawnow;
%
% hf.UserData.uihandles.hTextStatement.String='Done removing documentation'; drawnow;
% pause(0.5);

%%

hf.UserData.uihandles.hTextTitle.String='Done. GIBBON has been uninstalled.';

hf.UserData.uihandles.hTextStatement.String='Documentation integration removal is not supported yet. ';

hDoneButton = uicontrol('Style', 'pushbutton', 'String', 'Done!','Position',[W hf.Position(4)-W*3.5   5*W round(W/1.5)],'Callback',{@closeDone,{hf}},'FontSize',12);
hf.UserData.uihandles.hBrowseButton=hDoneButton;

end

%% Resizing

function resizeAll(~,~,inputCell)
hf=inputCell{1};
W=hf.UserData.W;

%Reposition axis
axisHeight=hf.Position(4)-7*W;
axisWidth=hf.Position(3)-2*W;
hf.UserData.uihandles.hAxis.Position=[round((hf.Position(3)-axisWidth)/2) round((hf.Position(4)-axisHeight))-6*W axisWidth axisHeight];

%Reposition gui items
if isfield(hf.UserData.uihandles,'hTextTitle')
    hf.UserData.uihandles.hTextTitle.Position=      [W hf.Position(4)-W     hf.Position(3)-W*2 round(W/1.5)];
end

if isfield(hf.UserData.uihandles,'hTextStatement')
    hf.UserData.uihandles.hTextStatement.Position=  [W hf.Position(4)-2*W   hf.Position(3)-W*2 round(W/1.5)];
end

if isfield(hf.UserData.uihandles,'hTextInfo1')
    hf.UserData.uihandles.hTextInfo1.Position=      [W hf.Position(4)-3*W hf.Position(3)-W*2 round(W/1.5)];
end

if isfield(hf.UserData.uihandles,'hTextInput1')
    hf.UserData.uihandles.hTextInput1.Position=     [W hf.Position(4)-3.5*W   hf.Position(3)-W*2 round(W/1.5)];
end

if isfield(hf.UserData.uihandles,'hBrowseButton')
    hf.UserData.uihandles.hBrowseButton.Position=     [W hf.Position(4)-W*4.5   5*W                round(W/1.5)];
end

if isfield(hf.UserData.uihandles,'hconfirmButton')
    hf.UserData.uihandles.hconfirmButton.Position=     [W hf.Position(4)-W*5.5   5*W                round(W/1.5)];
end

if isfield(hf.UserData.uihandles,'hDoneButton')
    hf.UserData.uihandles.hDoneButton.Position=     [W hf.Position(4)-W*3.5   5*W                round(W/1.5)];
end

end

%% Close/done

function closeDone(~,~,inputCell)

hf=inputCell{1};
close(hf);

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) 2017  Kevin Mattheus Moerman
%
% 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/>.


text:
function [outImg ] = crop2divisibleSize( inImg,factor)
% *************************************************************************
% Superresolution
% crop2divisibleSize
%
% Description:
% crops the image, so that it is divisible by "factor"
%
%
% Version 1.0
%
% Created by:   Armin Kappeler
% Date:         03/04/2015
%
% *************************************************************************

%check, if dividable by upscale factor
imsize_original = size(inImg);
imsize = floor(imsize_original/factor)*factor;

if (sum(imsize==imsize_original) ~= 2) %if needed, crop
    outImg = imcrop(inImg, [1 1 imsize(2)-1 imsize(1)-1]);
else
    outImg = inImg;
end

end

author: thomas-koehler

repo: SupER


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

    enumeration
        TwoBody('Analytic Two Body')
        Numerical('Numerical Integration');
    end

    properties
        name char = '';
    end

    methods
        function obj = BodyPropagationTypeEnum(name)
            obj.name = name;
        end
    end

    methods(Static)
        function listBoxStr = getListBoxStr()
            m = enumeration('BodyPropagationTypeEnum');
            [~,I] = sort({m.name});
            listBoxStr = {m(I).name};
        end

        function [ind, enum] = getIndForName(name)
            m = enumeration('BodyPropagationTypeEnum');
            [~,I] = sort({m.name});
            m = m(I);
            ind = find(ismember({m.name},name),1,'first');
            enum = m(ind);
        end

        function [enum, ind] = getEnumForListboxStr(nameStr)
            m = enumeration('BodyPropagationTypeEnum');
            ind = find(ismember({m.name},nameStr),1,'first');
            enum = m(ind);
        end
    end
end


text:
classdef MOEADM2M < ALGORITHM
% <multi> <real/integer>
% MOEA/D based on MOP to MOP
% K --- 10 --- Number of reference vectors

%------------------------------- Reference --------------------------------
% H. Liu, F. Gu, and Q. Zhang, Decomposition of a multiobjective
% optimization problem into a number of simple multiobjective subproblems,
% IEEE Transactions on Evolutionary Computation, 2014, 18(3): 450-455.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        function main(Algorithm,Problem)
            %% Parameter setting
            K = Algorithm.ParameterSet(10);

            %% Generate random population
            [W,K]      = UniformPoint(K,Problem.M);
            Problem.N  = ceil(Problem.N/K)*K;
            S          = Problem.N/K;
            Population = Problem.Initialization();
            Population = Associate(Population,W,S);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPoolLocal      = randi(S,S,K) + repmat(0:S:S*(K-1),S,1);
                MatingPoolGlobal     = randi(Problem.N,1,Problem.N);
                rnd                  = rand(S,K) < 0.7;
                MatingPoolLocal(rnd) = MatingPoolGlobal(rnd);
                Offspring  = Operator(Problem,Population,Population(MatingPoolLocal(:)));
                Population = Associate([Population,Offspring],W,S);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function options = selectSOSmodel(F,options,NonLinearParameterization,noRANK,IntegerData,UncertainData,obj)
switch options.sos.model
    case 0
        constraint_classes = constraintclass(F);
        noCOMPLICATING = ~any(ismember([7 8 9 10 12 13 14 15],constraint_classes));
        if noCOMPLICATING & ~NonLinearParameterization && noRANK && ~IntegerData && ~isa(obj,'logdet')
            options.sos.model = 1;
            if options.verbose>0;disp('Using kernel representation (options.sos.model=1).');end
        else
            if NonLinearParameterization
                if options.verbose>0;disp('Using image representation (options.sos.model=2). Nonlinear parameterization found');end
            elseif ~noRANK
                if options.verbose>0;disp('Using image representation (options.sos.model=2). SOS-rank constraint was found.');end
            elseif IntegerData
                if options.verbose>0;disp('Using image representation (options.sos.model=2). Integrality constraint was found.');end
            elseif UncertainData
                if options.verbose>0;disp('Using image representation (options.sos.model=2). Uncertain data was found.');end
            elseif isa(obj,'logdet')
                if options.verbose>0;disp('Using image representation (options.sos.model=2). Logdet objective was found.');end
            else
                if options.verbose>0;disp('Using image representation (options.sos.model=2). Integer data, KYPs or similar was found.');end
            end
            options.sos.model = 2;
        end
    case 1
        if NonLinearParameterization
            if options.verbose>0;disp('Switching to image model due to nonlinear parameterization (not supported in kernel model).');end
            options.sos.model = 2;
        end
        if ~noRANK
            if options.verbose>0;disp('Switching to image model due to SOS-rank constraints (not supported in kernel model).');end
            options.sos.model = 2;
        end
        if IntegerData
            if options.verbose>0;disp('Switching to image model due to integrality constraints (not supported in kernel model).');end
            options.sos.model = 2;
        end
    case 3
    otherwise
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function pass = test_cdr( )
% Test diskfun cdr() command

tol = 1000*chebfunpref().cheb2Prefs.chebfun2eps;

% Check it works:
f = diskfun( @(x,y,z) exp(-((x-.4).^2+(y-.9).^2)));
[C, D, R] = cdr( f );
pass(1) = norm( chebfun2(@(t,r) feval(f,t,r, 'polar'),[-pi,pi,-1,1])...
- C * D * R' ) < tol;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%GETIMHEIGHT Get featsize, function oudated
%
% IMHEIGHT = GETIMHEIGHT(A)
%
% Returns feature size and prints warning.
% This routine is outdated, it should be replaced by GETFEATSIZE or
% GETOBJSIZE. It is assumed to be GETFEATSIZE for the moment.

author: marianux

repo: ecg-kit

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


text:
% demo_triaxial     May 25, 2008

% this scriot demonstrates how to interact
% with the triaxial.m function which calculates
% the geodetic altitude of a spacecraft
% relative to a triaxial ellipsoid

% Orbital Mechanics with Matlab

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

clear all;

global flat

% conversion factor - degrees to radians

dtr = pi / 180.0;

% Earth gravitational constant (km**3/sec**2)

mu = 398600.4415;

% Earth equatorial radius (kilometers)

req = 6378.1363;

% Earth flattening factor

flat = 1.0 / 298.257;

% classical orbital elements

oev(1) = 8000.0d0;
oev(2) = 0.015d0;
oev(3) = 28.5d0 * dtr;
oev(4) = 120.0d0 * dtr;
oev(5) = 45.0d0 * dtr;
oev(6) = 30.0d0 * dtr;

% compute eci state vector

[rsc, vsc] = orb2eci(mu, oev);

% compute geodetic altitude

alt = triaxial(rsc);

clc; home;

fprintf('\n\nprogram demo_triaxial');

fprintf('\n\naltitude relative to a triaxial ellipsoid');
fprintf('\n-----------------------------------------');

fprintf('\n\naltitude =  %14.8f kilometers\n', alt);

oeprint1(mu, oev);


text:
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     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/>.

[filename,pathname,fileindex] = uigetfile('*.*','Please open a video file');
dname = uigetdir('c:/','Save images to');
dname = [dname '\'];
video_info = VideoReader([pathname filename]);

LEN = 50000;
nfmt = ['%0.' num2str(length(num2str(LEN))) 'd'];
ffmt = 'png';%file format

videoFReader = vision.VideoFileReader([pathname filename]);

Frames = 0;
% Frames = 41443;
while ~isDone(videoFReader)
    Frames = Frames+1;
    disp(['Processing ' num2str(Frames) ' frame']);
    video_image = step(videoFReader);
    resizeImage = imresize (video_image, [270, 480]);
    imwrite(resizeImage,strcat(dname,num2str(Frames,nfmt),['.' ffmt]),ffmt); %save each frame
end
disp('Finised');

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function cw=v_lpcpp2cw(pp)
%V_LPCPP2PZ LPC: Convert power spectrum polynomial in cos(w) to power spectrum zeros CW=(RP)
% pp is a polynomial such that |polyval(ra,e^jw)| = polyval(pp,cos(w))

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: v_lpcpp2cw.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nf,p1]=size(pp);
cw=zeros(nf,p1-1);
for k=1:nf
   cw(k,:)=roots(pp(k,:)).';
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
close all;
clearvars;
clc;

x = spx.ecg.synthetic.simple_multiple(500, 10);
figure;
plot(x);
xlabel('Time');
ylabel('Amplitude');
grid on;

author: indigits

repo: sparse-plex


text:
function test_suite = test_distancePointLine
%TESTDISTANCEPOINTEDGE  One-line description here, please.
%   output = test_distancePointLine(input)
%
%   Example
%   testDistancePointLine
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2012-10-24,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

test_suite = functiontests(localfunctions);

function testBasic(testCase) %#ok<*DEFNU>
point = [0 0];
line = [1 2 1 0];
testCase.assertEqual(2, distancePointLine(point, line), 'AbsTol', .01);

function testHorizontal(testCase)
% an horizontal line, with points all around
line = [2 2 2 0];
testCase.assertEqual(distancePointLine([1 1], line), 1, 'AbsTol', .01);
testCase.assertEqual(distancePointLine([2 1], line), 1, 'AbsTol', .01);
testCase.assertEqual(distancePointLine([3 1], line), 1, 'AbsTol', .01);
testCase.assertEqual(distancePointLine([4 1], line), 1, 'AbsTol', .01);
testCase.assertEqual(distancePointLine([5 1], line), 1, 'AbsTol', .01);
testCase.assertEqual(distancePointLine([5 2], line), 0, 'AbsTol', .01);
testCase.assertEqual(distancePointLine([5 3], line), 1, 'AbsTol', .01);

function testDiagonal(testCase)
% diagonal (slope 2)
line = [1 1 3 6];
testCase.assertEqual(0, distancePointLine([3 5], line));
testCase.assertEqual(sqrt(5), distancePointLine([3 0], line), 'AbsTol', .01);

function testSingleMulti(testCase)
point = [5 5];
lines = [0 0 3 0;10 0 0 3;10 10 -3 0;10 0 3 3];
exp = [5 5 5 5*sqrt(2)];
testCase.assertEqual(exp, distancePointLine(point, lines), 'AbsTol', .01);

function testSingleMultiWithInvalid(testCase)
point = [5 5];
lines = [0 0 3 0;10 0 0 3;10 10 -3 0;10 0 3 3;10 0 0 0];
exp = [5 5 5 5*sqrt(2) 5*sqrt(2)];
testCase.assertEqual(exp, distancePointLine(point, lines), 'AbsTol', .01);

function testMultiSingle(testCase)
line  = [10 10 5 0];
points = [10 10;15 10;20 10;10 0;30 10];
exp = [0; 0; 0; 10; 0];
testCase.assertEqual(exp, distancePointLine(points, line), 'AbsTol', .01);

function testMultiSingleInvalidLine(testCase)
line  = [15 10 0 0];
points = [10 10; 15 10; 20 10; 10 10; 30 10];
exp = [5; 0; 5; 5; 15];
testCase.assertEqual(exp, distancePointLine(points, line), 'AbsTol', .01);

function testMultiMulti(testCase)
lines  = [10 30 10 0; 20 30 0 10;20 40 -10 0;10 40 0 -10];
points = [14 33;15 38];
exp = [3 6 7 4;8 5 2 5];
testCase.assertEqual(exp, distancePointLine(points, lines), 'AbsTol', .01);

function testMultiMultiWithInvalid(testCase)
lines  = [10 30 10 0; 20 30 0 10;20 40 -10 0;10 40 0 -10;10 30 0 0];
points = [14 33;15 38];
exp = [3 6 7 4 5;8 5 2 5 hypot(8, 5)];
testCase.assertEqual(exp, distancePointLine(points, lines), 'AbsTol', .01);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
% Sample layer automatically
function [X, mu, sigma] = vargplvmSampleLatent(model, dim, X, startingPoint)
if nargin <4 || isempty(startingPoint)
     % This point will be initially drawn. Then, we will sample and alter
     % only one if its dimensions.
    startingPoint = 1;
end

if nargin < 4 || isempty(X)
    Xorig = model.vardist.means;
    N = size(Xorig,1);
    xmin = min(Xorig(:,dim));
    xmax = max(Xorig(:,dim));
    df = xmax - xmin;
    xmin = xmin - 4*df/N; % also catch some points before xmin
    xmax = xmax + 4*df/N;
    x = linspace(xmin,xmax, 3*N); % this is the series of changes made in a specific dimension
    X = repmat(Xorig(startingPoint,:), length(x),1); % Just select some initial point
    X(:,dim) = x';
end

if nargout > 2
    [mu sigma] = vargplvmPosteriorMeanVar(model, X);
else
    mu = vargplvmPosteriorMeanVar(model, X);
end

author: SheffieldML

repo: vargplvm


text:
function [KGain,XNEW,PT1,PT2,PT3,YNEW] = MSMTUPDT(Xinew1,Xinew2,Xinew3,Yinew1,Yinew2,Yinew3,Xbark,Ybark,Pbark,Wco,Wci,Wcin,XNEW,VTt);

PYbarkYbark1 = Wco*(Yinew1 - Ybark)*(Yinew1 - Ybark)';
PYbarkYbark2 = Wci*(Yinew2 - Ybark)*(Yinew2 - Ybark)';
PYbarkYbark3 = Wcin*(Yinew3 - Ybark)*(Yinew3 - Ybark)';

PXbarkYbark1 = Wco*(Xinew1 - Xbark)*(Yinew1 - Ybark)';
PXbarkYbark2 = Wci*(Xinew2 - Xbark)*(Yinew2 - Ybark)';
PXbarkYbark3 = Wcin*(Xinew3 - Xbark)*(Yinew3 - Ybark)';

DENUP = sum(sum(PYbarkYbark1) + sum(PYbarkYbark2) + sum(PYbarkYbark3));
NUMUP = sum(sum(PXbarkYbark1) + sum(PXbarkYbark2) + sum(PXbarkYbark3));

KGain = NUMUP/DENUP;

YNEW = (sin(XNEW)).^2 + exp(VTt);
XNEW = Xbark + KGain*(YNEW - Ybark);

PT1 = Pbark - KGain*PYbarkYbark1*KGain;
PT2 = Pbark - KGain*PYbarkYbark2*KGain;
PT3 = Pbark - KGain*PYbarkYbark3*KGain;


text:
function [ss_estimates,ss_estimates_constant,ss_estimates_exogenous,ss_estimates_contribution_exo]=ssestimates(ss_record,n,T,cband)

% function [ss_estimates]=ssestimates(ss_record,n,T,cband)
% calculates the point estimate (median), lower bound and upper bound of the steady-state from the posterior distribution
% inputs:  - cell 'ss_record': record of the gibbs sampler draws for the steady-state
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - integer 'T': number of sample time periods (defined p 7 of technical guide)
%          - scalar 'cband': confidence level for VAR coefficients
% outputs: - cell 'ss_estimates': lower bound, point estimates, and upper bound for the steady-state

% create first the cell that will contain the estimates
ss_estimates=cell(n,1);
ss_estimates_constant=cell(n,1);
ss_estimates_exogenous=cell(n,1);
ss_estimates_contribution_exo=cell(1,1);
% for each variable and each sample period, compute the median, lower and upper bound from the Gibbs sampler records
% consider variables in turn
for ii=1:n
   % consider sample periods in turn
   for jj=1:T
   % compute first the lower bound
   ss_estimates{ii,1}(1,jj)=quantile(ss_record{ii,1}(:,jj),(1-cband)/2);
   % then compute the median
   ss_estimates{ii,1}(2,jj)=quantile(ss_record{ii,1}(:,jj),0.5);
   % finally compute the upper bound
   ss_estimates{ii,1}(3,jj)=quantile(ss_record{ii,1}(:,jj),(1-(1-cband)/2));
% % %    if m>1
% % %    % same for only the constant
% % %    ss_estimates_constant{ii,1}(1,jj)=quantile(ss_record_constant{ii,1}(:,jj),(1-cband)/2);
% % %    % then compute the median
% % %    ss_estimates_constant{ii,1}(2,jj)=quantile(ss_record_constant{ii,1}(:,jj),0.5);
% % %    % finally compute the upper bound
% % %    ss_estimates_constant{ii,1}(3,jj)=quantile(ss_record_constant{ii,1}(:,jj),(1-(1-cband)/2));
% % %    % same for only the exogenous
% % %    ss_estimates_exogenous{ii,1}(1,jj)=quantile(ss_record_exogenous{ii,1}(:,jj),(1-cband)/2);
% % %    % then compute the median
% % %    ss_estimates_exogenous{ii,1}(2,jj)=quantile(ss_record_exogenous{ii,1}(:,jj),0.5);
% % %    % finally compute the upper bound
% % %    ss_estimates_exogenous{ii,1}(3,jj)=quantile(ss_record_exogenous{ii,1}(:,jj),(1-(1-cband)/2));
% % %    % compute contribution from only the exogenous part
% % %    ss_estimates_contribution_exo{ii,1}(1,jj)=(quantile(ss_record_constant{ii,1}(:,jj),0.5)-quantile(ss_record{ii,1}(:,jj),0.5));
% % %    else
% % %    ss_estimates_constant=[];
% % %    ss_estimates_exogenous=[];
% % %    end
   end
end

author: european-central-bank

repo: BEAR-toolbox


text:
function triangleDemo(varargin)
%TRIANGLEDEMO Demo file of geom2d lib: lines and circles of a triangle
%
%   Usage:
%   triangleDemo
%
%   The macro run automatically, and draw several liens and circles
%   associated with a basic triangle.
%
%   Example
%   triangleDemo
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-11-04,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

%% Triangle

% defines vertices
p1 = [2 4];
p2 = [18 6];
p3 = [4 16];

% concatenates vertices to form the polygon
triangle = [p1; p2; p3];

% draw the triangle
figure(1); clf;
hold on;
axis([0 20 0 20]);
axis equal;
drawPolygon(triangle, 'linewidth', 3);
drawPoint(triangle, 'marker', 'o', 'markersize', 10, 'linewidth', 2, ...
    'markerFaceColor', 'w');

%% altitudes

% create lines associated with each triangle edge
edge1 = createLine(p2, p3);
edge2 = createLine(p1, p3);
edge3 = createLine(p1, p2);

% altitudes of the triangle
alt1 = orthogonalLine(edge1, p1);
alt2 = orthogonalLine(edge2, p2);
alt3 = orthogonalLine(edge3, p3);

% compute also feet
foot1 = intersectLines(edge1, alt1);
foot2 = intersectLines(edge2, alt2);
foot3 = intersectLines(edge3, alt3);

% draw altitudes
drawLine(alt1, 'color', [0 0 .8]);
drawLine(alt2, 'color', [0 0 .8]);
drawLine(alt3, 'color', [0 0 .8]);

% draw feet
feet = [foot1; foot2 ;foot3];
drawPoint(feet, 'marker', 's', 'color', [0 0 .8], 'linewidth', 2, 'markerFaceColor', 'w');

orthoCenter = intersectLines(alt1, alt2);
drawPoint(orthoCenter, ...
    'marker', 'o', 'color', [0 0 .8], 'linewidth', 2, 'markerfacecolor', 'w');

%% Median rays and inscribed circle

% compute rays emanating from each vertex
ray1 = bisector(p2, p1, p3);
ray2 = bisector(p3, p2, p1);
ray3 = bisector(p1, p3, p2);

% draw rays (all in one call)
drawRay([ray1 ; ray2 ; ray3], 'color', [0 .5 0]);

% center of inscribed circle (assimilates rays to lines)
innerCircleCenter = intersectLines(ray1, ray2);

% radius of iscribed circle is computed as the distance to one of the sides
innerRadius = distancePointLine(innerCircleCenter, edge1);

% create the circle
innerCircle = [innerCircleCenter innerRadius];

% draw the inner circle
drawCircle(innerCircle, 'color', [0 .5 0], 'linewidth', 2);
drawPoint(innerCircleCenter, ...
    'color',  [0 .5 0], 'markerFaceColor', 'w', 'linewidth', 2);

%% Circumscribed circle

% edges midpoints
mid12 = midPoint(p1, p2);
mid13 = midPoint(p1, p3);
mid23 = midPoint(p2, p3);

% draw midpoints
midPoints = [mid12 ; mid13; mid23];
drawPoint(midPoints, 'marker', 's', 'color', 'r', 'linewidth', 2, 'markerFaceColor', 'w');

% perpendicular bisectors associated with each side
perp1 = orthogonalLine(edge1, mid23);
perp2 = orthogonalLine(edge2, mid13);
perp3 = orthogonalLine(edge3, mid12);

% draw perpendicular bisectors
drawLine([perp1 ; perp2 ; perp3], 'color', 'r');

% compute orthogonal circle
orthoCenter = intersectLines(perp1, perp2);
orthoRadius = distancePoints(orthoCenter, p1);

% draw Orthogonal center
drawCircle([orthoCenter orthoRadius], 'color', 'r', 'linewidth', 2);
drawPoint(orthoCenter, 'color', 'r', 'markerFaceColor', 'w', 'linewidth', 2);


text:
function om = add_legacy_cost(om, name, idx, varargin)
%ADD_LEGACY_COST  Adds a set of user costs to the model.
%
%   OM.ADD_LEGACY_COST(NAME, CP);
%   OM.ADD_LEGACY_COST(NAME, CP, VARSETS);
%   OM.ADD_LEGACY_COST(NAME, IDX_LIST, CP);
%   OM.ADD_LEGACY_COST(NAME, IDX_LIST, CP, VARSETS);
%
%   Adds a named block of user-defined costs to the model. Each set is
%   defined by the CP struct described below. All user-defined sets of
%   costs are combined together into a single set of cost parameters in
%   a single CP struct by BULD_COST_PARAMS. This full aggregate set of
%   cost parameters can be retreived from the model by GET_COST_PARAMS.
%
%   Examples:
%       cp1 = struct('N', N1, 'Cw', Cw1);
%       cp2 = struct('N', N2, 'Cw', Cw2, 'H', H, 'dd', dd, ...
%                     'rh', rh, 'kk', kk, 'mm', mm);
%       om.add_legacy_cost('usr1', cp1, {'Pg', 'Qg', 'z'});
%       om.add_legacy_cost('usr2', cp2, {'Vm', 'Pg', 'Qg', 'z'});
%
%       om.init_indexed_name('c', {2, 3});
%       for i = 1:2
%         for j = 1:3
%           om.add_legacy_cost('c', {i, j}, cp(i,j), ...);
%         end
%       end
%
%   Let x refer to the vector formed by combining the specified VARSETS,
%   and f_u(x, CP) be the cost at x corresponding to the cost parameters
%   contained in CP, where CP is a struct with the following fields:
%       N      - nw x nx sparse matrix (optional, identity matrix by default)
%       Cw     - nw x 1 vector
%       H      - nw x nw sparse matrix (optional, all zeros by default)
%       dd, mm - nw x 1 vectors (optional, all ones by default)
%       rh, kk - nw x 1 vectors (optional, all zeros by default)
%
%   These parameters are used as follows to compute f_u(x, CP)
%
%       R  = N*x - rh
%
%               /  kk(i),  R(i) < -kk(i)
%       K(i) = <   0,     -kk(i) <= R(i) <= kk(i)
%               \ -kk(i),  R(i) > kk(i)
%
%       RR = R + K
%
%       U(i) =  /  0, -kk(i) <= R(i) <= kk(i)
%               \  1, otherwise
%
%       DDL(i) = /  1, dd(i) = 1
%                \  0, otherwise
%
%       DDQ(i) = /  1, dd(i) = 2
%                \  0, otherwise
%
%       Dl = diag(mm) * diag(U) * diag(DDL)
%       Dq = diag(mm) * diag(U) * diag(DDQ)
%
%       w = (Dl + Dq * diag(RR)) * RR
%
%       f_u(x, CP) = 1/2 * w'*H*w + Cw'*w
%
%   See also OPT_MODEL, PARAMS_LEGACY_COST, EVAL_LEGACY_COST.

%   MATPOWER
%   Copyright (c) 2008-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if iscell(idx)
    cp = varargin{1};
    args = varargin(2:end);
else                            %% simple named set
    cp = idx;
    args = varargin;
    idx = {};
end

if isempty(args)
    varsets = {};
else
    varsets = args{1};
end

%% convert varsets from cell to struct array if necessary
varsets = om.varsets_cell2struct(varsets);
nv = om.varsets_len(varsets);   %% number of variables

if isfield(cp, 'N')
    [nw, nx] = size(cp.N);
else
    nw = length(cp.Cw);
    nx = nw;
    cp.N = speye(nw, nx);
end

%% check sizes
if nx ~= nv
    if nw == 0
        cp.N = sparse(nw, nx);
    else
        error('@opt_model/add_legacy_cost: number of columns in N (%d x %d) does not match\nnumber of variables (%d)\n', nw, nx, nv);
    end
end
if size(cp.Cw, 1) ~= nw
    error('@opt_model/add_legacy_cost: number of rows of Cw (%d x %d) and N (%d x %d) must match\n', size(cp.Cw), nw, nx);
end
if isfield(cp, 'H') && (size(cp.H, 1) ~= nw || size(cp.H, 2) ~= nw)
    error('@opt_model/add_legacy_cost: both dimensions of H (%d x %d) must match the number of rows in N (%d x %d)\n', size(cp.H), nw, nx);
end
if isfield(cp, 'dd') && size(cp.dd, 1) ~= nw
    error('@opt_model/add_legacy_cost: number of rows of dd (%d x %d) and N (%d x %d) must match\n', size(cp.dd), nw, nx);
end
if isfield(cp, 'rh') && size(cp.rh, 1) ~= nw
    error('@opt_model/add_legacy_cost: number of rows of rh (%d x %d) and N (%d x %d) must match\n', size(cp.rh), nw, nx);
end
if isfield(cp, 'kk') && size(cp.kk, 1) ~= nw
    error('@opt_model/add_legacy_cost: number of rows of kk (%d x %d) and N (%d x %d) must match\n', size(cp.kk), nw, nx);
end
if isfield(cp, 'mm') && size(cp.mm, 1) ~= nw
    error('@opt_model/add_legacy_cost: number of rows of mm (%d x %d) and N (%d x %d) must match\n', size(cp.mm), nw, nx);
end

%% add the legacy cost set
om.add_named_set('cost', name, idx, nw, cp, varsets);

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function meshfaces_test02 ( )

%*****************************************************************************80
%
%% MESHFACES_TEST02: A rectangle is subdivided into three.
%
%  Modified:
%
%    23 August 2014
%
%  Author:
%
%    Darren Engwirda
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MESHFACES_TEST02\n' );
  fprintf ( 1, '  A rectangle is subdivided into three.\n' );
  fprintf ( 1, '  The middle rectangle is thin.\n' );

  node = [0.0, 0.0; 1.0,0.0; 1.0,1.0; 0.0,1.0; 1.01,0.0; 1.01,1.0; 3.0,0.0; 3.0,1.0];
  edge = [1,2; 2,3; 3,4; 4,1; 2,5; 5,6; 6,3; 5,7; 7,8; 8,6];

  face{1} = [1,2,3,4];
  face{2} = [5,6,7,2];
  face{3} = [8,9,10,6];
%
%  Since we don't save the output, we will just see an image of the mesh.
%
  meshfaces ( node, edge, face );

  filename = 'test02.png'
  print ( '-dpng', filename );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  An image of the mesh was saved as "%s"\n', filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = cumsum( varargin )

%Disciplined convex/geometric programming information for SUM:
%   CUMSUM(X) and CUMSUM(X,DIM) are vectorized forms of addition. So
%   when CUMSUM is used in a DCP or DGP, elements in each subvector
%   must satisfy the corresponding combination rules for addition (see
%   PLUS). For example, suppose that X looks like this:
%      X = [ convex concave affine  ;
%            affine concave concave ]
%   Then CUMSUM(X,1) would be permittted, but CUMSUM(X,2) would not,
%   because the top row contains the sum of convex and concave terms, in
%   violation of the DCP ruleset. For DGPs, addition rules dictate that
%   the elements of X must be log-convex or log-affine.

persistent P
if isempty( P ),
    P.map      = cvx_remap( { 'constant' ; 'affine' ; 'convex' ; 'concave' } );
    P.funcs    = { @cumsum_1, @cumsum_2, @cumsum_2 };
    P.zero     = 0;
    P.reduce   = false;
    P.reverse  = true;
    P.dimarg   = 2;
    P.constant = 1;
    P.name     = 'cumsum';
end
y = cvx_reduce_op( P, varargin{:} );

function x = cumsum_1( x )
x = builtin( 'cumsum', x, 2 );

function x = cumsum_2( x )
s = x.size_;
if s(2) ~= 1,
    b = reshape( x.basis_, [], s(2) );
    b = cumsum( b, 2 );
    x = cvx( s, reshape( b, [], prod(s) ) );
end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

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


text:
% finputcheck() - check Matlab function {'key','value'} input argument pairs
%
% Usage: >> result = finputcheck( varargin, fieldlist );
%        >> [result varargin] = finputcheck( varargin, fieldlist, ...
%                                              callingfunc, mode, verbose );
% Input:
%   varargin  - Cell array 'varargin' argument from a function call using 'key',
%               'value' argument pairs. See Matlab function 'varargin'.
%               May also be a structure such as struct(varargin{:})
%   fieldlist - A 4-column cell array, one row per 'key'. The first
%               column contains the key string, the second its type(s),
%               the third the accepted value range, and the fourth the
%               default value.  Allowed types are 'boolean', 'integer',
%               'real', 'string', 'cell' or 'struct'.  For example,
%                       {'key1' 'string' { 'string1' 'string2' } 'defaultval_key1'}
%                       {'key2' {'real' 'integer'} { minint maxint } 'defaultval_key2'}
%  callingfunc - Calling function name for error messages. {default: none}.
%  mode        - ['ignore'|'error'] ignore keywords that are either not specified
%                in the fieldlist cell array or generate an error.
%                {default: 'error'}.
%  verbose     - ['verbose', 'quiet'] print information. Default: 'verbose'.
%
% Outputs:
%   result     - If no error, structure with 'key' as fields and 'value' as
%                content. If error this output contain the string error.
%   varargin   - residual varagin containing unrecognized input arguments.
%                Requires mode 'ignore' above.
%
% Note: In case of error, a string is returned containing the error message
%       instead of a structure.
%
% Example (insert the following at the beginning of your function):
%	result = finputcheck(varargin, ...
%               { 'title'         'string'   []       ''; ...
%                 'percent'       'real'     [0 1]    1 ; ...
%                 'elecamp'       'integer'  [1:10]   [] });
%   if isstr(result)
%       error(result);
%   end
%
% Note:
%   The 'title' argument should be a string. {no default value}
%   The 'percent' argument should be a real number between 0 and 1. {default: 1}
%   The 'elecamp' argument should be an integer between 1 and 10 (inclusive).
%
%   Now 'g.title' will contain the title arg (if any, else the default ''), etc.
%
% Author: Arnaud Delorme, CNL / Salk Institute, 10 July 2002

% Copyright (C) Arnaud Delorme, CNL / Salk Institute, 10 July 2002, 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 [g, varargnew] = finputcheck( vararg, fieldlist, callfunc, mode, verbose )

	if nargin < 2
		help finputcheck;
		return;
	end;
	if nargin < 3
		callfunc = '';
	else
		callfunc = [callfunc ' ' ];
	end;
    if nargin < 4
        mode = 'do not ignore';
    end;
    if nargin < 5
        verbose = 'verbose';
    end;
	NAME = 1;
	TYPE = 2;
	VALS = 3;
	DEF  = 4;
	SIZE = 5;

	varargnew = {};
	% create structure
	% ----------------
	if ~isempty(vararg)
        if isstruct(vararg)
            g = vararg;
        else
            for index=1:length(vararg)
                if iscell(vararg{index})
                    vararg{index} = {vararg{index}};
                end;
            end;
            try
                g = struct(vararg{:});
            catch
                vararg = removedup(vararg, verbose);
                try
                    g = struct(vararg{:});
                catch
                    g = [ callfunc 'error: bad ''key'', ''val'' sequence' ]; return;
                end;
            end;
        end;
	else
		g = [];
	end;

	for index = 1:size(fieldlist,NAME)
		% check if present
		% ----------------
		if ~isfield(g, fieldlist{index, NAME})
			g = setfield( g, fieldlist{index, NAME}, fieldlist{index, DEF});
		end;
		tmpval = getfield( g, {1}, fieldlist{index, NAME});

		% check type
		% ----------
        if ~iscell( fieldlist{index, TYPE} )
            res = fieldtest( fieldlist{index, NAME},  fieldlist{index, TYPE}, ...
                           fieldlist{index, VALS}, tmpval, callfunc );
            if isstr(res), g = res; return; end;
        else
            testres = 0;
            tmplist = fieldlist;
            for it = 1:length( fieldlist{index, TYPE} )
                if ~iscell(fieldlist{index, VALS})
                     res{it} = fieldtest(  fieldlist{index, NAME},  fieldlist{index, TYPE}{it}, ...
                                           fieldlist{index, VALS}, tmpval, callfunc );
                else res{it} = fieldtest(  fieldlist{index, NAME},  fieldlist{index, TYPE}{it}, ...
                                           fieldlist{index, VALS}{it}, tmpval, callfunc );
                end;
                if ~isstr(res{it}), testres = 1; end;
            end;
            if testres == 0,
                g = res{1};
                for tmpi = 2:length(res)
                    g = [ g 10 'or ' res{tmpi} ];
                end;
                return;
            end;
        end;
	end;

    % check if fields are defined
	% ---------------------------
	allfields = fieldnames(g);
	for index=1:length(allfields)
		if isempty(strmatch(allfields{index}, fieldlist(:, 1)', 'exact'))
			if ~strcmpi(mode, 'ignore')
				g = [ callfunc 'error: undefined argument ''' allfields{index} '''']; return;
			end;
			varargnew{end+1} = allfields{index};
			varargnew{end+1} = getfield(g, {1}, allfields{index});
		end;
	end;

function g = fieldtest( fieldname, fieldtype, fieldval, tmpval, callfunc );
	NAME = 1;
	TYPE = 2;
	VALS = 3;
	DEF  = 4;
	SIZE = 5;
    g = [];

    switch fieldtype
     case { 'integer' 'real' 'boolean' 'float' },
      if ~isnumeric(tmpval) && ~islogical(tmpval)
          g = [ callfunc 'error: argument ''' fieldname ''' must be numeric' ]; return;
      end;
      if strcmpi(fieldtype, 'boolean')
          if tmpval ~=0 && tmpval ~= 1
              g = [ callfunc 'error: argument ''' fieldname ''' must be 0 or 1' ]; return;
          end;
      else
          if strcmpi(fieldtype, 'integer')
              if ~isempty(fieldval)
                  if (any(isnan(tmpval(:))) && ~any(isnan(fieldval))) ...
                          && (~ismember(tmpval, fieldval))
                      g = [ callfunc 'error: wrong value for argument ''' fieldname '''' ]; return;
                  end;
              end;
          else % real or float
              if ~isempty(fieldval) && ~isempty(tmpval)
                  if any(tmpval < fieldval(1)) || any(tmpval > fieldval(2))
                      g = [ callfunc 'error: value out of range for argument ''' fieldname '''' ]; return;
                  end;
              end;
          end;
      end;

     case 'string'
      if ~isstr(tmpval)
          g = [ callfunc 'error: argument ''' fieldname ''' must be a string' ]; return;
      end;
      if ~isempty(fieldval)
          if isempty(strmatch(lower(tmpval), lower(fieldval), 'exact'))
              g = [ callfunc 'error: wrong value for argument ''' fieldname '''' ]; return;
          end;
      end;

     case 'cell'
      if ~iscell(tmpval)
          g = [ callfunc 'error: argument ''' fieldname ''' must be a cell array' ]; return;
      end;

     case 'struct'
      if ~isstruct(tmpval)
          g = [ callfunc 'error: argument ''' fieldname ''' must be a structure' ]; return;
      end;

     case 'function_handle'
      if ~isa(tmpval, 'function_handle')
          g = [ callfunc 'error: argument ''' fieldname ''' must be a function handle' ]; return;
      end;

     case '';
     otherwise, error([ 'finputcheck error: unrecognized type ''' fieldname '''' ]);
    end;

% remove duplicates in the list of parameters
% -------------------------------------------
function cella = removedup(cella, verbose)
% make sure if all the values passed to unique() are strings, if not, exist
%try
    [tmp indices] = unique_bc(cella(1:2:end));
    if length(tmp) ~= length(cella)/2
        myfprintf(verbose,'Note: duplicate ''key'', ''val'' parameter(s), keeping the last one(s)\n');
    end;
    cella = cella(sort(union(indices*2-1, indices*2)));
%catch
    % some elements of cella were not string
%    error('some ''key'' values are not string.');
%end;

function myfprintf(verbose, varargin)

if strcmpi(verbose, 'verbose')
    fprintf(varargin{:});
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% dipfit_nonlinear() - perform nonlinear dipole fit on one of the components
%                   to improve the initial dipole model. Only selected dipoles
%                   will be fitted.
%
% Usage:
%  >> EEGOUT = dipfit_nonlinear( EEGIN, optarg)
%
% Inputs:
%    ...
%
% Optional inputs are specified in key/value pairs and can be:
%    ...
%
% Output:
%    ...
%
% Author: Robert Oostenveld, SMI/FCDC, Nijmegen 2003
%         Thanks to Nicolas Robitaille for his help on the CTF MEG
%         implementation

% SMI, University Aalborg, Denmark http://www.smi.auc.dk/
% FC Donders Centre, University Nijmegen, the Netherlands http://www.fcdonders.kun.nl

% Copyright (C) 2003 Robert Oostenveld, SMI/FCDC roberto@smi.auc.dk
%
% 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 [EEGOUT] = dipfit_nonlinear( EEG, varargin )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% convert the optional arguments into a configuration structure that can be
% understood by FIELDTRIPs dipolefitting function
if nargin>2
  cfg = struct(varargin{:});
else
  help dipfit_nonlinear
  return
end

% specify the FieldTrip DIPOLEFITTING configuration
cfg.model      = 'moving';
cfg.gridsearch = 'no';
if ~isfield(cfg, 'nonlinear')
  % if this flag is set to 'no', only the dipole moment will be fitted
  cfg.nonlinear  = 'yes';
end
% add some additional settings from EEGLAB to the configuration
tmpchanlocs    = EEG.chanlocs;
cfg.channel    = { tmpchanlocs(EEG.dipfit.chansel).labels };
if isfield(EEG.dipfit, 'vol')
    cfg.vol        = EEG.dipfit.vol;
elseif isfield(EEG.dipfit, 'hdmfile')
    cfg.hdmfile    = EEG.dipfit.hdmfile;
else
    error('no head model in EEG.dipfit')
end

if isfield(EEG.dipfit, 'elecfile') & ~isempty(EEG.dipfit.elecfile)
    cfg.elecfile = EEG.dipfit.elecfile;
end
if isfield(EEG.dipfit, 'gradfile') & ~isempty(EEG.dipfit.gradfile)
    cfg.gradfile = EEG.dipfit.gradfile;
end

% set up the initial dipole model based on the one in the EEG structure
cfg.dip.pos = EEG.dipfit.model(cfg.component).posxyz;
cfg.dip.mom = EEG.dipfit.model(cfg.component).momxyz';
cfg.dip.mom = cfg.dip.mom(:);

% convert the EEGLAB data structure into a structure that looks as if it
% was computed using FIELDTRIPs componentanalysis function
comp = eeglab2fieldtrip(EEG, 'componentanalysis', 'dipfit');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Added code to handle CTF data with multipleSphere head model           %
%  This code is copy-pasted in dipfit_gridSearch, dipfit_nonlinear        %
%  The flag .isMultiSphere is used by dipplot                             %
%  Nicolas Robitaille, January 2007.                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Do some trick to force fieldtrip to use the multiple sphere model
if strcmpi(EEG.dipfit.coordformat, 'CTF')
   cfg = rmfield(cfg, 'channel');
   comp = rmfield(comp, 'elec');
   cfg.gradfile = EEG.dipfit.chanfile;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% END                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% fit the dipoles to the ICA component(s) of interest using FIELDTRIPs
% dipolefitting function

currentPath = pwd;
ptmp = which('ft_prepare_vol_sens');
ptmp = fileparts(ptmp);
if isempty(ptmp), error('Path to "forward" folder of Fieldtrip missing'); end;
cd(fullfile(ptmp, 'private'));
try,
    if exist('hlp_microcache','file')
        source = hlp_microcache('dipfit',@ft_dipolefitting,cfg, comp);
    else
        source = ft_dipolefitting(cfg, comp);
    end
catch,
    cd(currentPath);
    lasterr
    error(lasterr);
end;
cd(currentPath);

% reformat the output dipole sources into EEGLABs data structure
EEG.dipfit.model(cfg.component).posxyz  = source.dip.pos;
EEG.dipfit.model(cfg.component).momxyz  = reshape(source.dip.mom, 3, length(source.dip.mom)/3)';
EEG.dipfit.model(cfg.component).diffmap = source.Vmodel - source.Vdata;
EEG.dipfit.model(cfg.component).sourcepot = source.Vmodel;
EEG.dipfit.model(cfg.component).datapot   = source.Vdata;
EEG.dipfit.model(cfg.component).rv        = source.dip.rv;
%EEG.dipfit.model(cfg.component).rv = sum((source.Vdata - source.Vmodel).^2) / sum( source.Vdata.^2 );

EEGOUT = EEG;

author: PatternRecognition

repo: OpenBMI


text:
function filum_test ( )

%*****************************************************************************80
%
%% FILUM_TEST tests the FILUM library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2011
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FILUM_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the FILUM library.\n' );

  filum_test02 ( );
  filum_test03 ( );
  filum_test06 ( );
  filum_test085 ( );

  filum_test11 ( );
  filum_test13 ( );
  filum_test14 ( );
  filum_test15 ( );
  filum_test165 ( );
  filum_test17 ( );

  filum_test22 ( );
  filum_test225 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FILUM_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function [ y, next ] = p25_equil ( neqn, next )

%*****************************************************************************80
%
%% P25_EQUIL returns equilibrium solutions of problem p25.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, integer NEXT, the index of the previous
%    equilibrium, which should be 0 on first call.
%
%    Output, real Y(NEQN), the "next" equilibrium solution, if any.
%
%    Output, integer NEXT, the index of the current equilibrium,
%    or 0 if there are no more.
%
  if ( next == 0 )
    next = 1;
    y(1:neqn,1) = 0.0;
  else
    next = 0;
    y = [];
  end

  return
end


text:
classdef TestSequentialLaminateTestedWithNumerics < handle

    properties (Access = private)
        microFile = 'RVE_Square_Triangle_FineFine';
        fileOutputName = 'SeqLaminate';
    end

    properties (Access = protected)
        LaminateDirection
        Theta
        FractionVolume

        StiffTensor
        WeakTensor
        MaterialValues

        MixtureCh
        Rank2Ch
        NumericalCh
        SeqLamCh

        FiberDirection
    end

    methods (Access = public)

        function hasPassed = hasPassed(obj)
            ChNum   = obj.NumericalCh.getValue();
            ChSL    = obj.SeqLamCh.getValue();
            ChMix   = obj.MixtureCh.getValue();
            ChRank  = obj.Rank2Ch.getValue();
            firstCondition  = obj.relativeNorm(ChSL,ChNum)   < 1e-2;
            secondCondition = obj.relativeNorm(ChMix,ChNum)  < 1e-3;
            thirdCondition  = obj.relativeNorm(ChRank,ChNum) < 1e-10;
            hasPassed = firstCondition & secondCondition & thirdCondition;
        end

    end

    methods (Access = protected)

        function compute(obj)
            obj.init();
            obj.computeNumericallyChForLaminate();
            obj.loadFractionVolume()
            obj.computeWeakAndStiffTensorsFromNumericalHomogenizerData();
            obj.computeSequentialLaminateTensor();
            obj.computeRank2HomogenizerTensor();
            obj.computeMixtureTheoryTensor();
        end

    end

    methods (Access = private)

        function init(obj)
           obj.loadLaminateDirection()
           obj.loadFiberDirection()
        end

        function computeNumericallyChForLaminate(obj)
           d = obj.createNumericalHomogenizerDataBase();
           homog = NumericalHomogenizer(d);
           homog.compute();
           obj.NumericalCh    = obj.rotateCh(homog);
           obj.MaterialValues = homog.matValues;
           obj.FractionVolume = homog.cellVariables.volume;
        end

        function d = createNumericalHomogenizerDataBase(obj)
            nDB = NumericalHomogenizerDataBase(obj.microFile);
            d = nDB.dataBase;
            d.outFileName = obj.fileOutputName;
            d.hasToCaptureImage = false;
        end

        function loadFractionVolume(obj)
           obj.Theta = obj.FractionVolume;
        end

        function computeWeakAndStiffTensorsFromNumericalHomogenizerData(obj)
            E1  = obj.MaterialValues.E_plus;
            nu1 = obj.MaterialValues.nu_plus;
            E0  = obj.MaterialValues.E_minus;
            nu0 = obj.MaterialValues.nu_minus;
            obj.StiffTensor = IsotropicConstitutiveTensor(E1,nu1);
            obj.WeakTensor  = IsotropicConstitutiveTensor(E0,nu0);
        end

        function computeSequentialLaminateTensor(obj)
            C0 = obj.WeakTensor;
            C1 = obj.StiffTensor;
            dir{1} = obj.LaminateDirection;
            m1 = 1;
            SeqHomog = VoigtHomogPlaneStressHomogenizer(C0,C1,dir,m1,obj.Theta);
            obj.SeqLamCh  = SeqHomog.getPlaneStressHomogenizedTensor();
        end

        function computeRank2HomogenizerTensor(obj)
            C0 = obj.WeakTensor;
            C1 = obj.StiffTensor;
            dir{1} = obj.LaminateDirection;
            m1 = 1;
            SeqHomog = VoigtPlaneStressHomogHomogenizer(C0,C1,dir,m1,obj.Theta);
            obj.Rank2Ch  = SeqHomog.getPlaneStressHomogenizedTensor();
        end

        function computeMixtureTheoryTensor(obj)
            C1 = obj.StiffTensor;
            C0 = obj.WeakTensor;
            d = [0 0 1];
            dir = Vector3D;
            dir.setValue(d);
            dir.normalize();
            angle = -acos(dot(obj.FiberDirection.getValue(),[1 0 0]));
            vFrac = obj.Theta;
            homogenizer = MixtureTheoryHomogenizer(C1,C0,dir,angle,vFrac);
            obj.MixtureCh = homogenizer.Ch;
        end

        function Ch = rotateCh(obj,homog)
            dir = obj.FiberDirection;
            r = ChRotatorForFiberHomogenizer();
            Ch = r.rotate(dir,homog.cellVariables.Ch());
        end

    end

    methods (Access = private, Static)

        function relNorm = relativeNorm(A,B)
            relNorm = norm(A - B)/norm(B);
        end

    end

    methods (Abstract,Access = protected)
        loadLaminateDirection(obj)
        loadFiberDirection(obj)
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function b = r8to_vxm ( n, a, x )

%*****************************************************************************80
%
%% R8TO_VXM multiplies a vector by a R8TO matrix.
%
%  Discussion:
%
%    The R8TO storage format is used for a Toeplitz matrix, which is constant
%    along diagonals.  Thus, in an N by N Toeplitz matrix, there are at most
%    2*N-1 distinct entries.  The format stores the N elements of the first
%    row, followed by the N-1 elements of the first column (skipping the
%    entry in the first row).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, real A(2*N-1), the R8TO matrix.
%
%    Input, real X(N), the vector to be multiplied by A.
%
%    Output, real B(N), the product A' * X.
%
  for i = 1 : n

    b(i) = a(i:-1:1) * x(1:i)' + a(n+1:2*n-i) * x(i+1:n)';

  end

  return
end


text:
function images = loadMNISTImages(filename)
%loadMNISTImages returns a 28x28x[number of MNIST images] matrix containing
%the raw MNIST images

fp = fopen(filename, 'rb');
assert(fp ~= -1, ['Could not open ', filename, '']);

magic = fread(fp, 1, 'int32', 0, 'ieee-be');
assert(magic == 2051, ['Bad magic number in ', filename, '']);

numImages = fread(fp, 1, 'int32', 0, 'ieee-be');
numRows = fread(fp, 1, 'int32', 0, 'ieee-be');
numCols = fread(fp, 1, 'int32', 0, 'ieee-be');

images = fread(fp, inf, 'unsigned char');
images = reshape(images, numCols, numRows, numImages);
images = permute(images,[2 1 3]);

fclose(fp);

% Reshape to #pixels x #examples
images = reshape(images, size(images, 1) * size(images, 2), size(images, 3));
% Convert to double and rescale to [0,1]
images = double(images) / 255;

end

repo: MachineLearning


text:
function [err,PatchHandles] = ShowPlane (EqPlane,Opt)
%
%   [err,PatchHandles] = ShowPlane (EqPlane,Opt)
%
%Purpose:
%   This function displays a plane with equation ax + by + cz +d = 0
%   The plane is defined by one sqaure patch
%   The planes span the current axis limits of the figure.
%
%Input Parameters:
%   EqPlane is a Nx4 matrix containing the equation of the plane containing each
%       triplet in  Triplets. The plane passing by triplet i is speicifed in
%       EqPlane(i,:) the plane would be
%       EqPlane(i,1)x + EqPlane(i,2)y + EqPlane(i,3)z + EqPlane(i,4) = 0
%   Opt is the options structure
%     .Fig is a handle to the figure you want the planes displayed in,
%         default is the current figure.
%     .WriteIV if this string is not empty, the planes that are displayed
%         on the graph are written to an inventor format file
%     .units 'mm' or 'tesscon'. default is tesscon
%       if you specify mm, then the xyz coordinates are transformed to tesscon
%       before writing them out. (*319.7). The idea is to write all inventor
%       files in tesscon units. This option will only be used if WriteIV is not empty.
%     .OvrWrite (0/1) default is 0, flag for overwriting existing .iv file
%
%Output Parameters:
%   err : 0 No Problem
%       : 1 Mucho Problems
%
%   PatchHandles : the handle to the patches displayed on the figure
%
%More Info :
%   see Plane_Equation
%
%
%
%     Author : Ziad Saad
%     Date : Thu Oct 22 20:19:36 CDT 1998

%Define the function name for easy referencing
FuncName = 'ShowPlane';

%initailize return variables
err = 1;

%check on the size of input data
if (nargin == 1),
	Opt.Fig = [];
	Opt.OvrWrite = 0;
	Opt.WriteIV = '';
end

if (~isfield(Opt,'OvrWrite') | isempty(Opt.OvrWrite)), Opt.OvrWrite = 0; end

if (size(EqPlane,2) ~= 4),	err = ErrEval(FuncName,'Err_Bad size for EqPlane');	return;	end

Nplanes = size(EqPlane,1);
Nnodes = 4.* Nplanes;

%pop up a figure
if (~isfield(Opt,'Fig') | isempty(Opt.Fig)),
	Opt.Fig = gcf;
end

figure(Opt.Fig);

%get the axis roperties of the figure
Xlim = get(gca,'Xlim');
Ylim = get(gca,'ylim');
Zlim = get(gca,'Zlim');

inode = 0;
Node = zeros(Nnodes,3);
Pat = zeros(Nplanes,4);
ztmp = zeros(1,4);

for (i=1:1:Nplanes),
	%using the XY limits, find the corrspondign z values
	if (EqPlane(i,3) ~= 0),
		ztmp(1) = (-EqPlane(i,4) - EqPlane(i,1).*Xlim(1) - EqPlane(i,2).*Ylim(1)) ./ EqPlane(i,3);
		ztmp(2) = (-EqPlane(i,4) - EqPlane(i,1).*Xlim(2) - EqPlane(i,2).*Ylim(1)) ./ EqPlane(i,3);
		ztmp(3) = (-EqPlane(i,4) - EqPlane(i,1).*Xlim(2) - EqPlane(i,2).*Ylim(2)) ./ EqPlane(i,3);
		ztmp(4) = (-EqPlane(i,4) - EqPlane(i,1).*Xlim(1) - EqPlane(i,2).*Ylim(2)) ./ EqPlane(i,3);
		%form the four points on the plane
		Node(inode+1,:) = [Xlim(1) Ylim(1) ztmp(1)];
		Node(inode+2,:) = [Xlim(2) Ylim(1) ztmp(2)];
		Node(inode+3,:) = [Xlim(2) Ylim(2) ztmp(3)];
		Node(inode+4,:) = [Xlim(1) Ylim(2) ztmp(4)];

	elseif (EqPlane(i,2) ~= 0),
		ytmp(1) = (-EqPlane(i,4) - EqPlane(i,1).*Xlim(1) - EqPlane(i,3).*Zlim(1)) ./ EqPlane(i,2);
		ytmp(2) = (-EqPlane(i,4) - EqPlane(i,1).*Xlim(2) - EqPlane(i,3).*Zlim(1)) ./ EqPlane(i,2);
		ytmp(3) = (-EqPlane(i,4) - EqPlane(i,1).*Xlim(2) - EqPlane(i,3).*Zlim(2)) ./ EqPlane(i,2);
		ytmp(4) = (-EqPlane(i,4) - EqPlane(i,1).*Xlim(1) - EqPlane(i,3).*Zlim(2)) ./ EqPlane(i,2);
		%form the four points on the plane
		Node(inode+1,:) = [Xlim(1) ytmp(1) Zlim(1)];
		Node(inode+2,:) = [Xlim(2) ytmp(2) Zlim(1)];
		Node(inode+3,:) = [Xlim(2) ytmp(3) Zlim(2)];
		Node(inode+4,:) = [Xlim(1) ytmp(4) Zlim(2)];
	elseif (EqPlane(i,1) ~= 0),
		xtmp(1) = (-EqPlane(i,4) - EqPlane(i,2).*Ylim(1) - EqPlane(i,3).*Zlim(1)) ./ EqPlane(i,1);
		xtmp(2) = (-EqPlane(i,4) - EqPlane(i,2).*Ylim(2) - EqPlane(i,3).*Zlim(1)) ./ EqPlane(i,1);
		xtmp(3) = (-EqPlane(i,4) - EqPlane(i,2).*Ylim(2) - EqPlane(i,3).*Zlim(2)) ./ EqPlane(i,1);
		xtmp(4) = (-EqPlane(i,4) - EqPlane(i,2).*Ylim(1) - EqPlane(i,3).*Zlim(2)) ./ EqPlane(i,1);
		%form the four points on the plane
		Node(inode+1,:) = [xtmp(1) Ylim(1) Zlim(1)];
		Node(inode+2,:) = [xtmp(2) Ylim(2) Zlim(1)];
		Node(inode+3,:) = [xtmp(3) Ylim(2) Zlim(2)];
		Node(inode+4,:) = [xtmp(4) Ylim(1) Zlim(2)];
	end

	%verify that all points are on plane for debugging only
	%sum(EqPlane.*[Node(inode+1,:) 1])
	%sum(EqPlane.*[Node(inode+2,:) 1])
	%sum(EqPlane.*[Node(inode+3,:) 1])
	%sum(EqPlane.*[Node(inode+4,:) 1])

	%form the faceset connection, for the patch
	Pat(i,:) = [inode+1 inode+2 inode+3 inode+4];

	inode = inode + 4;

end

%Now display those patches
vc = 0:1./Nnodes:(1-1./Nnodes);
tcolor = [vc' (1-vc)' vc'];

PatchHandles = patch('vertices',Node,'faces',Pat,...
          'FaceVertexCData',tcolor,'FaceColor','flat');

%write to file ?
if (isfield(Opt,'WriteIV') & ~isempty(Opt.WriteIV)),
	fprintf (1,'Saving patches to iv file %s ...\n',Opt.WriteIV);
   Opt.OptIV.BaseCol = tcolor;
   if (~isfield(Opt,'units') | isempty(Opt.units)),	Opt.units = 'tesscon'; end
	Opt.OptIV.units = Opt.units;
	Opt.OptIV.OvrWrite = Opt.OvrWrite;
	Opt.OptIV.verbose = 0;
   [err] = WriteInv21Surf(Opt.WriteIV,Node,Pat,Opt.OptIV);  %exclude redundant last node
end

err = 0;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  In this script, we perform phase transition analysis
%  of Orthogonal matching pursuit.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
clear all;
clc;
rng('default');
% Create the directory for storing results
[status_code,message,message_id] = mkdir('bin');
target_file_path = 'bin/mmv_phase_transition_noiseless_s_16.mat';
N = 256;
S = 16;
pta = spx.pursuit.PhaseTransitionAnalysis(N);
% pta.NumTrials = 100;
dict_model = @(M, N) spx.dict.simple.gaussian_dict(M, N);
data_model = @(N, K) spx.data.synthetic.SparseSignalGenerator(N, K, S).gaussian;
recovery_solver = @(Phi, K, y) spx.pursuit.joint.CoSaMP(Phi, K).solve(y).Z;
pta.run(dict_model, data_model, recovery_solver);
pta.save_results(target_file_path);


text:
function mask = AngularMask(r, c)
%
%        mask = AngularMask(r, c)
%
%        This function creates a mask for a Angular/Spherical map
%
%        Input:
%           -r: rows of the image in the Angular/Spherical format
%           -c: columns of the image in the Angular/Spherical format
%        Output:
%           -mask: a mask where the Angular/Spherical is defined
%
%     Copyright (C) 2011-15  Francesco Banterle
%
%     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/>.
%

[X,Y] = meshgrid(1:c, 1:r);
X = X / c * 2 - 1;
Y = Y / r * 2 - 1;
R = sqrt(X.^2 + Y.^2);

tmpMask = ones(r, c);
tmpMask(R > 1) = 0;

mask = zeros(r, c, 3);

for i=1:3
    mask(:,:,i) = tmpMask;
end

end

author: banterle


text:
function [ y, next ] = p38_equil ( neqn, next )

%*****************************************************************************80
%
%% P38_EQUIL returns equilibrium solutions of problem p38.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, integer NEXT, the index of the previous
%    equilibrium, which should be 0 on first call.
%
%    Output, real Y(NEQN), the "next" equilibrium solution, if any.
%
%    Output, integer NEXT, the index of the current equilibrium,
%    or 0 if there are no more.
%
  if ( next == 0 )
    next = 1;
    y(1) = 0.0;
  elseif ( next == 1 )
    next = 2;
    y(1) = 1.0;
  else
    next = 0;
    y = [];
  end

  return
end


text:
function displayIm = shownsct( y )
% SHOWNSSC   Show nonsubsampled Contourlet transform coefficients.
%
%       shownsct(y)
% Input:
%	y:	a cell vector of length n+1, one for each layer of
%		subband images from NSCT, y{1} is the lowpass image
%
% NOTE:
%       It need further improvement later!!!!

% History:
%   08/08/2003  Created by Jianping Zhou.

% Level of decomposition.
clevels = length( y ) ;

% Show the subband images.
for i=1:clevels
    figure;
    if iscell( y{i} )
        % The number of directional subbands.
        csubband = length( y{i} ) ;
        if csubband > 7
            col = 4 ;
        else
            col = 2 ;
        end
        row = csubband / col ;
        for j = 1:csubband
            subplot( row, col, j ) ;
            imshow( uint8(y{i}{j}) );
            title( sprintf('NSSC coefficients: level %d', i) );
        end
    else
        imshow ( uint8(y{i}) ) ;
        title( sprintf('Nonsubsampled Contourlet coefficients level %d', i) );
    end
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
%% Extentions
% This script shows some of the extensions I have created for accessing
% (and writing) midi files.

%% Tracks
% The standard note matrix from the
% <https://www.jyu.fi/hum/laitokset/musiikki/en/research/coe/materials/miditoolbox/ Midi Toolbox>.
% does not have track information.  My extensions add an eighth column to
% the note matrix with track information.

% Read in a midi file an get the track information
nm = readmidi_java('202.kar',true);
fprintf('There are %d tracks with notes in this midi file.\n', ...
    length(unique(nm(:,8))));
% And write out tracks ...
writemidi_java(nm,'track_demo.mid');
% If we read this newly created midi file back in again, we should get the
% tracks...
nm_new = readmidi_java('track_demo.mid',true);
figure(1);
clf;
plot(nm(:,8) == nm_new(:,8),'b*');
title('Note track numbers check');
xlabel('note number');
ylim([-0.2 1.2]);
set(gca,'YTick',[0,1])
set(gca,'YTickLabel',{'different','same'})

%% Lyrics
% Some midi files have lyrics embedded in them.  If the file extension is
% .kar rather than .mid or .midi, it will almost certainly have lyrics.

% Get the lyrics in text form.  get_lyrics_as_bytes() will give you the raw
% bytes, which may be easier to use for more complex processing.  Some of
% the original lyric commands in the midi file may have unprintable
% characters.  These unprintable characters are converted into hex strings
% by get_lyrics_as_text().
[lyrics info] = get_lyrics_as_text('202.kar');

% Now go ahead and print out the first 20 lyric commands and their times in
% seconds
fprintf('********* FIRST 20 RAW LYRICS ********* \n\n');
for i=1:20
   fprintf('%g: |%s|\n',info(i,2),lyrics{i});
end
fprintf('\n\n');

% Now try to do something a little smarter so that we can see the lyrics
% as lines.  In this particular midi file, 0xd indicates and end of line
% and 0xa indicates a new stanza.
fprintf('********* LYRICS IN STANZAS ********* \n\n');
for i=1:length(lyrics)
    if strcmp(lyrics{i},'[0xd]')
        fprintf('\n');
    elseif strcmp(lyrics{i},'[0xa]')
        fprintf('\n');
    else
        fprintf('%s',lyrics{i});
    end
end

%% Instrumentation
% Program commands change the instrumentation for individual channels in
% tracks.  Each <http://www.midi.org/techspecs/gm1sound.php patch number>
% corresponds to a particular instrument.  My extensions allow you to read
% the patch commands out of an existing midi file and they allow you to
% specify a single patch per channel per track when writing midi files.
% You can listen to the original midi file - with the extension changed to
% .mid - <202.mid here> and the midi file with trumpet vocals
% <vocals_with_trumpet.mid here>.  The vocals start about 6 seconds in.

% first get the note matrix with track information
nm = readmidi_java('202.kar',true);

% now get the program changes
changes = get_program_changes('202.kar');

% let's find out what the patch was for the vocal channel
vocal_patch = changes(changes(:,5) == 3,:);

% it's a piccolo, number 73
fprintf('Vocal patch #: %d\n', vocal_patch(1));

% now let's create a new .mid file where the vocals are on a trumpet (57)
writemidi_seconds(nm,'vocals_with_trumpet.mid',[3 vocal_patch(4) 57]);

%% Tempos
% Tempo commands change the lengths of beats.  The note matrix returned by
% readmidi_java takes these tempo changes into account.

tempos = get_tempos('202.kar');
% There are quite a lot of these changes
figure(1);
clf;
plot(tempos(:,3), tempos(:,1),'*-');
title('Tempo changes in 202.kar');
xlabel('Time (s)');
ylabel('Beats per minute');
ylim([0 max(ylim)]);

% 2010-05-03 Christine Smit csmit@ee.columbia.edu

author: Sable

repo: mcbench-benchmarks


text:
function [result_tab, error_code, nb_warm_start, z]=SparseLP_old(c,polyhedron, zeps, max_warm_start)

% Compute the sparsest solution of a linear program (LP)
% min c'x
% s.t.  Sx=b
%       lb <= x <= ub
%

 n=size(c,1);error_code=0;

% Compute the optimal value of the LP
% [x0, err_init]=LPSolve(c,Aeq,beq,Ain,bin,lb,ub);
 [x0, err_init]=LPSolve(c,polyhedron);
 err_init
if err_init~=1
 init_result=length(find(abs(x0)>zeps));
 optvalue=x0'*c;

% Heuristic to get 0 if optimal value is 0.
 if optvalue==0 && max(polyhedron.Aeq*zeros(n,1)-polyhedron.beq)==0 && max(polyhedron.Ain*zeros(n,1)-polyhedron.bin)<=0 && max(polyhedron.lb-zeros(n,1))<=0 && max(zeros(n,1)-polyhedron.ub)<=0
	z=zeros(n,1);result_tab=[0]; error_code=0; nb_warm_start=0;
 else

% Generate final polyhedron with optimal value as a constraint.
 polyhedron=struct('Aeq',[polyhedron.Aeq;c'],'beq',[polyhedron.beq;optvalue],'Ain',polyhedron.Ain,'bin',polyhedron.bin,'lb',polyhedron.lb,'ub',polyhedron.ub);

%#################################################################################
% MAIN : Algorithm call with warm start if we didn't improve l1-norm solution or if we have precision problem

 nb_warm_start=-1;tic;
 [x, err_l1]=L1Solve(polyhedron); l1_result=length(find(abs(x)>zeps));
 err_l1
 if err_l1~=1
	result_tab=[init_result, l1_result];
 else
	result_tab=[init_result];
 end

% Parameters
 nb_itemax_r=50; nb_itemax_sla=2; update_param=0.5;rstop=sqrt(eps); move=1;

 while err_l1~=1 && nb_warm_start<max_warm_start && init_result>0 && ( nb_warm_start==-1 || nnz_z<=min(result_tab)) && move

  [z, nnz_z, error_code] = ThetaL0(polyhedron, @theta1, @dtheta1, nb_itemax_r, rstop, update_param, nb_itemax_sla, zeps, x);
% Update
	nb_warm_start=nb_warm_start+1;
	if nnz_z<min(result_tab)
		result_tab=[result_tab, nnz_z];x=z;
	elseif z==x
		move=0;
	end

 end
 if err_l1==1
	z=x0;nnz_z=-1;error_code=1;
 end
%################################################################################
 toc
 end
else
 result_tab=[]; error_code=1;nb_warm_start=0;z=x0;
end
end

%################################################################################
%# Declaration of functions thetas and their derivatives :

function y = theta1(t,r)
	y=abs(t)./(abs(t)+r);
end

function y = dtheta1(t,r)
	y=sign(t).*r./((abs(t)+r).^2);
end

function y = theta2(t,r)
	y=1-exp(-abs(t)/r);
end

function y = dtheta2(t,r)
	y=sign(t).*exp(-abs(t)./r)./r;
end

%######################################################################################################################
%######################################################################################################################
function [x, nnz_x, error_code] = ThetaL0(polyhedron, theta_func, dtheta_func, nb_itemax_r, rstop, update_param, nb_itemax_sla, zero_precision, init_point)

% Citation :  A Smoothing Method for Sparse Optimization over Polyhedral Sets, Modelling, Computation and Optimization in Information Systems and Management Sciences Advances in Intelligent Systems and Computing Volume 359, 2015, pp 369-379   Tangi Migot, Mounir Haddou, 2015

% error_code :
% -2 : maximum iteration in r
% -1 : maximum iteration in SLA
% 0 : success
% 1 : unbounded
% 2 : precision

 error_code=0;
% code d'erreur SLA : (init in 1 : in case no SLA are done)
 err_code=1;

% for debug :
minx=[];rtab=[];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialization

% if no initial point we take the solution of l1-norm problem
 if isempty(init_point)
	init_point=L1Solve(polyhedron);
 end

 sol=init_point;
 nnz_current=length(find(abs(sol)>zero_precision));nnzero=[];soltab=[];minx=[min(sol(sol>0))];
% initialization of parameter r
 i=1;r0=1;rold=r0;r=rold;rtab=[r0];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%  Main Loop

 while i<nb_itemax_r && theta_func(min(sol(sol>zero_precision)),rold)<nnz_current/(nnz_current+1) && r>rstop
% Scaling to put variables in [0,1] box
	scal=max(abs(sol))/2;z=sol/scal;

	spolyhedron=struct('Aeq',polyhedron.Aeq,'beq',polyhedron.beq/scal,'Ain',polyhedron.Ain,'bin',polyhedron.bin/scal,'lb',polyhedron.lb/scal,'ub',polyhedron.ub/scal);
% Concave program
	[z, err_code] = ThetaSLA(spolyhedron, z,r, nb_itemax_sla, theta_func, dtheta_func);

	    if err_code==1
	     error_code=1;i=nb_itemax_r; % unbounded
	    else
	     	sol=z*scal;soltab=[soltab;sol'];minx=[minx;min(sol(sol>0))];rtab=[rtab,r];
	     	nnzk=length(find(abs(sol)>zero_precision));nnzero=[nnzero, nnzk];
% Update
	     	i=i+1;rold=r;r=r*update_param;
	    end
	    if err_code==-1
	     error_code=-1; % max iteration SLA
	    end

 end
%%%%% End Main Loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Best solution find :
if err_code ~=1
 nnz_x=min(nnzero);index=find(nnzero==nnz_x);x=soltab(index(1),:)';
else
 nnz_x=length(find(abs(sol)>zero_precision));x=init_point;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Error code :
 if i==nb_itemax_r && error_code~=1
	error_code=-2; % max iteration in r
 end

 if ( (~isempty(polyhedron.Aeq) && norm(polyhedron.Aeq*x-polyhedron.beq,Inf)>sqrt(eps)) || (~isempty(polyhedron.Ain) && max(polyhedron.Ain*x-polyhedron.bin)>sqrt(eps)) || (~isempty(polyhedron.lb) && max(polyhedron.lb-x)>sqrt(eps)) || (~isempty(polyhedron.ub) && max(x-polyhedron.ub)>sqrt(eps))) && error_code~=1
	error_code=2; % precision
 end

end
%######################################################################################################################
%
% Solve concave program of minimization Theta with Successiv Linearization Algorithm
%
%######################################################################################################################
function [x, error_code] =ThetaSLA(polyhedron,init_point,r_param,ite_max, theta_func, dtheta_func)

% err_code :
% -1 : maximum iteration
%  0 : success
%  1 : unbounded linear program
%  2 : precision
 error_code=0;

 j=1;y=init_point;zero_precision=sqrt(eps);

 while j<ite_max && (j==1 || abs(c'*yold-c'*x)>zero_precision)
	% Objectiv theta function with scaling :

%	c=dtheta_func(y,r_param)/(theta_func(1,r_param))
	c=dtheta_func(y,r_param)*r_param;

	[x,bool_LP]=LPSolve(c,polyhedron);
	% update :
	yold=y;y=x;j=j+1;

	if bool_LP==1
		error_code=1;j=ite_max;
	end
 end

 if error_code~=1
	 if j==ite_max
		error_code=-1;
	 end

	 if (~isempty(polyhedron.Aeq) && norm(polyhedron.Aeq*x-polyhedron.beq,Inf)>zero_precision) || (~isempty(polyhedron.Ain) && max(polyhedron.Ain*x-polyhedron.bin)>zero_precision) || (~isempty(polyhedron.lb) &&max(polyhedron.lb-x)>sqrt(eps)) || (~isempty(polyhedron.ub) &&max(x-polyhedron.ub)>sqrt(eps))
		error_code=2;
	 end
 end
end
%######################################################################################################################
%
% Solve convex program of minimization l1-norm
%
%######################################################################################################################
function [x,error_code]=L1Solve(polyhedron)
%{
% here we use SeDuMi from CVX Matlab package
 n=size(polyhedron.lb,1);error_code=0;
% min c'x
% s.t.  Aeq x == beq
%       Ain x <= bin
%       lb <= x <= ub

 cvx_begin quiet
     cvx_precision best
%     cvx_solver sdpt3
     variable x(n)
     minimize(norm(x,1))
     subject to
         polyhedron.Aeq*x == polyhedron.beq;
         polyhedron.Ain*x <= polyhedron.bin;
         x >= polyhedron.lb;
         x <= polyhedron.ub;
 cvx_end

% error 1 : unbounded
 if cvx_optval<=-Inf || cvx_optval>=Inf
	error_code=1;
 end

% error 2 : precision
 if norm(polyhedron.Aeq*x-polyhedron.beq,Inf)>sqrt(eps) || max(polyhedron.Ain*x-polyhedron.bin)>sqrt(eps) || max(polyhedron.lb-x)>sqrt(eps) || max(x-polyhedron.ub)>sqrt(eps)
	error_code=2;
 end
%}
%######################################################################################################################
% Solveur GUROBI

% solver with Gurobi (reformulation comme un LP)
if isempty(polyhedron.Aeq)
	%On a pas besoin des égalités
	AEQ=[];
else
	m1=size(polyhedron.beq,1);n=size(polyhedron.Aeq,2);AEQ=[polyhedron.Aeq, zeros(m1,n)];
end
if isempty(polyhedron.Ain)
	n=size(polyhedron.lb,1);
	AIN=[eye(n), -eye(n); -eye(n), -eye(n)];
	BIN=[zeros(n,1); zeros(n,1)];
else
	m2=size(polyhedron.bin,1);n=size(polyhedron.Ain,2);
	AIN=[polyhedron.Ain, zeros(m2,n); eye(n), -eye(n); -eye(n), -eye(n)];
	BIN=[polyhedron.bin;zeros(n,1); zeros(n,1)];
end
if isempty(polyhedron.lb)
	LB=[];
else
	n=size(polyhedron.lb,1);LB=[polyhedron.lb; zeros(n,1)];
end
if isempty(polyhedron.ub)
	UB=[];
else
	n=size(polyhedron.ub,1);temp=max(abs(polyhedron.lb),abs(polyhedron.ub));UB=[polyhedron.ub; temp];
end

polyhedron_l1=struct('Aeq',AEQ,'beq',polyhedron.beq,'Ain',AIN,'bin',BIN,'lb',LB,'ub',UB);
c=[zeros(n,1);ones(n,1)];

[x,error_code]=LPSolve(c,polyhedron_l1);
x=x(1:n);

end
%######################################################################################################################
%
% Solve linear Program
%
%######################################################################################################################
function [x,error_code]=LPSolve(c,polyhedron)
%{
% Solveur SEDUMI
 n=size(c,1);error_code=0;
% min c'x
% s.t.  Aeq x == beq
%       Ain x <= bin
%       lb <= x <= ub

 cvx_begin quiet
    cvx_precision best
    cvx_solver sedumi
    variable x(n)
    minimize(c'*x)
    subject to
	polyhedron.Aeq*x == polyhedron.beq;
	polyhedron.Ain*x <= polyhedron.bin;
	x >= polyhedron.lb;
	x <= polyhedron.ub;
 cvx_end

% error 1 : unbounded
 if cvx_optval<=-Inf || cvx_optval>=Inf
	error_code=1;
 end

% error 2 : precision
 if norm(polyhedron.Aeq*x-polyhedron.beq,Inf)>sqrt(eps) || max(polyhedron.Ain*x-polyhedron.bin)>sqrt(eps) || max(polyhedron.lb-x)>sqrt(eps) || max(x-polyhedron.ub)>sqrt(eps)
	error_code=2;
 end
%}
%######################################################################################################################
% Solveur GUROBI

[x, fval, exitflag] = linprog_gurobi(c, polyhedron.Ain, polyhedron.bin, polyhedron.Aeq, polyhedron.beq, polyhedron.lb, polyhedron.ub);

if exitflag==1 || exitflag==0
 error_code=0;
elseif exitflag==-2 || exitflag==-3
 error_code=1;
end

% error 2 : precision
 if (~isempty(polyhedron.Aeq) && norm(polyhedron.Aeq*x-polyhedron.beq,Inf)>sqrt(eps)) || (~isempty(polyhedron.Ain) && max(polyhedron.Ain*x-polyhedron.bin)>sqrt(eps)) || (~isempty(polyhedron.lb) &&max(polyhedron.lb-x)>sqrt(eps)) || (~isempty(polyhedron.ub) &&max(x-polyhedron.ub)>sqrt(eps))
	error_code=2;
 end

end
function [x, fval, exitflag] = linprog_gurobi(f, A, b, Aeq, beq, lb, ub)
%LINPROG A linear programming example using the Gurobi MATLAB interface
%
%   This example is based on the linprog interface defined in the
%   MATLAB Optimization Toolbox. The Optimization Toolbox
%   is a registered trademark of The MathWorks, Inc.
%
%   x = LINPROG(f,A,b) solves the linear programming problem:
%
%   minimize     f'*x
%   subject to   A*x <= b
%
%
%   x = LINPROG(f,A,b,Aeq,beq) solves the problem:
%
%   minimize     f'*x
%   subject to     A*x <= b,
%                Aeq*x == beq.
%
%   x = LINPROG(f,A,b,Aeq,beq,lb,ub) solves the problem:
%
%   minimize     f'*x
%   subject to     A*x <= b,
%                Aeq*x == beq,
%          lb <=     x <= ub.
%
%   You can set lb(j) = -inf, if x(j) has no lower bound,
%   and ub(j) = inf, if x(j) has no upper bound.
%
%   [x, fval] = LINPROG(f, A, b) returns the objective value
%   at the solution. That is, fval = f'*x.
%
%   [x, fval, exitflag] = LINPROG(f, A, b) returns an exitflag
%   containing the status of the optimization. The values for
%   exitflag and corresponding status codes are:
%      1 - OPTIMAL,
%      0 - ITERATION_LIMIT,
%     -2 - INFEASIBLE,
%     -3 - UNBOUNDED.
%

if nargin < 3
    error('linprog(f, A, b)')
end

if nargin > 7
    error('linprog(f, A, b, Aeq, beq, lb, ub)');
end

if ~isempty(A)
    n = size(A, 2);
elseif nargin > 4 && ~isempty(Aeq)
    n = size(Aeq, 2);
else
    error('No linear constraints specified')
end

if ~issparse(A)
    A = sparse(A);
end

if nargin > 3 && ~issparse(Aeq)
    Aeq = sparse(Aeq);
end

model.obj = f;

if nargin < 4
    model.A = A;
    model.rhs = b;
    model.sense = '<';
else
    model.A = [A; Aeq];
    model.rhs = [b; beq];
    model.sense = [repmat('<', size(A,1), 1); repmat('=', size(Aeq,1), 1)];
end

if nargin < 6
    model.lb = -inf(n,1);
else
    model.lb = lb;
end

if nargin == 7
   model.ub = ub;
end

params.outputflag = 0;
result = gurobi(model, params);

if strcmp(result.status, 'OPTIMAL')
    exitflag = 1;
elseif strcmp(result.status, 'ITERATION_LIMIT')
    exitflag = 0;
elseif strcmp(result.status, 'INF_OR_UNBD')
    params.dualreductions = 0;
    result = gurobi(model, params);
    if strcmp(result.status, 'INFEASIBLE')
        exitflag = -2;
    elseif strcmp(result.status, 'UNBOUNDED')
        exitflag = -3;
    else
        exitflag = nan;
    end
elseif strcmp(result.status, 'INFEASIBLE')
    exitflag = -2;
elseif strcmp(result.status, 'UNBOUNDED')
    exitflag = -3;
else
    exitflag = nan;
end

if isfield(result, 'x')
    x = result.x;
else
    x = nan(n,1);
end

if isfield(result, 'objval')
    fval = result.objval;
else
    fval = nan;
end
end

author: opencobra

repo: cobratoolbox


text:
function [mWalkout, fR95, fProb, PHI, R] = calc_Schusterwalk(mCatalog)
% function  [mWalkout, fR95, fProb, PHI, R] = calc_Schusterwalk(mCatalog)
% ------------------------------------------------------------------------
% Calculate random walkout (Schuster's method). See Rydelek & Sacks, Nature
% 1989.
%
% Incoming variables:
% mCatalog : EQ catalog
%
% Outgoing variables:
% mWalkout : Random walkout phases
% fR95     : Critical radius at 95% level
% fProb    : Probability to obtain vector of length >= R
% PHI      : Sum of phase angles
% R        : Possible maximum length of vector
%
% J. Woessner, woessner@seismo.ifg.ethz.ch
% last update: 13.03.03

% Get hours and minutes, calculate minutes
vH = mCatalog(:,8);
vMin = mCatalog(:,9);
vTime = vH*60+vMin;

% Calculate degrees and radians
vThetaDegree = vTime*180/720;
vThetaRad = vTime*pi/720;

mWalkout(1,1)=0;
mWalkout(1,2)=0;
for i=1:length(mCatalog(:,8))
    mWalkout(i+1,1)=mWalkout(i,1)+sin(vThetaRad(i));
    mWalkout(i+1,2)=mWalkout(i,2)+cos(vThetaRad(i));
end

A=sum(sin(vThetaRad));
B=sum(cos(vThetaRad));
R=sqrt(A^2+B^2);
PHI=atan(B/A);
% length(mCatalog(:,8))
% if mWalkout(length(mCatalog(:,8)),2)<0
%     PHI=PHI+pi;
% end

N=length(mCatalog(:,8));
% Critical radius at 95% confidence level
fR95=1.73*sqrt(N); % Rydelek & Sacks, Nature 1989
% Probability to obtain vector of length >= R from a random walkout
fProb=exp(-R^2/N);

author: CelsoReyes


text:
function handles = setSimGraphEdits(hObject, handles)
% setSimGraphEdits Enables or disables the text edits depending on the
% current similarity graph type

switch handles.SimGraphType
    case 1
        enableNeighbors = 'off';
        enableEps = 'off';
        enableSigma = 'on';
    case {2,3}
        enableNeighbors = 'on';
        enableEps = 'off';
        enableSigma = 'on';
    case 4
        enableNeighbors = 'off';
        enableEps = 'on';
        enableSigma = 'off';
end

set(handles.edtSimGraphNeighbors, 'Enable', enableNeighbors);
set(handles.edtSimGraphEps, 'Enable', enableEps);
set(handles.edtSimGraphSigma, 'Enable', enableSigma);

guidata(hObject, handles);

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval