text:
function testFinalModelsAerts_HN_Old(pathExperiments,nExp,fSetNames)
% - paramCells: One cell of parameters for each fSet
% - nameCells: One cell of names for each fSet

startpath = pwd;
nFset = numel(fSetNames);
for exp = 1:nExp
    cd(fullfile(pathExperiments,['Experiment',num2str(exp)])), pathExperiment = pwd;
    load('training'), load('testing')
    nameOutcomes = fieldnames(training); nOutcomes = numel(nameOutcomes);
    cd('AertsSign'), pathFinalModels = pwd;
    for o = 1:nOutcomes
        for f = 1:nFset
            results = []; results = struct;
            cd(fullfile(pathFinalModels,nameOutcomes{o},fSetNames{f}))
            load('finalModel'), load('coeff'), load('response'), order = 4;
            results.model.Name = finalModel.Name; results.model.coeff = coeff; results.model.order = order;
            results.trainData.data = finalModel.Data; results.trainData.response = response; results.trainData.outcome = training.(nameOutcomes{o}).outcome;

            % TESTING MODEL
            outcome = testing.(nameOutcomes{o}).outcome;
            resp = zeros(numel(outcome),1);
            data = zeros(numel(outcome),order);
            for n = 1:order
                data(:,n) = testing.(nameOutcomes{o}).sign.(fSetNames{f})(:,n);
                resp = resp + coeff(n)*data(:,n);
            end
            resp = resp + coeff(end);
            results.testData.data = data; results.testData.response = resp; results.testData.outcome = outcome;
            [AUC,sensitivity,specificity,accuracy] = calcPerformMetrics(resp,outcome,0);
            results.AUC = AUC;
            results.Sensitivity = sensitivity;
            results.Specificity = specificity;
            results.Accuracy = accuracy;
            cd(pathExperiment), save(['testResultsAerts_',fSetNames{f},'_',nameOutcomes{o}],'results')
        end
    end
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
% eeg_lat2point() - convert latencies in time units relative to the
%                   time locking event of an eeglab() data epoch to
%                   latencies in data points (assuming concatenated epochs).
% Usage:
%       >> [newlat] = eeg_lat2point( lat_array, epoch_array,...
%                                 srate, timelimits, timeunit);
% Inputs:
%   lat_array   - latency array in 'timeunit' units (see below)
%   epoch_array - epoch number for each latency
%   srate       - data sampling rate in Hz
%   timelimits  - [min max] epoch timelimits in 'timeunit' units (see below)
%   timeunit    - time unit relative to seconds. Default is 1 = seconds.
%
% Outputs:
%   newlat      - converted latency values in points assuming concatenated
%                 data epochs (see eeglab() event structure)
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2 Mai 2002
%
% See also: eeg_point2lat(), eeglab()

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

function newlat = eeg_lat2point( lat_array, epoch_array, srate, timewin, timeunit);

if nargin <4
    help eeg_lat2point;
    return;
end;
if nargin <5
	timeunit = 1;
end;

if length(lat_array) ~= length(epoch_array)
	if length(epoch_array)~= 1
		disp('eeg_lat2point: latency and epochs must have the same length'); return;
	else
		epoch_array = ones(1,length(lat_array))*epoch_array;
	end;
end;
if length(timewin) ~= 2
    disp('eeg_lat2point: timelimits must have length 2'); return;
end;
if iscell(epoch_array)
	epoch_array = [ epoch_array{:} ];
end;
if iscell(lat_array)
	lat_array = [ lat_array{:} ];
end

timewin = timewin*timeunit;
pnts = (timewin(2)-timewin(1))*srate+1;
newlat  = (lat_array*timeunit-timewin(1))*srate+1 + (epoch_array-1)*pnts;

author: PatternRecognition

repo: OpenBMI


text:
function rgb = skull

% This returns a predefined color as [red green blue] values
%   red               = [255   0   0]/255;
%   green             = [  0 192   0]/255;
%   blue              = [  0   0 255]/255;
%   magenta           = [255 255   0]/255;
%   cyan              = [  0 255 255]/255;
%   yellow            = [255 255   0]/255;
%   white             = [255 255 255]/255;
%   black             = [  0   0   0]/255;
%
%   skull             = [140  85  85]/255
%   cortex            = [255 213 119]/255;
%   cortex_light      = [199 194 169]/255;
%   cortex_dark       = [100  97  85]/255;
%   skin              = [249 223 192]/255;
%   skin_light        = [249 223 192]/255;
%   skin_medium_light = [225 194 158]/255;
%   skin_medium       = [188 142 106]/255;
%   skin_medium_dark  = [155 102	65]/255;
%   skin_dark         = [ 91  71  61]/255;

rgb = [140 85 85]/255;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% CUBEMAXWELLSADDLE solves Maxwell type equations in a cube using linear order element.
% This is a special case of div u = g being nozero.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

clear; close all;

%% Defacult setting
[node,elem] = cubemesh([-1,1,-1,1,-1,1],1);
%%
pde.J = @(p) [sin(p(:,1)).*cos(p(:,2)).*sin(p(:,3)), ...
              cos(p(:,1)).*sin(p(:,2)).*sin(p(:,3)), ...
              2*cos(p(:,1)).*cos(p(:,2)).*cos(p(:,3))];
pde.exactu = @(p)[0*p(:,1), 0*p(:,2), ...
    cos(p(:,1)).*cos(p(:,2)).*cos(p(:,3))];
pde.g_D = pde.exactu;
pde.curlu = @(p) [-cos(p(:,1)).*sin(p(:,2)).*cos(p(:,3)), ...
              sin(p(:,1)).*cos(p(:,2)).*cos(p(:,3)), 0*p(:,3)];
pde.g = @(p) -cos(p(:,1)).*cos(p(:,2)).*sin(p(:,3));
pde.mu = 1;

%%
bdFlag = setboundary3(node,elem,'Dirichlet');
option.printlevel = 0;
% option.solver = 'direct';
option.solver = 'mg';
option.solver = 'diag';

%% Parameters
maxIt = 3;
N = zeros(maxIt,1);
h = zeros(maxIt,1);
energyErr = zeros(maxIt,1);
L2Err = zeros(maxIt,1);
uIuhErr = zeros(maxIt,1);

%% Finite Element Method
for k = 1:maxIt
    [node,elem,bdFlag] = uniformrefine3(node,elem,bdFlag);
    [soln,eqn,info] = Maxwell1saddle(node,elem,bdFlag,pde,option);
    u = soln.u;
    fprintf('\n\n # of DoFs = %d \n',length(u));
    % compute error
    uI = edgeinterpolate1(pde.exactu,node,eqn.edge);
    energyErr(k) = getHcurlerror3ND1(node,elem,pde.curlu,u);
    L2Err(k) = getL2error3ND1(node,elem,pde.exactu,u);
    uIuhErr(k) = sqrt((u-uI)'*(eqn.A)*(u-uI));
%     L2Err(k) = sqrt((u-uI)'*(eqn.M)*(u-uI));
    fprintf('\n ||curl(u-u_h)|| is %g \n',energyErr(k))
    N(k) = length(u);
    h(k) = 1./(size(node,1)^(1/3)-1);
end

%% Plot convergence rates
figure(1);
showrateh3(h,energyErr,1,'k-+','|| curl (u-u_h) ||',...
           h,uIuhErr,1,'r-+','|| curl (u_I-u_h) ||',...
           h,L2Err,1,'b-+','|| u-u_h||');

repo: ifem


text:
%TESTDATASIZE of datafiles and convert to dataset
%
%	 B = TESTDATASIZE(A,STRING)
%  I = TESTDATASIZE(A,STRING,FALSE)
%  I = TESTDATASIZE(N)
%
% INPUT
%  A         datafile or dataset
%  STRING    'data' (default) or 'features' or 'objects'
%  N         Given data size to be tested
%
% OUTPUT
%  B         Converted dataset
%  I         TRUE:  conversion possible
%            FALSE: conversion not possible
%
% DESCRIPTION
% Depending on the value of PRMEMORY and the size of the datafile A, it is
% converted to a dataset, otherwise an error is generated.
% In case the third parameter is FALSE or the first is a scalar just a test
% is executed. In case of no output arguments an error is generated if
% conversion in impossible.
%
% The parameter STRING controls the type of comparison:
%
% 'data'          PROD(SIZE(A)) < PRMEMORY
% 'objects'       SIZE(A,1).^2  < PRMEMORY
% 'features'      SIZE(A,2).^2  < PRMEMORY
%
% SEE ALSO (PRTools Guide)
% DATASETS, DATAFILES, PRMEMORY

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function b = testdatasize(a,type,flag)

	if nargin < 3,
		flag = 1;
	end

	if nargin < 2
		type = 'data';
  end

  b = true;
	if nargin == 1 & isdouble(a) & numel(a) == 1
    % estimated data size in a
    if a > prmemory
      if nargout == 0
        error(['Dataset too large for memory.' newline ...
          'Size (Mega elements) is ' int2str(round(a/1000000)) ', memory is ' int2str(prmemory/1000000) newline ...
          'Possible solutions:' newline ...
          '- decrease data size' newline ...
          '- increase PRMEMORY, see prmemory' newline ...
          '- consider batch processing, see setbatch']);
      else
        b = false;
      end
    end
    return
  end

	if isdataset(a) | isdouble(a)
		if flag
			b = a;
		end
		return
	end

  % Now we have a datafile
	a = setfeatsize(a,0); % featsize of datafiles is unreliable
	switch type
		case 'data'
			if prod(size(a)) > prmemory
				if flag
					error(['Dataset too large for memory.' newline ...
          'Size is ' int2str(prod(size(a))) ', memory is ' int2str(prmemory) newline ...
          'Possible solutions:' newline ...
          '- decrease data size' newline ...
          '- increase PRMEMORY, see prmemory' newline ...
          '- consider batch processing, see setbatch']);
				else
					b = false;
				end
			end
		case 'objects'
			if size(a,1).^2 > prmemory
				if flag
					error(['Number of objects too large for memory.' newline ...
          'Size is ' int2str(size(a,1).^2) ', memory is ' int2str(prmemory) newline ...
          'Possible solutions:' newline ...
          '- decrease data size' newline ...
          '- increase PRMEMORY, see prmemory' newline ...
          '- consider batch processing, see setbatch']);
				else
					b = false;
				end
			end
		case 'features'
			if size(a,2).^2 > prmemory
				if flag
					error(['Number of features too large for memory.' newline ...
          'Size is ' int2str(size(a,2).^2) ', memory is ' int2str(prmemory) newline ...
          'Possible solutions:' newline ...
          '- decrease data size' newline ...
          '- increase PRMEMORY, see prmemory' newline ...
          '- consider batch processing, see setbatch']);
				else
					b = false;
				end
			end
		otherwise
			error('Unknown test requested')
	end
	if nargout > 0
		if flag
			b = prdataset(a);
		end
	end

return

author: marianux

repo: ecg-kit

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


text:
function VBA_spm_defaults
% Set the defaults which are used by SPM
%__________________________________________________________________________
%
% If you want to customise some defaults for your installation, do not
% modify this file directly, but create a file named spm_my_defaults.m
% instead, accessible from MATLAB search path; e.g., it can be saved in
% MATLAB Startup Folder: userhome/Documents/MATLAB.
%
% Example: create the following file to change the image file extension:
% ----------- file /home/karl/Documents/MATLAB/spm_my_defauls.m -----------
% global defaults
% defaults.images.format = 'img';
%--------------------------------------------------------------------------
%
% spm_defaults should not be called directly in any script or function
% (apart from SPM internals).
% To load the defaults, use spm('Defaults',modality).
% To get/set the defaults, use spm_get_defaults.
%
%                 ** This file should not be edited **
%__________________________________________________________________________
% Copyright (C) 1994-2014 Wellcome Trust Centre for Neuroimaging

% SPM
% $Id: spm_defaults.m 6221 2014-09-30 20:52:29Z guillaume $

global defaults

% Command Line Usage default
%==========================================================================
defaults.cmdline  = 0;

% User Interface defaults
%==========================================================================
defaults.ui.monitor = NaN;
defaults.ui.colour  = [0.58 0.77 0.57];
defaults.ui.fs      = 14;  % unused
defaults.ui.print   = 'ps';
defaults.renderer   = 'zbuffer';

% File format specific
%==========================================================================
% Note that defaults.analyze.flip is no longer used.  Specifying the
% left/right handedness of the voxel indices is now done entirely by
% spm_flip_analyze_images.m
defaults.images.format  = 'nii'; % options: 'img', 'nii'
defaults.mat.format     = '-v6'; % options: '-mat', '-v6', '-v7.0', '-v7.3'

% Toolboxes defaults
%==========================================================================
defaults.tbx.dir = { fullfile(spm('Dir'),'toolbox') };

% DICOM Import defaults
%==========================================================================
defaults.dicom.root = 'flat'; % Folder hierarchy

% fMRI design defaults
%==========================================================================
defaults.stats.fmri.t   = 16;
defaults.stats.fmri.t0  = 8;
defaults.stats.fmri.hpf = 128;
defaults.stats.fmri.cvi = 'AR(1)';

% Mask defaults
%==========================================================================
defaults.mask.thresh    = 0.8;

% Stats defaults
%==========================================================================
defaults.stats.maxmem      = 2^26;
defaults.stats.maxres      = 64;
defaults.stats.resmem      = false;
defaults.stats.fmri.ufp    = 0.001;  % Upper tail F-probability
defaults.stats.pet.ufp     = 0.05;
defaults.stats.eeg.ufp     = 0.05;
defaults.stats.topoFDR     = 1;
defaults.stats.rft.nonstat = 0;
defaults.stats.results.volume.distmin =  8;
defaults.stats.results.volume.nbmax   =  3;
defaults.stats.results.svc.distmin    =  4;
defaults.stats.results.svc.nbmax      = 16;
defaults.stats.results.mipmat         = {fullfile(spm('dir'),'MIP.mat')};

% Filename prefix defaults
%==========================================================================
defaults.slicetiming.prefix     = 'a';
defaults.realign.write.prefix   = 'r';
defaults.coreg.write.prefix     = 'r';
defaults.unwarp.write.prefix    = 'u';
defaults.normalise.write.prefix = 'w';
defaults.smooth.prefix          = 's';

% Realignment defaults
%==========================================================================
defaults.realign.estimate.quality = 0.9;
defaults.realign.estimate.interp  = 2;
defaults.realign.estimate.wrap    = [0 0 0];
defaults.realign.estimate.sep     = 4;
defaults.realign.estimate.fwhm    = 5;
defaults.realign.estimate.rtm     = 1;
defaults.realign.write.mask       = 1;
defaults.realign.write.interp     = 4;
defaults.realign.write.wrap       = [0 0 0];
defaults.realign.write.which      = [2 1];

% Unwarp defaults
%==========================================================================
defaults.unwarp.estimate.rtm      = 0;
defaults.unwarp.estimate.fwhm     = 4;
defaults.unwarp.estimate.basfcn   = [12 12];
defaults.unwarp.estimate.regorder = 1;
defaults.unwarp.estimate.regwgt   = 1e5;
defaults.unwarp.estimate.foe      = [4 5];
defaults.unwarp.estimate.soe      = [];
defaults.unwarp.estimate.rem      = 1;
defaults.unwarp.estimate.jm       = 0;
defaults.unwarp.estimate.noi      = 5;
defaults.unwarp.estimate.expround = 'Average';
%
% Unwarp uses defaults.realign.write
% defaults for writing.
%

% Coregistration defaults
%==========================================================================
defaults.coreg.estimate.cost_fun = 'nmi';
defaults.coreg.estimate.sep      = [4 2];
defaults.coreg.estimate.tol      = [0.02 0.02 0.02 0.001 0.001 0.001 0.01 0.01 0.01 0.001 0.001 0.001];
defaults.coreg.estimate.fwhm     = [7 7];
defaults.coreg.write.interp      = 4;
defaults.coreg.write.wrap        = [0 0 0];
defaults.coreg.write.mask        = 0;

% Spatial Normalisation defaults
%==========================================================================
defaults.normalise.write.preserve   = 0;
defaults.normalise.write.bb         = [[-78 -112 -70];[78 76 85]];
defaults.normalise.write.vox        = [2 2 2];
defaults.normalise.write.interp     = 4;

% Old Spatial Normalisation defaults
%==========================================================================
defaults.old.normalise.estimate.smosrc  = 8;
defaults.old.normalise.estimate.smoref  = 0;
defaults.old.normalise.estimate.regtype = 'mni';
defaults.old.normalise.estimate.weight  = '';
defaults.old.normalise.estimate.cutoff  = 25;
defaults.old.normalise.estimate.nits    = 16;
defaults.old.normalise.estimate.reg     = 1;
defaults.old.normalise.write.preserve   = 0;
defaults.old.normalise.write.bb         = [[-78 -112 -70];[78 76 85]];
defaults.old.normalise.write.vox        = [2 2 2];
defaults.old.normalise.write.interp     = 1;
defaults.old.normalise.write.wrap       = [0 0 0];
defaults.old.normalise.write.prefix     = 'w';

% VBM Preprocessing defaults
%==========================================================================
defaults.old.preproc.tpm     = cellstr(char(...
    fullfile(spm('Dir'),'toolbox','OldSeg','grey.nii'),...
    fullfile(spm('Dir'),'toolbox','OldSeg','white.nii'),...
    fullfile(spm('Dir'),'toolbox','OldSeg','csf.nii'))); % Prior probability maps
defaults.old.preproc.ngaus          = [2 2 2 4]'; % Gaussians per class
defaults.old.preproc.warpreg        = 1;          % Warping Regularisation
defaults.old.preproc.warpco         = 25;         % Warp Frequency Cutoff
defaults.old.preproc.biasreg        = 0.0001;     % Bias regularisation
defaults.old.preproc.biasfwhm       = 60;         % Bias FWHM
defaults.old.preproc.regtype        = 'mni';      % Affine Regularisation
defaults.old.preproc.fudge          = 5;          % Fudge factor, can not be set in batch GUI
defaults.old.preproc.samp           = 3;          % Sampling distance
defaults.old.preproc.output.GM      = [0 0 1];
defaults.old.preproc.output.WM      = [0 0 1];
defaults.old.preproc.output.CSF     = [0 0 0];
defaults.old.preproc.output.biascor = 1;
defaults.old.preproc.output.cleanup = 0;

% Smooth defaults
%==========================================================================
defaults.smooth.fwhm = [8 8 8];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-Prevent users from making direct calls to spm_defaults
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
persistent runOnce
try
    if ~isdeployed && isempty(runOnce)
        d = dbstack;
        if isempty(intersect({'spm','spm_get_defaults'},{d.name}))
            fprintf(['Direct calls to spm_defauts are deprecated.\n' ...
                'Please use spm(''Defaults'',modality) ' ...
                'or spm_get_defaults instead.\n']);
            runOnce = 1;
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-Execute user-specified defaults files
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
def = defaults;
user_defaults = 'spm_my_defaults.m';
if exist(user_defaults,'file')
    if isdeployed && exist(fullfile(spm('Dir'),user_defaults),'file')
        user_defaults_file = cellstr(fullfile(spm('Dir'),user_defaults));
    else
        user_defaults_file = which(user_defaults,'-ALL');
    end
    for i=1:numel(user_defaults_file)
        try
            spm('run', user_defaults_file{i});
        catch
            lr = lasterror;
            warning(lr.message);
        end
    end
    if spm_check_version('matlab','8.0') >= 0, my_isequaln = @isequaln;
    else my_isequaln = @isequalwithequalnans; end
    if ~my_isequaln(def,defaults)
        fprintf('Defaults settings have been modified by file(s):\n');
        for i=1:numel(user_defaults_file)
            fprintf('  %s\n',user_defaults_file{i});
        end
        fn0 = fieldnames(def);
        mf = fn0(~cellfun(@(x) my_isequaln(def.(x),defaults.(x)),fn0));
        if ~isempty(mf)
            fprintf('Modified fields: ');
            for i=1:numel(mf)
                fprintf('%s ',mf{i});
            end
            fprintf('\n');
        end
    end
end

author: MBB-team

repo: VBA-toolbox


text:
function out = feval(f, x, y)
%FEVAL   Evaluate a SEPARABLEAPPROX at one or more points.
%   FEVAL(F, X, Y) evaluates the SEPARABLEAPPROX F and the point(s) in (X, Y),
%   where X and Y are doubles.
%
%   FEVAL(F, X) evaluates the SEPARABLEAPPROX F at the point(s) in X, where X
%   is a double, interpreted as F(real(X), imag(X)).
%
% See also SUBSREF.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Empty check:
if ( isempty(f) )
    out = [];
    return
end

% Get the low rank representation for f.
[cols, D, rows] = cdr(f);

if ( strcmpi(x, ':') && strcmpi(y, ':') )    % f(:, :)
    % Simply return the SEPARABLEAPPROX:
    out = f;

elseif ( strcmpi(x, ':') && isnumeric( y ) ) % f(:, y)
    % Make evaluation points a vector.
    y = y(:);
    % Evaluate (returns a column chebfun):
    out = feval( cols, y ) * D * rows.';
    % Simplify:
    out = simplify( out, [], 'globaltol' );

elseif ( isnumeric( x ) && strcmpi(y, ':') ) % f(x, :)
    % Make evaluation points a vector.
    x = x( : );
    % Evaluate (returns a row chebfun):
    out = cols * D * feval( rows, x ).';
    % Simplify:
    out = simplify( out, [], 'globaltol' );

elseif ( isnumeric(x) && isnumeric(y) )  % f(x, y)
    if ( ndims(x) >= 3 && isequal(size(x), size(y)) )
        % x and y are tensors. Call from CHEBFUN3.
        sizeX = size(x);
        x = x(:);
        y = y(:);
        out = feval(f, x, y);
        %% RESHAPE FOR OUTPUT:
        out = reshape(out, sizeX);
        return
    end

    takeTranspose = 0;

    % If the evaluation points are derived from meshgrid, then there is a
    % fast way to evaluate a separableApprox. Check for this property.
    if ( min(size(x)) > 1 && all(size(x) == size(y)) && numel(size(x)) == 2 )
        % Check to see if the input is a meshgrid:
        if ( max(max(abs(bsxfun(@minus, x, x(1,:))))) <= 10*eps  && ...
                max(max(abs(bsxfun(@minus, y, y(:,1))))) <= 10*eps )
            % This allows someone to do:
            % [xx,yy] = meshgrid(linspace(-1,1));
            % f(xx,yy)

            x = x(1,:);
            y = y(:,1);

        elseif ( max(max(abs(bsxfun(@minus, y, y(1,:))))) <= 10*eps && ...
                max(max(abs(bsxfun(@minus, x, x(:,1))))) <= 10*eps )
            % This allows someone to do:
            % [yy,xx] = meshgrid(linspace(-1,1));
            % f(xx,yy)

            x = x(:,1);
            y = y(1,:);
            takeTranspose = 1;
        else
            % Evaluate at matrices, but they're not from meshgrid:
            [m,n] = size( x );
            out = zeros( m, n );
            % Unroll the loop that is the longest
            if m > n
                for ii = 1:n
                    out(:,ii) = dot(feval(cols, ...
                        y(:,ii))*D,feval(rows, x(:,ii)),2);
                end
            else
                for jj = 1:m
                    out(jj,:) = dot(feval(cols, ...
                        y(jj,:).')*D,feval(rows, x(jj,:).'),2);
                end
            end
            return
        end
    else
    end

    % Evaluate:
    if ( isvector(x) && ~isscalar(x) && all(size(x) == size(y)) )
        % Determine whether inputs are pure vectors.
        out = feval(cols, y(:)) .* feval(rows, x(:)) * diag(D);
        if ( size(x,1) == 1 )
            takeTranspose = 1;
        end
    else
        out = feval(cols, y(:)) * D * feval(rows, x(:)).';
    end

    % Take transpose:
    if ( takeTranspose )
        out = transpose( out );
    end

else
    error('CHEBFUN:SEPARABLEAPPROX:feval:inputs', ...
        'Unrecognized arguments for evaluation.');

end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% Display Start message
disp('Running Convolutional Codes BER Simulation set')

disp('By default, this script does not run the Monte Carlo simulations.')
disp('It simply loads results from a previous run.')
sim_mode = input('Do you want to (R)un simulations or (L)oad results from disk? [R/L]','s');
if isempty(sim_mode)
    sim_mode = 'l';
end
sim_mode = lower(sim_mode);

%% Plot theoretical curves
[h_fig, h_lines] = ConvCode_BER_Curves;
%% Run Monte Carlo simulations
disp('By default, this script does not run the Monte Carlo simulations.')
disp('It simply loads results from a previous run.  If you are interested in')
disp('running the simulations, edit run_me.m and uncomment the following')
disp('line:  CC_BER = ConvCode_Simulate;')

SNRs = 0:0.5:7;
total_time = tic;
% Pre-allocate
CC_BER = zeros(6,length(SNRs));
figure(h_fig)
drawnow
if strcmp(sim_mode,'r')
    hold on
    % Create place-holder plots
    simLines = semilogy(SNRs, CC_BER,'*-');
    CC_BER = ConvCode_Simulate(SNRs, simLines);
else
    % Load from Disk (NO SIMULATIONS ARE RUN)
    load CC_BER  % Loads simulation results from disk instead of running them
    hold on
    simLines = semilogy(SNRs, CC_BER,'*-');
end
toc(total_time)

%% Plot simulation results
% hold on
% simLines = semilogy(0:0.5:7, CC_BER);

%% Add Legend
% Do MATLAB graphics magic to create concise legend
% See "Controling Legends" in MATLAB doc

% group lines together
simGrp = hggroup('DisplayName','Simulation');
theoGrp = hggroup('DisplayName','Theoretical Bound');
set(simLines,'Parent',simGrp)
set(h_lines,'Parent',theoGrp)
set(get(get(simGrp,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','on'); % Include this hggroup in the legend
set(get(get(theoGrp,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','on'); % Include this hggroup in the legend
legend show

figure1 = h_fig;
annotation(figure1,'textarrow',[0.3489 0.4267],[0.3911 0.4381],...
    'TextEdgeColor','none',...
    'TextBackgroundColor',[1 1 1],...
    'String',{'Soft-decision','Decoding'});

% Create textarrow
annotation(figure1,'textarrow',[0.3911 0.4956],[0.2974 0.3603],...
    'TextEdgeColor','none',...
    'TextBackgroundColor',[1 1 1],...
    'String',{'3-bit quantization'});

% Create textbox
annotation(figure1,'textbox',[0.5985 0.5804 0.1164 0.04873],...
    'String',{'Hard-decision','Decoding'},...
    'HorizontalAlignment','center',...
    'LineStyle','none',...
    'BackgroundColor',[1 1 1]);

author: Sable

repo: mcbench-benchmarks


text:
function value = i4_width ( i )

%*****************************************************************************80
%
%% I4_WIDTH returns the "width" of an integer.
%
%  Example:
%
%        I  VALUE
%    -----  -------
%    -1234    5
%     -123    4
%      -12    3
%       -1    2
%        0    1
%        1    1
%       12    2
%      123    3
%     1234    4
%    12345    5
%
%  Discussion:
%
%    The width of an integer is the number of characters necessary to print it.
%
%    The width of an integer can be useful when setting the appropriate output
%    format for a vector or array of values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the number whose width is desired.
%
%    Output, integer VALUE, the number of characters necessary to represent
%    the integer in base 10, including a negative sign if necessary.
%
  if ( 0 <= i )
    value = i4_log_10 ( i ) + 1;
  else
    value = i4_log_10 ( i ) + 2;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function lambda = kershaw_eigenvalues ( )

%*****************************************************************************80
%
%% KERSHAW_EIGENVALUES returns the eigenvalues of the KERSHAW matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    David Kershaw,
%    The Incomplete Cholesky-Conjugate Gradient Method for the Iterative
%    Solution of Systems of Linear Equations,
%    Journal of Computational Physics,
%    Volume 26, Number 1, January 1978, pages 43-65.
%
%  Parameters:
%
%    Output, real LAMBDA(4,1), the eigenvalues of the matrix.
%
  lambda(1:4,1) = [ ...
    5.828427124746192, ...
    5.828427124746188, ...
    0.171572875253810, ...
    0.171572875253810 ];

  return
end


text:
%  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., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
function amtriant
%AMTRIANT Unit test for the function AMTRIANG.

%	O. Lemoine - February 1996.

N=256; t0=149; T=50;
sig=amtriang(N,t0,T);
if abs(sig(t0)-1)>sqrt(eps),				% sig(t0)=1
 error('amtriang test 1 failed ');
end
[tm,T1]=loctime(sig);
if abs(T-T1)>=sqrt(1/N),				% width
 error('amtriang test 3 failed ');
end
dist=1:min([N-t0,t0-1]);
if any(abs(sig(t0-dist)-sig(t0+dist))>sqrt(eps))~=0, 	% symmetry
 error('amtriang test 3 failed ');
end

N=12; t0=5; T=7;
sig=amtriang(N,t0,T);
if abs(sig(t0)-1)>sqrt(eps),				% sig(t0)=1
 error('amtriang test 4 failed ');
end
[tm,T1]=loctime(sig);
if abs(T-T1)>=sqrt(1/N),				% width
 error('amtriang test 5 failed ');
end
dist=1:min([N-t0,t0-1]);
if any(abs(sig(t0-dist)-sig(t0+dist))>sqrt(eps))~=0, 	% symmetry
 error('amtriang test 6 failed ');
end

N=535; t0=354; T=101;
sig=amtriang(N,t0,T);
if abs(sig(t0)-1)>sqrt(eps),				% sig(t0)=1
 error('amtriang test 7 failed ');
end
[tm,T1]=loctime(sig);
if abs(T-T1)>=sqrt(1/N),				% width
 error('amtriang test 8 failed ');
end
dist=1:min([N-t0,t0-1]);
if any(abs(sig(t0-dist)-sig(t0+dist))>sqrt(eps))~=0, 	% symmetry
 error('amtriang test 9 failed ');
end

author: HeLiangHIT


text:
function [ out ] = MotorImagery_online( eeg, online )
%MOTORIMAGERY_ONLINE Summary of this function goes here
%   Detailed explanation goes here
if ~isfield(online,'train'); warning('training module is not exist'); end
if ~isfield(online,'apply'); warning('applying module is not exist'); end

if ~isfield(online,'option'); disp('applying module is not exist');
    opt={
        'device', 'BrainVision'
        'paradigm', 'MotorImagery'
        'Feedback','off'
        'host', 'JohnKim-PC'
        'port','51244'
        }; % default setting
    opt=opt_cellToStruct(online.option{:});
else
    opt=opt_cellToStruct(online.option{:});
end

in=eeg;
if isfield(online,'train')
    for k=1:length(online.train)
        myFun=str2func(online.train{k});
        switch online.train{k}
            case 'func_csp'
                [out CSP_W, CSP_D]= feval(myFun, in, online.train{k,2},0);
            case 'func_featureExtraction'
                out=feval(myFun, in, online.train{k,2},0);
            case 'classifier_trainClassifier'
                CF_PARAM=feval(myFun, in, online.train{k,2},0);
            otherwise
                out=feval(myFun, in, online.train{k,2},0);
        end
        in=out;
    end
end

switch lower(opt.device)
    case 'brainvision'
        H = bv_open(opt.host);

end

while (true)
    dat=bv_read(H);
    in=dat';
    if ~isempty(dat)
        if isfield(online,'apply')
            for k=1:length(online.apply)
                myFun=str2func(online.apply{k});
                switch online.apply{k}
                    case 'func_projection'
                        [out]=feval(myFun, in, CSP_W);
                    case 'func_featureExtraction'
                        [out]=feval(myFun, in, online.apply{k,2},0);
                    case 'classifier_applyClassifier'
                        [out]=feval(myFun, in, CF_PARAM,0);
                    otherwise
                        out=feval(myFun, in, online.apply{k,2},0);
                end
                in=out;
            end
            if strcmp(opt.Feedback,'on')
            end
        end
    end
end

end


text:
function t_opf_tspopf_pdipm(quiet)
%T_OPF_TSPOPF_PDIPM  Tests for PDIPM-based optimal power flow.

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

if nargin < 1
    quiet = 0;
end

num_tests = 159;

t_begin(num_tests, quiet);

[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;

casefile = 't_case9_opf';
if quiet
    verbose = 0;
else
    verbose = 0;
end

t0 = 'PDIPMOPF : ';
mpopt = mpoption('opf.violation', 1e-6, 'pdipm.gradtol', 1e-8, ...
        'pdipm.comptol', 1e-8, 'pdipm.costtol', 1e-9);
mpopt = mpoption(mpopt, 'out.all', 0, 'verbose', verbose, 'opf.ac.solver', 'PDIPM');

%% turn off some warnings
if have_feature('matlab')
    sing_matrix_warn_id = 'MATLAB:nearlySingularMatrix';
    s2 = warning('query', sing_matrix_warn_id);
    warning('off', sing_matrix_warn_id);
end

if have_feature('pdipmopf')
    %% set up indices
    ib_data     = [1:BUS_AREA BASE_KV:VMIN];
    ib_voltage  = [VM VA];
    ib_lam      = [LAM_P LAM_Q];
    ib_mu       = [MU_VMAX MU_VMIN];
    ig_data     = [GEN_BUS QMAX QMIN MBASE:APF];
    ig_disp     = [PG QG VG];
    ig_mu       = (MU_PMAX:MU_QMIN);
    ibr_data    = (1:ANGMAX);
    ibr_flow    = (PF:QT);
    ibr_mu      = [MU_SF MU_ST];
    ibr_angmu   = [MU_ANGMIN MU_ANGMAX];

    %% get solved AC OPF case from MAT-file
    load soln9_opf;     %% defines bus_soln, gen_soln, branch_soln, f_soln

    %% run OPF
    t = t0;
    [baseMVA, bus, gen, gencost, branch, f, success, et] = runopf(casefile, mpopt);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %% run with automatic conversion of single-block pwl to linear costs
    t = [t0 '(single-block PWL) : '];
    mpc = loadcase(casefile);
    mpc.gencost(2, NCOST) = 2;
    [r, success] = runopf(mpc, mpopt);
    [f, bus, gen, branch] = deal(r.f, r.bus, r.gen, r.branch);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);
    xr = [r.var.val.Va;r.var.val.Vm;r.var.val.Pg;r.var.val.Qg;0;r.var.val.y];
    t_is(r.x, xr, 8, [t 'raw x returned from OPF']);

%     %% get solved AC OPF case from MAT-file
%     load soln9_opf_Plim;       %% defines bus_soln, gen_soln, branch_soln, f_soln
%
%     %% run OPF with active power line limits
%     t = [t0 '(P line lim) : '];
%     mpopt1 = mpoption(mpopt, 'opf.flow_lim', 'P');
%     [baseMVA, bus, gen, gencost, branch, f, success, et] = runopf(casefile, mpopt1);
%     t_ok(success, [t 'success']);
%     t_is(f, f_soln, 3, [t 'f']);
%     t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
%     t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
%     t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
%     t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
%     t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
%     t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
%     t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
%     t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
%     t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
%     t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);
%
%     %% run OPF with active power line limits
%     t = [t0 '(P^2 line lim) : '];
%     mpopt1 = mpoption(mpopt, 'opf.flow_lim', '2');
%     [baseMVA, bus, gen, gencost, branch, f, success, et] = runopf(casefile, mpopt1);
%     t_ok(success, [t 'success']);
%     t_is(f, f_soln, 3, [t 'f']);
%     t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
%     t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
%     t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
%     t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
%     t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
%     t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
%     t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
%     t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
%     t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
%     t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %%-----  test OPF with quadratic gen costs moved to generalized costs  -----
    mpc = loadcase(casefile);
    mpc.gencost = [
        2   1500    0   3   0.11    5   0;
        2   2000    0   3   0.085   1.2 0;
        2   3000    0   3   0.1225  1   0;
    ];
    [baseMVA, bus_soln, gen_soln, gencost, branch_soln, f_soln, success, et] = runopf(mpc, mpopt);
    branch_soln = branch_soln(:,1:MU_ST);

    A = sparse(0,0);
    l = [];
    u = [];
    nb = size(mpc.bus, 1);      % number of buses
    ng = size(mpc.gen, 1);      % number of gens
    thbas = 1;                thend    = thbas+nb-1;
    vbas     = thend+1;       vend     = vbas+nb-1;
    pgbas    = vend+1;        pgend    = pgbas+ng-1;
    qgbas    = pgend+1;       qgend    = qgbas+ng-1;
    nxyz = 2*nb + 2*ng;
    N = sparse((1:ng)', (pgbas:pgend)', mpc.baseMVA * ones(ng,1), ng, nxyz);
    fparm = [ 1    0   0 1;
              1 -100 100 1;
              1  -10  10 1 ];
    H = 2 * spdiags(mpc.gencost(:, 5), 0, ng, ng);
    Cw = mpc.gencost(:, 6);
    mpc.gencost(:, 5:7) = 0;

    %% run OPF with quadratic gen costs moved to generalized costs
    t = [t0 'w/quadratic generalized gen cost : '];
    [r, success] = opf(mpc, A, l, u, mpopt, N, fparm, H, Cw);
    [f, bus, gen, branch] = deal(r.f, r.bus, r.gen, r.branch);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);
    t_is(r.cost.usr, f, 12, [t 'user cost']);

    %%-----  run OPF with legacy costs and deadzone  -----
    load soln9_opf;
    mpc = loadcase(casefile);
    mpc.N = sparse((1:nb)', (vbas:vend)', ones(nb,1), nb, nxyz);
    mpc.fparm = ones(nb,1) * [ 2 1.08 0.02 1e8 ];
    mpc.Cw = ones(nb, 1);
    t = [t0 'w/legacy cost, in deadzone : '];
    r = runopf(mpc, mpopt);
    [f, bus, gen, branch] = deal(r.f, r.bus, r.gen, r.branch);
    t_ok(r.success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);
    t_is(r.cost.usr, 0, 12, [t 'user cost']);

    t = [t0 'w/legacy cost, not in deadzone : '];
    mpc.fparm = ones(nb,1) * [ 2 1.08 0.01 1e8 ];
    r = runopf(mpc, mpopt);
    [f, bus, gen, branch] = deal(r.f, r.bus, r.gen, r.branch);
    t_ok(r.success, [t 'success']);
    t_is(f, 9009.0890, 3, [t 'f']);
    t_is([min(bus(:, VM)) mean(bus(:, VM)) max(bus(:, VM))], ...
        [1.066624, 1.083980, 1.091698], 5, [t 'bus voltage']);
    t_is(r.cost.usr, 1673.065465, 5, [t 'user cost']);

    %%-----  run OPF with extra linear user constraints & costs  -----
    %% single new z variable constrained to be greater than or equal to
    %% deviation from 1 pu voltage at bus 1, linear cost on this z
    %% get solved AC OPF case from MAT-file
    load soln9_opf_extras1;   %% defines bus_soln, gen_soln, branch_soln, f_soln
    A = sparse([1;1;2;2],[10;25;10;25],[-1;1;1;1],2,25);
    u = [Inf; Inf];
    l = [-1; 1];

    N = sparse(1, 25, 1, 1, 25);    %% new z variable only
    fparm = [1 0 0 1];              %% w = r = z
    H = sparse(1,1);                %% no quadratic term
    Cw = 100;

    t = [t0 'w/extra constraints & costs 1 : '];
    [r, success] = opf(casefile, A, l, u, mpopt, N, fparm, H, Cw);
    [f, bus, gen, branch] = deal(r.f, r.bus, r.gen, r.branch);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);
    t_is(r.var.val.z, 0.025419, 6, [t 'user variable']);
    t_is(r.cost.usr, 2.5419, 4, [t 'user cost']);

    %%-----  test OPF with capability curves  -----
    mpc = loadcase('t_case9_opfv2');
    %% remove angle diff limits
    mpc.branch(1, ANGMAX) = 360;
    mpc.branch(9, ANGMIN) = -360;

    %% get solved AC OPF case from MAT-file
    load soln9_opf_PQcap;   %% defines bus_soln, gen_soln, branch_soln, f_soln

    %% run OPF with capability curves
    t = [t0 'w/capability curves : '];
    [baseMVA, bus, gen, gencost, branch, f, success, et] = runopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %%-----  test OPF with angle difference limits  -----
    mpc = loadcase('t_case9_opfv2');
    %% remove capability curves
    mpc.gen(2:3, [PC1, PC2, QC1MIN, QC1MAX, QC2MIN, QC2MAX]) = zeros(2,6);

    %% get solved AC OPF case from MAT-file
    load soln9_opf_ang;   %% defines bus_soln, gen_soln, branch_soln, f_soln

    %% run OPF with angle difference limits
    t = [t0 'w/angle difference limits : '];
    [baseMVA, bus, gen, gencost, branch, f, success, et] = runopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  1, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);
    t_is(branch(:,ibr_angmu ), branch_soln(:,ibr_angmu ),  2, [t 'branch angle mu']);

    %%-----  test OPF with ignored angle difference limits  -----
    %% get solved AC OPF case from MAT-file
    load soln9_opf;   %% defines bus_soln, gen_soln, branch_soln, f_soln

    %% run OPF with ignored angle difference limits
    t = [t0 'w/ignored angle difference limits : '];
    mpopt1 = mpoption(mpopt, 'opf.ignore_angle_lim', 1);
    [baseMVA, bus, gen, gencost, branch, f, success, et] = runopf(mpc, mpopt1);
    %% ang limits are not in this solution data, so let's remove them
    branch(1, ANGMAX) = 360;
    branch(9, ANGMIN) = -360;
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %% angle bounded above by 0, unbounded below
    %% for issue/18
    t = [t0 'w/angle difference limit = 0 : '];
    mpc = loadcase(casefile);
    b = 5;
    mpc.branch(b, ANGMAX) = 0;
    r = runopf(mpc, mpopt);
    t_ok(success, [t 'success']);
    diff = r.bus(r.branch(b, F_BUS), VA) - r.bus(r.branch(b, T_BUS), VA);
    t_is(diff, 0, 5, [t 'angle diff']);

    %%-----  OPF with ref bus not = bus 1, ref angle not = 0  -----
    t = [t0 'ref bus ~= 1, ref ang ~= 0 : '];
    mpc = loadcase(casefile);
    mpc.bus([1;3], BUS_TYPE) = [PV; REF];   %% swap bus types
    bus_soln([1;3], BUS_TYPE) = bus_soln([3;1], BUS_TYPE);   %% swap bus types
    mpc.bus(3, VA) = 3.3014277;
    r = runopf(mpc, mpopt);
    [success, f, bus, gen, branch] = deal(r.success, r.f, r.bus, r.gen, r.branch);
    t_ok(success, [t 'success']);
    t_is(f, f_soln, 3, [t 'f']);
    t_is(   bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %%-----  test OPF with opf.use_vg  -----
    %% get solved AC OPF case from MAT-file
    load soln9_opf_vg;  %% defines bus_soln, gen_soln, branch_soln, f_soln
                        %%    and bus_soln1, gen_soln1, branch_soln1, f_soln1

    %% run with opf.use_vg = 1
    t = [t0 'w/opf.use_vg = 1 : '];
    mpc = loadcase(casefile);
    mpc.gen = mpc.gen([1 2 1 3], :);
    mpc.gencost = mpc.gencost([1 2 1 3], :);
    mpc.gen([1 3], [PMAX PMIN]) = mpc.gen([1 3], [PMAX PMIN]) / 2;
    mpc.gen(3, [QMIN, QMAX]) = 0;   %% no reactive capability for gen 3
    mpc.gencost([1 3], COST:end) = mpc.gencost([1 3], COST:end) / 2;
    mpc.gen(1, VG) = 1.05;
    mpc.gen(3, VG) = 1.06;
    mpopt1 = mpoption(mpopt, 'opf.use_vg', 1);
    r = runopf(mpc, mpopt1);
    t_ok(r.success, [t 'success']);
    t_is(r.f, f_soln, 3, [t 'f']);
    t_is(   r.bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
    t_is(   r.bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   r.bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   r.bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   r.gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
    t_is(   r.gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   r.gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(r.branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
    t_is(r.branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(r.branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

    %% run with opf.use_vg = 0.9
    t = [t0 'w/opf.use_vg = 0.9 : '];
    mpopt1 = mpoption(mpopt, 'opf.use_vg', 0.9);
    r = runopf(mpc, mpopt1);
    t_ok(r.success, [t 'success']);
    t_is(r.f, f_soln1, 3, [t 'f']);
    t_is(   r.bus(:,ib_data   ),    bus_soln1(:,ib_data   ), 10, [t 'bus data']);
    t_is(   r.bus(:,ib_voltage),    bus_soln1(:,ib_voltage),  3, [t 'bus voltage']);
    t_is(   r.bus(:,ib_lam    ),    bus_soln1(:,ib_lam    ),  3, [t 'bus lambda']);
    t_is(   r.bus(:,ib_mu     ),    bus_soln1(:,ib_mu     ),  2, [t 'bus mu']);
    t_is(   r.gen(:,ig_data   ),    gen_soln1(:,ig_data   ), 10, [t 'gen data']);
    t_is(   r.gen(:,ig_disp   ),    gen_soln1(:,ig_disp   ),  3, [t 'gen dispatch']);
    t_is(   r.gen(:,ig_mu     ),    gen_soln1(:,ig_mu     ),  3, [t 'gen mu']);
    t_is(r.branch(:,ibr_data  ), branch_soln1(:,ibr_data  ), 10, [t 'branch data']);
    t_is(r.branch(:,ibr_flow  ), branch_soln1(:,ibr_flow  ),  3, [t 'branch flow']);
    t_is(r.branch(:,ibr_mu    ), branch_soln1(:,ibr_mu    ),  2, [t 'branch mu']);

    t = [t0 'hi-deg polynomial Pg costs : '];
    mpc = loadcase(casefile);
    mpc.gencost = [
        2   1500    0   6   1e-6/5  0   0   0   0   0;
        2   3000    0   5   1e-4/4  0   0   0   0   0;
        2   2000    0   3   1/2     0   0   0   0   0;
    ];
    r = runopf(mpc, mpopt);
    [f, bus, gen, branch] = deal(r.f, r.bus, r.gen, r.branch);
    t_ok(r.success, [t 'success']);
    t_is(f, 11899.4652, 4, [t 'f']);
    t_is(gen(:, PG), [100.703628; 88.719864; 128.679485], 5, [t 'Pg']);
    t_is(gen(:, QG), [8.151114; -27.105856; 12.471647], 5, [t 'Qg']);
    t_is([min(bus(:, VM)) mean(bus(:, VM)) max(bus(:, VM))], ...
        [1.059191 1.079404 1.1], 5, [t 'bus voltage']);

    t = [t0 'hi-deg polynomial Qg costs : '];
    mpc = loadcase(casefile);
    mpc.gencost = [
        2   1500    0   6   1e-6/5  0   0   0   0   0;
        2   3000    0   5   1e-4/4  0   0   0   0   0;
        2   2000    0   3   1/2     0   0   0   0   0;
        2   0       0   6   1e-6/5  0   0   0   0   0;
        2   0       0   3   1/2     0   0   0   0   0;
        2   0       0   5   1e-4/4  0   0   0   0   0;
    ];
    r = runopf(mpc, mpopt);
    [f, bus, gen, branch] = deal(r.f, r.bus, r.gen, r.branch);
    t_ok(r.success, [t 'success']);
    t_is(f, 11998.8345, 4, [t 'f']);
    t_is(gen(:, PG), [100.853092; 89.258007; 128.518423], 5, [t 'Pg']);
    t_is(gen(:, QG), [-20.834352; -7.966093; 34.591995], 5, [t 'Qg']);
    t_is([min(bus(:, VM)) mean(bus(:, VM)) max(bus(:, VM))], ...
        [1.007729 1.050776 1.1], 5, [t 'bus voltage']);

%     %% OPF with user-defined nonlinear constraints
%     t = [t0 'w/nonlin eq constraint : '];
%     mpc = loadcase('case30');
%     mpc.user_constraints.nle = {
%         {'Pg_usr', 1, 'opf_nle_fcn1', 'opf_nle_hess1', {'Pg'}, {}}
%     };
%     r = runopf(mpc, mpopt);
%     t_ok(r.success, [t 'success']);
%     t_is(r.gen(1, PG) * r.gen(2, PG) / 100, r.gen(6, PG), 8, t);
%     t_is(r.gen(6, PG), 20.751163, 5, t);
%     %% check for indexing error, issue #77
%     t_is(r.raw.pimul(61), 0, 6, [t 'issue #77']);
%     t_is(r.mu.nln.u(61), 0, 6, [t 'issue #77']);

    %% OPF with all buses isolated
    t = [t0 'all buses isolated : '];
    mpc.bus(:, BUS_TYPE) = NONE;
    try
        r = runopf(mpc, mpopt);
        t_is(r.success, 0, 12, [t 'success = 0']);
    catch
        t_ok(0, [t 'unexpected fatal error']);
    end

    %% OPF with no branch limits
    t = [t0 'w/no branch limits : '];
    mpc = loadcase(casefile);
    mpc.branch(:, RATE_A) = 0;
    r = runopf(mpc, mpopt);
    t_ok(r.success, [t 'success']);
    t_is(r.f, 5496.128635, 4, [t 'f']);
    t_is(r.gen(:, PG), [90; 220.463932; 10], 5, [t 'Pg']);
    t_is([min(r.bus(:, VM)) mean(r.bus(:, VM)) max(r.bus(:, VM))], ...
        [1.070692 1.090449 1.1], 5, [t 'bus voltage']);
else
    t_skip(num_tests, [t0 'not available']);
end

if have_feature('matlab')
    warning(s2.state, sing_matrix_warn_id);
end

t_end;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function v = cat(dim,varargin)
% implements cat for S2Grid
%

v = cat(dim,vector3d(varargin{1}),varargin{2:end});

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
%SOLVELIN Solve linear system Ax=b.
%
%   [ X, FLAG, T_SOLVE ] = SOLVELIN( A, B, TYPE, X0, VARARGIN ) Solves
%   the linear sparse system Ax = b with solver of TYPE (backslash,
%   mumps, gmres, bicgstab, or amg). X0 is an optional initial guess
%   for the iterative solver types (gmres/bicgstab/amg) and T_SOLVE is
%   the total time. FLAG returns 0 for success and ~0 otherwise.
%
%   See also SOLVESTAT, SOLVETIME

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
report_this_filefun(mfilename('fullpath'));

bv2 = [];
bv3 = [] ;
me = [];
def = {'150'};
ni2 = inputdlg('Number of events in each window?','Input',1,def);
l = ni2{:};
ni = str2double(l);

think

for i = 1:ni/10:length(newt2)-ni
    % [bv magco,  stan] =  bvalcalc(newt2(i:i+ni,:));
    [bv magco stan ] =  bvalca2(newt2(i:i+ni,:));

    bv2 = [bv2 ; magco newt2(i,3)];
end

% Find out of figure already exists
%
[existFlag,figNumber]=figure_exists('Mc with time',1);
newdepWindowFlag=~existFlag;
bdep= figNumber;

% Set up the window

if newdepWindowFlag
    Mcfig = figure_w_normalized_uicontrolunits( ...
        'Name','Mc with time',...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'NextPlot','replace', ...
        'backingstore','on',...
        'Visible','on');

    uicontrol('Units','normal',...
        'Position',[.0 .85 .08 .06],'String','Info ',...
         'Callback','infoz(1)');
    matdraw
end

hold on
figure_w_normalized_uicontrolunits(Mcfig)
hold on
delete(gca)
delete(gca)
axis off

rect = [0.15 0.30 0.7 0.45];
axes('position',rect)
pl = plot(bv2(:,2),bv2(:,1),'^r');
set(pl,'LineWidth',1.5,'MarkerSize',10,...
    'MarkerFaceColor','y','MarkerEdgeColor','r')
hold on
pl = plot(bv2(:,2),bv2(:,1),'b')
set(pl,'LineWidth',1.0)

grid
set(gca,'Color',color_bg)
set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.2)
%
ylabel('Mc')
%set(gca,'Xlim',[t0b teb]);

xlabel('Time')
tist = [  name ' - b(t), ni = ' num2str(ni) ];
title(tist)
done

nu = [];
for i = 1:length(bv2)
    l = newt2.Date >= bv2(i,2) &  newt2.Magnitude >= bv2(i,1);
    nu = [nu length(newt2(l,1))];
end

figure
plot(bv2(:,2),nu,'o')
hold on
plot(bv2(:,2),nu)

figure
plot(bv2(:,1),nu,'o')
hold on
plot(bv2(:,1),nu)

author: CelsoReyes


text:
function [mm2] = km22mm2(km2)
% Convert area from square kilometers to square millimeters.
% Chad A. Greene 2012
mm2 = km2*1000000000000;


text:
% Compile MEX files for Synchrosqueezing toolbox
%
%---------------------------------------------------------------------------------
%    Synchrosqueezing Toolbox
%    Authors: Eugene Brevdo (http://www.math.princeton.edu/~ebrevdo/)
%---------------------------------------------------------------------------------
mex curve_ext.c
mex synsq_cwt_squeeze_mex.c
mex diff_mex.c


text:
function mb = trimeshMeanBreadth(vertices, faces)
%TRIMESHMEANBREADTH Mean breadth of a triangular mesh.
%
%   MB = trimeshMeanBreadth(VERTICES, FACES)
%   Computes the mean breadth (proporitonal to the integral of mean
%   curvature) of a triangular mesh.
%
%   Example
%     [V, F] = createCube;
%     F2 = triangulateFaces(F);
%     MB = trimeshMeanBreadth(V, F2)
%     MB =
%         1.5000
%
%   See also
%   meshes3d, trimeshSurfaceArea, trimeshEdgeFaces, polyhedronMeanBreadth
%
%   References
%   Stoyan D., Kendall W.S., Mecke J. (1995) "Stochastic Geometry and its
%       Applications", John Wiley and Sons, p. 26
%   Ohser, J., Muescklich, F. (2000) "Statistical Analysis of
%       Microstructures in Materials Sciences", John Wiley and Sons, p.352

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2015-08-19, using Matlab 8.5.0.197613 (R2015a)
% Copyright 2015-2022 INRA - Cepia Software Platform

%% Check input validity

if size(faces, 2) ~= 3
    error('meshes3d:trimeshMeanBreadth:NonTriangularMesh', ...
        'Requires a triangular mesh as input');
end

%% Compute edge and edgeFaces arrays
% Uses the same code as in trimeshEdgeFaces

% compute vertex indices of each edge (in increasing index order)
edges = sort([faces(:,[1 2]) ; faces(:,[2 3]) ; faces(:,[3 1])], 2);

% create an array to keep indices of faces "creating" each edge
nFaces = size(faces, 1);
edgeFaceInds = repmat( (1:nFaces)', 3, 1);

% sort edges, keeping indices
[edges, ia, ib] = unique(edges, 'rows'); %#ok<ASGLU>
nEdges = size(edges, 1);

% allocate memory for result
edgeFaces = zeros(nEdges, 2);

% iterate over edges, to identify incident faces
for iEdge = 1:nEdges
    inds = find(ib == iEdge);
    edgeFaces(iEdge, 1:length(inds)) = edgeFaceInds(inds);
end

%% Compute dihedral angle for each edge

% compute normal of each face
normals = meshFaceNormals(vertices, faces);

% allocate memory for resulting angles
alpha = zeros(nEdges, 1);

% iterate over edges
for iEdge = 1:nEdges
    % indices of adjacent faces
    indFace1 = edgeFaces(iEdge, 1);
    indFace2 = edgeFaces(iEdge, 2);

    % normal vector of adjacent faces
    normal1 = normals(indFace1, :);
    normal2 = normals(indFace2, :);

    % compute dihedral angle of two vectors
    alpha(iEdge) = vectorAngle3d(normal1, normal2);
end

%% Compute mean breadth
% integrate the dihedral angles weighted by the length of each edge

% compute length of each edge
lengths = meshEdgeLength(vertices, edges);

% compute product of length by angles
mb = sum(alpha .* lengths) / (4*pi);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
      function h = specgram(s, T, varargin)
         %specgram - plots spectrogram of traces
         %  h = TraceSpectra.specgram(traces) Generates a spectrogram from the
         %  trace(s), overwriting the current figure.  The return value is a
         %  handle to the spectrogram, and is optional. The spectrograms will be
         %  created in the same shape as the passed traces.  ie, if W is a 2x3
         %  matrix of traces, then specgram(TraceSpectra,W) will generate a 2x3
         %  plot of spectra.
         %
         %  Many additional behaviors can be modified through the passing of
         %  additional parameters, as listed further below.  These parameters are
         %  always passed in pairs, as such:
         %
         %  TraceSpectra.specgram(traces,'PARAM1',VALUE1,...,'PARAMn',VALUEn)
         %    Any number of these parameters may be passed to specgram.
         %
         %  specgram(..., 'axis', AXIS_HANDLE)
         %    Specify the axis AXIS_HANDLE within which the spectrogram will be
         %    generated.  The boundary of the axis becomes the boundary for the
         %    entire spectra plot.  For a matrix of traces, this area is
         %    subdivided into NxM subplots, where N and M are the size of the
         %    trace matrix.
         %
         %  specgram(..., 'xunit', XUNIT)
         %    Spedifies the x-unit scale to be used with the spectrogram.  The
         %    default unit is 'seconds'.
         %    valid xunits:
         %     'seconds','minutes','hours','days','doy' (day of year),and 'date'
         %
         %  specgram(..., 'colormap', ALTERNATEMAP)
         %    Instead of using the default colormap, any colormap may be used.  An
         %    alternate way of setting the global map is by using the SETMAP
         %    function.  ALTERNATEMAP will either be a name (eg. grayscale) or an
         %    Nx3 numeric. Type HELP GRAPH3D to see additional useful colormaps.
         %
         %
         %  specgram(..., 'colorbar', COLORBAR_OPTION)
         %    Generates a spectrogram from the trace and uses a specific map
         %    valid COLORBAR_OPTION values: 'horiz' (default),'vert','none',
         %      'HORIZ' places a single colorbar below all plots
         %      'VERT' places a single colorbar to the right of all plots
         %      'NONE' supresses the colorbar placement
         %
         %  specgram(..., 'yscale', YSCALE)
         %    Choosing 'log' Allows the y-axis to be generated on a log-frequency
         %    scale, with uneven vertical cell spacing.  The default value is
         %    'normal', and provides the standard spectrogram view.
         %    valid yscales: 'normal', 'log' (see NOTE below)
         %
         %    NOTE: In order to use the log scale, UIMAGESC needs to be available on
         %    the matlab path.  This routine was created by Frederic Moisy, and may
         %    be downloaded from the maltabcentral fileexchange (File ID: 11368).
         %    If this routine is not found,then the original spectrogram will be
         %    created.
         %
         %   h = specgram(..., 'fontsize', FONTSIZE)
         %   Specify the font size for a spectrogram.  The default font size is 8.
         %
         %  specgram2(..., 'innerLabels', SHOWINNERLABELS)
         %    Suppress the labling of the inside graphs by setting SHOWINNERLABELS
         %    to false.  If this is false, then the frequency label only shows on
         %    the leftmost spectrograms, and the X-unit label only shows on the
         %    bottommost spectrograms.
         %
         %  The following plots a trace using an alternate mapping, an xunit of
         %  'hours', and with the y-axis plotted using a log scale.
         %  ex. specgram(TraceSpectra, trace,...
         %      'colormap', alternateMap,'xunit','h','yscale','log')
         %
         %
         %  Example:
         %    % create an arbitrary subplot, and then plot multiple spectra
         %    a = subplot(3,2,1);
         %    specgram(TraceSpectra,waves,'axis',a); % waves is an NxM trace
         %
         %
         %  See also SPECTRALOBJECT/SPECGRAM2, trace/PLOT, trace/FILLGAPS

         % Thanks to Jason Amundson for providing the way to do log scales

         % currFontSize = 8; % replaced in the parser
         %enforce input arguments.
         if ~isa(T,'TraceData')
            try
               T = SeismicTrace(T);
            catch er
               error('Second input argument should be a TraceData, not a %s',class(T));
            end
         end
         p = parseSpecgramInputs(s, varargin);

         %% search for relevent property pairs passed as parameters
         myaxis = p.Results.axis; % area to be used
         % position not used...
         alternateMap = p.Results.colormap;
         xChoice = p.Results.xunit; %time units for plot
         currFontSize = p.Results.fontsize; % used for all laels within plot
         colorbarpref = p.Results.colorbar; %position of colorbar
         suppressLabels = p.Results.innerlabels; %show only outside
         useXlabel = p.Results.useXlabel; % no x labels
         useYlabel = p.Results.useYlabel; % no y labels

         logscale = strcmpi(p.Results.yscale,'log');

         %% figure out exactly WHERE to plot the spectrogram(s)
         %find out area(axis) in which the spectrograms will be plotted
         clabel= 'Relative Amplitude  (dB)';

         if myaxis == 0,
            clf;
            myaxis = gca;
            get(gca,'position');
         else
            get(myaxis,'position');
         end

         %% If there are multiple traces...
         % subdivide the axis and loop through specgram2 with individual traces.

         if numel(T) > 1
            %create the colorbar if desired
            TraceSpectra.createcolorbar(s, colorbarpref, clabel, currFontSize)
            h = TraceSpectra.subdivide_axes(myaxis,size(T));
            remainingproperties = TraceSpectra.buildParameterList(p.Unmatched);
            for n=1:numel(h)
               keepYlabel =  ~suppressLabels || (n <= size(h,1));
               keepXlabel = ~suppressLabels || (mod(n,size(h,2))==0);
               specgram(s,T(n),...
                  'xunit',xChoice,...
                  'axis',h(n),...
                  'fontsize',currFontSize,...
                  'useXlabel',keepXlabel,...
                  'useYlabel',keepYlabel,...
                  'colorbar','none',...
                  remainingproperties{:});
            end
            return
         end

         axes(myaxis);

         if T.hasnan % only 1 trace guaranteed at this point...
            warning('TraceSpectra:specgram:nanValue',...
               ['This trace has at least one NaN value, which will blank',...
               'the related spectrogram segment. ',...
               'Remove NaN values by either splitting up the ',...
               'trace into non-NaN sections or by using TraceData/fillgaps',...
               ' to replace the NaN values.']);
         end

         [xunit, xfactor] = TraceSpectra.parse_xunit(xChoice);

         switch lower(xunit)
            case 'date'
               % we need the actual times...
               Xvalues = get(T,'timevector');

            case 'day of year'
               startvec = datevec(get(T,'start'));
               Xvalues = get(T,'timevector');
               dec31 = datenum(startvec(1)-1,12,31); % 12/31/xxxx of previous year in Matlab format
               Xvalues = Xvalues - dec31;
               xunit = [xunit, ' (', datestr(startvec,'yyyy'),')'];

            otherwise,
               dl= 1:T.nsamples(); %dl : DataLength
               Xvalues = dl .* T.period ./ xfactor;
         end

         %%  once was function specgram(d, NYQ, nfft, noverlap, freqmax, dBlims)

         nx = length(T.data);
         window = hanning(s.nfft);
         nwind = length(window);
         if nx < nwind    % zero-pad x if it has length less than the window length
            T.data(nwind) = 0;
            nx = nwind;
         end

         ncol = fix( (nx - s.over) / (nwind - s.over) );

         %added "floor" below
         colindex = 1 + floor(0:(ncol-1))*(nwind- s.over);
         rowindex = (1:nwind)';
         if length(T.data)<(nwind+colindex(ncol)-1)
            T.data(nwind+colindex(ncol)-1) = 0;   % zero-pad x
         end

         y = zeros(nwind,ncol);

         % put x into columns of y with the proper offset
         % should be able to do this with fancy indexing!
         A_ = colindex(  ones(nwind, 1) ,: )    ;
         B_ = rowindex(:, ones(1, ncol)    )    ;
         y(:) = T.data(fix(A_ + B_ -1));
         clear A_ B_

         for k = 1:ncol;         %  remove the mean from each column of y
            y(:,k) = y(:,k)-mean(y(:,k));
         end

         % Apply the window to the array of offset signal segments.
         y = window(:,ones(1,ncol)).*y;

         % USE FFT
         % now fft y which does the columns
         y = fft(y,s.nfft);
         if ~any(any(imag(T.data)))    % x purely real
            if rem(s.nfft,2),    % nfft odd
               select = 1:(s.nfft+1)/2;
            else
               select = 1:s.nfft/2+1;
            end
            y = y(select,:);
         else
            select = 1:s.nfft;
         end
         f = (select - 1)' * T.samplerate / s.nfft;

         NF = s.nfft/2 + 1;
         nf1=round(f(1) / T.nyquist * NF);                     %frequency window
         if nf1==0, nf1=1; end
         nf2=NF;

         y = 20*log10(abs(y(nf1:nf2,:)));

         F = f(f <= s.freqmax);

         if F(1)==0,
            F(1)=0.001;
         end

         if logscale
            t = linspace(Xvalues(1), Xvalues(end),ncol); % Replaced by TCB - 01/18/2012
            try
               h = uimagesc(t, log10(F), y(nf1:length(F),:), s.dBlims);
            catch exception
               if strcmp(exception.identifier, 'MATLAB:UndefinedFunction')
                  warning('Spectralobject:specgram:uimageNotInstalled',...
                     ['Cannot plot with log spacing because uimage, uimagesc ',...
                     ' not installed or not visible in matlab path.']);
                  h = imagesc(Xvalues,F,y(nf1:length(F),:),s.dBlims);
                  logscale = false;
               else
                  rethrow(exception)
               end
            end
         else
            h = imagesc(Xvalues, F, y(nf1:length(F),:), s.dBlims);
         end
         set(gca, 'fontsize', currFontSize);
         if strcmpi(xunit, 'date')
            datetick('x', 'keepticks');
         end
         titlename = [T.station '-' T.channel '  from:' T.start];
         title (titlename);
         axis xy;
         colormap(alternateMap);
         shading flat
         axis tight;
         if useYlabel
            if ~logscale
               ylabel ('Frequency (Hz)')
            else
               ylabel ('Log Frequency (log Hz)')
            end

         end
         if useXlabel
            xlabel(['Time - ',xunit]);
         end

         %create the colorbar if desired
         TraceSpectra.createcolorbar(s,colorbarpref, clabel, currFontSize);

         %% added a series of functions that help with argument parsing.
         % These were ported from my trace/plot function.

      end %specgram

author: geoscience-community-codes

repo: GISMO


text:
function [muqOrig SqOrig] = modelPriorKernGrad(dynModel)
% MODELPRIORKERNGRAD
% see vargpTimeDynamicsPriorKernGrad. This function here is just a wrapper.
% VARGPLVM

fhandle = str2func([dynModel.type 'PriorKernGrad']);
[muqOrig SqOrig] = fhandle(dynModel);

% if isfield(model, 'paramGroups')
%  g = g*model.paramGroups;
% end

author: SheffieldML

repo: vargplvm


text:
function varargout = randomAngle3d(varargin)
%RANDOMANGLE3D Return a 3D angle uniformly distributed on unit sphere.
%
%   usage
%   [THETA PHI] = randomAngle3d
%   Generate an angle unformly distributed on the surface of the unit
%   sphere.
%
%   "Mathematical" convention is used: theta is the colatitude (angle with
%   vertical axis, 0 for north pole, +pi for south pole, pi/2 for points at
%   equator) with z=0.
%   phi is the same as matlab cart2sph: angle from Ox axis, counted
%   positively counter-clockwise.
%
%   [THETA PHI] = randomAngle3d(N)
%   generates N random angles (N is a scalar). The result is a N-by-2
%   array.
%
%   Example:
%     % Draw some points on the surface of a sphere
%     figure;
%     drawSphere; hold on;
%     drawPoint3d(pts, '.');
%     axis equal;
%
%   See also:
%   angles3d, sph2cart2, cart2sph2
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% created the 18/02/2005.
% Copyright INRA - Cepia Software platform

%   HISTORY
%   2007-01-04 change angle order, update doc
%   2011-06-27 fix bug in input parsing, add doc

N = 1;
if ~isempty(varargin)
    N = varargin{1};
end

phi = 2*pi*rand(N, 1);
theta = asin(2*rand(N, 1)-1) + pi/2;

if nargout<2
    var = [theta phi];
    varargout{1} = var;
else
    varargout{1} = theta;
    varargout{2} = phi;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
%EVALEXPR Evaluates an expression in specified points.
%
%   [ VEVAL ] = EVALEXPR( S_EXPR, XP, PROB, SOLNUM, EVAL_TYPE, TOL )
%   Evaluates the expression S_EXPR in the points specified in XP (in
%   global x/y/z-coordinates). PROB is a valid finite element problem
%   struct. The EVAL_TYPE flag (default 1) toggles vectorized
%   evaluation for dof variables with P0-P2 fem shape functions, or
%   allows faster evalution via intermediate nodal linear
%   interpolation. TOL is a tolerance used by ismembertol/deduplicate
%   to determined which if any evaluation points in XP are aligned
%   with grid points.
%
%       Input       Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       s_expr      string                 String expression to evaluate
%       xp          [n_sdim,n_xp]          Coordinates of evaluation points
%       prob        struct                 Problem definition struct
%       solnum      scalar  {n_sols}       Solution number/time to evaluate
%       eval_type   scalar  {1}            Evaluation type
%                                             0 - Standard/full evaluation
%                                             1 - Vectorized evaluation for dep. variables
%                                             2 - Linear interpolation via nodes (fast)
%                                             3 - Like 2 and allow invalid (nan) values
%       tol         scalar  {1e-5}         Tolerance for evalation/grid points
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       vEval       [n_xp,1]               Output vector of evaluated values

% Copyright 2013-2022 Precise Simulation, Ltd.


text:
function imgs = modcrop(imgs, modulo)
if size(imgs,3)==1
    sz = size(imgs);
    sz = sz - mod(sz, modulo);
    imgs = imgs(1:sz(1), 1:sz(2));
else
    tmpsz = size(imgs);
    sz = tmpsz(1:2);
    sz = sz - mod(sz, modulo);
    imgs = imgs(1:sz(1), 1:sz(2),:);
end

author: huangzehao

repo: caffe-vdsr


text:
function blobs_size = caffe_get_blobs_size(net, blob_names)
%
% This file is part of the code that implements the following ICCV2015 accepted paper:
% title: "Object detection via a multi-region & semantic segmentation-aware CNN model"
% authors: Spyros Gidaris, Nikos Komodakis
% institution: Universite Paris Est, Ecole des Ponts ParisTech
% Technical report: http://arxiv.org/abs/1505.01749
% code: https://github.com/gidariss/mrcnn-object-detection
%
%
% AUTORIGHTS
% --------------------------------------------------------
% Copyright (c) 2015 Spyros Gidaris
%
% "Object detection via a multi-region & semantic segmentation-aware CNN model"
% Technical report: http://arxiv.org/abs/1505.01749
% Licensed under The MIT License [see LICENSE for details]
% ---------------------------------------------------------

num_blobs  = length(blob_names);
blobs_size = cell(num_blobs,1);
for i = 1:num_blobs
    size_this     = net.blobs(blob_names{i}).shape;
    blobs_size{i} = ones(1,4);
    blobs_size{i}((end-length(size_this)+1):end) = size_this;
end
end

author: gidariss

repo: LocNet

save_path: github-repos/MATLAB/gidariss-LocNet


text:
function gbasis=groebner(polyset,ord,varnames,tol)
% GROEBNER - calculate the reduced Groebner basis of a set of polynomials
%  usage: gbasis = groebner(polyset,ord)
%         gbasis = groebner(polyset,ord,varnames)
%         gbasis = groebner(polyset,ord,varnames,tol)
%
% INPUTS: polyset, ord, varnames, tol
%  polyset is a cell array of multivariate polynomial coefficients
%    the cell elements can be all strings or all arrays.
%    If the polynomials are specified with arrays, each row of the array
%    represents a single term.  The first element is the coefficient and
%    the others are the degrees of each unknown in the term.
%    e.g. [-1,2,1;1,0,1] represents -x1^2*x2+x2
%    If the polynomials are specified as strings, they must use variable
%    names 'x1','x2',... (unless varnames is provided) and be valid inputs
%    for petschel.str2poly.m
%  ord is a string specifying the monomial ordering:
%    'lex': lexicographical, order by highest power of most significant
%      indeterminate, e.g. x1>x2^2, x1*x2>x1.
%    'grlex': graded lex, order by total degree then lex,
%      e.g. x1^2*x2 > x1*x2^2 > x1^2
%    'grevlex': graded reverse lex, order by total degree then by lowest
%      power of least significant indeterminate, e.g. x1*x2^2 > x1*x2*x3
%    (for all cases have x1>x2>...>xn)
%    note: revlex is not a well-ordering because 1>x1>x1^2>...
%  varnames (optional) is a cell array of variable names if polyset is a
%    cell array of strings and the variable names are not 'x1', 'x2', etc.
%    If specified, there must always be at least 2 variables.
%  tol (optional) is the zero tolerance for coefficients, otherwise
%    catastrophic cancellation may occur.  Default value is 0
%
% OUTPUTS: gbasis
%  gbasis is the reduced Groebner basis of the set of polynomials, in the
%    same format (strings or matrices) as polyset{1}.
%
% ALGORITHM:
%  A modified Buchberger's algorithm is used.  In Buchberger's algorithm,
%  at each step calculate Sij = (Lij/ai)*Pi - (Lij/aj)*Pj, for all pairs of
%  polynomials Pi, Pj, where ai and aj are the leading terms of Pi and Pj,
%  and Lij is the least common multiple of ai and aj; then reduce Sij with
%  respect to the polynomials and add it to the set (note that Sij always
%  reduces to zero if ai and aj have no variables in common).  The
%  algorithm terminates when no new polynomials can be added to the set.
%  This function reduces the polynomials at each step in order to find the
%  reduced Groebner basis.
%
% EXAMPLE: simplify equations x^2+2xy^2=0, xy+2y^3=1
%  A1=[1,2,0;2,1,2]; % 1*x^2*y^0 + 2*x^1*y^2
%  A2=[-1,0,0;1,1,1;2,0,3]; % -1*x^0*y^0 + 1*x^1*y^1 + 2*x^0*y^3
%  y=groebner({A1,A2},'lex');
%  petschel.poly2str(y) % returns {'x2^3-0.5', 'x1'}
%
% EXAMPLE: same, but with equations specified as strings:
%  groebner({'x^2+2*x*y^2','x*y+2*y^3-1'},'lex',{'x','y'})
%  % returns {'y^3-0.5','x'}
%
% EXAMPLE: show that a set of equations is inconsistent
%  groebner({'x + y', 'x^2 - 1', 'y^2 - 2*x'}, 'lex', {'x', 'y'})
%  % returns {'1'}
%
%
% KNOWN ISSUES
%
%  0. GENERAL: while a significant effort has been made to ensure that the
%    program works as described, unidentified bugs may remain and the results
%    should be regarded with caution.  In particular, the suite of tests
%    performed was by no means comprehensive.
%
%  1. ACCURACY: the algorithm by default works in floating-point which
%    ultimately cannot perform exact arithmetic, so unexpected silent errors may
%    occur even when the coefficients are specified as integers.  A workaround
%    is to use other data types but this has its own issues (see below).
%
%  1a. EXAMPLE (roundoff - linearly independent polynomials)
%    eps % returns 2.2204e-16 (32-bit PC windows machine epsilon)
%    groebner({'x1+x2','x1+1.000000000000001*x2'},'lex') % ok, returns {'x2','x1'}
%    groebner({'x1+x2','x1+1.0000000000000001*x2'},'lex') % not ok, returns {'x2+x1'}
%
%  1b. EXAMPLE (roundoff - linearly dependent polynomials)
%    2/sqrt(2)-sqrt(2) % expect 0, returns -2.2204e-016
%    petschel.poly2str(groebner({[1,1,0;sqrt(2),0,1],[1/sqrt(2),1,0;1,0,1]},'lex')) % ok, returns {'1.41421*x2+x1'}
%    petschel.poly2str(groebner({[1,1,0;sqrt(2),0,1],[sqrt(2),1,0;2,0,1]},'lex')) % not ok, returns {'x2','x1'}
%    petschel.poly2str(groebner({[1,1,0;sqrt(2),0,1],[sqrt(2),1,0;2,0,1]},'lex',{},1e-14)) % ok, returns {'1.41421*x2+x1'}
%
%  1c. EXAMPLE (courtesy of Christophe Lauwerys - catastrophic cancellation & tolerance bug)
%    groebner({'t^3+x+y','t^2+0.5*x^2-x-z^2','t^2+y-z^2'},'lex',{'t','x','y','z'}) % ans{1} contains Inf
%
%  1d. EXAMPLE (from Mathematica documentation - catastrophic cancellation because of large intermediate coefficients)
%    % expect 3 polynomials returns, one of degree 8 in z; instead can get '1'
%    groebner({'x^2 + y^2 + z^2 - 1', 'x*y - z + 2', 'z^2 - 2*x + 3*y'}, 'lex', {'x', 'y','z'})
%
%  2. NON-DOUBLE INPUTS: use of non-double coefficients is not supported.  While
%    it may be tempting to introduce symbolic coefficients, it is not clear how
%    they should be consistently treated - for example 'A*x' could reduce to '0'
%    or 'x' or both, depending upon the assumptions on the coefficient A.
%
%  3. EFFICIENCY: The reduction algorithms aren't terribly efficient.  One of
%    the major bottlenecks at present is the calculation of the leading term,
%    which currently involves a call to sortrows each time.  It could be more
%    efficient to maintain an ordered list or even a priority queue (e.g.
%    implemented as a heap).  Also lattice reduction algorithms such as LLL
%    might speed things up.
%
%
% SEE ALSO:
%  petschel.poly2str, petschel.str2poly, petschel.polynsolve

% Author: Ben Petschel 19/6/2009
%
% Change history:
%  19/6/2009 - first release (ord='grlex' and ord='grevlex' not
%    implemented yet)
%  22/6/2009 - implemented ord='grlex' and ord='grevlex'
%    reduction algorithm now reduces lower-order terms as well
%  23/6/2009 - removed ord='revlex' because it is not a well-ordering
%  17/7/2009 - fixed bug in handling polynomials with >=3 variables
%  20/3/2010 - changed polynomial representation from multidimensional
%    arrays to 2d arrays (more memory efficient with large number of vars)
%  11/10/2010 - bugfix for handling error tolerance

if (numel(polyset)>0) && ischar(polyset{1})

    charout = true;

    if nargin<3

        polyset = petschel.str2poly(polyset);

    else

        polyset = petschel.str2poly(polyset,varnames);

    end

else

    charout = false;

end

if nargin<4

    tol = 0;

end

% determine the number of degrees
nd=max(cellfun(@(x)size(x,2)-1,polyset)); % error if = 0

if nd == 0

    error('number of variables must be >= 1');

end

gbasis = fullreduce(polyset,ord,tol,nd);

oldgbasis = {};

while ~isequal(oldgbasis,gbasis)

    oldgbasis = gbasis;

    Qset = SPoly(gbasis,ord,tol,nd);

    gbasis = fullreduce([gbasis,Qset],ord,tol,nd);

end

if charout

    if nargin<3

        gbasis = petschel.poly2str(gbasis);

    else

        gbasis = petschel.poly2str(gbasis,varnames);

    end

end

end % main function groebner(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Qset = SPoly(Pset,ord,tol,nd)
% calculates all reduced S-polynomials of Pset (Pset must be reduced)
%   Sij = (Lij/ai)*Pi - (Lij/aj)*Pj,
% for all pairs of polynomials Pi, Pj, where ai and aj are the leading
% terms of Pi and Pj, and Lij is the least common multiple of ai and aj
% then reduces Sij wrt Pset

Qset = {};

for i=1:numel(Pset)-1

    for j=i+1:numel(Pset)

        [c1,d1] = leadterm(Pset{i},ord,tol,nd); %#ok<*ASGLU>

        [c2,d2] = leadterm(Pset{j},ord,tol,nd);

        if any(and(d1>0,d2>0))
            % leading terms have a variable in common, so S-poly is nontrivial
            L = max(d1,d2); % LCM of leading terms

            S = addpolys(multiplyterm(Pset{i},L-d1),multconst(multiplyterm(Pset{j},L-d2),-1));

            S = reduceset(S,Pset,ord,tol,nd);

            %if any(S(:)~=0),
            if abs(leadterm(S,ord,tol,nd))>tol

                Qset{end+1}=S; % add S to the list

            end

        end

    end

end

end % helper function SPoly(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Qset = fullreduce(Pset,ord,tol,nd)
% reduces all polynomials wrt each other
Qset = Pset;

oldQset = {};

while ~isequal(oldQset,Qset)

    oldQset = Qset;

    for i=1:numel(Qset)-1

        for j=i+1:numel(Qset)

            Qset{i}=reduce(Qset{i},Qset{j},ord,tol,nd);

            Qset{j}=reduce(Qset{j},Qset{i},ord,tol,nd);

        end

    end

end

% keep only the nonzero results (to within zero tolerance tol)
Qset = Qset(cellfun(@(x)abs(leadterm(x,ord,tol,nd))>tol,Qset));

if numel(Qset)==1

    % special case, a single equation can slip through without being reduced
    c = leadterm(Qset{1},ord,tol,nd);

    % coefficients occupy first column of the array
    Qset{1} = multconst(Qset{1},1/c); % know c~=0 because have already removed zero eqns

end

end % helper function fullreduce(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q = reduceset(P,Pset,ord,tol,nd)
% reduces P wrt all polynomials in Pset
Q = P;

oldQ = [];

while ~isequal(oldQ,Q)

    oldQ = Q;

    for i=1:numel(Pset)

        Q=reduce(Q,Pset{i},ord,tol,nd);

    end

end

end % helper function reduceset(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q1=reduce(P1,P2,ord,tol,nd)
% reduces a polynomial P1 by subtracting multiples of P2 so that leading
% term coefficient is 1 and leading term of P2 does not divide that of P1

Q2=P2;

[c2,d2] = leadterm(Q2,ord,tol,nd);

if abs(c2)<=tol

    % ignore terms less than tol
    keepgoing = false;

    Q1 = P1;

else

    keepgoing = true;

    Q2 = multconst(Q2,1/c2);

    Q1 = 0; % will successively add terms to Q
    % remove the leading terms, in order to reduce wrt lower terms
    remain = P1;

end

while keepgoing

    [c1,d1] = leadterm(remain,ord,tol,nd);

    if abs(c1)<=tol

        % ignore terms less than tol
        keepgoing = false;

    else
        %reduce remain by Q2, if possible, or remove leading term
        if d1>=d2
            % can reduce Q1 by Q2, so subtract a multiple of Q2 from Q1
            remain = addpolys(remain,multconst(multiplyterm(Q2,d1-d2),-c1));

        else
            % cannot reduce any further wrt leading term, but try lower terms
            lead = multiplyterm(c1,d1);

            remain = addpolys(remain,multconst(lead,-1));

            Q1 = addpolys(Q1,lead); % add irreducible terms to Q1

        end

    end

end

% reduce leading coefficient, if possible
c1 = leadterm(Q1,ord,tol,nd);

if abs(c1)>tol

    Q1 = multconst(Q1,1/c1);

end

end % helper function reduce(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q=multconst(P,c)
% returns c*P where c is a constant

Q = P;

Q(:,1)=Q(:,1)*c;

end % helper function multconst(...)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q=addpolys(P1,P2)
% adds two polynomials

s1 = size(P1,2);

s2 = size(P2,2);

if s1>s2

    Q = [P1;P2,zeros(size(P2,1),s1-s2)];

else

    Q = [P1,zeros(size(P1,1),s2-s1);P2];

end

Q = sortrows(Q,2:size(Q,2));

i=1;

while i<size(Q,1)
    % merge rows representing terms with the same exponents
    if Q(i,1)==0

        Q(i,:)=[];

    elseif all(Q(i,2:end)==Q(i+1,2:end))

        Q(i,1)=Q(i,1)+Q(i+1,1);

        Q(i+1,:)=[];

        if Q(i,1)==0

            % remove terms that add to zero
            Q(i,:)=[];

        end

    else

        i=i+1;

    end

end

end % helper function addpolys

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Q=multiplyterm(P,d)
% multiplies a polynomial by a term (no coefficient)

% padding a matrix at the end is same as multiplying by x
d=d(:).';

sP = size(P,2);

sd = size(d,2)+1; % d doesn't include coefficient

if sP<sd

    Q=[P,zeros(size(P,1),sd-sP)];

else

    Q=P;

    d=[(d(:)).',zeros(1,sP-sd)];

end

Q=Q+repmat([0,d],size(Q,1),1);

end % helper function multiplyterm

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [coeff,deg]=leadterm(P,ord,tol,nd)
% returns the leading term of P
% deg is a row vector of length size(P,2)-1
% coeff is the leading coefficient, deg(i) is the degree of xi
% ignores terms with coefficients <= tol
% nd is the minimum number of degrees (in order to handle cases where deg=0)

d=size(P,2);

if all(abs(P(:,1))<=tol)
    % P=0
    coeff=0;
    deg=zeros(1,d-1);
else
    switch ord
        case 'lex'
            P=sortrows(P,2:d);

        case 'revlex'
            error('revlex is not a well-ordering');

        case 'grlex'
            P = [P,sum(P(:,2:d),2)];
            P = sortrows(P,[d+1,2:d]);
            P = P(:,1:d);

        case 'grevlex'
            P = [P,sum(P(:,2:d),2)];
            P = sortrows(P,[d+1,-(d:-1:2)]);
            P = P(:,1:d);

        otherwise
            error('ordering %s not recognized',ord);

    end % switch ord

    ind = find(abs(P(:,1))>tol,1,'last');

    if isempty(ind)

        coeff = 0;

        deg = zeros(1,d-1);

    else

        coeff = P(ind,1);

        deg = P(ind,2:end);

    end

end % if all(P(:)==0)

if nd>d-1

    deg = [deg,zeros(1,nd-d+1)];

end

end % helper function leadterm(...)

author: jmaih


text:
function a = r8cvv_rinc ( mn, a, m, roff, i, dai )

%*****************************************************************************80
%
%% R8CVV_RINC increments row I in an R8CVV.
%
%  Discussion:
%
%    An R8CVV is a "vector of vectors" of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 November 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer MN, the size of the cell array.
%
%    Input, real A(MN), the cell array.
%
%    Input, integer M, the number of rows in the array.
%
%    Input, integer ROFF(M+1), the row offsets.
%
%    Input, integer I, the row.
%    1 <= I <= M.
%
%    Input, real DAI(NR(I)), the increment for A(I,*).
%
%    Output, real A(MN), the cell array.
%
  k1 = roff(i) + 1;
  k2 = roff(i+1);
  nv = k2 + 1 - k1;
  a(k1:k2) = a(k1:k2) + dai(1:nv);

  return
end