text:
function [F,L,KL] = spm_vb_Fn(Y,block)
% Compute voxel-wise contributions to model evidence
% FORMAT [F,L,KL] = spm_vb_Fn(Y,block)
%
% Y          - [T x N] time series
% block      - data structure (see spm_vb_glmar)
%
% F          - [N x 1] vector where nth entry is unique contribution to
%              model evidence from voxel n
% L          - [N x 1] Average Likelihood
% KL.w       - [N x 1] KL w - unique contribution
% KL.a       - [N x 1] KL a - unique contribution
% KL.lam     - [N x 1] KL Lambda
% KL.alpha   - Scalar
% KL.beta    - Scalar
%__________________________________________________________________________
% Copyright (C) 2005-2014 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_vb_Fn.m 6079 2014-06-30 18:25:37Z spm $

T  = block.T;
p  = block.p;
k  = block.k;
N  = block.N;
X  = block.X;

C2 = block.C2;

Bk = kron(diag(block.mean_alpha),block.Dw);
B  = block.Hw*Bk*block.Hw';

if p>0
    if ~strcmp(block.priors.A,'Discrete')
        Jk = kron(diag(block.mean_beta),block.Da);
        J  = block.Ha*Jk*block.Ha';
    end
end

%-Get KL-alpha
%--------------------------------------------------------------------------
KL.alpha = 0;
for j=1:k
    KL.alpha = KL.alpha + spm_kl_gamma(block.b_alpha(j),block.c_alpha(j),block.b_alpha_prior(j),block.c_alpha_prior(j));
end

% Get KL-beta
%--------------------------------------------------------------------------
KL.beta = 0;
if p > 0
    if strcmp(block.priors.A,'Discrete')
        for j=1:p
            for s=1:block.priors.S
                KL.beta = KL.beta + spm_kl_gamma(block.b_beta(j,s),block.c_beta(j,s),block.b_beta_prior(j,s),block.c_beta_prior(j,s));
            end
        end
    else
        for j = 1:p
            KL.beta = KL.beta + spm_kl_gamma(block.b_beta(j),block.c_beta(j),block.b_beta_prior(j),block.c_beta_prior(j));
        end
    end
end

% Get average Likelihood, KL-Lambda and terms for KL-w, KL-a
%--------------------------------------------------------------------------
for n=1:N
    subblock_n = [(n-1)*k+1:n*k];

    if p>0
        Gn = spm_vb_get_Gn(Y,block,n);
    else
        wc = block.w_cov{n};
        en = (Y(:,n)-X*block.w_mean(subblock_n,1));
        Gn = trace(wc*block.XTX)+en'*en;
    end

    L(n)   = -0.5*block.mean_lambda(n)*Gn;
    L(n)   = L(n) + 0.5*(T-p)*(psi(block.c_lambda(n)) + log(block.b_lambda(n)));
    L(n)   = L(n)-0.5*block.C2/N;

    KL.lam(n) = spm_kl_gamma(block.b_lambda(n),block.c_lambda(n),block.b_lambda_prior(n),block.c_lambda_prior(n));

    KL_w1  = -0.5*sum(log(block.mean_alpha))-0.5*log(det(block.w_cov{n}));
    KL_w1  = KL_w1-0.5*k*block.log_det_Dw/N;
    KL_w2  = 0.5*trace(B(subblock_n,subblock_n)*block.w_cov{n});

    subblock_ni = [1:N*k];
    subblock_ni(subblock_n) = [];
    Bnn         = B(subblock_n,subblock_n);
    Bni         = B(subblock_n,subblock_ni);
    Bin         = B(subblock_ni,subblock_n);

    w_mean  = block.w_mean(subblock_n,1);
    KL_w3   = 0.5*w_mean'*Bnn*w_mean+0.5*w_mean'*Bni*block.w_mean(subblock_ni,1);
    KL_w3   = KL_w3-0.5*k;
    KL.w(n) = KL_w1+KL_w2+KL_w3;

    if p>0
        asubblock_n = [(n-1)*p+1:n*p];
        a_mean      = block.a_mean(asubblock_n,1);
        if strcmp(block.priors.A,'Discrete')
            ibeta_n = diag(block.priors.gamma(n,:)*(1./block.mean_beta'));
            a_n     = block.priors.gamma(n,:)*block.as';
            KL.a(n) = spm_kl_normal(a_mean,block.a_cov{n},a_n,ibeta_n);
%             a_mean
%             sqrt(block.a_cov{n})
%             a_n
%             sqrt(ibeta_n)
        else
            asubblock_ni = [1:N*p];
            asubblock_ni(asubblock_n) = [];
            Jnn = J(asubblock_n,asubblock_n);
            Jni = J(asubblock_n,asubblock_ni);
            KL_a1 = -0.5*sum(log(block.mean_beta)) - 0.5*log(det(block.a_cov{n}));
            KL_a1 = KL_a1-0.5*p*block.log_det_Da/N;
            KL_a2 = 0.5*trace(Jnn*block.a_cov{n});
            KL_a3 = 0.5*a_mean'*Jnn*a_mean+0.5*a_mean'*Jni*block.a_mean(asubblock_ni,1);
            KL_a3 = KL_a3-0.5*p;
            KL.a(n) = KL_a1 + KL_a2 + KL_a3;
        end
    else
        KL.a(n) = 0;
    end
end

F = L -KL.w  -KL.lam -KL.a -KL.alpha/N -KL.beta/N;


text:
 function y = reshapee(x, varargin)
%|function y = reshapee(x, varargin)
%|
%| reshape function that allows possibly one null argument, and all
%| other arguments can be vectors, unlike matlab that requires scalars.
%| example: reshape(rand(2*3*5,7), [2 3], [], 7) will become [2 3 5 7]
%|
%| in
%|	x		[(*dim)]
%|	varargin	dimensions
%|
%| out
%|	y		[dim]
%|
%| Copyright 2004-8-22, Jeff Fessler, University of Michigan

if nargin == 1 && streq(x, 'test'), reshapee_test, return, end
if nargin < 2, ir_usage, end

dim_i = size(x);

ndim = 0;
edim = [];
dim_o = [];
for ii=1:length(varargin)
	arg = varargin{ii};
	ndim = ndim + length(arg);
	if isempty(arg)
		if ~isempty(edim)
			fail 'only one empty dim allowed'
		end
		edim = 1+length(dim_o);
		dim_o = [dim_o, 1]; % trick: place holder
	else
		dim_o = [dim_o, arg];
	end
end

if ~isempty(edim) % fill in empty dim if present
	if prod(dim_o) <= 0, fail('nonpositive dims?'), end
	dim_o(edim) = prod(dim_i) / prod(dim_o);
	if round(dim_o(edim)) ~= dim_o(edim)
		pr dim_i
		pr dim_o
		fail('bad dim')
	end
end

y = reshape(x, dim_o);

% self test
function reshapee_test
dim = [2 3 5 7];
x = 1:prod(dim);
y = reshapee(x, dim(1:2), [], dim(4));
z = reshape(x, dim);
jf_equal(y,z)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [ds_all, bs_all, targets] = bboxpred_data(name)
% Collect training data for bounding box prediction.
%   [ds, bs, targets] = bboxpred_data(name)
%
% Return values
%   ds_all    Predicted bounding boxes (clipped to the image)
%             One cell percomponent
%   bs_all    All filter bounding boxes (unclipped)
%             One cell percomponent
%   targets   Ground-truth bounding boxes (clipped)
%             One cell percomponent
%
% Argument
%   name      Object class

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2011-2012 Ross Girshick
% Copyright (C) 2008, 2009, 2010 Pedro Felzenszwalb, Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

conf = voc_config();

try
  load([conf.paths.model_dir name '_bboxdata']);
catch
  % load final model for class
  load([conf.paths.model_dir name '_final']);
  % get training data
  pos = pascal_data(model.class, model.year);

  numpos = length(pos);
  model.interval = conf.training.interval_fg;
  pixels = model.minsize * model.sbin / 2;
  minsize = prod(pixels);
  nrules = length(model.rules{model.start});
  parb = cell(1,numpos);
  part = cell(1,numpos);

  % compute latent filter locations and record target bounding boxes
  parfor i = 1:numpos
    pard{i} = cell(1,nrules);
    parb{i} = cell(1,nrules);
    part{i} = cell(1,nrules);
    fprintf('%s %s: bboxdata: %d/%d\n', procid(), name, i, numpos);
    bbox = pos(i).boxes;
    % skip small examples
    if (bbox(3)-bbox(1)+1)*(bbox(4)-bbox(2)+1) < minsize
      continue;
    end
    % get example
    im = imreadx(pos(i));
    [im, bbox] = croppos(im, bbox);
    [pyra, model_dp] = gdetect_pos_prepare(im, model, bbox, 0.7);
    [ds, bs] = gdetect_pos(pyra, model_dp, 1, ...
                            1, 0.7, [], 0.5);
    if ~isempty(ds)
      % component index
      c = ds(1,end-1);
      bs = reduceboxes(model, bs);
      ds = clipboxes(im, ds);
      pard{i}{c} = [pard{i}{c}; ds(:,1:end-2)];
      parb{i}{c} = [parb{i}{c}; bs(:,1:end-2)];
      part{i}{c} = [part{i}{c}; bbox];
    end
  end
  ds_all = cell(1,nrules);
  bs_all = cell(1,nrules);
  targets = cell(1,nrules);
  for i = 1:numpos
    for c = 1:nrules
      ds_all{c} = [ds_all{c}; pard{i}{c}];
      bs_all{c} = [bs_all{c}; parb{i}{c}];
      targets{c} = [targets{c}; part{i}{c}];
    end
  end
  save([conf.paths.model_dir name '_bboxdata'], ...
       'ds_all', 'bs_all', 'targets');
end

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
classdef TestImreadmulti
    %TestImreadmulti

    methods (Static)
        function test_read_pages
            filename = which('mri.tif');
            if isempty(filename)
                error('mexopencv:testskip', 'missing data');
            end

            imgs = cv.imreadmulti(filename);
            validateattributes(imgs, {'cell'}, {'vector', 'nonempty'});
            cellfun(@(im) validateattributes(im, {'uint8'}, ...
                {'nonempty', 'ndims',2}), imgs);
        end

        function test_compare
            % TIFF file with 27 grayscale images
            filename = which('mri.tif');
            if isempty(filename)
                error('mexopencv:testskip', 'missing data');
            end

            % number of pages
            info = imfinfo(filename);
            num_pages = numel(info);

            % read all pages
            imgs1 = cv.imreadmulti(filename, 'AnyColor',true, 'AnyDepth',true);
            validateattributes(imgs1, {'cell'}, {'vector', 'numel',num_pages});

            % compare against imread from MATLAB
            imgs2 = cell(1,num_pages);
            for i=1:num_pages
                % (convert indexed images to grayscale)
                [X,map] = imread(filename, 'Info',info, 'Index',i);
                imgs2{i} = ind2gray(X, map);
            end
            assert(isequal(imgs1(:), imgs2(:)), 'Images are not equal');
        end

        function test_single_page
            filename = fullfile(mexopencv.root(),'test','fruits.jpg');
            imgs = cv.imreadmulti(filename);
            validateattributes(imgs, {'cell'}, {'vector', 'numel',1});
        end

        function test_error_non_existant_file
            try
                cv.imreadmulti('non_existant_image.tif');
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
        function test_error_argnum
            try
                cv.imreadmulti();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function X = circle_graph_layout(G,radius)
% CIRCLE_LAYOUT Layout the vertices of a graph on a circle
%
% X = circle_layout(G) generates a layout of graph with vertices uniformly
% placed on the circle.  This function does no interesting layout and just
% places the vertices around the circle in order.
%
% X = circle_layout(G,radius) places the vertices with a radius other than
% 1.0.
%
% Example:
%   G = cycle_graph(6);
%   X = circle_graph_layout(G);
%   gplot(G,X);

% David F. Gleich
% Copyright, Stanford University, 2008

%% History
%  2008-09-25: Initial coding
%%

if ~exist('radius','var') || isempty(radius), radius = 1.0; end;
pi = 3.14159;
n = num_vertices(G);
X = zeros(n,2);
X(:,1) = radius*cos( (0:n-1)'*2*pi/n );
X(:,2) = radius*sin( (0:n-1)'*2*pi/n );

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [hpe] = TW2hpe(TW)
% Convert power from terawatts to electric horsepower.
% Chad A. Greene 2012
hpe = TW*1340482573.7265;

author: Sable

repo: mcbench-benchmarks


text:
function [translatedValue] = mrViewGetDictionary(valueToTranslate)
%
%
% USAGE:
%   This function is called to translate a parameter value to the internal
%   key assigned to that value. This internal key will be used for the
%   other dictionaries containing Help and struct location information.
%
% INPUT:
%   valueToTranslate - The parameter passed in to this dictionary
%
% OUTPUT:
%   translatedValue - The internal key used by all of the
%   viewGetDictionaries
%
% EXAMPLE:
%   a = mrViewGetDictionary('Anatomy Data Size')
%   a =
%       anatdatasize
%
%
GLOBAL viewGetDict;

if (isempty(viewGetDict))
    %If this is the first run and viewGetDict hasn't been initialized,
    %let's do that.
    viewGetDict = containers.Map;
    viewGetDict = mrViewGetDictionarySetup(viewGetDict); %Initializing it
end

if viewGetDict.isKey(valueToTranslate)
    translatedValue = viewGetDict(valueToTranslate);
else
    translateValue = ''; %Return empty
    warning('vista:ViewGetDict','The value supplied to mrViewGetDictionary is invalid. Please try again.';
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = frame2oell(varargin)
% VL_FRAMES2OELL   Convert a geometric frame to an oriented ellipse
%   EFRAME = VL_FRAME2OELL(FRAME) converts the generic FRAME to an
%   oriented ellipses EFRAME. FRAME and EFRAME can be matrices, with
%   one frame per column.
%
%   A frame is either a point, a disc, an oriented disc, an ellipse,
%   or an oriented ellipse. These are represented respectively by 2,
%   3, 4, 5 and 6 parameters each, as described in VL_PLOTFRAME().  An
%   oriented ellipse is the most general geometric frame; hence, there
%   is no loss of information in this conversion.
%
%   If FRAME is an oriented disc or ellipse, then the conversion is
%   immediate. If, however, FRAME is not oriented (it is either a
%   point or an unoriented disc or ellipse), then an orientation must
%   be assigned. The orientation is chosen in such a way that the
%   affine transformation that maps the standard oriented frame into
%   the output EFRAME does not rotate the Y axis. If frames represent
%   detected visual features, this convention corresponds to assume
%   that features are upright.
%
%   If FRAME is a point, then the output is an ellipse with null area.
%
%   See: feature frames,
%   VL_PLOTFRAME(), VL_HELP().
[varargout{1:nargout}] = vl_frame2oell(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function [EigenvectorsDiscrete,EigenVectors]=discretisation(EigenVectors)
%
% EigenvectorsDiscrete=discretisation(EigenVectors)
%
% Input: EigenVectors = continuous Ncut vector, size = ndata x nbEigenvectors
% Output EigenvectorsDiscrete = discrete Ncut vector, size = ndata x nbEigenvectors
%
% Timothee Cour, Stella Yu, Jianbo Shi, 2004

[n,k]=size(EigenVectors);

vm = sqrt(sum(EigenVectors.*EigenVectors,2));
EigenVectors = EigenVectors./repmat(vm,1,k);

R=zeros(k);
R(:,1)=EigenVectors(1+round(rand(1)*(n-1)),:)';
c=zeros(n,1);
for j=2:k
    c=c+abs(EigenVectors*R(:,j-1));
    [minimum,i]=min(c);
    R(:,j)=EigenVectors(i,:)';
end

lastObjectiveValue=0;
exitLoop=0;
nbIterationsDiscretisation = 0;
nbIterationsDiscretisationMax = 20;%voir
while exitLoop== 0
    nbIterationsDiscretisation = nbIterationsDiscretisation + 1 ;
    EigenvectorsDiscrete = discretisationEigenVectorData(EigenVectors*R);
    [U,S,V] = svd(EigenvectorsDiscrete'*EigenVectors,0);
    NcutValue=2*(n-trace(S));

    if abs(NcutValue-lastObjectiveValue) < eps | nbIterationsDiscretisation > nbIterationsDiscretisationMax
        exitLoop=1;
    else
        lastObjectiveValue = NcutValue;
        R=V*U';
    end
end

author: jwyang


text:
% STRINGS
%
% Files
%   strendswith   - Determines whether a string ends with a specified pattern
%   strjoin       - Creates a string by joining multiple terms
%   strsplit      - Splits a string into multiple terms
%   strsplit_re   - Splits a string into multiple terms with regex delimiter
%   strstartswith - Determines whether a string starts with a specified pattern
%   strgsub       - Replaces all substrings matching a specified pattern


text:
function [ftAllNew,transMdl] = ftTrans_pca(ftAll,maSrc,target,maLabeled,param)
%A baseline to other domain adaptation methods, use PCA or KPCA to extract
% features.
%
% Application scope:
%	+ Any scenario since no domain adaptation is done.
%
% ftAll:	All samples in all domains. n-by-m matrix, n is the number of
%	samples, m is the dimension of features.
% maSrc,target,maLabeled:	useless.
%
% param: Struct of hyper-parameters, please see the first cell of this
%	program ("default parameters") for details. You can set parameter p to
%	x by setting param.p = x. For parameters that are not set, default
%	values will be used.
%
% ftAllNew:	All samples in the learned subspace.
% transMdl:	A struct containing the model, transMdl.W is the projection
%	matrix.
%
% The ftProc_(k)pca_tr function is a part of the PRTools toolbox.
% if kerName is set in param, use KPCA; else, use original linear PCA.
%
% Copyright 2016 Ke YAN, Tsinghua Univ. http://yanke23.com , xjed09@gmail.com

%% default parameters
pcaCoef = 0; % see ftProc_pca_tr
kerName = 'poly'; % see ftProc_kpca_tr
kerSigma = 1;
degree = 2;

defParam

%% compute and project the samples
if ~isfield(param,'kerName')
	[ftAllNew,pcaModel] = ftProc_pca_tr(ftAll,[],struct('pcaCoef',pcaCoef));
else
	[ftAllNew,pcaModel] = ftProc_kpca_tr(ftAll,[],struct('pcaCoef',pcaCoef,...
		'kernel',kerName,'kerSigma',kerSigma,'degree',degree));
end

transMdl = pcaModel;

end

author: viggin

repo: domain-adaptation-toolbox


text:
function [model, deletedSEEDRxns] = deleteSeedGapfilledReactions(model, biomassReaction)
% Part of the DEMETER pipeline. Deletes reactions gapfilled by the Model
% SEED pipeline that are no longer needed after the reconstruction was
% refined.
%
% INPUT
% model             COBRA model structure
% biomassReaction   Biomass reaction abbreviation
%
% OUTPUT
% model             COBRA model structure
% deletedSEEDRxns   deleted gapfilled reactions
%
% .. Authors:
% Almut Heinken and Stefania Magnusdottir, 2016-2019

global CBT_LP_SOLVER
if isempty(CBT_LP_SOLVER)
    initCobraToolbox
end
solver = CBT_LP_SOLVER;

tol = 1e-8;

gfRxns = model.rxns(find(strcmp(model.grRules,'')));
gfRxns = union(gfRxns, model.rxns(strncmp('Unknown', model.grRules, length('Unknown'))));
gfRxns = union(gfRxns, model.rxns(strncmp('0000000.0.peg', model.grRules, length('0000000.0.peg'))));
gfRxns = union(gfRxns, model.rxns(strncmp('AUTOCOMPLETION', model.grRules, length('AUTOCOMPLETION'))));
gfRxns = union(gfRxns, model.rxns(strncmp('INITIALGAPFILLING', model.grRules, length('INITIALGAPFILLING'))));

% not consider exchange, demand, and sink reactions
gfRxns(find(strncmp(gfRxns,'EX_',3)),:)=[];
gfRxns(find(strncmp(gfRxns,'DM_',3)),:)=[];
gfRxns(find(strncmp(gfRxns,'sink_',5)),:)=[];

% only consider reactions that were already in draft reconstructions
translateRxns = readtable('ReactionTranslationTable.txt', 'Delimiter', '\t');
translateRxns=table2cell(translateRxns);
gfRxns=intersect(gfRxns,translateRxns(:,2));

noDelete = {biomassReaction; 'O2t'; 'O2t5i'; 'CO2t'; 'MG2abc'; 'H2Ot'; 'L2A6ODs'; 'G5SADs'; 'AOBUTDs'; 'ACGAMK'; 'DM_atp_c_'; 'FE3abc'; 'SO4t2'; 'TSULabc'; 'PDHbr'; 'PDHc'; 'INDOLEt2r'; 'SUCCt2r'; 'PHEt2r'; 'TRPt2r'};
checkDelete = setdiff(gfRxns, noDelete);

transp = {  % delete exchange reaction with following transporters
            'AHCYSts', 'EX_ahcys(e)'
            'ALAHISabc', 'EX_alahis(e)'
            'ALALEUabc', 'EX_alaleu(e)'
            'CBIabc', 'EX_cbi(e)'
            'CGLYabc', 'EX_cgly(e)'
            'CHORt', 'EX_chor(e)'
            'CMPt2', 'EX_cmp(e)'
            'COAt', 'EX_coa(e)'
            'Cut1', 'EX_cu2(e)'
            'CYTDt4', 'EX_cytd(e)'
            'DATPt', 'EX_datp(e)'
            'DDCAt', 'EX_ddca(e)'
            'DGSNt2', 'EX_dgsn(e)'
            'DGTPt', 'EX_dgtp(e)'
            'DPCOAt', 'EX_dpcoa(e)'
            'DTMPt', 'EX_dtmp(e)'
            'DTTPti', 'EX_dttp(e)'
            'ETHAt2', 'EX_etha(e)'
            'G6Pt6_2', 'EX_g6p(e)'
            'GLYC3Pt', 'EX_glyc3p[e]'
            'HDCEAtr', 'EX_hdcea(e)'
            'MALTHPabc', 'EX_malthp(e)'
            'MALTHXabc', 'EX_malthx(e)'
            'MK7te', 'EX_mqn7(e)'
            'NADPt', 'EX_nadp(e)'
            'OCDCAt2', 'EX_ocdca(e)'
%             'PTRCt2', 'EX_ptrc(e)'
            'SALCt2', 'EX_salc(e)'
            'THMMPt2', 'EX_thmmp(e)'
            'TTDCAt2', 'EX_ttdca(e)'
            'UREAt', 'EX_urea(e)'
            };
cnt = 1;
deletedSEEDRxns = {};
model = changeRxnBounds(model, model.rxns(strmatch('EX_', model.rxns)), -1000, 'l');
model = changeRxnBounds(model, 'EX_o2(e)', 0, 'l');

% minFlux = zeros(size(model.rxns));
% maxFlux = zeros(size(model.rxns));
% % reaction vector is not sorted, so cannot run fastFVA on all reactions at
% % once
% for n = 1:length(model.rxns)
%     model = changeObjective(model, model.rxns{n});
%     sol = optimizeCbModel(model, 'max');
%     maxFlux(n) = sol.f;
%     sol = optimizeCbModel(model, 'min');
%     minFlux(n) = sol.f;
% end

checkDelete=intersect(model.rxns,checkDelete,'stable');

% perform flux variability analysis
currentDir=pwd;
try
    [minFlux, maxFlux, ~, ~] = fastFVA(model, 0, 'max', 'ibm_cplex', ...
        checkDelete, 'S');
catch
    warning('fastFVA could not run, so fluxVariability is instead used. Consider installing fastFVA for shorter computation times.');
    cd(currentDir)
    [minFlux, maxFlux] = fluxVariability(model, 0, 'max', checkDelete);
end

model = changeObjective(model, biomassReaction);

% first check for and remove all blocked reactions
for i = 1:size(checkDelete, 1)
    if ~isempty(find(ismember(transp(:, 1), checkDelete{i, 1})))
        if abs(minFlux(i, 1)) < tol && abs(maxFlux(i, 1)) < tol
            modelTest = changeRxnBounds(model, checkDelete{i, 1}, 0, 'b');
            FBA = optimizeCbModel(modelTest, 'max');
            if FBA.f > tol
                deletedSEEDRxns{cnt, 1} = transp{find(strcmp(transp(:,1),checkDelete{i, 1})), 1};
                deletedSEEDRxns{cnt + 1, 1} = transp{find(strcmp(transp(:,1),checkDelete{i, 1})), 2};
                cnt = cnt + 2;
            end
        end
    else
        if abs(minFlux(i, 1)) < tol && abs(maxFlux(i, 1)) < tol
            modelTest = changeRxnBounds(model, checkDelete{i, 1}, 0, 'b');
            FBA = optimizeCbModel(modelTest, 'max');
            if FBA.f > tol
                deletedSEEDRxns{cnt, 1} = checkDelete{i, 1};
                cnt = cnt + 1;
            end
        end
    end
end

% then the reactions that are not blocked but may be safely deleted
% load complex medium
constraints = readtable('ComplexMedium.txt', 'Delimiter', 'tab');
constraints=table2cell(constraints);
constraints=cellstr(string(constraints));

% apply complex medium
modelTest = useDiet(model,constraints);
for i = 1:size(checkDelete, 1)
    if ~isempty(find(ismember(transp(:, 1), checkDelete{i, 1})))  % if deleted reaction is in transporter list
        trListInd = find(ismember(transp(:, 1), checkDelete{i, 1}));
        trInd = find(ismember(modelTest.rxns, transp{trListInd, 1}));
        exInd = find(ismember(modelTest.rxns, transp{trListInd, 2}));
        savedTrLB = modelTest.lb(trInd);
        savedTrUB = modelTest.ub(trInd);
        savedExLB = modelTest.lb(exInd);
        savedExUB = modelTest.ub(exInd);
        modelTest = changeRxnBounds(modelTest, transp{trListInd, 1}, 0, 'b');
        modelTest = changeRxnBounds(modelTest, transp{trListInd, 2}, 0, 'b');
        FBA = optimizeCbModel(modelTest, 'max');
        if FBA.f > 1e-6
            deletedSEEDRxns{cnt, 1} = transp{trListInd, 1};
            deletedSEEDRxns{cnt + 1, 1} = transp{trListInd, 2};
            cnt = cnt + 2;
            modelTest = removeRxns(modelTest, {transp{trListInd, 1}, transp{trListInd, 2}});
        else
            modelTest.lb(trInd) = savedTrLB;
            modelTest.ub(trInd) = savedTrUB;
            modelTest.lb(exInd) = savedExLB;
            modelTest.ub(exInd) = savedExUB;
        end
    else
        rxnInd = find(ismember(modelTest.rxns, checkDelete{i, 1}));
        savedRxnLB = modelTest.lb(rxnInd);
        savedRxnUB = modelTest.ub(rxnInd);
        modelTest = changeRxnBounds(modelTest, checkDelete{i, 1}, 0, 'b');
        FBA = optimizeCbModel(modelTest, 'max');
        if FBA.f > tol
            deletedSEEDRxns{cnt, 1} = checkDelete{i, 1};
            cnt = cnt + 1;
            modelTest = removeRxns(modelTest, checkDelete{i, 1});
        else
            modelTest.lb(rxnInd) = savedRxnLB;
            modelTest.ub(rxnInd) = savedRxnUB;
        end
    end
end
model = removeRxns(model, deletedSEEDRxns);

% change back to unlimited medium
% list exchange reactions
exchanges = model.rxns(strncmp('EX_', model.rxns, 3));
% open all exchanges
model = changeRxnBounds(model, exchanges, -1000, 'l');
model = changeRxnBounds(model, exchanges, 1000, 'u');

end

author: opencobra

repo: cobratoolbox


text:
%GSP_DEMO Tutorial on the GSPBox
%
%   In this demo, we are going to show the basic operations of the GSPBox.
%   To lauch the toolbox, just go into the repository where the GSPBox was
%   extracted and type:
%
%           gsp_start;
%
%   A banner will popup telling you that everything happens correctly. To
%   speedup some processing, you might want to compile some mexfile. Refer
%   to |gsp_make| for more informations. However, if the compilation is not
%   working on your computer, keep quiet, everything should still work and
%   most of the routine are implemented only in matlab.
%
%   Most likely, the first thing you would like to do is to create a graph.
%   To do so, you only need the adjacendy or the weight matrix $W$. Once
%   you have it, you can construct a graph using::
%
%           G = gsp_graph(W);
%
%   This function will create a full structure ready to be used with the
%   toolbox. To know a bit more about what is in this structure, you can
%   refer to the help of the function |gsp_graph_default_parameters|.
%
%   The GSPBox contains also a list of graph generators. To see a full list
%   of these graphs, type:::
%
%           help graphs
%
%   For this demo, we will use the graph |gsp_logo|. You can load it
%   using:::
%
%           G = gsp_logo
%
%   Here observe the attribute of the structure *G*.
%
%   * *G.W*: Weight matrix
%   * *G.A*: Adacency matrix
%   * *G.N*: Number of nodes
%   * *G.type*: Type of graph
%   * *G.directed*: 1 if the graph is directed, 0 if not
%   * *G.lap_type*: Laplacian type
%   * *G.d*: Degree vector
%   * *G.Ne*: Number of edges
%   * *G.coords*: Coordinates of the vertices
%   * *G.plotting*: Plotting parameters
%
%   In the folder 'plotting', the GSPBox contains some plotting routine.
%   For instance, we can plot a graph using::
%
%           gsp_plot_graph(G);
%
%   .. figure::
%
%      GSP graph
%
%      This figure shows the result of the command 'gsp_plot_graph(G)'
%
%   Wonderful! Isn't it? Now, let us start to analyse this graph. To compute
%   graph Fourier transform or exact graph filtering, you need to
%   precompute the Fourier basis of the graph. This operation could be
%   relatively long since it involves a full diagonalization of the
%   Laplacian. Don't worry, you do not need to perform this operation to
%   filter signals on graph. The fourier basis is computed by::
%
%           G = gsp_compute_fourier_basis(G);
%
%   The function |gsp_compute_fourier_basis| add two new fields to the
%   structure *G*:
%
%   * *G.U*: The eigenvectors of the Fourier basis
%   * *G.e*: The eigenvalues
%
%   The fourier eigenvectors does look like a sinusoide on the graph. Let's
%   plot the second and the third ones. (The first one is constant!)::
%
%           gsp_plot_signal(G,G.U(:,2));
%           title('Second eigenvector')
%           subplot(212)
%           gsp_plot_signal(G,G.U(:,3));
%           title('Third eigenvector')
%
%   .. figure::
%
%      Eigenvectors
%
%
%
%   Now, we are going to show a basic filtering operation. Filters are usually
%   defined in the spectral domain. To define the following filter
%
%   ..   h(x) = 1/(1+tau*x),
%
%   .. math:: h(x) =\frac{1}{1+\tau x},
%
%   just write in Matlab::
%
%           tau = 1;
%           h = @(x) 1./(1+tau*x);
%
%   Hint: You can define filterbank using cell array!
%
%   Let's display this filter::
%
%           gsp_plot_filter(G,h);
%
%   .. figure::
%
%      Low pass filter $h$
%
%      The filter $h$ is plotted along all the spectrum of the graph.
%      The black cross are the eigenvalues of the Laplacian. They are the
%      points where the continuous filter will be evaluated to create a
%      discrete filter.
%
%   To apply the filter to a given signal, you only need to run a single
%   function::
%
%           % Create a signal
%           f = zeros(G.N,1);
%           f(G.info.idx_g) = -1;
%           f(G.info.idx_s) = 1;
%           f(G.info.idx_p) = -0.5;
%           f = f + 0.3*randn(G.N,1);
%           % Remove the noise
%           f2 = gsp_filter(G,h,f);
%
%   `gsp_filter` is actually a shortcut to |gsp_filter_analysis|.
%   `gsp_filter_analysis` performs the analysis operator associated to a
%   filterbank. See the |gsp_demo_wavelet| for more information.
%
%   Finnaly, we display the result of this low pass filtering on the graph::
%
%           figure;
%           subplot(211)
%           gsp_plot_signal(G,f);
%           title('Signal with noise')
%           subplot(212)
%           gsp_plot_signal(G,f2);
%           title('Signal denoised');
%
%   .. figure::
%
%      Result of filtering
%
%      The noise is largely removed thanks to the filter. However, some
%      energy is diffused between the letters. This is the typical
%      behaviour of a low pass filter.
%
%   Enjoy the GSPBOX !
%

% Author: Nathanael Perraudin
% Date : 14 August 2014

clear;
close all;

G = gsp_logo;

% display the graph
figure;
gsp_plot_graph(G);

%% Compute the Fourier basis

G = gsp_compute_fourier_basis(G);

% Display an eigenvector
figure;
subplot(211)
gsp_plot_signal(G,G.U(:,2));
title('Second eigenvector')
subplot(212)
gsp_plot_signal(G,G.U(:,3));
title('Third eigenvector')

%% Create a signal

f = zeros(G.N,1);
f(G.info.idx_g) = -1;
f(G.info.idx_s) = 1;
f(G.info.idx_p) = -0.5;
f = f + 0.3*randn(G.N,1);

%% Define a low pass filter
tau = 1;
h = @(x) 1./(1+tau*x);

figure
gsp_plot_filter(G,h);
title('Filter h')

%% Perform the filtering operation
f2 = gsp_filter(G,h,f);
% f2 = gsp_filter_analysis(G,h,f);

% Display the result
figure;
subplot(211)
gsp_plot_signal(G,f);
title('Signal with noise')
subplot(212)
gsp_plot_signal(G,f2);
title('Signal denoised');

repo: gspbox


text:
%POLYGON Polygon class
%
% A general class for manipulating polygons and vectors of polygons.
%
% Methods::
%  plot           Plot polygon
%  area           Area of polygon
%  moments        Moments of polygon
%  centroid       Centroid of polygon
%  perimeter      Perimter of polygon
%  transform      Transform polygon
%  inside         Test if points are inside polygon
%  intersection   Intersection of two polygons
%  difference     Difference of two polygons
%  union          Union of two polygons
%  xor            Exclusive or of two polygons
%  display        print the polygon in human readable form
%  char           convert the polgyon to human readable string
%
% Properties::
%  vertices   List of polygon vertices, one per column
%  extent     Bounding box [minx maxx; miny maxy]
%  n          Number of vertices
%
% Notes::
% - This is reference class object
% - Polygon objects can be used in vectors and arrays
%
% Acknowledgement::
%
% The methods: inside, intersection, difference, union, and xor are based on code
% written by:
%  Kirill K. Pankratov, kirill@plume.mit.edu,
%  http://puddle.mit.edu/~glenn/kirill/saga.html
% and require a licence.  However the author does not respond to email regarding
% the licence, so use with care, and modify with acknowledgement.

% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

% TODO
%  split the code in two.  Simple polygon functions in Polgon class, subclass with
%  Pankratov code to Polygon2.
%  add method to detect empty polygon, overload isempty

classdef Polygon < handle

    properties
        vertices
        extent
    end

    properties (Dependent=true)
        n
        x
        y
    end

    methods

        function p = Polygon(v, wh)
            %Polygon.Polygon Polygon class constructor
            %
            % P = Polygon(V) is a polygon with vertices given by V, one column per
            % vertex.
            %
            % P = Polygon(C, WH) is a rectangle centred at C with dimensions
            % WH=[WIDTH, HEIGHT].

            if nargin == 0
                p.n = 0;
                p.vertices = [];
                return;
            end

            if nargin < 2
                if numrows(v) ~= 2
                    error('vertices must have two rows');
                end
                p.vertices = v;
            end
            if nargin == 2
                if length(v) ~= 2
                    error('first argument must be polygon centre');
                end
                if length(wh) ~= 2
                    error('second arugment must be width height');
                end

                p.vertices = [
                    v(1)-wh(1)/2  v(1)+wh(1)/2  v(1)+wh(1)/2  v(1)-wh(1)/2
                    v(2)-wh(2)/2  v(2)-wh(2)/2  v(2)+wh(2)/2 v(2)+wh(2)/2 ];
            end

            % compute the extent
            p.extent(1,1) = min(p.x);
            p.extent(1,2) = max(p.x);
            p.extent(2,1) = min(p.y);
            p.extent(2,2) = max(p.y);
        end

        function r = get.n(p)
            r = numcols(p.vertices);
        end

        function r = get.x(p)
            r = p.vertices(1,:)';
        end

        function r = get.y(p)
            r = p.vertices(2,:)';
        end

        function r = set.n(p)
            error('cant set property');
        end
        function r = set.x(p)
            error('cant set property');
        end
        function r = set.y(p)
            error('cant set property');
        end

        function s = char(p)
            %Polygon.char String representation
            %
            % S = P.char() is a compact representation of the polgyon in human
            % readable form.
            s = sprintf('%d vertices', p.n);
        end

        function display(p)
            %Polygon.display Display polygon
            %
            % P.display() displays the polygon in a compact human readable form.
            %
            % See also Polygon.char.
            loose = strcmp( get(0, 'FormatSpacing'), 'loose');
            if loose
                disp(' ');
            end
            disp([inputname(1), ' = '])
            if loose
                disp(' ');
            end
            disp(char(p))
            if loose
                disp(' ');
            end
        end

        function plot(plist, varargin)
            %Polygon.plot Draw polygon
            %
            % P.plot() draws the polygon P in the current plot.
            %
            % P.plot(LS) as above but pass the arguments LS to plot.
            %
            % Notes::
            % - The polygon is added to the current plot.

            opt.fill = [];
            [opt,args] = tb_optparse(opt, varargin);

            ish = ishold
            hold all

            for p=plist
                % for every polygon in the list

                % get the vertices
                X = p.vertices(1,:)';
                Y = p.vertices(2,:)';

                while true

                    % look for NaNs which indicate disjoint vertex sets
                    k = find(isnan(X));

                    if length(k) > 0
                        % if a NaN chop out the segment before and after
                        k = k(1);
                        x = X(1:k-1);
                        y = Y(1:k-1);
                        X = X(k+1:end);
                        Y = Y(k+1:end);
                    else
                        x = X; y = Y;
                    end

                    % close the polygon
                    x = [x; x(1)];
                    y = [y; y(1)];

                    if opt.fill
                        patch(x, y, opt.fill);
                    else
                        plot(x, y, args{:});
                    end

                    if length(k) == 0
                        break;
                    end
                end
            end
            hold(ish)
        end

        function a = area(p)
            %Polygon.area Area of polygon
            %
            % A = P.area() is the area of the polygon.
            %
            % See also Polygon.moments.
            a = p.moments(0, 0);
        end

        function m = moments(p, mp, mq)
            %Polygon.moments Moments of polygon
            %
            % A = P.moments(p, q) is the pq'th moment of the polygon.
            %
            % See also Polygon.area, Polygon.centroid, mpq_poly.
            m = mpq_poly(p.vertices, mp, mq);
        end

        function q = transform(p, T)
            %Polygon.transform Transform polygon vertices
            %
            % P2 = P.transform(T) is a new Polygon object whose vertices have
            % been transformed by the SE(2) homgoeneous transformation T (3x3).
            if length(T) == 3
                T = se2(T);
            end
            q = Polygon( homtrans(T, p.vertices) );
        end

        function f = inside(p, points)
            %Polygon.inside Test if points are inside polygon
            %
            % IN = P.inside(P) tests if points given by columns of P (2xN) are inside
            % the polygon.  The corresponding elements of IN (1xN) are either true or
            % false.
            IN = inpolygon(points(1,:), points(2,:), p.x, p.y)
        end

        function c = centroid(p)
            %Polygon.centroid Centroid of polygon
            %
            % X = P.centroid() is the centroid of the polygon.
            %
            % See also Polygon.moments.
            xc = p.moments(1,1) / p.area();
        end

        function r = perimeter(p)
            %Polygon.perimeter Perimeter of polygon
            %
            % L = P.perimeter() is the perimeter of the polygon.
            p = sum(sqrt(diff(p.x).^2+diff(p.y).^2));
        end

        function f = intersect(p, plist)
            %Polygon.intersect Intersection of polygon with list of polygons
            %
            % I = P.intersect(PLIST) indicates whether or not the Polygon P
            % intersects with
            %
            % i(j) = 1 if p intersects polylist(j), else 0.

            % Based on ISINTPL
            %  Copyright (c) 1995 by Kirill K. Pankratov,
            %       kirill@plume.mit.edu.
            %       06/20/95, 08/25/95
            f = [];
            for q=plist
                f = [f isintpl(p.x, p.y, q.x, q.y)];
            end
        end

        function f = intersect_line(p, l)
            %Polygon.intersect_line Intersection of polygon and line segment
            %
            % I = P.intersect_line(L) is the intersection points of a polygon P with
            % the line segment L=[x1 x2; y1 y2].  I (2xN) has one column per
            % intersection, each column is [x y]'.

            f = [];
            % find intersections
            for i=1:p.n
                in = mod(i, p.n)+1;

                xv = [p.x(i); p.x(in)];
                yv = [p.y(i); p.y(in)];
                intsec = iscross(xv, yv, l(1,:)', l(2,:)');
                if intsec
                    [x,y] = intsecl(xv, yv, l(1,:)', l(2,:)');
                    f = [f [x;y]];
                end
            end
        end

        function r = difference(p, q)
            %Polygon.difference Difference of polygons
            %
            % D = P.difference(Q) is polygon P minus polygon Q.
            %
            % Notes::
            % - If polygons P and Q are not intersecting, returns
            %   coordinates of P.
            % - If the result D is not simply connected or consists of
            %   several polygons, resulting vertex list will contain NaNs.

            % POLYDIFF  Difference of 2 polygons.
            %       [XO,YO] = POLYDIFF(X1,Y1,X2,Y2) Calculates polygon(s) P
            %   of difference of polygons P1 and P1 with coordinates
            %   X1, Y1 and X2, Y2.
            %   The resulting polygon(s) is a set of all points which belong
            %   to P1 but not to to P2: P = P1 & ~P2.
            %   The input polygons must be non-self-intersecting and
            %   simply connected.
            %
            %       If polygons P1, P2 are not intersecting, returns
            %   coordinates of the first polygon X1, X2.
            %   If the result P is not simply connected or consists of several
            %   polygons, resulting boundary consists of concatenated
            %   coordinates of these polygons, separated by NaN.
            %  Copyright (c) 1995 by Kirill K. Pankratov,
            %       kirill@plume.mit.edu.
            %       06/25/95

            % Call POLYBOOL with flag=3
            [xo,yo,ind] = polybool(p.x, p.y, q.x, q.y, 3);

            r = Polygon([xo(:) yo(:)]');
        end

        function r = intersection(p, q)
            %Polygon.intersection Intersection of polygons
            %
            % I = P.intersection(Q) is a Polygon representing the
            % intersection of polygons P and Q.
            %
            % Notes::
            % - If these polygons are not intersecting, returns empty polygon.
            % - If intersection consist of several disjoint polygons
            %   (for non-convex P or Q) then vertices of I is the concatenation
            %   of the vertices of these polygons.

            % POLYINTS  Intersection of 2 polygons.
            %   [XO,YO] = POLYINTS(X1,Y1,X2,Y2) Calculates polygon(s)
            %   if intersection of polygons with coordinates X1, Y1
            %   and X2, Y2.
            %       The resulting polygon(s) is a set of all points which
            %   belong to both P1 and P2: P = P1 & P2.
            %   These polygons must be non-self-intersecting and
            %   simply connected.
            %
            %   If these polygons are not intersecting, returns empty.
            %   If intersection consist of several disjoint polygons
            %   (for non-convex P1 or P2) output vectors XO, YO consist
            %   of concatenated cooddinates of these polygons,
            %  Copyright (c) 1995 by Kirill K. Pankratov,
            %       kirill@plume.mit.edu.
            %       06/25/95

            % Call POLYBOOL with flag=1
            [xo,yo,ind] = polybool(p.x, p.y, q.x, q.y, 1);

            r = Polygon([xo(:) yo(:)]');
        end

        function r = union(p, q)
            %Polygon.union Union of polygons
            %
            % I = P.union(Q) is a polygon representing the
            % union of polygons P and Q.
            %
            % Notes::
            % - If these polygons are not intersecting, returns a polygon with
            %   vertices of both polygons separated by NaNs.
            % - If the result P is not simply connected (such as a polygon
            %   with a "hole") the resulting contour consist of counter-
            %   clockwise "outer boundary" and one or more clock-wise
            %   "inner boundaries" around "holes".

            % POLYUNI  Union of 2 polygons.
            %       [XO,YO] = POLYINT(X1,Y1,X2,Y2) Calculates polygon(s) P
            %       which is (are) union of polygons P1 and P2 with coordinates
            %   X1, Y1 and X2, Y2.
            %   The resulting polygon(s) is a set of all points which belong
            %   either to P1 or to P2: P = P1 | P2.
            %   The input polygons must be non-self-intersecting and
            %   simply connected.
            %
            %       If polygons P1, P2 are not intersecting, returns
            %   coordinates of the both polygons separated by NaN.
            %   If both P1 and P2 are convex, their boundaries can have no
            %   more than 2 intersections. The result is also a convex
            %   polygon.
            %       If the result P is not simply connected (such as a polygon
            %   with a "hole") the resulting contour consist of counter-
            %   clockwise "outer boundary" and one or more clock-wise
            %   "inner boundaries" around "holes".
            %  Copyright (c) 1995 by Kirill K. Pankratov,
            %       kirill@plume.mit.edu.
            %       06/25/95

            % Call POLYBOOL with flag=2 ..........
            [xo,yo,ind] = polybool(p.x, p.y, q.x, q.y, 2);

            r = Polygon([xo(:) yo(:)]');
        end

        function r = xor(p, q)
            %Polygon.xor Exclusive or of polygons
            %
            % I = P.union(Q) is a polygon representing the
            % exclusive-or of polygons P and Q.
            %
            % Notes::
            % - If these polygons are not intersecting, returns a polygon with
            %   vertices of both polygons separated by NaNs.
            % - If the result P is not simply connected (such as a polygon
            %   with a "hole") the resulting contour consist of counter-
            %   clockwise "outer boundary" and one or more clock-wise
            %   "inner boundaries" around "holes".

            % POLYXOR  Exclusive OR of 2 polygons.
            %       [XO,YO] = POLYXOR(X1,Y1,X2,Y2) Calculates polygon(s) P
            %   of difference of polygons P1 and P1 with coordinates
            %   X1, Y1 and X2, Y2.
            %   The resulting polygon(s) is a set of all points which belong
            %   either to P1 or to P2 but not to both:
            %   P = (P1 & ~P2) | (P2 & ~P1).
            %   The input polygons must be non-self-intersecting and
            %   simply connected.
            %
            %       If polygons P1, P2 are not intersecting, returns
            %   coordinates of the both polygons separated by NaN.
            %   If the result P is not simply connected or consists of several
            %   polygons, resulting boundary consists of concatenated
            %   coordinates of these polygons, separated by NaN.

            %  Copyright (c) 1995 by Kirill K. Pankratov,
            %       kirill@plume.mit.edu.
            %       06/25/95

            % Call POLYBOOL twice with flag=3
            [xx,yy,ind] = polybool(p.x, p.y, q.x, q.y, 3);

            xo = [xx; NaN]; yo = [yy; NaN];
            [xx,yy,ind] = polybool(q.x, q.y, p.x, p.y, 3);

            xo = [xo; xx]; yo = [yo; yy];

            r = Polygon([xo(:) yo(:)]');

        end
    end % methods
end % classdef

function [is,in,un] = interval(x1,x2)

    % Intersection and union of 2 intervals.
    %   [IS,IN,UN] = INTERVAL(X1,X2) calculates pair-wise
    %   intersection IN and union UN of N pairs of
    %   intervals with coordinates X1 and X2 (both are
    %   2 by N vectors). Returns 1 by N boolean vector IS
    %   equal to 1 if intervals have non-empty intersection
    %   and 0 if they don't.

    %  Copyright (c) 1995 by Kirill K. Pankratov,
    %       kirill@plume.mit.edu.
    %       08/24/95

    % Handle input ...........................
    if nargin==0, help interval, return, end
    if nargin==1
        un = x1;
    else
        un = [x1; x2];
    end

    [in,un] = sort(un);     % Sort both intervals together
    un = un(1:2,:)-1;
    is = sum(floor(un/2));  % Check for [0 0 1 1] or [1 1 0 0]
    is = (is==1);
    ii = find(in(2,:)==in(3,:));
    is(ii) = .5*ones(size(ii));

    % Extract intersection and union from sorted coordinates
    if nargout>1
        un = in([1 4],:);
        in = in(2:3,:);
        in(:,~is) = flipud(in(:,~is));
    end
end

function [is,S] = iscross(x1,y1,x2,y2,tol)

    % ISCROSS  Finds whether pairs of lines cross each other
    %   [IS,S] = ISCROSS(X1,Y1,X2,Y2) where arguments X1, Y1,
    %   X2, Y2 are all 2 by N matrices are coordinates of
    %   ends of the pairs of line segments.
    %   Returns vector  IS (1 by N)  consisting of ones if
    %   corresponding pairs cross each other, zeros if they
    %   don't and .5 if an end of one line segment lies on
    %   another segment.
    %   Also returns a matrix  S (4 by N) with each row
    %   consisting of cross products (double areas of
    %   corresponding triangles) built on the following points:
    %   (X2(1,:),Y2(1,:)),(X1(1,:),Y1(1,:)),(X2(2,:),Y2(2,:)),
    %   (X2(1,:),Y2(1,:)),(X1(2,:),Y1(2,:)),(X2(2,:),Y2(2,:))
    %   (X1(1,:),Y1(1,:)),(X2(1,:),Y2(1,:)),(X1(2,:),Y1(2,:))
    %   (X1(1,:),Y1(1,:)),(X2(2,:),Y2(2,:)),(X1(2,:),Y1(2,:))
    %   The signs of these 4 areas can be used to determine
    %   whether these lines and their continuations cross each
    %   other.
    %   [IS,S] = ISCROSS(X1,Y1,X2,Y2,TOL) uses tolerance TOL
    %   for detecting the crossings (default is 0).

    %  Copyright (c) 1995 by Kirill K. Pankratov
    %       kirill@plume.mit.edu
    %       08/14/94, 05/18/95, 08/25/95

    % Defaults and parameters .......................
    tol_dflt = 0; % Tolerance for area calculation
    is_chk = 1;   % Check input arguments

    % Handle input ..................................
    if nargin==0, help iscross, return, end
    if nargin<4           % Check if all 4 entered
        error('  Not enough input arguments')
    end
    if nargin<5, tol = tol_dflt; end
    if tol < 0, is_chk = 0; tol = 0; end

    % Check the format of arguments .................
    if is_chk
        [x1,y1,x2,y2] = linechk(x1,y1,x2,y2);
    end

    len = size(x1,2);
    o2 = ones(2,1);

    % Find if the ranges of pairs of segments intersect
    [isx,S,A] = interval(x1,x2);
    scx = diff(A);
    [isy,S,A] = interval(y1,y2);
    scy = diff(A);
    is = isx & isy;

    % If S values are not needed, extract only those pairs
    % which have intersecting ranges ..............
    if nargout < 2
        isx = find(is);  % Indices of pairs to be checked
        % further
        x1 = x1(:,isx);
        x2 = x2(:,isx);
        y1 = y1(:,isx);
        y2 = y2(:,isx);
        is = is(isx);
        if isempty(is), is = zeros(1,len); return, end
        scx = scx(isx);
        scy = scy(isx);
    end

    % Rescale by ranges ...........................
    x1 = x1.*scx(o2,:);
    x2 = x2.*scx(o2,:);
    y1 = y1.*scy(o2,:);
    y2 = y2.*scy(o2,:);

    % Calculate areas .............................
    S = zeros(4,length(scx));
    S(1,:) = (x2(1,:)-x1(1,:)).*(y2(2,:)-y1(1,:));
    S(1,:) = S(1,:)-(x2(2,:)-x1(1,:)).*(y2(1,:)-y1(1,:));

    S(2,:) = (x2(1,:)-x1(2,:)).*(y2(2,:)-y1(2,:));
    S(2,:) = S(2,:)-(x2(2,:)-x1(2,:)).*(y2(1,:)-y1(2,:));

    S(3,:) = (x1(1,:)-x2(1,:)).*(y1(2,:)-y2(1,:));
    S(3,:) = S(3,:)-(x1(2,:)-x2(1,:)).*(y1(1,:)-y2(1,:));

    S(4,:) = (x1(1,:)-x2(2,:)).*(y1(2,:)-y2(2,:));
    S(4,:) = S(4,:)-(x1(2,:)-x2(2,:)).*(y1(1,:)-y2(2,:));

    % Find if they cross each other ...............
    is = (S(1,:).*S(2,:)<=0)&(S(3,:).*S(4,:)<=0);

    % Find very close to intersection
    isy = min(abs(S));
    ii = find(isy<=tol & is);
    is(ii) = .5*ones(size(ii));

    % Output
    if nargout < 2
        isy = zeros(1,len);
        isy(isx) = is;
        is = isy;

    else
        isy = scx.*scy;
        ii = find(~isy);
        isy(ii) = ones(size(ii));
        S = S./isy(ones(4,1),:);

    end

end

function [xo,yo,ind] = polybool(x1,y1,x2,y2,flag)

    %   [XO,YO] = POLYBOOL(X1,Y1,X2,Y2,FLAG)
    %   calulates results of Boolean operations on
    %   a pair of polygons.
    %   FLAG Specifies the type of the operation:
    %   1 - Intersection (P1 & P2)
    %   2 - Union (P1 | P2)
    %   3 - Difference (P1 & ~P2)

    %  Copyright (c) 1995 by Kirill K. Pankratov,
    %       kirill@plume.mit.edu.
    %       06/25/95, 09/07/95

    %   This program calls the following functions:
    %   AREA, ISINTPL, ISCROSS, INTSECL.

    % Algorithm:
    %  1. Check boundary contour directions (area).
    %     For intersection and union make all
    %     counter-clockwise. For difference make the second
    %    contour clock-wise.
    %  2. Calculate matrix of intersections (function ISINTPL).
    %     Quick exit if no intersections.
    %  3. For intersecting segments calculate intersection
    %     coordinates (function INTSECL).
    %  4. Sort intersections along both contours.
    %  5. Calculate sign of cross-product between intersectiong
    %     segments. This will give which contour goes "in" and
    %     "out" at intersections.
    %
    %  6. Start with first intersection:
    %     Determine direction to go ("in" for intersection,
    %     "out" for union).
    %     Move until next intersection, switch polygons at each
    %     intersection until coming to the initial point.
    %     If not all intersections are encountered, the
    %     resulting polygon is disjoint. Separate output
    %     coordinates by NaN and repeat procedure until all
    %     intersections are counted.

    % Default for flag
    flag_dflt = 1; % 1- intersec., 2-union, 3 - diff.

    % Handle input
    if nargin==0, help polybool, return, end
    if nargin < 4
        error(' Not enough input arguments')
    end
    if nargin<5, flag = flag_dflt; end

    x1 = x1(:); y1 = y1(:);
    x2 = x2(:); y2 = y2(:);
    l1 = length(x1);
    l2 = length(x2);

    % Check areas and reverse if negative
    nn1 = area(x1,y1);
    if nn1<0, x1 = flipud(x1); y1 = flipud(y1); end
    nn2 = area(x2,y2);
    if (nn2<0 & flag<3) | (nn2>0 & flag==3)
        x2 = flipud(x2); y2 = flipud(y2);
    end

    % If both polygons are identical ........
    if l1==l2
        if all(x1==x2) & all(y1==y2)
            if flag<3, xo = x1; yo = y1; ind = 1:l1;
            else, xo = []; yo = []; ind = []; end
            return
        end
    end

    % Calculate matrix of intersections .....
    [is,C] = isintpl(x1,y1,x2,y2);
    is = any(any(C));

    % Quick exit if no intersections ........
    if ~is
        if flag==1       % Intersection
            xo=[]; yo = [];
        elseif flag==2   % Union
            xo = [x1; nan; x2];
            yo = [y1; nan; y2];
        elseif flag==3   % Difference
            xo = x1; yo = y1;
        end
        return
    end

    % Mark intersections with unique numbers
    i1 = find(C);
    ni = length(i1);
    C(i1) = 1:ni;

    % Close polygon contours
    x1 = [x1; x1(1)]; y1 = [y1; y1(1)];
    x2 = [x2; x2(1)]; y2 = [y2; y2(1)];
    l1 = length(x1);  l2 = length(x2);

    % Calculate intersections themselves
    [i1,i2,id] = find(C);
    xs1 = [x1(i1) x1(i1+1)]'; ys1 = [y1(i1) y1(i1+1)]';
    xs2 = [x2(i2) x2(i2+1)]'; ys2 = [y2(i2) y2(i2+1)]';

    % Call INTSECL ............................
    [xint,yint] = intsecl(xs1,ys1,xs2,ys2);

    % For sements belonging to the same line
    % find interval of intersection ...........
    ii = find(xint==inf);
    if ~isempty(ii)
        [is,inx] = interval(xs1(:,ii),xs2(:,ii));
        [is,iny] = interval(ys1(:,ii),ys2(:,ii));
        xint(ii) = mean(inx);
        yint(ii) = mean(iny);
    end

    % Coordinate differences of intersecting segments
    xs1 = diff(xs1); ys1 = diff(ys1);
    xs2 = diff(xs2); ys2 = diff(ys2);

    % Calculate cross-products
    cp = xs1.*ys2-xs2.*ys1;
    cp = cp>0;
    if flag==2, cp=~cp; end % Reverse if union
    cp(ii) = 2*ones(size(ii));

    % Sort intersections along the contours
    ind = (xint-x1(i1)').^2+(yint-y1(i1)').^2;
    ind = ind./(xs1.^2+ys1.^2);
    cnd = min(ind(ind>0));
    ind = ind+i1'+i2'/(ni+1)*cnd*0;
    [xo,ii] = sort(ind);
    xs1 = id(ii);
    [xo,ind] = sort(xs1);
    ind = rem(ind,ni)+1;
    xs1 = xs1(ind);

    ind = (xint-x2(i2)').^2+(yint-y2(i2)').^2;
    ind = ind./(xs2.^2+ys2.^2);
    cnd = min(ind(ind>0));
    [xo,ii] = sort(i2'+ind+i1'/(ni+1)*cnd*0);
    xs2 = id(ii);
    [xo,ind] = sort(xs2);
    ind = rem(ind,ni)+1;
    xs2 = xs2(ind);

    % Combine coordinates in one vector
    x1 = [x1; x2]; y1 = [y1; y2];

    % Find max. possible length of a chain
    xo = find(any(C'));
    xo = diff([xo xo(1)+l1]);
    mlen(1) = max(xo);
    xo = find(any(C));
    xo = diff([xo xo(1)+l2]);
    mlen(2) = max(xo);

    % Check if multiple intersections in one segment
    xo = diff([i1 i2]);
    is_1 = ~all(all(xo));

    % Begin counting intersections *********************

    % Initialization ..................
    int = zeros(size(xint));
    nn = 1;   % First intersection
    nn1 = i1(nn); nn2 = i2(nn);
    b = cp(nn);
    is2 = b==2;
    xo = []; yo = []; ind = [];
    closed = 0;

    % Proceed until all intersections are counted
    while ~closed  % begin counting `````````````````````0

        % If contour closes, find new starting point
        if int(nn) & ~all(int)
            ii = find(int);
            C(id(ii)) = zeros(size(ii));
            nn = min(find(~int));  % Next intersection
            nn1 = i1(nn);
            nn2 = i2(nn);
            xo = [xo; nan];        % Separate by NaN
            yo = [yo; nan];
            ind = [ind; nan];
            % Choose direction ......
            b = cp(nn);
        end

        % Add current intersection ......
        xo = [xo; xint(nn)];
        yo = [yo; yint(nn)];
        ind = [ind; 0];
        int(nn) = 1;
        closed = all(int);

        % Find next segment
        % Indices for next intersection
        if ~b, nn = xs1(nn);
        else,  nn = xs2(nn);
        end
        if ~b, pt0 = nn1; else,  pt0 = nn2; end

        nn1 = i1(nn);
        nn2 = i2(nn);

        if b, pt = nn2; else, pt = nn1; end

        if b, pt0 = pt0+l1; pt = pt+l1; end
        ii = (pt0+1:pt);

        % Go through the beginning ..............
        cnd = pt<pt0 | (pt==pt0 & is_1 & flag>1);
        if cnd
            if ~b,  ii = [pt0+1:l1 1:pt];
            else,   ii = [pt0+1:l1+l2 l1+1:pt];
            end
        end
        len = length(ii);
        cnd = b & len>mlen(2);
        cnd = cnd | (~b & len>mlen(1));
        if is2 | cnd, ii=[]; end

        % Add new segment
        xo = [xo; x1(ii)];
        yo = [yo; y1(ii)];
        ind = [ind; ii'];

        % Switch direction
        if cp(nn)==2, b = ~b; is2 = 1;
        else, b = cp(nn); is2 = 0;
        end

    end    % End while (all intersections) '''''''''''''''0

    % Remove coincident successive points
    ii = find(~diff(xo) & ~diff(yo));
    xo(ii) = []; yo(ii) = []; ind(ii) = [];

    % Remove points which are
    ii = find(isnan(xo));
    if ~isempty(ii)
        i2 = ones(size(xo));
        ii = [ii; length(xo)+1];

        i1 = find(diff(ii)==3);
        i1 = ii(i1);
        i1 = [i1; i1+1; i1+2];
        i2(i1) = zeros(size(i1));

        i1 = find(diff(ii)==2);
        i1 = ii(i1);
        i1 = [i1; i1+1];
        i2(i1) = zeros(size(i1));

        xo = xo(i2); yo = yo(i2); ind = ind(i2);
    end
end

function  [xo,yo] = intsecpl(xv,yv,xl,yl,trace)

    % INTSECPL Intersection of a polygon and a line.
    %   [XI,YI] = INTSECPL(XV,YV,XL,YL) calculates
    %   intersections XI, YI of a polygon with vertices XV,
    %   YV and a line specified by pairs of end coordinates
    %   XL = [XL0 XL1], YL = [YL0 YL1]. Line is assumed to
    %   continue beyond the range of end points.
    %   INTSECPL(XV,YV,[A B]) uses another specification for
    %   a line: Y = A*X+B.
    %
    %   If a line does not intersect polygon, returns empty
    %   XI, YI.
    %   For convex polygon maximum number of intersections is
    %   2, for non-convex polygons multiple intersections are
    %   possible.
    %
    %   INTSECPL(XV,YV,XL,YL)  by itself or
    %   [XI,YI] = INTSECPL(XV,YV,XL,YL,1) plots polygon,
    %   a line segment and intersection segment(s)
    %   (part(s) of the same line inside the polygon).

    %  Copyright (c) 1995 by Kirill K. Pankratov,
    %       kirill@plume.mit.edu.
    %       06/25/95, 08/27/95, 09/27/95

    %  Calls ISCROSS, INTSECL programs.

    % Defaults and parameters .................................
    tol = 1e-14;  % Tolerance
    marg = tol;   % Margins for polygon frame
    is_ab = 0;    % Default A*X+B  mode

    % Handle input ............................................
    if nargin==0, help intsecpl, return, end
    if nargin < 3
        error(' Not enough input arguments')
    end
    if nargin<5, trace = 0; end
    if nargin==4  % Check if 4-th arg is trace
        if max(size(yl))==1, trace = yl;  is_ab = 1; end
    end
    if nargin==3, is_ab = 1; end
    trace = trace | nargin<2;
    if length(xv)~=length(yv)
        error(' Vectors X, Y must have the same size')
    end

    % Auxillary ...........
    xv = [xv(:); xv(1)];
    yv = [yv(:); yv(1)];
    ii = find(abs(diff(xv))<tol & abs(diff(yv))<tol);
    xv(ii) = []; yv(ii) = [];
    nv = length(xv);
    ov = ones(nv-1,1);

    % Polygon frame
    lim = [min(xv)-marg max(xv)+marg min(yv)-marg max(yv)+marg];
    % Estimate for diameter
    d = sqrt((lim(2)-lim(1)).^2+(lim(4)-lim(3)).^2);

    % Form line segment depending on how line is specified
    if is_ab       % A*X+B  mode ...................

        xl = xl(:);
        if length(xl)<2
            error(' Line is specified by at least two parameters')
        end

        a = xl(1); b = xl(2);
        xl = [lim(1)-1 lim(2)+1];
        yl = a*xl+b;

    else          % [X1 X2],[Y1 Y2]  mode ..........

        x0 = (xl(1)+xl(2))/2;
        y0 = (yl(1)+yl(2))/2;
        dx = xl(2)-x0;
        dy = yl(2)-y0;
        dl = max(abs([lim(1:2)-x0 lim(3:4)-y0]));
        d = max(d,dl);
        dl = sqrt(dx.^2+dy.^2);
        dl = max(d,d/dl);
        dx = dx*dl; dy = dy*dl;
        xl = [x0-dx x0+dx];
        yl = [y0-dy y0+dy];

    end

    % Find intersecting segments ..............................
    is = iscross([xv(1:nv-1) xv(2:nv)]',[yv(1:nv-1) yv(2:nv)]',...
        xl(ov,:)',yl(ov,:)',0);

    % Quick exit if no intersections .........................
    if ~any(is)
        if trace
            % Intersection with polygon frame
            [xl,yl] = intsecpl(lim([1 2 2 1]),lim([3 3 4 4]),xl,yl);
            %plot(xv,yv,xl,yl) % Plotting itself
        end
        xo = []; yo = [];
        return
    end

    % For segments touching the line (is==.5) find whether pairs of
    % successive segments are on the same side
    ii = find(is==.5)';
    if ii~=[]
        xo = [ii-1 ii+1];
        xo = xo+(nv-1)*(xo<1);
        yo = iscross([xv(xo(:,1)) xv(xo(:,2))]',[yv(xo(:,1)) yv(xo(:,2))]',...
            xl,yl,tol);
        ii = ii(find(yo==1));
        is(ii) = zeros(size(ii));
    end

    % Calculate intersection coordinates ......................
    ii = find(is);
    oi = ones(size(ii));
    [xo,yo] = intsecl([xv(ii) xv(ii+1)]',[yv(ii) yv(ii+1)]',...
        xl(oi,:)',yl(oi,:)');

    dx = find(~finite(xo));
    xo(dx) = []; yo(dx) = [];
    ii(dx) = []; oi(dx) = [];

    % Sort intersections along the line ..........
    xo = xo(:); yo = yo(:);
    if any(diff(xo))
        [xo,ii] = sort(xo);
        yo = yo(ii);

    else

        [yo,ii] = sort(yo);
        xo = xo(ii);
    end

    % Exclude repeated points (degenerate cases)
    % ///////// It seems that this is not needed, figure this
    % later ///////////////////
    if length(ii)>1 & 0  % Do not execute
        xx = [xo yo];
        yy = diff(xx)';
        ii = [1 find(any(abs(yy)>tol))+1];
        xo = xx(ii,1); yo = xx(ii,2);
        oi = ones(size(xo));
    end

    % Plotting ................................................
    if trace
        oi(3:2:length(oi)) = oi(3:2:length(oi))+1;
        oi = cumsum(oi);
        len = max(oi);
        xp = nan*ones(len,1); yp = xp;
        xp(oi) = xo;
        yp(oi) = yo;

        % Intersection with polygon frame
        [xl,yl] = intsecpl(lim([1 2 2 1]),lim([3 3 4 4]),xl,yl);

        plot(xv,yv,xl,yl,xp,yp) % Plotting itself
    end
end

function [is,S] = isintpl(x1,y1,x2,y2)

    % ISINTPL Check for intersection of polygons.
    %   [IS,S] = ISINTPL(X1,Y1,X2,Y2) Accepts coordinates
    %   X1,Y1 and X2, Y2 of two polygons. Returns scalar
    %   IS equal to 1 if these polygons intersect each other
    %   and 0 if they do not.
    %   Also returns Boolean matrix S of the size length(X1)
    %   by length(X2) so that {ij} element of which is 1 if
    %   line segments i to i+1 of the first polygon and
    %   j to j+1 of the second polygon intersect each other,
    %   0 if they don't and .5 if they "touch" each other.

    %  Copyright (c) 1995 by Kirill K. Pankratov,
    %       kirill@plume.mit.edu.
    %       06/20/95, 08/25/95

    % Handle input ...................................
    if nargin==0, help isintpl, return, end
    if nargin<4
        error(' Not enough input arguments ')
    end

    % Make column vectors and check sizes ............
    x1 = x1(:); y1 = y1(:);
    x2 = x2(:); y2 = y2(:);
    l1 = length(x1);
    l2 = length(x2);
    if length(y1)~=l1 | length(y2)~=l2
        error('(X1,Y1) and (X2,Y2) must pair-wise have the same length.')
    end

    % Quick exit if empty
    if l1<1 | l2<1, is = []; S = []; return, end

    % Check if their ranges are intersecting .........
    lim1 = [min(x1) max(x1)]';
    lim2 = [min(x2) max(x2)]';
    isx = interval(lim1,lim2);  % X-ranges
    lim1 = [min(y1) max(y1)]';
    lim2 = [min(y2) max(y2)]';
    isy = interval(lim1,lim2);  % Y-ranges
    is = isx & isy;
    S = zeros(l2,l1);

    if ~is, return, end  % Early exit if disparate limits

    % Indexing .......................................
    [i11,i12] = meshgrid(1:l1,1:l2);
    [i21,i22] = meshgrid([2:l1 1],[2:l2 1]);
    i11 = i11(:); i12 = i12(:);
    i21 = i21(:); i22 = i22(:);

    % Calculate matrix of intersections ..............
    S(:) = iscross([x1(i11) x1(i21)]',[y1(i11) y1(i21)]',...
        [x2(i12) x2(i22)]',[y2(i12) y2(i22)]')';

    S = S';
    is  = any(any(S));
end

function [xo,yo] = intsecl(x1,y1,x2,y2,tol)

    % INTSECL Intersection coordinates of two line segments.
    %       [XI,YI] = INTSECL(X1,Y1,X2,Y2) where all 4
    %   arguments are 2 by N matrices with coordinates
    %   of ends of N pairs of line segments (so that
    %       the command such as PLOT(X1,Y1,X2,Y2) will plot
    %       these pairs of lines).
    %       Returns 1 by N vectors XI and YI consisting of
    %       coordinates of intersection points of each of N
    %   pairs of lines.
    %
    %   Special cases:
    %   When a line segment is degenerate into a point
    %   and does not lie on line through the other segment
    %   of a pair returns XI=NaN while YI has the following
    %   values: 1 - when the first segment in a pair is
    %   degenerate, 2 - second segment, 0 - both segments
    %   are degenerate.
    %   When a pair of line segments is parallel, returns
    %   XI = Inf while YI is 1 for coincident lines,
    %   0 - for parallel non-coincident ones.
    %   INTSECL(X1,Y1,X2,Y2,TOL) also specifies tolerance
    %   in detecting coincident points in different line
    %   segments.

    %  Copyright (c) 1995 by Kirill K. Pankratov
    %       kirill@plume.mit.edu
    %       04/15/94, 08/14/94, 05/10/95, 08/23/95

    % Defaults and parameters .........................
    tol_dflt = 0; % Tolerance for coincident points
    is_chk = 1;   % Check input arguments

    % Handle input ....................................
    if nargin==0, help intsecl, return, end
    if nargin<4           % Check if all 4 entered
        error('  Not enough input arguments')
    end
    if nargin<5, tol = tol_dflt; end
    if tol < 0, is_chk = 0; tol = 0; end

    % Check the format of arguments .......
    if is_chk
        [x1,y1,x2,y2] = linechk(x1,y1,x2,y2);
    end

    % Auxillary
    o2 = ones(2,1);
    i_pt1 = []; i_pt2 = []; i_pt12 = [];

    % Make first points origins ...........
    xo = x1(1,:);
    yo = y1(1,:);
    x2 = x2-xo(o2,:);
    y2 = y2-yo(o2,:);

    % Differences of first segments .......
    a = x1(2,:)-x1(1,:);
    b = y1(2,:)-y1(1,:);
    s = sqrt(a.^2+b.^2);  % Lengths of first segments
    i_pt1 = find(~s);
    s(i_pt1) = ones(size(i_pt1));
    rr = rand(size(i_pt1));
    a(i_pt1) = cos(rr);
    b(i_pt1) = sin(rr);

    % Normalize by length .................
    a = a./s; b = b./s;

    % Rotate coordinates of the second segment
    tmp = x2.*a(o2,:)+y2.*b(o2,:);
    y2 = -x2.*b(o2,:)+y2.*a(o2,:);
    x2 = tmp;

    % Calculate differences in second segments
    s = x2(2,:)-x2(1,:);
    tmp = y2(2,:)-y2(1,:);
    cc = tmp(i_pt1);

    % Find some degenerate cases .......................

    % Find zeros in differences
    izy2 = find(~tmp);
    tmp(izy2) = ones(size(izy2));

    % Find degenerate and parallel segments
    bool = ~s(izy2);
    i_par = izy2(~bool);
    i_pt2 = izy2(bool);

    bool = abs(y2(1,i_pt2))<=tol;
    i_pt2_off = i_pt2(~bool);
    i_pt2_on = i_pt2(bool);

    if ~isempty(i_par)
        bool = abs(y2(1,i_par))<=tol;
        i_par_off = i_par(~bool);
        i_par_on = i_par(bool);
    end

    % Calculate intercept with rotated x-axis ..........
    tmp = s./tmp;   % Slope
    tmp = x2(1,:)-y2(1,:).*tmp;

    % Rotate and translate back to original coordinates
    xo = tmp.*a+xo;
    yo = tmp.*b+yo;

    % Mark special cases ...................................
    % First segments are degenerate to points
    if ~isempty(i_pt1)
        bool = ~s(i_pt1) & ~cc;
        i_pt12 = i_pt1(bool);
        i_pt1 = i_pt1(~bool);

        bool = abs(tmp(i_pt1))<=tol;
        i_pt1_on = i_pt1(bool);
        i_pt1_off = i_pt1(~bool);

        xo(i_pt1_on) = x1(1,i_pt1_on);
        yo(i_pt1_on) = y1(1,i_pt1_on);

        oo = ones(size(i_pt1_off));
        xo(i_pt1_off) = nan*oo;
        yo(i_pt1_off) = oo;
    end

    % Second segments are degenerate to points ...
    if ~isempty(i_pt2)
        oo = ones(size(i_pt2_off));
        xo(i_pt2_off) = nan*oo;
        yo(i_pt2_off) = 2*oo;
    end

    % Both segments are degenerate ...............
    if ~isempty(i_pt12)
        bool = x1(i_pt12)==xo(i_pt12);
        i_pt12_on = i_pt12(bool);
        i_pt12_off = i_pt12(~bool);

        xo(i_pt12_on) = x1(1,i_pt12_on);
        yo(i_pt12_on) = y1(1,i_pt12_on);

        oo = ones(size(i_pt12_off));
        xo(i_pt12_off) = nan*oo;
        yo(i_pt12_off) = 0*oo;
    end

    % Parallel segments .........................
    if ~isempty(i_par)
        oo = ones(size(i_par_on));
        xo(i_par_on) = inf*oo;
        yo(i_par_on) = oo;

        oo = ones(size(i_par_off));
        xo(i_par_off) = inf*oo;
        yo(i_par_off) = 0*oo;
    end

end

function [x1,y1,x2,y2] = linechk(x1,y1,x2,y2)

    % LINECHK Input checking for line segments.

    %  Copyright (c) 1995 by Kirill K. Pankratov
    %       kirill@plume.mit.edu
    %       08/22/95,

    % String for transposing
    str = ['x1=x1'';'; 'y1=y1'';'; 'x2=x2'';'; 'y2=y2'';'];

    % Sizes
    sz = [size(x1); size(y1); size(x2); size(y2)]';
    psz = prod(sz);

    % Check x1, y1
    if psz(1)~=psz(2)
        error('  Arguments  x1 and y1 must have the same size')
    end

    % Check x2, y2
    if psz(3)~=psz(3)
        error('  Arguments  x2 and y2 must have the same size')
    end

    % Check if any arguments are less than 2 by 1
    if any(max(sz)<2)
        error('  Arguments  x1, y1, x2, y2 must be at least 2 by 1 vectors')
    end

    % Check if no size is equal to 2
    if any(all(sz~=2))
        error('  Arguments  x1, y1, x2, y2 must be 2 by 1 vectors')
    end

    % Find aruments to be transposed .............................
    ii = find(sz(1,:)~=2);
    for jj = 1:length(ii)
        eval(str(ii(jj),:)); % Transpose if neccessary
    end
    sz(:,ii) = flipud(sz(:,ii));

    % If vectors, extend to 2 by n matrices ......................
    n = max(sz(2,:));
    on = ones(1,n);
    if sz(2,1)<n, x1 = x1(:,on); end
    if sz(2,2)<n, y1 = y1(:,on); end
    if sz(2,3)<n, x2 = x2(:,on); end
    if sz(2,4)<n, y2 = y2(:,on); end
end

author: Allopart

repo: rbpf-gmapping


text:
function [flag, u, v, t] = rayTriangleIntersection (o, d, p0, p1, p2)
% Ray/triangle intersection using the algorithm proposed by Möller and Trumbore (1997).
%
% Input:
%    o : origin.
%    d : direction.
%    p0, p1, p2: vertices of the triangle.
% Output:
%    flag: (0) Reject, (1) Intersect.
%    u,v: barycentric coordinates.
%    t: distance from the ray origin.
% Author:
%    Jesus Mena

    epsilon = 0.00001;

    e1 = p1-p0;
    e2 = p2-p0;
    q  = cross(d,e2);
    a  = dot(e1,q); % determinant of the matrix M

    if (a>-epsilon && a<epsilon)
        % the vector is parallel to the plane (the intersection is at infinity)
        [flag, u, v, t] = deal(0,0,0,0);
        return;
    end;

    f = 1/a;
    s = o-p0;
    u = f*dot(s,q);

    if (u<0.0)
        % the intersection is outside of the triangle
        [flag, u, v, t] = deal(0,0,0,0);
        return;
    end;

    r = cross(s,e1);
    v = f*dot(d,r);

    if (v<0.0 || u+v>1.0)
        % the intersection is outside of the triangle
        [flag, u, v, t] = deal(0,0,0,0);
        return;
    end;

    t = f*dot(e2,r); % verified!
    flag = 1;
    return;
end


text:
%  myKnownRegimesFilter - filtering procedure for state-space models with
%  known switching dates
%
%  ::
%
%
%    [LogLik,Incr,retcode,Filters]=myKnownRegimesFilter(syst,y,U,z,options,regs)
%
%  Args:
%
%     - **syst** [struct]: structure provided by dsge.filter
%
%     - **y** [matrix]: matrix of data provided by dsge.filter
%
%     - **U** [matrix]: matrix of trends provided by dsge.filter
%
%     - **z** [matrix]: matrix of deterministic terms provided by dsge.filter
%
%     - **options** [struct]: options provided by dsge.filter
%
%     - **regs** [vector]: history of regimes, must be of the same length as
%        **y**
%
%  Returns:
%     :
%
%     - **LogLik** [numeric]: value of the log likelihood
%
%     - **Incr** [vector]: contributions to the likelihood in each period
%
%     - **retcode** [numeric]: flag equal to 0 if there is no problem
%
%     - **Filters** [struct]: structure containing all the filtering
%       information
%
%  Note:
%
%     - If the filter is run on a constant parameter model, then the last
%       input argument need not be specified.
%
%     - If the function is passed through a rise/dsge object, then it should
%       be called as
%       ff=filter(m,'kf_user_algo',{@myKnownRegimesFilter,regs}). That is,
%       RISE will provide all input arguments expect the last one that needs
%       to be provided by the user.
%
%     - It is the responsibility of the user to make sure that the solution
%       of the model is consistent with agents expecting no transition from
%       one regime to another. This occurs when the original model is
%       backward-looking and/or the transition matrix is diagonal.
%
%     See also: myConstantParamFilter
%

author: jmaih


text:
function layer = vl_nnconvt(varargin)
%VL_NNCONVT Additional options for vl_nnconvt (CNN deconvolution)
%   Y = Layer.vl_nnconvt(X, F, B) computes the convolution-transpose of the
%   image X with the filter bank F and biases B. See help vl_nnconvt for
%   more details.
%
%   This method overloads MatConvNet's vl_nnconvt function for Layer
%   objects, so that instead of executing vl_nnconvt, a new Layer object is
%   returned. X, F and B can be other Layers, including Params, or
%   constants.
%
%   The overloaded method accepts the same options as Layer.vl_nnconv.

  layer = vl_nnconv(varargin{:}, 'transpose', true) ;

end

repo: MFT


text:
% std_erpimageplot() - Commandline function to plot cluster ERPimage or channel erpimage.
%
% Usage:
%   >> [STUDY] = std_erpimageplot(STUDY, ALLEEG, key1, val1, key2, val2);
%   >> [STUDY data times freqs pgroup pcond pinter] = ...
%                std_erpimageplot(STUDY, ALLEEG ...);
% Inputs:
%   STUDY      - EEGLAB STUDY set comprising some or all of the EEG datasets in ALLEEG.
%   ALLEEG     - global EEGLAB vector of EEG structures for the datasets in the STUDY.
%                Note: ALLEEG for a STUDY set is typically created using load_ALLEEG().
%
% Additional help:
% Inputs and output of this function are strictly identical to the std_erspplot().
% See the help message of this function for more information.
%
% See also: std_erspplot()
%
% Authors: Arnaud Delorme, UCSD/CERCO, August, 2011-

% Copyright (C) Arnaud Delorme, arno@ucsd.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 [STUDY, allitc, alltimes, allfreqs, pgroup, pcond, pinter, events] = std_erpimageplot(STUDY, ALLEEG, varargin)

if nargin < 2
    help std_erpimageplot;
    return;
end;

events = [];
STUDY = pop_erpimparams(STUDY, 'default');
if strcmpi(STUDY.etc.erpimparams.concatenate, 'off')

    % use std_erspplot for stats when concatenate is off
    [STUDY allitc alltimes allfreqs pgroup pcond pinter events] = std_erspplot(STUDY, ALLEEG, 'datatype', 'erpim', varargin{:});

else
    params = STUDY.etc.erpimparams;
    [ opt moreparams ] = finputcheck( varargin, { ...
        'design'      'integer' [] STUDY.currentdesign;
        'topotime'    'real'    [] params.topotime;
        'topotrial'  'real'     [] params.topotrial;
        'timerange'   'real'    [] params.timerange;
        'trialrange'  'real'    [] params.trialrange;
        'colorlimits' 'real'    [] params.colorlimits; % ERPimage
        'statistics'  'string'  [] params.statistics;
        'groupstats'  'string'  [] params.groupstats;
        'condstats'   'string'  [] params.condstats;
        'threshold'   'real'    [] params.threshold;
        'naccu'       'integer' [] params.naccu;
        'mcorrect'    'string'  [] params.mcorrect;
        'erpimageopt' 'cell'    [] params.erpimageopt;
        'concatenate'  'string' { 'on','off' }  params.concatenate;
        'channels'    'cell'    []              {};
        'clusters'    'integer' []              [];
        'comps'       {'integer','string'}  []              []; % for backward compatibility
        'plotsubjects' 'string' { 'on','off' }  'off';
        'plotmode'    'string' { 'normal','condensed','none' }  'normal';
        'subject'     'string'  []              '' }, 'std_erpimageplot', 'ignore');

    if ~isempty(opt.topotime) && ~isempty(opt.topotrial)
        error('Cannot plot topography when ERP-image is in trial concatenation mode');
    end;
    if ~isempty(opt.trialrange)
        error('Cannot select trial range when ERP-image is in trial concatenation mode');
    end;
    if strcmpi(opt.groupstats, 'on') || strcmpi(opt.condstats, 'on')
        disp('Warning: cannot perform statistics when ERP-image is in trial concatenation mode');
    end;

    % options
    if ~isempty(opt.colorlimits), options = { 'caxis' opt.colorlimits opt.erpimageopt{:} };
    else                          options = { 'cbar' 'on' opt.erpimageopt{:}  };
    end;

    if ~isempty(opt.channels)
        [STUDY allerpimage alltimes alltrials tmp events] = std_readersp(STUDY, ALLEEG, 'channels', opt.channels, 'infotype', 'erpim', 'subject', opt.subject, ...
            'concatenate', 'on', 'timerange', opt.timerange, 'design', opt.design);

        % get figure title
        % ----------------
        locs = eeg_mergelocs(ALLEEG.chanlocs);
        locs = locs(std_chaninds(STUDY, opt.channels));
        allconditions = STUDY.design(opt.design).variable(1).value;
        allgroups     = STUDY.design(opt.design).variable(2).value;
        alltitles = std_figtitle('condnames', allconditions, 'cond2names', allgroups, 'chanlabels', { locs.labels }, ...
            'subject', opt.subject, 'valsunit', 'ms', 'datatype', 'ERPIM');

        figure;
        for iCond = 1:length(allconditions)
            for iGroup = 1:length(allgroups)
                tmpevents = events{iCond, iGroup};
                if isempty(tmpevents), tmpevents = zeros(1, size(allerpimage{iCond, iGroup},2)); end;
                subplot(length(allconditions), length(allgroups), (iCond-1)*length(allgroups) + iGroup);

                % use color scale for last plot
                if ~isempty(opt.colorlimits) && iCond == length(allconditions) && iGroup == length(allgroups)
                    options = { options{:} 'cbar' 'on' };
                end;
                erpimage(allerpimage{iCond, iGroup}, tmpevents, alltimes, alltitles{iCond, iGroup}, params.smoothing, params.nlines, options{:});
            end;
        end;

    else
        for cInd = 1:length(opt.clusters)
            [STUDY allerpimage alltimes alltrials tmp events] = std_readersp(STUDY, ALLEEG, 'clusters', opt.clusters(cInd), 'infotype', 'erpim', 'subject', opt.subject, ...
                'concatenate', 'on', 'timerange', opt.timerange, 'design', opt.design);
            % get figure title
            % ----------------
            locs = eeg_mergelocs(ALLEEG.chanlocs);
            locs = locs(std_chaninds(STUDY, opt.channels));
            allconditions = STUDY.design(opt.design).variable(1).value;
            allgroups     = STUDY.design(opt.design).variable(2).value;
            alltitles = std_figtitle('condnames', allconditions, 'cond2names', allgroups, 'clustname', STUDY.cluster(opt.clusters(cInd)).name, ...
                'subject', opt.subject, 'valsunit', 'ms', 'datatype', 'ERPIM');

            figure;
            for iCond = 1:length(allconditions)
                for iGroup = 1:length(allgroups)
                    tmpevents = events{iCond, iGroup};
                    if isempty(tmpevents), tmpevents = zeros(1, size(allerpimage{iCond, iGroup},2)); end;
                    subplot(length(allconditions), length(allgroups), (iCond-1)*length(allgroups) + iGroup);

                    % use color scale for last plot
                    if ~isempty(opt.colorlimits) && iCond == length(allconditions) && iGroup == length(allgroups)
                        options = { options{:} 'cbar' 'on' };
                    end;
                    erpimage(allerpimage{iCond, iGroup}, tmpevents, alltimes, alltitles{iCond, iGroup}, params.smoothing, params.nlines, options{:});
                end;
            end;
        end;
    end;
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ x, obj_val,exitflag ] = QCQP_LB2( H_wave,y_wave,N,l,u,obj_prev)
%UNTITLED9 Summary of this function goes here
%   Detailed explanation goes here
Q = 2*H_wave*H_wave.';
f = -2*H_wave*y_wave;
c = y_wave.'*y_wave;
A = zeros(N,2*N);
for ii = 1:N
    A(ii,ii) = cos((l(ii)+u(ii))/2);
    A(ii,ii+N) = sin((l(ii)+u(ii))/2);         %Linear Constraint
end

function [y,grady] = quadobj(x,Q,f,c)         %Objective Function
    y = 1/2*x.'*Q*x+f.'*x+c;
    if nargout>1
        grady = Q*x + f;
    end
end

function [y_con,yeq,grady,gradyeq] = quadconstr(x)      %Nonlinear Constraint
    for nn = 1:N
        y_con(nn) = x(nn)^2+x(nn+N)^2 - 1;
    end
    y_con(N+1) = 1/2*x.'*Q*x+f.'*x+c-obj_prev;
    yeq = [];
    if nargout>2
       gradyeq = zeros(length(x),N+1);
       for nn = 1:N
           grady(nn,nn) = 2*x(nn);
           grady(nn+N,nn) = 2*x(nn+N);
       end
       grady(:,N+1) = Q*x + f;
    end
    gradyeq = [];
end

Acon = -A;

for ii = 1:N
    b(ii) = -cos((u(ii)-l(ii))/2);
end

% options = optimoptions(@fmincon,'Algorithm','sqp',...
%     'SpecifyObjectiveGradient',true,'SpecifyConstraintGradient',true,...
%     'Display','off','OptimalityTolerance',1e-6,'MaxIterations',400);

options = optimoptions('fmincon','Algorithm','sqp', 'GradObj','on');

fun = @(x)quadobj(x,Q,f,c);
nonlconstr = @(x)quadconstr(x);
x1 = ones(2*N,1);
[x,ignore,exitflag] = fmincon(fun,x1,Acon,b,[],[],[],[],nonlconstr,options);
obj_val = objval_func(x,H_wave,y_wave);
end

author: yuanhao-cui

repo: Must-Reading-on-ISAC

path: github-repos/MATLAB/yuanhao-cui-Must-Reading-on-ISAC/Must-Reading-on-ISAC-34cd6615c52ebca121428a979e756c608b195040/Codes/Fan2018TSP/Codes for DFRC Waveform Design/Constant Modulus/QCQP_LB2.m


text:
function [node,elem]=readasc(fname)
%
% [node,elem]=readasc(fname)
%
% read FreeSurfer ASC mesh format
%
% author: Qianqian Fang <q.fang at neu.edu>
% date: 2009/04/02
%
% input:
%      fname: name of the asc file
%
% output:
%      node: node positions of the mesh
%      elem: element list of the mesh
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

node=[];
elem=[];
fid=fopen(fname,'rt');
if(fid==-1)
        error(['can not read file ' fname]);
end

line=fgetl(fid); % the first line is #!ascii ....
dim=fscanf(fid,'%d',2);
node=fscanf(fid,'%f',[4,dim(1)])';
elem=fscanf(fid,'%f',inf);
fclose(fid);

if(length(elem)==4*dim(2))
    elem=reshape(elem,[4,dim(2)])';
elseif(length(elem)==8*dim(2))
    elem=reshape(elem,[8,dim(2)])';
end

if(~any(node(:,end)))
	node=node(:,1:end-1);
end
if(~any(elem(:,end)))
        elem=elem(:,1:end-1);
end

elem=elem+1;

author: fangq


text:
function value = p00_param_default ( test )

%*****************************************************************************80
%
%% P00_PARAM_DEFAULT defaults the parameters of any test problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer TEST, the test problem index.
%
%    Output, integer VALUE, the number of problem parameters.
%
  if ( test == 1 )
    value = 0;
  elseif ( test == 2 )
    value = 0;
  elseif ( test == 3 )
    value = 0;
  elseif ( test == 4 )
    value = 0;
  elseif ( test == 5 )
    value = 0;
  elseif ( test == 6 )
    value = 0;
  elseif ( test == 7 )
    value = 0;
  elseif ( test == 8 )
    value = 0;
  elseif ( test == 9 )
    value = 0;
  elseif ( test == 10 )
    value = 0;
  elseif ( test == 11 )
    value = 0;
  elseif ( test == 12 )
    value = 0;
  elseif ( test == 13 )
    value = 0;
  elseif ( test == 14 )
    value = 0;
  elseif ( test == 15 )
    value = p15_param ( 'DEFAULT' );
  elseif ( test == 16 )
    value = p16_param ( 'DEFAULT' );
  elseif ( test == 17 )
    value = p17_param ( 'DEFAULT' );
  elseif ( test == 18 )
    value = p18_param ( 'DEFAULT' );
  elseif ( test == 19 )
    value = p19_param ( 'DEFAULT' );
  elseif ( test == 20 )
    value = p20_param ( 'DEFAULT' );
  elseif ( test == 21 )
    value = 0;
  elseif ( test == 22 )
    value = 0;
  elseif ( test == 23 )
    value = 0;
  elseif ( test == 24 )
    value = 0;
  elseif ( test == 25 )
    value = 0;
  elseif ( test == 26 )
    value = p26_param ( 'DEFAULT' );
  elseif ( test == 27 )
    value = 0;
  elseif ( test == 28 )
    value = 0;
  elseif ( test == 29 )
    value = 0;
  elseif ( test == 30 )
    value = 0;
  elseif ( test == 31 )
    value = p31_param ( 'DEFAULT' );
  elseif ( test == 32 )
    value = p32_param ( 'DEFAULT' );
  elseif ( test == 33 )
    value = p33_param ( 'DEFAULT' );
  elseif ( test == 34 )
    value = p34_param ( 'DEFAULT' );
  elseif ( test == 35 )
    value = p35_param ( 'DEFAULT' );
  elseif ( test == 36 )
    value = 0;
  elseif ( test == 37 )
    value = p37_param ( 'DEFAULT' );
  elseif ( test == 38 )
    value = p38_param ( 'DEFAULT' );
  elseif ( test == 39 )
    value = p39_param ( 'DEFAULT' );
  elseif ( test == 40 )
    value = p40_param ( 'DEFAULT' );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_PARAM_DEFAULT - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem index TEST = %d\n', test );
    error ( 'P00_PARAM_DEFAULT - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function s=mergestruct(s1,s2)
%
% s=mergestruct(s1,s2)
%
% merge two struct objects into one
%
% authors:Qianqian Fang (q.fang <at> neu.edu)
% date: 2012/12/22
%
% input:
%      s1,s2: a struct object, s1 and s2 can not be arrays
%
% output:
%      s: the merged struct object. fields in s1 and s2 will be combined in s.
%
% license:
%     BSD or GPL version 3, see LICENSE_{BSD,GPLv3}.txt files for details
%
% -- this function is part of jsonlab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

if(~isstruct(s1) || ~isstruct(s2))
    error('input parameters contain non-struct');
end
if(length(s1)>1 || length(s2)>1)
    error('can not merge struct arrays');
end
fn=fieldnames(s2);
s=s1;
for i=1:length(fn)
    s=setfield(s,fn{i},getfield(s2,fn{i}));
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
%% Copyright (C) 2019 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod  @@sym ezmesh (@var{z})
%% @defmethodx @@sym ezmesh (@var{f1}, @var{f2}, @var{f3})
%% @defmethodx @@sym ezmesh (@dots{}, @var{dom})
%% @defmethodx @@sym ezmesh (@dots{}, @var{N})
%% Simple 3D wireframe mesh plots of symbolic expressions.
%%
%% Example 3D surface mesh plot:
%% @example
%% @group
%% syms x y
%% z = sin(2*x)*sin(y)
%%   @result{} z = (sym) sin(2⋅x)⋅sin(y)
%% ezmesh(z)                                    % doctest: +SKIP
%% @end group
%% @end example
%%
%% Example parametric mesh of a Möbius strip:
%% @example
%% @group
%% syms u v
%% x = (1+v*cos(u/2))*cos(u)
%%   @result{} x = (sym)
%%       ⎛     ⎛u⎞    ⎞
%%       ⎜v⋅cos⎜─⎟ + 1⎟⋅cos(u)
%%       ⎝     ⎝2⎠    ⎠
%% y = (1+v*cos(u/2))*sin(u);
%% z = v*sin(u/2);
%%
%% ezmesh(x, y, z, [0 2*pi -0.5 0.5], 32)       % doctest: +SKIP
%% axis equal
%% @end group
%% @end example
%%
%% See help for the (non-symbolic) @code{ezmesh}, which this
%% routine calls after trying to convert sym inputs to
%% anonymous functions.
%%
%% @seealso{ezmesh, @@sym/ezsurf, @@sym/ezplot, @@sym/function_handle}
%% @end defmethod

function varargout = ezmesh(varargin)

  % first input is handle, shift
  if (ishandle(varargin{1}))
    firstpotsym = 2;
  else
    firstpotsym = 1;
  end

  maxnumsym = 3;
  firstsym = [];

  for i = firstpotsym:nargin
    if (isa(varargin{i}, 'sym'))
      if (i < firstpotsym + maxnumsym)
        % one of the fcns to plot, covert to handle fcn

        % Each is function of one var, and its the same var for all
        thissym = symvar(varargin{i});
        assert(length(thissym) <= 2, ...
          'ezmesh: parameterized: functions should have at most two inputs');
        if (isempty(thissym))
          % a number, create a constant function in a dummy variable
          % (0*t works around some Octave oddity on 3.8 and hg Dec 2014)
          thisf = inline(sprintf('%g + 0*t', double(varargin{i})), 't');
          %thisf = @(t) 0*t + double(varargin{i});  % no
        else
          % check variables match (sanity check)
          if (isempty(firstsym))
            firstsym = thissym;
          else
            assert(all(logical(thissym == firstsym)), ...
              'ezmesh: all functions must be in terms of the same variables');
          end
          thisf = function_handle(varargin{i});
        end

        varargin{i} = thisf;

      else
        % plot ranges, etc, convert syms to doubles
        varargin{i} = double(varargin{i});
      end
    end
  end

  h = ezmesh(varargin{:});

  if (nargout)
    varargout{1} = h;
  end

end

%!error <all functions must be in terms of the same variables>
%! syms u v t
%! ezmesh(u*v, 2*u*v, 3*v*t)

%!error <functions should have at most two inputs>
%! syms u v t
%! ezmesh(u*v, 2*u*v, u*v*t)

repo: octsympy


text:
function [subs,vals] = find(t)
%FIND Find subscripts of nonzero elements in a sparse tensor.
%
%   [SUBS,VALS] = FIND(T) returns the subscripts and corresponding
%   values of the nonzero elements of T.
%
%   Note that unlike the standard MATLAB find function for an array,
%   find does not return linear indices. Instead, it returns an M x N
%   array where M is the number of nonzero values and N = ndims(T).
%   Thus, I(k,:) specifies the subscript of value V(k).
%
%   See also SPTENSOR, FIND.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

subs = t.subs;
vals = t.vals;

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt