text:
function model = changeCOBRAConstraints(model, constraintID, varargin)
% Modify an existing COBRA constraint by providing new settings for the
% constraint to update existing settings.
% USAGE:
%    model = changeCOBRAConstraints(model, constraintID, varargin)
%
% INPUTS:
%    model:             model structure
%    constraintID:      The ID of the constraint (or the index in the
%                       ctrs field)
%
% OPTIONAL INPUTS:
%    varargin:      all elements (name, d, dsense and c can be modified as indicated below.
%                   * idList:           cell array of ids either from either the rxns, or the evars vectors. Can also be a a double vector of indices (in which case indices in evars have to be set off by numel(rxns).
%                   * c:                the elements of the C matrix. If idList is empty, this has to be a vector of length #rxns + #evars
%                   * dsense:           the constraint sense ('L': <= ,'G': >=, 'E': =), or a vector for multiple constraints
%                   * d:                The right hand side of the C*v <= d constraint (or a vector, for multiple simultaneous addition)
%                   * name:             The new, descriptive name of the constraint.
%
% OUTPUT:
%    model:         constrained model
%
% EXAMPLE:
%    Modify the constraint 'A_and_B_Lower_10' to read A + B + C <= 10
%    model = addCOBRAConstraints(model, 'A_and_B_Lower_10', 'idList', {'A','B','C'}, 'c', [1,1,1]);
%
% NOTE:
%    If c is provided, the whole constraint coefficients will be reset and
%    existing coefficients will be removed!
%
% Author: Thomas Pfau, Oct 2018

parser = inputParser();
parser.addRequired('model',@isstruct);
parser.addRequired('constraintID',@(x) ischar(x) && any(ismember(model.ctrs,x)));
parser.addParameter('d',[],@isnumeric);
parser.addParameter('c',[],@(x) isnumeric(x));
parser.addParameter('dsense','', @ischar );
parser.addParameter('idList',{},@(x) iscell(x) );
parser.addParameter('name',[],@(x) ischar(x) );
parser.parse(model,constraintID,varargin{:});

coefs = columnVector(parser.Results.c)';
d = parser.Results.d;
dsense = parser.Results.dsense;
idList = parser.Results.idList;
name = parser.Results.name;
c = parser.Results.c;

% get some model properties
[~,nRxns] = size(model.S);
if isfield(model,'evars')
    [nCtrs,nVars] = size(model.D);
else
    nVars = 0;
end

if ischar(constraintID)
    constraintID = find(ismember(model.ctrs,constraintID));
end

% check idList to determine new C/D values.
if ~isempty(c)
    if isempty(idList)
        if size(coefs,2) ~= nVars + nRxns
            error('If no idList is provided, the c vector has to contain one element for each reaction and each variable in the model');
        end
        newC = c(1:nRxns);
        newD = c((nRxns+1):end);
    else
        if size(c,2) ~= numel(idList)
            error('If idList is provided, the c vector has to contain one element for each element in idList');
        end
        % get the positions of the provided ids in rxns/evars
        [pos,pres] = getIDPositions(model,idList,'rxns');
        if any(~pres)
            missingreactions = idList(~pres);
            error('The following ids were not found in the model:\n%s\nNo Constraint was added',strjoin(missingreactions,', '));
        end
        % create the respective rows.
        newC = sparse(1,nRxns);
        newD = sparse(1,nVars);
        % get the positions of the variables
        vars = pos > nRxns;
        varPos = pos(vars) - nRxns;
        varCoefs = c(vars);
        rxns = pos <= nRxns;
        rxnPos = pos(rxns);
        rxnCoefs = c(rxns);
        newC(rxnPos) = rxnCoefs;
        newD(varPos) =varCoefs;
    end

    % set the rows if applicable
    model.C(constraintID,:) = newC;
    if isfield(model, 'D')
        model.D(constraintID,:) = newD;
    end
end

% update d
if ~isempty(d)
    model.d(constraintID) = d;
end

% update dsense
if ~isempty(dsense)
    model.dsense(constraintID) = dsense;
end

% update name:
if ~isempty(name)
    if ~isfield(model,'ctrNames')
        model = createEmptyFields(model,'ctrNames');
    end
    model.ctrNames{constraintID} = name;
end

author: opencobra

repo: cobratoolbox


text:
function [c, ceq, cGrad, ceqGrad] = pathConstraint(x)
% [c, ceq, cGrad, ceqGrad] = pathConstraint(x)
%
% This function implements a simple path constraint to keep the knee joint
% of the robot from hyer-extending.
%

q1 = x(1,:);
q2 = x(2,:);
q4 = x(4,:);
q5 = x(5,:);

c = [...
    q1-q2;    %Stance knee joint limit
    q5-q4];   %Swing knee joint limit

ceq = [];

if nargout == 4 %Analytic gradients
    % Gradients with respect to:
    % [t,q1,q2,q3,q4,q5,dq1,dq2,dq3,dq4,dq5,u1,u2,u3,u4,u5] = 1+5+5+5
    nCst = 2;   %stance leg ; swing leg
    nGrad = 16;  %time, angles, rates, torques
    nTime = size(x,2);
    cGrad = zeros(nCst,nGrad,nTime);
    cGrad(1,3,:) = -1;  % cst stance wrt q2
    cGrad(1,2,:) = 1; % cst stance wrt q1
    cGrad(2,5,:) = -1;  % cst swing wrt q4
    cGrad(2,6,:) = 1; % cst swing wrt q5

    ceqGrad = [];
end

end

author: MatthewPeterKelly

repo: OptimTraj


text:
function value = c8_div ( z1, z2 )

%*****************************************************************************80
%
%% C8_DIV divides two C8's.
%
%  Discussion:
%
%    A C8 is a complex value.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 February 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, complex Z1, Z2, the arguments.
%
%    Output, complex VALUE, the function value.
%
  a = c8_real ( z1 );
  b = c8_imag ( z1 );
  c = c8_real ( z2 );
  d = c8_imag ( z2 );

  e = c * c + d * d;

  f = ( a * c + b * d ) / e;
  g = ( b * c - a * d ) / e;

  value = f + g * i;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Script demonstrating usage of the ccmod function.
%
% Author: Brendt Wohlberg <brendt@lanl.gov>  Modified: 2015-04-09
%
% This file is part of the SPORCO library. Details of the copyright
% and user license can be found in the 'Copyright' and 'License' files
% distributed with the library.

% Training images
S0 = zeros(512, 512, 5);
S0(:,:,1) = single(stdimage('lena.grey')) / 255;
S0(:,:,2) = single(stdimage('barbara.grey')) / 255;
S0(:,:,3) = single(stdimage('kiel.grey')) / 255;
S0(:,:,4) = single(rgb2gray(stdimage('mandrill'))) / 255;
tmp = single(stdimage('man.grey')) / 255;
S0(:,:,5) = tmp(101:612, 101:612);

%Reduce images size to speed up demo script
tmp = zeros(256, 256, 5);
for k = 1:size(S0,3),
  tmp(:,:,k) = imresize(S0(:,:,k), 0.5);
end
S0 = tmp;

% Filter input images and compute highpass images
npd = 16;
fltlmbd = 5;
[Sl, Sh] = lowpass(S0, fltlmbd, npd);

% Load dictionary
load([sporco_path '/Data/ConvDict.mat']);
dmap = containers.Map(ConvDict.Label, ConvDict.Dict);
D0 = dmap('12x12x36');

% Set up cbpdn parameters
lambda = 0.1;
opt = [];
opt.Verbose = 1;
opt.MaxMainIter = 200;
opt.AutoRho = 1;
opt.AutoRhoPeriod = 1;
opt.RelaxParam = 1.8;

% Compute sparse representation on current dictionary
[X, optinf] = cbpdn(D0, Sh, lambda, opt);

% Set up ccmod parameters
opt = [];
opt.Verbose = 1;
opt.MaxMainIter = 500;
opt.sigma = size(Sh,3);
opt.AutoSigma = 1;
opt.AutoSigmaPeriod = 1;
opt.RelaxParam = 1.8;
opt.AuxVarObj = 1;

% Update dictionary for training set S
[D1, optinf1] = ccmod(X, Sh, size(D0), opt);

% Plot functional value and residuals
figure;
subplot(1,3,1);
plot(optinf1.itstat(:,2));
xlabel('Iterations');
ylabel('Functional value');
subplot(1,3,2);
semilogy(optinf1.itstat(:,4));
xlabel('Iterations');
ylabel('Primal residual');
subplot(1,3,3);
semilogy(optinf1.itstat(:,5));
xlabel('Iterations');
ylabel('Dual residual');

% Update dictionary with new filter sizes for training set S
dsz = [repmat([12 12]', [1 24]) repmat([8 8]', [1 12])];
[D2, optinf2] = ccmod(X, Sh, dsz, opt);

% Display dictionaries
figure;
subplot(1,3,1);
imdisp(tiledict(D0));
title('D0');
subplot(1,3,2);
imdisp(tiledict(D1));
title('D1');
subplot(1,3,3);
imdisp(tiledict(D2, dsz));
title('D2');

% Plot functional value evolution
figure;
plot(optinf1.itstat(:,2), 'r');
hold on;
plot(optinf2.itstat(:,2), 'b');
hold off;
xlabel('Iterations');
ylabel('Functional value');
legend('D2', 'D3');

author: thfylsty


text:
function [bar] = mmHg2bar(mmHg)
% Convert pressure from millimeters of mercury at 0 degrees C to bar
% Chad Greene 2012
bar = mmHg*0.00133322;

author: Sable

repo: mcbench-benchmarks


text:
function [res]=tt_smv(sttm, vec)
%TT-matrix with sparse factors by TT-vector multiplication
%   [res]=TT_SMV(STTM,VEC) Multiplies the TT-matrix stored in the TT1.0
%   format but with "sparse" cores by a vector VEC
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

d=size(vec,1);
res=cell(d,1);

n=sttm{d+2}(1);
% m=size(vec{1},1);
rm1=sttm{d+1}(1);
rv1=size(vec{1},2);
%mat{1} is (n x m x rm1 , vec{1} is mxrv1)
%            n x rm1 x rv1
res{1}=sttm{1}*vec{1};
res{1}=reshape(permute(reshape(res{1},[n,rm1,rv1]),[1,3,2]),[n,rv1*rm1]);

if (d>1)
    n=sttm{d+2}(d);
    % m=size(mat{d},2);
    rm1=sttm{d+1}(d-1);
    rv1=size(vec{d},2);
    %mat{1} is (n x rm1 x rv1 , vec{1} is mxrv1)
    res{d}=sttm{d}*vec{d};
    res{d}=reshape(permute(reshape(res{d},[n,rm1,rv1]),[1,3,2]),[n,rv1*rm1]);
end;

for i=2:d-1
    n=sttm{d+2}(i);
    m=size(vec{i},1);
    rm1=sttm{d+1}(i-1);
    rm2=sttm{d+1}(i);
    rv1=size(vec{i},2);
    rv2=size(vec{i},3);
    %mat{i} is n x m x rm1 x rm2, vec{i} is m x rv1 x rv2
    %n x (rv1*rv2)*rm1*rm2
    %n x rm2 x rm1 x m, n x rm2 x rm1 x (rv1*rv2)
    %n x rm2 x rmx1 x rv1 x rv2
    %want: n x rv1*rm1 x rv2*rm2
    res{i}=sttm{i}*reshape(vec{i},[m,rv1*rv2]);
    res{i}=reshape(res{i},[n,rm2,rm1,rv1,rv2]);
    res{i}=permute(res{i},[1,4,3,5,2]);
    res{i}=reshape(res{i},[n,rv1*rm1,rv2*rm2]);
    %res{i}=permute(reshape(permute(mat{i},[1,3,4,2])*reshape(vec{i},[m,rv1*rv2]),[n,rv1,rv2,rm1,rm2]),[1,2,4,3,5]);
end

end

author: oseledets

repo: TT-Toolbox


text:
function varargout = mpvar(var,n,m,opt)
% function P = mpvar(cstr,N,M,opt);
%
% DESCRIPTION
%   Create a polynomial matrix or vector variable
%
% INPUTS
%   cstr: Character string to be used in creating the coefficient vector.
%   N,M: row and column dimensions of polynomial matrix.
%   opt: If N==M, then set opt = 's' to generate a symmetric
%        matrix variable.
%
% OUTPUTS
%   P: polynomial matrix
%
% SYNTAX
%   P = mpvar('c',N)
%       Creates an NxN polynomial matrix with entries c_i_j.
%   P = mpvar('c',N,M)
%   P = mpvar('c',[N,M])
%       Creates an NxM polynomial matrix with entries c_i_j.
%   P = mpvar('c',N,N,'s')
%       Creates an NxN symmetric polynomial matrix with entries c_i_j.
%   P = mpvar('c',[N,1])
%   P = mpvar('c',[1,N])
%       Creates an Nx1 or 1xN polynomial vector with entries c_i if
%       N>1.  If N=1 then this creates a pvar named c.
%   mpvar(cstr,N,M)
%      Equivalent to calling eval([cstr '=mpvar(cstr,N,M);']).
%
% EXAMPLE
%   P = mpvar('p',[2,3])
%
% See also pvar

% 11/10/2002 PP    Initial Coding
% 11/10/2002 PJS   Minor modifications for speed and allow
%                  calling sequence with no outputs
% 12/09/2009 PJS   Modified symmetric call for speed
% 11/07/2010 PJS   Added syntax for [N,M] and for vectors

% Argument checking
error(nargchk(2,4,nargin))
error(nargchk(0,1,nargout))
if nargin==2
    m = [];
    opt = [];
elseif nargin==3
    if ischar(m)
        opt = m;
        m = [];
    else
        opt = [];
    end
end
if isempty(m)
    if isscalar(n)
        m = n;
    else
        m = n(2);
        n = n(1);
    end
end

if n~=m && strcmp(opt,'s')
    error('Symmetric option only valid for square matrices');
end

% Coefficient matrix
nt = n*m;
coefficient = speye(nt);

% Degree matrix
degmat = speye(nt);

% Variable names
stridx = int2str( (1:max([n m]))' );
stridx = strjust( stridx ,'left');
stridx = cellstr(stridx);
varname = cell(nt,1);
if n==1 && m==1
    % Scalar
    varname{1} = var;
elseif n>1 && m>1
    % Matrix
    for i1=1:n ;
        tempvar = [var '_' stridx{i1} '_'];
        for i2 =1:m;
            varname{i1+(i2-1)*n} = [tempvar stridx{i2}];
        end
    end
else
    % Vector
    for i1=1:nt
        varname{i1} = [var '_' stridx{i1}];
    end
end

% Create matrix dimension
matdim = [n m];
chkval = 0; % skip validity check
P = polynomial(coefficient,degmat,varname,matdim,chkval);
if strcmp(opt,'s')
    % Single indices into lower/upper parts of a symmetric matrix
    % sorted so that single index for (i,j) is aligned with (j,i)
    M = reshape(1:n^2,[n n]);
    Ml = tril(M,-1);
    Ms = Ml+Ml';
    [junk,idx]=sort(Ms(:));
    idx=reshape(idx(n+1:end),[2,n*(n-1)/2]);
    lidx = idx(1,:);
    uidx = idx(2,:);

    %Create symmetric matrix of pvars
    % [Working with degmat/varname is faster than indexing into P]
    d2 = degmat;
    d2(:,uidx) = d2(:,lidx)+d2(:,uidx);
    d2(:,lidx) = [];
    v2 = varname;
    v2(lidx) = [];
    P = polynomial(coefficient,d2,v2,matdim,chkval);
    P = combine(P);
end

if nargout==0
    assignin('caller', var, P);
else
    varargout{1} = P;
end

author: yu-jiang

repo: radpbook

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


text:
function o = horzcat(varargin)
% Horizontal concatenation of file_array objects
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

%
% $Id: horzcat.m 1143 2008-02-07 19:33:33Z spm $

o    = cat(2,varargin{:});
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% MEMMAPDATA - create a memory-mapped data class
%
% Usage:
%   >> data_class = memmapdata(data);
%
% Inputs:
%   data         - input data or data file
%
% Outputs:
%   data_class    - output dataset class
%
% Author: Arnaud Delorme, SCCN, INC, UCSD, Nov. 2008

% Copyright (C) 2008 Arnaud Delorme, SCCN, INC, UCSD
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function dataout = memmapdata(data, datadims);

    dataout.fileformat = 'normal';
    if ischar(data)
        if length(data) > 3
            if strcmpi('.dat', data(end-3:end))
                dataout.fileformat = 'transposed';
            end
        end

        % check that the file is not empty
        % --------------------------------
        a = dir(data);
        if isempty(a)
            error([ 'Data file ''' data '''not found' ]);
        elseif a(1).bytes == 0
            error([ 'Empty data file ''' data '''' ]);
        end

        if ~strcmpi(dataout.fileformat, 'transposed')
            dataout.data = memmapfile(data, 'writable', false, 'format', { 'single' datadims 'x' });
        else
            dataout.data = memmapfile(data, 'writable', false, 'format', { 'single' [ datadims(2:end) datadims(1) ] 'x' });
        end
        dataout = class(dataout,'memmapdata');
    else
        dataout = data;
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function pageFinish(api)

gui = localCreatePage();

api.setWizardTitle('Import Data');
api.setWizardDescription('Select Method');

api.setLeavePageCallback(@leavePage);
api.setGotoPageCallback(@gotoPage);

set(gui.hVarname     ,'Callback',          @localExportBehavior);
set(gui.hFinishGroup ,'SelectionChangeFcn',@localExportBehavior);

  function nextPage = leavePage

    nextPage = [];

  end

  function gotoPage

    api.Progress.enableNext(false);
    api.Progress.enableFinish(true);

    localUpdateGUI();

  end

  function localExportBehavior(varargin)

    varName = regexprep(get(gui.hVarname,'String'),'\s','');
    api.Export.setWorkspaceName(varName);

    doScript = get(gui.hRadio(1),'Value');
    api.Export.setGenerateScriptFile(doScript);

    localUpdateGUI();

  end

  function localUpdateGUI()

    doScript = api.Export.getGenerateScriptFile();
    varname  = api.Export.getWorkspaceName();

    val_alias = {'off','on'};
    set(gui.hVarname,'Visible',val_alias{~doScript+1})
    set(gui.hVarname,'String',varname);

    set(gui.hRadio(1),'Value',doScript)
    set(gui.hRadio(2),'Value',~doScript)

    data = api.getDataTransformed();

    set(gui.hTitle,'String',[ 'Summary of ' class(data) ' data to be imported:']);
    set(gui.hPreview,'String',cellstr(char(data)));

  end

  function gui = localCreatePage()
    % final import wizard page
    page = api.hPanel;

    h = api.Spacings.PageHeight;
    w = api.Spacings.PageWidth;
    m = api.Spacings.Margin;
    fs = api.Spacings.FontSize;

    summarytitle = uicontrol(...
      'Parent',page,...
      'String','Data summary:',...
      'Position',[1 h-20-m w 20],...
      'Style','text',...
      'FontSize',fs,...
      'HorizontalAlignment','left');

    preview = uicontrol(...
      'Parent',page,...
      'BackgroundColor',[1 1 1],...
      'FontName','monospaced',...
      'FontSize',fs,...
      'HorizontalAlignment','left',...
      'Max',2,...
      'Position',[1 h-210 w-2 175],...
      'String',blanks(0),...
      'Style','edit');

    finish_exp = uibuttongroup('title','Import to',...
      'Parent',page,...
      'FontSize',fs,...
      'units','pixels','position',[1 h-270 w-2 50]);

    radio_exp(2) = uicontrol(...
      'Parent',finish_exp,...
      'HorizontalAlignment','left',...
      'Position',[160 10 150 20 ],...
      'String',' workspace variable',...
      'FontSize',fs,...
      'Value',0,...
      'Style','radio');

    radio_exp(1) = uicontrol(...
      'Parent',finish_exp,...
      'HorizontalAlignment','left',...
      'Position',[10 10 150 20 ],...
      'String',' script (m-file)',...
      'FontSize',fs,...
      'Value',1,...
      'Style','radio');

    workspace(1) = uicontrol(...
      'Parent',page,...
      'BackgroundColor',[1 1 1],...
      'FontName','monospaced',...
      'HorizontalAlignment','left',...
      'FontSize',fs,...
      'Position',[320 15  150 22 ],...
      'String',blanks(0),...
      'Style','edit');

    gui.hFinishGroup = finish_exp;
    gui.hRadio       = radio_exp;
    gui.hVarname     = workspace;
    gui.hTitle       = summarytitle;
    gui.hPreview     = preview;

  end

end

author: mtex-toolbox

repo: mtex

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


text:
%% BASIC COLLISION EVENT (collisionEvent.m) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This class is designed as a child class to a generic "event", used to
% define the characteristics of a collision event.

% Author: James A. Douthwaite

classdef collisionEvent < interactionEvent

%% INITIALISE THE COLLISION EVENT SPECIFIC PARAMETERS
    properties
        % Collision unique properties
    end
%%  CLASS METHODS
    methods
        % CONSTRUCTION METHOD FOR THE COLLISION EVENT
        function obj = collisionEvent(time,objectA,objectB,summaryInfo,enum)
            % INPUTS:
            % time          - Detection time (s)
            % objectA       - Class for the reference object
            % objectIDB     - Class for the second object
            % summaryInfo   - Event description string
            % enum          - Provided alternative enum (nullified?)
            % OUTPUTS:
            % obj           - The collision object

            % INITIALISE THE GENERIC EVENT SUPERCLASS
            obj = obj@interactionEvent(time,'COLLISION',objectA,objectB,summaryInfo);

            % ASSIGN AN EVENT NUMERATION
%             if ~exist('enum','var') || ~isa(enum,'uint8')
            if nargin < 5 || ~isa(enum,'uint8')
                obj.type = eventType.collision;
            else
                obj.type = enum;
            end
        end
    end
end

repo: OpenMAS


text:
function  readframe = klb_get_readframe_fcn(filename,varargin)

[dim,numThreads] = myparse(varargin,'dim',3,'numThreads',[]);
readframe = @(f) klbreadframe(f);

  function im = klbreadframe(f)
    im = readKLBslice(filename,f,dim,numThreads);
  end
end

author: kristinbranson

repo: JAABA


text:
function boxes3d(varargin)
%BOXES3D Description of functions operating on 3D boxes.
%
%   A box defined by its coordinate extents:
%   BOX = [XMIN XMAX YMIN YMAX ZMIN ZMAX].
%
%   Example
%   % Draw a polyhedron together with its bounding box
%   [n e f]= createIcosahedron;
%   drawPolyhedron(n, f);
%   hold on;
%   drawBox3d(point3dBounds(n))
%
%
%   See also
%   boundingBox3d, box3dVolume, drawBox3d
%   intersectBoxes3d, mergeBoxes3d, randomPointInBox3d
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-07-26,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function test_ft_spikedensity()

% MEM 1gb
% WALLTIME 00:10:00
% DEPENDENCY ft_spikedensity

data = [];
nTrials = 100;
nSamples = 1001;
R = 0.51;
% create artificial data.
for iTrial = 1:nTrials
    data.trial{iTrial}(1,:) = round(R*rand(1,nSamples));
    data.trial{iTrial}(2,:) = (R*rand(1,nSamples));
end

data.time(1:nTrials) = {linspace(0,1,nSamples)};
data.fsample = 1000;
data.label{1} = 'spk1';
data.label{2} = 'eeg1';
data.hdr = [];
data.cfg.trl = [];

%%
cfgSdf.timwin        = [-0.02 0.02];
cfgSdf.winfunc        = 'gausswin';
cfgSdf.latency       = [0 3];
cfgSdf.keeptrials = 'yes';
cfgSdf.spikechannel = 1
[sdf sdfdata] = ft_spikedensity(cfgSdf,data);
%% show that we can also enter a spike input
spike = ft_checkdata(data, 'datatype', 'spike', 'feedback', 'yes');
%%
[sdf2 sdfdata] = ft_spikedensity(cfgSdf,spike);

any(sdf2.trial(:)-sdf.trial(:))>0
%%
nTrials = length(data.trial);
for iTrial = 1 : 5
    figure, plot(sdfdata.time{iTrial}, sdfdata.trial{iTrial}), hold on,

    spks = find(data.trial{iTrial}(1,:));
    xx = [data.time{iTrial}(1,spks); data.time{iTrial}(1,spks)];
    yy = [ones(1,length(spks))*max(sdfdata.trial{iTrial});ones(1,length(spks))*min(sdfdata.trial{iTrial})];
    hold on
    plot(xx,yy,'r'),
end
nanmean(sdf.avg) % expect 20 hz
close all
%%
% do the same for rectwin and alphawin
cfgSdf.winfunc = 'rectwin';

tic
[sdf sdfdata] = ft_spikedensity(cfgSdf,data);
toc
%%
nTrials = length(data.trial);
for iTrial = 1 : 5
    figure, plot(sdfdata.time{iTrial}, sdfdata.trial{iTrial}), hold on,

    spks = find(data.trial{iTrial}(1,:));
    xx = [data.time{iTrial}(1,spks); data.time{iTrial}(1,spks)];
    yy = [ones(1,length(spks))*max(sdfdata.trial{iTrial});ones(1,length(spks))*min(sdfdata.trial{iTrial})];
    hold on
    plot(xx,yy,'r'),
end
close all

% note that the firing rate matches to 1000/41ms ~ 25 hz

%% for alphawin
cfgSdf.winfunc = 'alphawin';

[sdf sdfdata] = ft_spikedensity(cfgSdf,data);
nTrials = length(data.trial);
for iTrial = 1 : 5
    figure, plot(sdfdata.time{iTrial}, sdfdata.trial{iTrial}), hold on,

    spks = find(data.trial{iTrial}(1,:));
    xx = [data.time{iTrial}(1,spks); data.time{iTrial}(1,spks)];
    yy = [ones(1,length(spks))*max(sdfdata.trial{iTrial});ones(1,length(spks))*min(sdfdata.trial{iTrial})];
    hold on
    plot(xx,yy,'r'),
end
close all

%% now check if we use different timwinsif we see that back with a gaussian
cfgSdf = [];
cfgSdf.timwin = [0 0.02];
cfgSdf.winfunc = 'gausswin';
cfgSdf.latency = [0 1];

[sdf sdfdata] = ft_spikedensity(cfgSdf,data);
nTrials = length(data.trial);
for iTrial = 1 : 5
    figure, plot(sdfdata.time{iTrial}, sdfdata.trial{iTrial}), hold on,

    spks = find(data.trial{iTrial}(1,:));
    xx = [data.time{iTrial}(1,spks); data.time{iTrial}(1,spks)];
    yy = [ones(1,length(spks))*max(sdfdata.trial{iTrial});ones(1,length(spks))*min(sdfdata.trial{iTrial})];
    hold on
    plot(xx,yy,'r'),
end
close all

%% try if this function also works with different latencies as promised in the help

clear
nTrials = 100;
nSamples = 1001;
R = 0.51;
% create artificial data.
data.fsample = 1000;

for iTrial = 1:nTrials
    % first put the latency as an integer of the sample frequency
    latency = 0 + round(20*(rand-0.5))/1000; % is going to include some trials, and exlude some
    latencyEnd = 1 + round(20*(rand-0.5))/1000;
    timeaxis = latency:(1/data.fsample):latencyEnd;
    data.time{iTrial}  = timeaxis;
    data.trial(iTrial) = {round(R*rand(1,length(timeaxis)))};
    latencies(iTrial,:) = minmax(timeaxis);
end
close all

data.label{1} = 'spk1';
data.hdr = [];
data.cfg.trl = [];
%% set the defaults
cfgSdf = [];
cfgSdf.latency = 'minperiod'
%cfgSdf.checkspikechan = 'no'
cfgSdf.timwin = [-0.25 0.25]
cfgSdf.winfunc = 'rectwin'
[sdf sdfdata] = ft_spikedensity(cfgSdf,data); % warning here!
figure, plot(sdf.time,sdf.avg)

%% allow variable trial length, check the dof
cfgSdf = [];
cfgSdf.keeptrials = 'yes';
cfgSdf.timwin = [-0.05 0.05];
cfgSdf.winfunc = 'rectwin';
cfgSdf.latency = 'maxperiod';
cfgSdf.vartriallen = 'yes';
[sdf sdfdata] = ft_spikedensity(cfgSdf,data);
figure, plot(sdf.time,sdf.avg)
any(cellfun(@length,data.trial)-cellfun(@length,sdfdata.trial))
% note the discrepancy of one sample, why?
%%
nTrials = length(data.trial);
for iTrial = 1 : 5
    figure, plot(sdfdata.time{iTrial}, sdfdata.trial{iTrial}), hold on,

    spks = find(data.trial{iTrial});
    xx = [data.time{iTrial}(spks); data.time{iTrial}(spks)];
    yy = [ones(1,length(spks))*max(sdfdata.trial{iTrial});ones(1,length(spks))*min(sdfdata.trial{iTrial})];
    hold on
    plot(xx,yy,'r'),
end
close all

%sdfdata.trial has maximum length if it fits the window, otherwise it has less length
figure, plot(sdf.time,sdf.avg)
pause(1)
close all

%% do not allow variable trial length
cfgSdf = [];
cfgSdf.timwin = [-0.001 0.001];
cfgSdf.winfunc = 'rectwin';
%cfgSdf.winfuncopt = 0.0;
cfgSdf.latency = [0 1];
cfgSdf.vartriallen = 'yes';
sum(latencies(:,1)>0| latencies(:,2)<1)
[sdf sdfdata] = ft_spikedensity(cfgSdf,data);
figure, plot(sdf.time,sdf.avg)
%% do not allow variable trial length
cfgSdf = [];
cfgSdf.timwin = [-0.02 0.02];
cfgSdf.winfunc = 'rectwin';
cfgSdf.latency = [0 1];
cfgSdf.vartriallen = 'no';
[sdf sdfdata] = ft_spikedensity(cfgSdf,data);
latencies(sdf.cfg.trials,:)
figure, plot(sdf.time,sdf.avg)
pause(1)
close all

%% check all options one for one
cfgSdf = [];
[sdf sdfdata] = ft_spikedensity(cfgSdf,data);
%%
cfgSdf.keeptrials = 'no';
[sdf sdfdata] = ft_spikedensity(cfgSdf,data);
cfgSdf = [];
%%
cfgSdf.latency = 'prestim'
[sdf sdfdata] = ft_spikedensity(cfgSdf,data);
cfgSdf = [];
%%
cfgSdf.winfunc = 'gauss'
[sdf sdfdata] = ft_spikedensity(cfgSdf,data);
cfgSdf = [];


text:
function [ a, det ] = dpodi ( a, lda, n, job )

%*****************************************************************************80
%
%% DPODI computes the determinant and inverse of a certain matrix.
%
%  Discussion:
%
%    The matrix is real symmetric positive definite.
%    DPODI uses the factors computed by DPOCO, DPOFA or DQRDC.
%
%    A division by zero will occur if the input factor contains
%    a zero on the diagonal and the inverse is requested.
%    It will not occur if the subroutines are called correctly
%    and if DPOCO or DPOFA has set INFO == 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 June 2005
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Dongarra, Moler, Bunch and Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real A(LDA,N), the output A from DPOCO or DPOFA, or the output
%    X from DQRDC.
%
%    Input, integer LDA, the leading dimension of the array A.
%
%    Input, integer N, the order of the matrix A.
%
%    Input, integer JOB, specifies the task.
%    11, both determinant and inverse.
%    01, inverse only.
%    10, determinant only.
%
%    Output, real A(LDA,N), if DPOCO or DPOFA was used to factor A then
%    DPODI produces the upper half of inverse(A).  If DQRDC was used to
%    decompose X then DPODI produces the upper half of inverse(X'*X)
%    where X' is the transpose.  Elements of A below the diagonal are
%    unchanged.  If the units digit of JOB is zero, A is unchanged.
%
%    Output, real DET(2), the determinant of A or of X'*X
%    if requested.
%      determinant = DET(1) * 10.0**DET(2)
%    with 1.0 <= DET(1) < 10.0 or DET(1) == 0.0.
%

%
%  Compute the determinant.
%
  if ( job / 10 ~= 0 )

    det(1) = 1.0;
    det(2) = 0.0;
    s = 10.0;

    for i = 1 : n

      det(1) = a(i,i) * a(i,i) * det(1);

      if ( det(1) == 0.0 )
        break
      end

      while ( det(1) < 1.0 )
        det(1) = s * det(1);
        det(2) = det(2) - 1.0;
      end

      while ( s <= det(1) )
        det(1) = det(1) / s;
        det(2) = det(2) + 1.0;
      end

    end

  end
%
%  Compute inverse(R).
%
  if ( mod ( job, 10 ) ~= 0 )

    for k = 1 : n

      a(k,k) = 1.0 / a(k,k);
      t = -a(k,k);
      a(1:k-1,k) = dscal ( k-1, t, a(1:k-1,k), 1 );

      for j = k+1 : n
        t = a(k,j);
        a(k,j) = 0.0;
        a(1:k,j) = daxpy ( k, t, a(1:k,k), 1, a(1:k,j), 1 );
      end

    end
%
%  Form inverse(R) * (inverse(R))'.
%
    for j = 1 : n
      for k = 1 : j-1
        t = a(k,j);
        a(1:k,k) = daxpy ( k, t, a(1:k,j), 1, a(1:k,k), 1 );
      end
      t = a(j,j);
      a(1:j,j) = dscal ( j, t, a(1:j,j), 1 );
    end

  end

  return
end


text:
%% An example of MATLAB GUI to control KUKA iiwa manipualtor using KST toolbox

%% To run the program:
% 1- Start the MatlabToolboxServer application from the teach pendant of
% the robot
% 2- Run the following script

% Copyright: Mohammad SAFEEA

% Add KST path to MATLAB
cDir = pwd;
cDir=getTheKSTDirectory(cDir);
addpath(cDir);

% Run the interface
a0_interface

repo: KST-Kuka-Sunrise-Toolbox


text:
function h = prtPlotUtilScatter(cX)
% Internal function PRT function

nPlotDimensions = size(cX,2);
if nPlotDimensions < 1
    warning('prt:prtPlotUtilScatter:NoPlotDimensionality','No plot dimensions requested.');
    return
end
if nPlotDimensions > 3
    error('prt:prtPlotUtilScatter:plotDimensionality','The number of requested plot dimensions (%d) is greater than 3.',nPlotDimensions);
end

switch nPlotDimensions
    case 1
        h = scatter(cX, ones(size(cX)));
    case 2
        h = scatter(cX(:,1), cX(:,2));
    case 3
        h = scatter3(cX(:,1), cX(:,2), cX(:,3));
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function isUpdated = bst_update(AskConfirm)
% BST_UPDATE:  Download and install the latest version of Brainstorm.
%
% USAGE:  isUpdated = bst_update(AskConfirm)
%         isUpdated = bst_update()
%
% INPUT:
%    - AskConfirm: {0,1}, If 1, ask user confirmation before proceeding to update

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2009-2022

% Java imports
import org.brainstorm.icon.*;
% Parse inputs
if (nargin == 0) || isempty(AskConfirm)
    AskConfirm = 0;
end
isUpdated = 0;

% === ASKING CONFIRMATION ===
if AskConfirm
    res = java_dialog('confirm', ['Download latest Brainstorm update ?' 10 10 ...
                                  'To turn off automatic updates, edit software preferences.' 10 10], 'Update');
    if ~res
        return
    end
end

% === DOWNLOAD NEW VERSION ===
% Get update zip file
urlUpdate  = 'http://neuroimage.usc.edu/bst/getupdate.php?c=UbsM09&src=1';
[installDir, bstDir] = fileparts(fileparts(fileparts(fileparts(mfilename('fullpath')))));
zipFile    = fullfile(installDir, 'brainstorm_update.zip');

% Check permissions
if ~file_attrib(installDir, 'w') || ~file_attrib(fullfile(installDir, bstDir), 'w')
    strMsg = 'Error: Installation folder is read-only...';
    if ispc
        strMsg = [strMsg 10 10 ...
                  'On some Windows computers, the user folders Documents and Downloads' 10 ...
                  'and the system folder C:\Programs\ are seen as read-only by Matlab.' 10 ...
                  'Try with admin privileges (right-click on Matlab > Run as Administrator)' 10 ...
                  'and if you still cannot update, move the brainstorm3 folder somewhere else.'];
    end
    disp(['BST> ' strrep(strMsg, char(10), [char(10) 'BST> '])]);
    if AskConfirm
        java_dialog('msgbox', [strMsg 10 10], 'Update');
    end
    return
end

% Download file
errMsg = gui_brainstorm('DownloadFile', urlUpdate, zipFile, 'Brainstorm update');
% If file was not downloaded correctly
if ~isempty(errMsg)
    disp('BST> Update: Unable to download updates.');
    disp(['BST>     ' strrep(errMsg, char(10), [10 'BST>     '])]);
    if AskConfirm
        java_dialog('msgbox', ['Could not download new packages:' 10 errMsg 10 10 ...
                               'Software was not updated.' 10 10], 'Update');
    end
    return
end

% === STOP BRAINSTORM ===
if isappdata(0, 'BrainstormRunning')
    bst_exit();
end

% === DELETE THE PREVIOUS INSTALLATION ===
% New progress bar
jProgressBar = javax.swing.JProgressBar();
jProgressBar.setIndeterminate(1);
jProgressBar.setStringPainted(0);
jLabel = javax.swing.JLabel('Installing new version...');
jPanel = gui_river();
jPanel.add('p hfill vfill', jLabel);
jPanel.add('p hfill', jProgressBar);
jDialog = javax.swing.JFrame('Brainstorm update');
jDialog.setAlwaysOnTop(1);
jDialog.setResizable(0);
jDialog.setPreferredSize(java.awt.Dimension(350, 130));
jDialog.setDefaultCloseOperation(javax.swing.JFrame.DISPOSE_ON_CLOSE);
jDialog.getContentPane().add(jPanel);
jDialog.pack();
jDialog.setLocationRelativeTo([]);
jDialog.setVisible(1);
jDialog.getContentPane().repaint();
jDialog.setIconImage(IconLoader.ICON_APP.getImage());
disp('BST> Update: Removing previous installation...');

% Go to zip folder (to make sure we are not in a folder we are deleting)
cd(installDir);
% Try the folders separately
warning('off', 'MATLAB:RMDIR:RemovedFromPath');
try
    rmdir(fullfile(installDir, bstDir, 'toolbox'), 's');
end
try
    rmdir(fullfile(installDir, bstDir, 'external'), 's');
end
try
    rmdir(fullfile(installDir, bstDir, 'bin'), 's');
end
try
    rmdir(fullfile(installDir, bstDir, 'defaults', 'anatomy', 'ICBM152'), 's');
end
try
    rmdir(fullfile(installDir, bstDir, 'defaults', 'eeg', 'ICBM152'), 's');
end
try
    rmdir(fullfile(installDir, bstDir, 'defaults', 'anatomy', 'Colin27'), 's');
end
try
    rmdir(fullfile(installDir, bstDir, 'defaults', 'eeg', 'Colin27'), 's');
end
try
    rmdir(fullfile(installDir, bstDir, 'defaults', 'anatomy', 'MNI_Colin27'), 's');
end
try
    rmdir(fullfile(installDir, bstDir, 'defaults', 'eeg', 'MNI_Colin27'), 's');
end
try
    rmdir(fullfile(installDir, bstDir, 'defaults', 'eeg', 'NotAligned'), 's');
end
warning('on', 'MATLAB:RMDIR:RemovedFromPath');

% === UNZIP FILE ===
disp('BST> Update: Unzipping...');
% Unzip update file
unzip(zipFile);
% Delete temporary update file
delete(zipFile);
% If the current install folder is not named 'brainstorm3': Copy files to destination folder
if ~strcmp(bstDir, 'brainstorm3')
    % Replace all old files with new files
    copyfile(fullfile(installDir, 'brainstorm3'), fullfile(installDir, bstDir), 'f');
    % Delete downloaded files
    rmdir(fullfile(installDir, 'brainstorm3'), 's');
end
% Add some folders to the path again
addpath(fullfile(installDir, bstDir, 'toolbox', 'misc'));
addpath(fullfile(installDir, bstDir, 'toolbox', 'core'));
addpath(fullfile(installDir, bstDir, 'toolbox', 'io'));
addpath(fullfile(installDir, bstDir, 'toolbox', 'gui'));

% === DISPLAY RELEASE NOTES ===
% Close waitbar
jDialog.setVisible(0);
jDialog.dispose();
% Display the latest updates
bst_mutex('create', 'ReleaseNotes');
jFrame = view_text({fullfile(installDir, bstDir, 'doc', 'updates.txt'), ...
                    fullfile(installDir, bstDir, 'doc', 'updates_2020.txt')}, 'Release notes', 1);
java_setcb(jFrame, 'WindowClosingCallback', @(h,ev)bst_mutex('release', 'ReleaseNotes'));
bst_mutex('waitfor', 'ReleaseNotes');

% === RESET ENVIRONMENT ===
% Delete the brainstorm.jar to force downloading a new one when Brainstorm restarts
jarFile = fullfile(installDir, bstDir, 'java', 'brainstorm.jar');
if exist(jarFile, 'file')
    delete(jarFile);
    % The brainstorm.jar file could not be deleted because it was still in use: delete it at next startup
    if exist(jarFile, 'file')
        % Create file to indicate that brainstorm.jar should be deleted
        UpdateFile = fullfile(installDir, bstDir, 'java', 'outdated_jar.txt');
        fid = fopen(UpdateFile, 'w');
        fwrite(fid, ['delete(''' jarFile ''');']);
        fclose(fid);
    end
end
% Update is finished
disp('BST> Update: Done.');

% === RESTART MATLAB ===
h = msgbox(['Brainstorm was updated successfully.' 10 10 ...
            'Matlab will now close.' 10 ...
            'Restart Matlab and run brainstorm.m to finish the installation.' 10 10], 'Update');
waitfor(h);
exit;

author: brainstorm-tools


text:
function imgs = readspefile(inputfile)
% imgs = readspefile(inputfile)
% Reads SPE file. This is copied from CSSTORM.m

fid = fopen(inputfile, 'r', 'l');
header = fread(fid,2050,'*uint16');

% parsing .spe header
x_dim = double(header(22));
y_dim = double(header(329));
z_dim = double(header(724));
img_mode = header(55);
if img_mode == 0
    imgs = fread(fid, inf, 'float32');
elseif img_mode == 1
    imgs = fread(fid, inf, 'uint32');
elseif img_mode == 2
    imgs = fread(fid, inf, 'int16');
elseif img_mode == 3
    imgs = fread(fid, inf, 'uint16');
end

fclose(fid);

imgs = reshape(imgs, [x_dim y_dim z_dim]);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% Exercise 4.31: Design of a cantilever beam (GP)
% Boyd & Vandenberghe "Convex Optimization"
% Almir Mutapcic - 01/30/06
% Updated to use GP mode 02/08/06
% (a figure is generated)
%
% We have a segmented cantilever beam with N segments. Each segment
% has a unit length and variable width and height (rectangular profile).
% The goal is minimize the total volume of the beam, over all segment
% widths w_i and heights h_i, subject to constraints on aspect ratios,
% maximum allowable stress in the material, vertical deflection y, etc.
%
% The problem can be posed as a geometric program (posynomial form)
%     minimize    sum( w_i* h_i)
%         s.t.    w_min <= w_i <= w_max,       for all i = 1,...,N
%                 h_min <= h_i <= h_max
%                 S_min <= h_i/w_i <= S_max
%                 6*i*F/(w_i*h_i^2) <= sigma_max
%                 6*F/(E*w_i*h_i^3) == d_i
%                 (2*i - 1)*d_i + v_(i+1) <= v_i
%                 (i - 1/3)*d_i + v_(i+1) + y_(i+1) <= y_i
%                 y_1 <= y_max
%
% with variables w_i, h_i, d_i, (i = 1,...,N) and v_i, y_i (i = 1,...,N+1).
% (Consult the book for other definitions and a recursive formulation of
% this problem.)

% optimization variables
N = 8;

% constants
wmin = .1; wmax = 100;
hmin = .1; hmax = 6;
Smin = 1/5; Smax = 5;
sigma_max = 1;
ymax = 10;
E = 1; F = 1;

cvx_begin gp
  % optimization variables
  variables w(N) h(N) v(N+1) y(N+1);

  % objective is the total volume of the beam
  % obj = sum of (widths*heights*lengths) over each section
  % (recall that the length of each segment is set to be 1)
  minimize( w'*h )
  subject to
    % non-recursive formulation
    d = 6*F*ones(N,1)./(E*ones(N,1).*w.*h.^3);
    for i = 1:N
      (2*i-1)*d(i) + v(i+1) <= v(i); %#ok
      (i-1/3)*d(i) + v(i+1) + y(i+1) <= y(i); %#ok
    end

    % constraint set
    wmin <= w    <= wmax; %#ok
    hmin <= h    <= hmax; %#ok
    Smin <= h./w <= Smax; %#ok
    6*F*[1:N]'./(w.*(h.^2)) <= sigma_max; %#ok
    y(1) <= ymax; %#ok
cvx_end

% display results
disp('The optimal widths and heights are: ');
w, h
fprintf('The optimal minimum volume of the beam is %3.4f.\n', sum(w.*h));

% plot the 3D model of the optimal cantilever beam
figure, clf
cantilever_beam_plot([h; w])


text:
function [reconstruction]=gsp_pyramid_synthesis_old(coarsest_approximation,prediction_errors,Gs,varargin)
%GSP_PYRAMID_SYNTHESIS Synthesizes a signal from its graph pyramid transform coefficients
%   Usage:  reconstruction=gsp_pyramid_synthesis(coarsest_approximation,prediction_errors,Gs);
%           reconstruction=gsp_pyramid_synthesis(coarsest_approximation,prediction_errors,Gs,param);
%
%   Input parameters:
%         coarsest_approximation    : The coarsest approximation of the original signal.
%         prediction_errors         : Cell array with the prediction errors at each level.
%         Gs                        : A multiresolution sequence of graph structures, including the idx parameters tracking the subsampling pattern.
%   Output parameters:
%         reconstruction            : The synthesized signal.
%   Additional parameters:
%         param.use_exact           : To use exact graph spectral filtering instead of the Chebyshev approximation.
%         param.order               : Degree of the Chebyshev approximation (default=30).
%         param.least_squares       : Set to 1 to use the least squares synthesis (default=0)
%         param.h_filters           : The filters used in the analysis operator. These are required for least squares synthesis, but not for the direct synthesis method
%
%   'gsp_pyramid_synthesis(coarsest_approximation,prediction_errors,Gs)'
%   synthesizes a signal from its graph pyramid transform coefficients.
%
%   See also:
%
%   Demos:
%
%   References:

%   AUTHOR : David I Shuman.
%   TESTING:
%   REFERENCE:

% Read input parameters
if nargin>3
    param=varargin{1};
else
    param=0;
end

num_levels=length(prediction_errors);

% Compute the pyramid transform
current_coarse_approximation=coarsest_approximation;

for i=num_levels:-1:1
    if isfield(param, 'least_squares')
        if param.least_squares
            if ~isfield(param, 'h_filters')
                error('h-filter not provided');
            else
                param.h_filter=param.h_filters{i};
            end
        end
    end
    current_coarse_approximation=gsp_pyramid_synthesis_single_interpolation_old(current_coarse_approximation,prediction_errors{i},Gs{i},Gs{i+1}.idx,param);
end
reconstruction=current_coarse_approximation;

repo: gspbox


text:
% Jason Joseph Rebello
% Carnegie Mellon University (Jan 2012 - May 2013)
% MS in Electrical & Computer Engineering
% Principal Component Analysis

% This program uses Principal Component Analysis to reduce the number
% of features used in face recognition. This program allows you to set K
% if you know the number of Principal components needed or calculates K
% based on how much variance you would like to preserve in the images.
% Note: K has been commented out since I am calculating K based on
% variance.

clear all;
clc;
close all;

fprintf('Principal Component Analysis used for Face Recognition\n\n');
t = cputime;

%% Initializing Variable

fprintf('Initializing Variables');
%K = 100; % Number of reduced features (principal components)
dispFaces = 12; % How many faces to display
variance = 90; % Enter a number between 1 and 99 both inclusive based on
                % how much variance to preserve.
topeig = 12; % Used to display the top 10 eigen vectors found
fprintf('...done\n\n');

%% Loading and Visualizing data

fprintf('Loading and Visualizing data');
load('faces.mat');
displayData(X(1:dispFaces, :)); % function taken from Machine Learning course
title('Original faces')         % by Prof Andrew Ng.
fprintf('...done\n\n');
pause(1);

%% Normalizing features

fprintf('Normalizing features');
[X mu stddev]  = normalizeFeatures(X);
fprintf('...done\n\n');

%% Perform PCA to get eigen vectors and eigen values

fprintf('Performing PCA & displaying top %d eigen vectors found', topeig);
[U,S] = performPCA(X);
displayData(U(:, 1:topeig)');
title('Top Eigen Vecotrs found');
fprintf('...done\n\n');
pause(2);

%% Calculate K based on variance

fprintf('Find best value of K based on Variance');
K = findK(S, variance);
fprintf('...done\n\n');

%% Get data with reduced features

fprintf('Displaying data with reduced features');
reducedData = reducedFeatures(X, U, K);
fprintf('...done\n\n');

%% Recover Original Data

fprintf('Recovering Original Data from reduced features');
XRecovered = recoverData(reducedData, U, K);
fprintf('...done\n\n');

%% Display Original and Recovered Data

fprintf('Displaying Original and Recovered Data');
show(X, XRecovered, dispFaces);
fprintf('...done\n\n');

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


text:
function v = r8vec_house_column ( n, a_vec, k )

%*****************************************************************************80
%
%% R8VEC_HOUSE_COLUMN defines a Householder premultiplier that "packs" a column.
%
%  Discussion:
%
%    The routine returns a vector V that defines a Householder
%    premultiplier matrix H(V) that zeros out the subdiagonal entries of
%    column K of the matrix A.
%
%       H(V) = I - 2 * v * v'
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix A.
%
%    Input, real A_VEC(N,1), column K of the matrix A.
%
%    Input, integer K, the column of the matrix to be modified.
%
%    Output, real V(N,1), a vector of unit L2 norm which defines an
%    orthogonal Householder premultiplier matrix H with the property
%    that the K-th column of H*A is zero below the diagonal.
%

%
%  Destroy all row vectors!
%
  a_vec = a_vec(:);

  v(1:n,1) = 0.0;

  if ( k < 1 || n <= k )
    return
  end

  s = sqrt ( sum ( a_vec(k:n,1).^2 ) );

  if ( s == 0.0 )
    return
  end

  if ( a_vec(k,1) < 0.0 )
    v(k,1) = a_vec(k,1) - abs ( s ) ;
  else
    v(k,1) = a_vec(k,1) + abs ( s );
  end
  v(k+1:n,1) = a_vec(k+1:n,1);
%
%  Normalize V.
%
  s = sqrt ( sum ( v(k:n,1).^2 ) );

  v(k:n,1) = v(k:n,1) / s;

  return
end


text:
% A patch group is a matrix pg of size
%
%     px x px x ch x pt x np
%
% where px is the spatial  size of the patch
%       pt is the temporal size of the patch
%       ch is the number of channels
%       np is the number of patches in the group
%
% This function transforms the patch group in to
% an RGB image with all the patches in the group,
% of size
%
%     pt*(px+sw)-sw x np*(px+sw)-sw x ch
%
% The temporal slices of the patch are stacked vertically,
% with separators of sw pixels. Each patch is then stacked
% horizontally, also with separators of sw pixels.
%
% This function was written to use from patch_group.
%
function pp = build_patch_image(pp, force_one_row, sc, sw)

px = size(pp,1);
pt = size(pp,4);
np = size(pp,5);

if nargin < 2,
	force_one_row = false;
end

if nargin < 3,
	sc = 255;
end

if nargin < 4,
	sw = 1;
end

if (size(sc,2) ~= 1) && (size(sc,2) ~= 3),
	error('The separator colormap should have either 1 or 3 channels')
end

ch = max(size(pp,3),size(sc,2));

if (max(size(sc)) == 1),

	pp = permute(pp, [2 1 4 5 3]);
	pp = cat(2, pp, 0*sc*ones(px, sw, pt, np, ch));
	pp = reshape(pp, [px, pt*(px + sw), np, ch]);

	pp = permute(pp, [2 1 3 4]);
	pp = cat(2, pp, sc*ones(pt*(px + sw), sw, np, ch));
	pp = reshape(pp, [pt*(px + sw), np*(px + sw), ch]);

	pp = pp(1:pt*(px + sw) - sw, 1:np*(px + sw) - sw, :);

	% if resulting image is too long and narrow, split it in more rows
	if (force_one_row == false) && (10*pt < np)

		% factorize np into m*n minimizing the perimeter m+n
		m = floor(sqrt(np));
		while mod(np,m) ~= 0,
			m = m - 1;
		end
		n = np/m;

		if (10*m*pt < n)
			m = floor(sqrt(np));
			n = ceil(np/m);
		end

		% separator between rows of patches has to be wider than separator between
		% the frames of a patch
		swr = sw * (1 + (pt > 1));
		ppp = sc*ones(m*(pt*(px + sw) - sw + swr) - swr, n*(px + sw) - sw, ch);
		for i = 0:m-1,

			start_pp  = i*n*(px + sw) + 1;
			until_pp  = min(i*n*(px + sw) + n*(px + sw) - sw, size(pp,2));
			range_ppp = i*(pt*(px + sw) - sw + swr) + [1:pt*(px + sw) - sw];
			ppp(range_ppp,1:1+until_pp - start_pp,:) = pp(:,start_pp:until_pp,:);

		end

		pp = ppp;

	end

else

	m = 1;

	% if resulting image is too long and narrow, split it in more rows
	if (force_one_row == false) && (10*pt < np)

		% factorize np into m*n minimizing the perimeter m+n
		m = floor(sqrt(np));
		while mod(np,m) ~= 0,
			m = m - 1;
		end

	end

	n = np/m;

	if size(sc,2) < ch, sc = repmat (sc, [1  1 ch]);
	else                sc = reshape(sc, [np 1 ch]);
	end

	if size(pp,3) < ch, pp = repmat(pp,[1 1 ch]);
	end

	ppp = zeros(m*(pt*(px + sw) + sw),n*(px + 2*sw), ch);
	for i = 1:np,

		scp = sc(i,:,:);

		p = pp(:,:,:,:,i);
		p = permute(p, [2 1 4 5 3]);
		sep = zeros([px, sw, pt, 1, ch]);
		p = cat(2, p, repmat(reshape(scp,[1, 1, 1, 1, ch]), px, sw, pt));
		p = squeeze(reshape(p, [px, pt*(px + sw), 1, ch]));

		% usando cat, agregar bordes derecho, izq, y superior
		p = permute(p, [2 1 3]);
		p = cat(1, repmat(scp, sw, px), p);
		p = cat(2, repmat(scp, pt*(px+sw) + sw, sw), p);
		p = cat(2, p, repmat(scp, pt*(px+sw) + sw, sw));

		[ni,mi] = ind2sub([n,m],i);
		ppp((mi-1)*size(p,1) + [1:size(p,1)],...
		    (ni-1)*size(p,2) + [1:size(p,2)],:) = p;

	end

	pp = ppp;

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% Max Jaderberg 2011

function result = run_dot_m( file_to_run, arguments )
%RUN_DOT_M Runs the given .m file with the argument struct given
%   For exmaple run_dot_m('/path/to/function.m', args);
%   args is a struct containing the arguments. function.m must take only
%   one parameter, the argument structure

    [dir, func_name, ext] = fileparts(file_to_run);

    if ~size(ext)
        result = 'Error: Need to give path to .m file';
        return
    end

    if ~strcmp(ext, '.m')
        result = 'Error: Need to give path to .m file';
        return
    end

%     Add function path to current path
    addpath(dir);
    if isstruct(arguments)
        result = feval(func_name, arguments);
    else
        result = feval(func_name);
    end
end

author: jaderberg

repo: python-matlab-bridge


text:
classdef WebSocketServer < handle
    %WEBSOCKETSERVER WebSocketServer is an ABSTRACT class that allows
    %MATLAB to start a java-websocket server instance. It then becomes
    %possible to send messages to any client that connect to it.
    %
    %   In order to make a valid implementation of the class, some methods
    %   must be defined in the superclass:
    %    onOpen(obj,conn,message)
    %    onTextMessage(obj,conn,message)
    %    onBinaryMessage(obj,conn,bytearray)
    %    onError(obj,conn,message)
    %    onClose((obj,conn,message)
    %   The "callback" behaviour of the server can be defined there. If
    %   the server needs to perform actions that are not responses to a
    %   client-caused event, these actions must be performed outside of
    %   these callback methods.

    properties (SetAccess = private)
        Port % Server port
        Secure % True if the server is a secure websocket server (wss)
        Status % Server status
        Connections % Stores active connections' hash and address as well as their java websocket object
        ServerObj % Java-WebSocket server object
    end

    properties (Access = private)
        KeyStore % Location of the keystore
        StorePassword % Keystore password
        KeyPassword % Key password
    end

    methods
        function obj = WebSocketServer(port,keyStore,storePassword,keyPassword)
            % Constructor
            obj.Port = port;
            if nargin>1
                obj.Secure = true;
                obj.KeyStore = keyStore;
                obj.StorePassword = storePassword;
                obj.KeyPassword = keyPassword;
            end
            % Start server
            obj.start();
        end

        function conns = get.Connections(obj)
            % Get current connections as a struct, listing HashCode,
            % Address and Port, use the struct2table method on the returned
            % struct for a better display
            N = size(obj.ServerObj.connections);
            connsObj = obj.ServerObj.connections.toArray;
            codes = num2cell(int32(zeros(N,1)));
            adds = repmat(cellstr(''),N,1);
            ports = codes;
            for n = 1:N
                codes{n} = int32(connsObj(n).hashCode());
                adds{n} = char(connsObj(n).getRemoteSocketAddress.getHostName());
                ports{n} = int32(connsObj(n).getRemoteSocketAddress.getPort());
            end
            conns = struct('HashCode',codes,'Address',adds,'Port',ports);
        end

        function start(obj)
            % Start the WebSocket server
            if obj.Status; error('The server is already running'); end
            % Create the java server object in with specified port
            if obj.Secure
                import io.github.jebej.matlabwebsocket.*;
%                 import io.github.jebej.matlabwebsocket.MatlabWebSocketSSLServer;
                obj.ServerObj = handle(MatlabWebSocketSSLServer(obj.Port,obj.KeyStore,obj.StorePassword,obj.KeyPassword),'CallbackProperties');
            else
%                 import io.github.jebej.matlabwebsocket.MatlabWebSocketServer;
                import io.github.jebej.matlabwebsocket.*;
                obj.ServerObj = handle(MatlabWebSocketServer(obj.Port),'CallbackProperties');
            end
            % Set callbacks
            set(obj.ServerObj,'OpenCallback',@(~,e)obj.openCallback(e));
            set(obj.ServerObj,'TextMessageCallback',@(~,e)obj.textMessageCallback(e));
            set(obj.ServerObj,'BinaryMessageCallback',@(~,e)obj.binaryMessageCallback(e));
            set(obj.ServerObj,'ErrorCallback',@(~,e)obj.errorCallback(e));
            set(obj.ServerObj,'CloseCallback',@(~,e)obj.closeCallback(e));
            % Start the server
            obj.ServerObj.start();
            obj.Status = true;
        end

        function stop(obj,timeout)
            % Stop the server with a timeout to close connections
            if ~obj.Status; error('The server is not running!'); end
            if nargin<2; timeout=5000; end;
            obj.ServerObj.stop(int32(timeout));
            % Explicitely delete the server object
            delete(obj.ServerObj); obj.ServerObj=[];
            obj.Status = false;
        end

        function delete(obj)
            % Destructor
            if obj.Status
                % Stop the server if it is running
                obj.stop();
            end
        end

        function conn = getConnection(obj,hashCode)
            % Get a WebSocketConnection to the client identified by the
            % HashCode
            if ~obj.Status; error('The server is not running!'); end
            try
                conn = WebSocketConnection(obj.ServerObj.getConnection(hashCode));
            catch err
                error(char(err.ExceptionObject.getMessage));
            end
        end

        function sendTo(obj,hashCode,message)
            % Directly send a message to a particular client, as identified
            % by its HashCode
            if ~obj.Status; error('The server is not running!'); end
            if ~isa(message,'char') && ~isa(message,'int8');
                error('You can only send character arrays or int8 arrays!');
            end
            try
                obj.ServerObj.sendTo(hashCode,message);
            catch err
                error(char(err.ExceptionObject.getMessage));
            end
        end

        function sendToAll(obj,message)
            % Send a message to all connected clients
            if ~obj.Status; error('The server is not running!'); end
            if ~isa(message,'char') && ~isa(message,'int8');
                error('You can only send character arrays or int8 arrays!');
            end
            obj.ServerObj.sendToAll(message);
        end

        function close(obj,hashCode)
            % Directly close connection to a particular client, as
            % identified by its HashCode
            if ~obj.Status; error('The server is not running!'); end
            try
                obj.ServerObj.close(hashCode);
            catch err
                error(char(err.ExceptionObject.getMessage));
            end
        end

        function closeAll(obj)
            % Close connection to all connected clients
            if ~obj.Status; error('The server is not running!'); end
            obj.ServerObj.closeAll();
        end
    end

    % Implement these methods in a subclass.
    methods (Abstract, Access = protected)
        onOpen(obj,conn,message)
        onTextMessage(obj,conn,message)
        onBinaryMessage(obj,conn,bytearray)
        onError(obj,conn,message)
        onClose(obj,conn,message)
    end

    % Private methods triggered by the callbacks defined above. This is
    % where the reactive behaviour of the server is defined.
    methods (Access = private)
        function openCallback(obj,e)
            % Define behavior in an onOpen method of a subclass
            obj.onOpen(WebSocketConnection(e.conn),char(e.message));
        end

        function textMessageCallback(obj,e)
            % Define behavior in an onTextMessage method of a subclass
            obj.onTextMessage(WebSocketConnection(e.conn),char(e.message));
        end

        function binaryMessageCallback(obj,e)
            % Define behavior in an onBinaryMessage method of a subclass
            obj.onBinaryMessage(WebSocketConnection(e.conn),e.blob.array);
        end

        function errorCallback(obj,e)
            % Define behavior in an onError method of a subclass
            obj.onError(WebSocketConnection(e.conn),char(e.message));
        end

        function closeCallback(obj,e)
            % Define behavior in an onClose method of a subclass
            obj.onClose(WebSocketConnection(e.conn),char(e.message));
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [It, Ix, Iy] = partial_deriv(images, uv_prev, interpolation_method, deriv_filter, b)

%PARTIAL_DERIV   Spatio-temporal derivatives
%   P = PARTIAL_DERIV(IMAGES, INIT) computes the spatio-temporal derivatives
%
%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2007-11-30 $
%   $Revision $
%
% Copyright 2007-2010, Brown University, Providence, RI. USA
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

if nargin == 2
    interpolation_method = 'cubic';
end;

if nargin == 4
    h = deriv_filter;
else
    h = [1 -8 0 8 -1]/12; % used in Wedel etal "improved TV L1"
end;

if nargin < 5
    b = 0.5;    %blending ratio
end;

H   = size(images, 1);
W   = size(images, 2);

[x,y]   = meshgrid(1:W,1:H);
x2      = x + uv_prev(:,:,1);
y2      = y + uv_prev(:,:,2);

% Record out of boundary pixels
B = (x2>W) | (x2<1) | (y2>H) | (y2<1);

if size(images, 4) ~= 1
    B = repmat(B, [1 1 size(images, 3)]);
end;

if size(images, 4) == 1
    img1 = images(:,:,1);
    img2 = images(:,:,2);
else
    img1 = images(:,:,:,1);
    img2 = images(:,:,:,2);
end;

if strcmp(interpolation_method, 'bi-cubic')

    %h = [-0.5, 0, 0.5]; % consistent with bi-cubic interpolation

    % Bicubic interpolation
    if nargout == 1

        if size(images, 4) == 1
            % gray-level
            warpIm = interp2_bicubic(images(:,:,2),x2,y2, h);
        else
            % color
            warpIm  = zeros(size(images(:,:,:,1)));
            for j = 1:size(images,3)
                warpIm(:,:,j) = interp2_bicubic(images(:,:,j, 2),x2,y2, h);
            end;
        end;

    elseif nargout == 3

        if size(images, 4) == 1
            % gray-level
            [warpIm Ix Iy] = interp2_bicubic(images(:,:,2),x2,y2, h);
        else
            % color
            warpIm  = zeros(size(images(:,:,:,1)));
            Ix      = warpIm;
            Iy      = warpIm;
            for j = 1:size(images,3)
                [warpIm(:,:,j) Ix(:,:,j) Iy(:,:,j)] = interp2_bicubic(images(:,:,j,2),x2,y2, h);
            end;
        end;

    else
        error('partial_deriv: number of output wrong!');
    end;

    indx        = isnan(warpIm);
    if size(images, 4) == 1
        It          = warpIm - images(:,:,1);
    else
        It          = warpIm - images(:,:,:,1);
    end;

    % Disable those out-of-boundary pixels in warping
    It(indx)    = 0;
    if nargout == 3

        % Temporal average
        I1x = imfilter(img1, h,  'corr', 'symmetric', 'same');  %
        I1y = imfilter(img1, h', 'corr', 'symmetric', 'same');

        % b = 0.6;    % recommended in Wedal etal "improved TV-L1" 2008
        % b = 0.5;
%         b = 1;

        Ix  = b*Ix+(1-b)*I1x;
        Iy  = b*Iy+(1-b)*I1y;

        Ix(indx) = 0;
        Iy(indx) = 0;
    end;

elseif strcmp(interpolation_method, 'bi-linear') || strcmp(interpolation_method, 'cubic')

    if strcmp(interpolation_method, 'bi-linear')
        method = 'linear';
    else
        method = 'cubic';
    end;

    % Matlab built-in
    if size(images, 4) == 1
        % Gray-level
        warpIm  = interp2(x,y,images(:,:,2),x2,y2,method);
        It      = warpIm - images(:,:,1);

    else
        % Color
        warpIm  = zeros(size(images(:,:,:,1)));
        for j = 1:size(images,3)
            warpIm(:,:,j) = interp2(x,y,images(:,:,j,2),x2,y2,method, NaN);
        end;
        It      = warpIm - images(:,:,:,1);

    end;

        % spline-based bicubic interpolation code with incorrect derivatives (below)
%         if size(images, 4) == 1
%             % gray-level
%             warpIm = interp2_bicubic(images(:,:,2),x2,y2, h);
%             It      = warpIm - images(:,:,1);
%         else
%             % color
%             warpIm  = zeros(size(images(:,:,:,1)));
%             for j = 1:size(images,3)
%                 warpIm(:,:,j) = interp2_bicubic(images(:,:,j, 2),x2,y2, h);
%             end;
%             It      = warpIm - images(:,:,:,1);
%         end;

    if nargout == 3

        % First compute derivative then warp
        I2x = imfilter(img2, h,  'corr', 'symmetric', 'same');
        I2y = imfilter(img2, h', 'corr', 'symmetric', 'same');

        if size(images, 4) == 1
            % Gray-level
            Ix  = interp2(x,y,I2x,x2,y2,method);
            Iy  = interp2(x,y,I2y,x2,y2,method);

%             % spline-based bicubic interpolation code with incorrect derivatives
%             Ix  = interp2_bicubic(I2x,x2,y2, h);
%             Iy  = interp2_bicubic(I2y,x2,y2, h);

        else
            % Color
            Ix  = zeros(size(images(:,:,:,1)));
            Iy  = Ix;

            for j = 1:size(images,3)
                Ix(:,:,j)  = interp2(x,y,I2x(:,:,j),x2,y2,method);
                Iy(:,:,j)  = interp2(x,y,I2y(:,:,j),x2,y2,method);
            end;

        end;

        % warp then derivative
%         if size(images, 4) == 1
%             tmp           = images(:,:,1);
%         else
%             tmp           = images(:,:,:,1);
%         end;
%         warpIm(B) = tmp(B);
%         Ix        = imfilter(warpIm, h,  'corr', 'symmetric', 'same');  %
%         Iy        = imfilter(warpIm, h', 'corr', 'symmetric', 'same');
        % end of warp then derivative

        % Temporal average
        I1x = imfilter(img1, h,  'corr', 'symmetric', 'same');  %
        I1y = imfilter(img1, h', 'corr', 'symmetric', 'same');

        % b = 0.6;    % recommended in Wedal etal "improved TV-L1" 2008
        % b = 0.5;

        Ix  = b*Ix+(1-b)*I1x;
        Iy  = b*Iy+(1-b)*I1y;

    end;

    % Disable those out-of-boundary pixels in warping
    It(B)   = 0;
    if nargout == 3
        Ix(B) = 0;
        Iy(B) = 0;
    end;

else
    error('partial_deriv: unknown interpolation method!');
end;


text:
function tf = cvx_use_sparse( sz, nz, isr )

%CVX_USE_SPARSE   Sparse/dense matrix efficiency test.
%   CVX_USE_SPARSE(SZ,NZ,ISR) is an internal function used by CVX to determine
%   if a matrix of size SZ with NZ nonzeros would require less memory to store
%   in sparse form or dense form. The matrix is assumed to be real if ISR=1 or
%   ISR is omitted, and complex otherwise.

if nargin == 1,
    ss = size( sz );
    if length( ss ) > 2,
        tf = false;
        return
    end
    isr = isreal( sz );
    if issparse( sz ),
        nz = nzmax( sz );
    else
        nz = nnz( sz );
    end
    sz = ss;
elseif length( sz ) > 2,
    tf = false;
    return
elseif nargin < 3,
    isr = true;
end
if isr,
    tf = 1 + ( 1 - 2 * sz( 1 ) ) * sz( 2 ) + 3 * nz < 0;
else
    tf = 1 + ( 1 - 4 * sz( 1 ) ) * sz( 2 ) + 5 * nz < 0;
end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function send_mail_message(id,subject,message,attachment)
%% SEND_MAIL_MESSAGE send email to gmail once calculation is done
% Example
% send_mail_message('its.neeraj','Simulation finished','This is the message area','results.doc')

% Pradyumna
% June 2008
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Your gmail ID and password
%(from which email ID you would like to send the mail)
mail = 'GMailId@gmail.com';    %Your GMail email address
password = 'GmailPassword';          %Your GMail password
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin == 1
    message = subject;
    subject = '';
elseif nargin == 2
    message = '';
    attachment = '';
elseif nargin == 3
    attachment = '';
end

% Send Mail ID
emailto = strcat(id,'@gmail.com');
%% Set up Gmail SMTP service.
% Then this code will set up the preferences properly:
setpref('Internet','E_mail',mail);
setpref('Internet','SMTP_Server','smtp.gmail.com');
setpref('Internet','SMTP_Username',mail);
setpref('Internet','SMTP_Password',password);

% Gmail server.
props = java.lang.System.getProperties;
props.setProperty('mail.smtp.auth','true');
props.setProperty('mail.smtp.socketFactory.class', 'javax.net.ssl.SSLSocketFactory');
props.setProperty('mail.smtp.socketFactory.port','465');

%% Send the email
if strcmp(mail,'GmailId@gmail.com')
    disp('Please provide your own gmail.')
    disp('You can do that by modifying the first two lines of the code')
    disp('after the comments.')
end

sendmail(emailto,subject,message,attachment)

end

repo: mcbench-benchmarks


text:
function [ mK ] = CreateConvMtx1D( vK, numElements, convShape )
% ----------------------------------------------------------------------------------------------- %
% [ mK ] = CreateConvMtx1D( vK, numElements, convShape )
% Generates a Convolution Matrix for 1D Kernel (The Vector vK) with
% support for different convolution shapes (Full / Same / Valid). The
% matrix is build such that for a signal 'vS' with 'numElements = size(vS
% ,1)' the following are equivalent: 'mK * vS' and conv(vS, vK,
% convShapeString);
% Input:
%   - vK                -   Input 1D Convolution Kernel.
%                           Structure: Vector.
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - numElements       -   Number of Elements.
%                           Number of elements of the vector to be
%                           convolved with the matrix. Basically set the
%                           number of columns of the Convolution Matrix.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, ...}.
%   - convShape         -   Convolution Shape.
%                           The shape of the convolution which the output
%                           convolution matrix should represent. The
%                           options should match MATLAB's conv2() function
%                           - Full / Same / Valid.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3}.
% Output:
%   - mK                -   Convolution Matrix.
%                           The output convolution matrix. The product of
%                           'mK' and a vector 'vS' ('mK * vS') is the
%                           convolution between 'vK' and 'vS' with the
%                           corresponding convolution shape.
%                           Structure: Matrix (Sparse).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References:
%   1.  MATLAB's 'convmtx()' - https://www.mathworks.com/help/signal/ref/convmtx.html.
% Remarks:
%   1.  The output matrix is sparse data type in order to make the
%       multiplication by vectors to more efficient.
%   2.  In case the same convolution is applied on many vectors, stacking
%       them into a matrix (Each signal as a vector) and applying
%       convolution on each column by matrix multiplication might be more
%       efficient than applying classic convolution per column.
% TODO:
%   1.
%   Release Notes:
%   -   1.1.000     19/07/2021  Royi Avital
%       *   Updated to use modern MATLAB arguments validation.
%   -   1.0.000     20/01/2019  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

arguments
    vK (:, :) {mustBeNumeric, mustBeVector}
    numElements (1, 1) {mustBeNumeric, mustBeReal, mustBePositive, mustBeInteger}
    convShape (1, 1) {mustBeNumeric, mustBeMember(convShape, [1, 2, 3])} = 1
end

CONVOLUTION_SHAPE_FULL         = 1;
CONVOLUTION_SHAPE_SAME         = 2;
CONVOLUTION_SHAPE_VALID        = 3;

kernelLength    = length(vK);

switch(convShape)
    case(CONVOLUTION_SHAPE_FULL)
        rowIdxFirst = 1;
        rowIdxLast  = numElements + kernelLength - 1;
        outputSize  = numElements + kernelLength - 1;
    case(CONVOLUTION_SHAPE_SAME)
        rowIdxFirst = 1 + floor(kernelLength / 2);
        rowIdxLast  = rowIdxFirst + numElements - 1;
        outputSize  = numElements;
    case(CONVOLUTION_SHAPE_VALID)
        rowIdxFirst = kernelLength;
        rowIdxLast  = (numElements + kernelLength - 1) - kernelLength + 1;
        outputSize  = numElements - kernelLength + 1;
end

mtxIdx = 0;

% The sparse matrix constructor ignores values of zero yet the Row / Column
% indices must be valid indices (Positive integers). Hence 'vI' and 'vJ'
% are initialized to 1 yet for invalid indices 'vV' will be 0 hence it has
% no effect.
vI = ones(numElements * kernelLength, 1);
vJ = ones(numElements * kernelLength, 1);
vV = zeros(numElements * kernelLength, 1);

for jj = 1:numElements
    for ii = 1:kernelLength
        if((ii + jj - 1 >= rowIdxFirst) && (ii + jj - 1 <= rowIdxLast))
            % Valid otuput matrix row index
            mtxIdx = mtxIdx + 1;
            vI(mtxIdx) = ii + jj - rowIdxFirst;
            vJ(mtxIdx) = jj;
            vV(mtxIdx) = vK(ii);
        end
    end
end

mK = sparse(vI, vJ, vV, outputSize, numElements);

end

author: RoyiAvital

repo: StackExchangeCodes


text:
function cordic_test007 ( )

%*****************************************************************************80
%
%% TEST007 demonstrates the use of EXP_CORDIC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 June 2007
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST007:\n' );
  fprintf ( 1, '  EXP_CORDIC computes the exponential function\n' );
  fprintf ( 1, '  using the CORDIC algorithm.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '      X    N         Exp(X)      Exp(X)          Difference\n' );
  fprintf ( 1, ...
    '                     Tabulated   CORDIC\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, theta, t1 ] = exp_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '\n' );

    for n = 0 : 5 : 25

      t2 = exp_cordic ( theta, n );
      d = t1 - t2;

      fprintf ( 1, '  %12f  %4d  %16.8e  %16.8e  %9e\n', theta, n, t1, t2, d );

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function linplus_test1565 ( )

%*****************************************************************************80
%
%% TEST1565 tests R8BUT_INDICATOR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  n = 6
  mu = 2

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST1565\n' );
  fprintf ( 1, '  R8BUT_INDICATOR sets up a R8BUT indicator matrix.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N     = %d\n', n );
  fprintf ( 1, '  Upper bandwidth MU = %d\n', mu );
%
%  Set the matrix.
%
  a = r8but_indicator ( n, mu );

  r8but_print ( n, mu, a, '  The indicator matrix:' );

  return
end


text:
function y1 = p00_rk_step ( test, neqn, order, t0, y0, t1 )

%*****************************************************************************80
%
%% P00_RK_STEP takes a single Runge-Kutta step from (T0,Y0) to (T1,Y1).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer TEST, the problem number.
%
%    Input, integer NEQN, the number of equations.
%
%    Input, integer ORDER, the order of the Runge-Kutta method to be
%    employed.  Legal values are 1 through 5.
%
%    Input, real T0, Y0(NEQN), the arguments of the derivative
%    function.
%
%    Input, real T1, the point at which an estimate of the solution
%    is desired.
%
%    Output, real Y1(NEQN), the estimated solution at T1.
%
  dt = t1 - t0;

  if ( order == 1 )

    yp0 = p00_fun ( test, neqn, t0, y0 );

    y1(1:neqn,1) = y0(1:neqn,1) + dt * yp0(1:neqn,1);

  elseif ( order == 2 )

    yp0 = p00_fun ( test, neqn, t0, y0 );

    yk1(1:neqn,1) = y0(1:neqn,1) + dt * yp0(1:neqn,1);

    tk1 = t0 + dt;

    ypk1 = p00_fun ( test, neqn, tk1, yk1 );

    y1(1:neqn,1) = y0(1:neqn,1) + dt * ( yp0(1:neqn,1) + ypk1(1:neqn,1) ) / 2.0;

  elseif ( order == 3 )

    yp0 = p00_fun ( test, neqn, t0, y0 );

    yk1(1:neqn,1) = y0(1:neqn,1) + 0.5 * dt * yp0(1:neqn,1);

    tk1 = t0 + 0.5 * dt;
    ypk1 = p00_fun ( test, neqn, tk1, yk1 );

    yk2(1:neqn,1) = y0(1:neqn,1) + dt * ( 2.0 * ypk1(1:neqn,1) - yp0(1:neqn,1) );

    tk2 = t0 + dt;
    ypk2 = p00_fun ( test, neqn, tk2, yk2 );

    y1(1:neqn,1) = y0(1:neqn,1) + ( dt / 6.0 ) ...
      * ( yp0(1:neqn,1) + 4.0 * ypk1(1:neqn,1) + ypk2(1:neqn,1) );

  elseif ( order == 4 )

    yp0 = p00_fun ( test, neqn, t0, y0 );

    yk1(1:neqn,1) = y0(1:neqn,1) + 0.5 * dt * yp0(1:neqn,1);

    tk1 = t0 + 0.5 * dt;

    ypk1 = p00_fun ( test, neqn, tk1, yk1 );

    yk2(1:neqn,1) = y0(1:neqn,1) + 0.5 * dt * ypk1(1:neqn,1);

    tk2 = t0 + 0.5 * dt;
    ypk2 = p00_fun ( test, neqn, tk2, yk2 );

    tk3 = t0 + dt;

    yk3(1:neqn,1) = y0(1:neqn,1) + dt * ypk2(1:neqn,1);

    ypk3 = p00_fun ( test, neqn, tk3, yk3 );

    y1(1:neqn,1) = y0(1:neqn,1) + ( dt / 6.0 ) * ( ...
                  yp0(1:neqn,1) ...
          + 2.0 * ypk1(1:neqn,1) ...
          + 2.0 * ypk2(1:neqn,1) ...
      +           ypk3(1:neqn,1) );

  elseif ( order == 5 )

    yp0 = p00_fun ( test, neqn, t0, y0 );

    yk1(1:neqn,1) = y0(1:neqn,1) + 0.25 * dt * yp0(1:neqn,1);

    tk1 = t0 + 0.25 * dt;

    ypk1 = p00_fun ( test, neqn, tk1, yk1 );

    yk2(1:neqn,1) = y0(1:neqn,1) + dt * ( ...
        3.0 * yp0(1:neqn,1) ...
      + 9.0 * ypk1(1:neqn,1) ) / 32.0;

    tk2 = t0 + 3.0 * dt / 8.0;

    ypk2 = p00_fun ( test, neqn, tk2, yk2 );

    yk3(1:neqn,1) = y0(1:neqn,1) + dt * ( ...
        1932.0 * yp0(1:neqn,1) ...
      - 7200.0 * ypk1(1:neqn,1) ...
      + 7296.0 * ypk2(1:neqn,1) ) / 2197.0;

    tk3 = t0 + 12.0 * dt / 13.0;

    ypk3 = p00_fun ( test, neqn, tk3, yk3 );

    yk4(1:neqn,1) = y0(1:neqn,1) + dt * ( ...
      + (  439.0 /  216.0 ) * yp0(1:neqn,1) ...
      -      8.0            * ypk1(1:neqn,1) ...
      + ( 3680.0 /  513.0 ) * ypk2(1:neqn,1) ...
      - (  845.0 / 4104.0 ) * ypk3(1:neqn,1) );

    tk4 = t0 + dt;

    ypk4 = p00_fun ( test, neqn, tk4, yk4 );

    yk5(1:neqn,1) = y0(1:neqn,1) + dt * ( ...
      - (    8.0 /   27.0 ) * yp0(1:neqn,1) ...
      + (    2.0          ) * ypk1(1:neqn,1) ...
      - ( 3544.0 / 2565.0 ) * ypk2(1:neqn,1) ...
      + ( 1859.0 / 4104.0 ) * ypk3(1:neqn,1) ...
      - (   11.0 /   40.0 ) * ypk4(1:neqn,1) );

    tk5 = t0 + 0.5 * dt;

    ypk5 = p00_fun ( test, neqn, tk5, yk5 );

    y1(1:neqn,1)  = y0(1:neqn,1) + dt * ( ...
        (    16.0 / 135.0   ) * yp0(1:neqn,1) ...
      + (  6656.0 / 12825.0 ) * ypk2(1:neqn,1) ...
      + ( 28561.0 / 56430.0 ) * ypk3(1:neqn,1) ...
      - (     9.0 / 50.0    ) * ypk4(1:neqn,1) ...
      + (     2.0 / 55.0    ) * ypk5(1:neqn,1) );

  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_RK_STEP - Fatal error!\n' );
    fprintf ( 1, '  Unavailable Runge Kutta order = %d\n', order );
    error ( 'P00_RK_STEP - Fatal error!' );
  end

  return
end


text:
function cyclic_reduction_test ( )

%*****************************************************************************80
%
%% CYCLIC_REDUCTION_TEST tests the CYCLIC_REDUCTION library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 May 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CYCLIC_REDUCTION_TEST\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the CYCLIC_REDUCTION library.\n' );

  cyclic_reduction_test01 ( );
  cyclic_reduction_test02 ( );
  cyclic_reduction_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CYCLIC_REDUCTION_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
classdef TestPolarToCart
    %TestPolarToCart

    methods (Static)
        function test_radians
            theta = linspace(0, pi, 10);
            rho = rand(size(theta));
            [x, y] = cv.polarToCart(rho, theta);
            validateattributes(x, {class(theta)}, {'size',size(theta)});
            validateattributes(y, {class(theta)}, {'size',size(theta)});

            theta = single(theta);
            rho = single(rho);
            [x, y] = cv.polarToCart(rho, theta);
            validateattributes(x, {class(theta)}, {'size',size(theta)});
            validateattributes(y, {class(theta)}, {'size',size(theta)});
        end

        function test_degrees
            theta = linspace(0, 180, 10);
            rho = rand(size(theta));
            [x, y] = cv.polarToCart(rho, theta, 'Degrees',true);
            validateattributes(x, {class(theta)}, {'size',size(theta)});
            validateattributes(y, {class(theta)}, {'size',size(theta)});
        end

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

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [abaqusData]=importAbaqusDat(varargin)

% function [abaqusData]=importAbaqusDat(fileName,optionStruct)
%-------------------------------------------------------------------------
%
% Form:
%         NODE FOOT-  COOR1          COOR2          COOR3
%         NOTE
%
%         1     -5.0000000E+00 -5.0000000E+00 -5.0000000E+00
%         2     -5.0000000E+00 -1.1555553E+00 -5.0000000E+00
%         ....
%         n     -5.0000000E+00 -1.1555553E+00 -5.0000000E+00
%
%         MAXIMUM          6.533       6.533       2.000
%         AT NODE                13          20          49
%
%         MINIMUM         -5.000      -5.000      -5.000
%         AT NODE                 1           1           111

% Form:
%         ELEMENT  PT FOOT-       E11         E22         E33         E12         E13         E23
%         NOTE
%
%         1   1         0.1427      0.1427     -0.3566      3.0358E-18  1.9120E-17  1.7250E-18
%         1   2         0.1427      0.1427     -0.3566      2.0817E-17  3.9274E-17 -2.5351E-17
%         ....
%         n   8         0.1427      0.1427     -0.3566      6.5052E-17  3.6108E-17 -4.4747E-17
%
%         MAXIMUM                0.1427      0.1427     -0.3566      9.8446E-17  2.0358E-16  2.1077E-16
%         ELEMENT                     6           1           3          24           7          15
%
%         MINIMUM                0.1427      0.1427     -0.3566     -8.8471E-17 -1.6717E-16 -1.7420E-16
%         ELEMENT                     8           3          10           2          12          11
%
%-------------------------------------------------------------------------

%% Parse input

switch nargin
    case 1
        fileName=varargin{1};
        optionStruct=[];
    case 2
        fileName=varargin{1};
        optionStruct=varargin{2};
end

defaultOptionStruct.nodeData.dataTypeStr='N O D E   O U T P U T';
defaultOptionStruct.nodeData.startTarget='    NODE FOOT-';
defaultOptionStruct.elementData.dataTypeStr='E L E M E N T   O U T P U T';
defaultOptionStruct.elementData.startTarget='    ELEMENT  PT';
defaultOptionStruct.stopTarget=' MAXIMUM';
defaultOptionStruct.getElementData=1;
defaultOptionStruct.getNodeData=1;

[optionStruct]=structComplete(optionStruct,defaultOptionStruct,1); %Complement provided with default if missing or empty

%%

T=txtfile2cell(fileName);

for qStep=1:1:1

    incrementLines = find(~cellfun(@isempty,regexp(T,'INCREMENT.*SUMMARY')));
    incrementLines(end+1)=numel(T);

    for qIncrement=1:1:numel(incrementLines)-1

        %Get increment data
        s=T{incrementLines(qIncrement)};
        %Scans this type of line: '                                INCREMENT     1 SUMMARY'
        abaqusData.STEP(qStep).INCREMENT(qIncrement).id=str2double(cell2mat(regexp(s,'\d+\.?\d*|-\d+\.?\d*|\.?\d+|-\.?\d+','match')));

        %Scans this type of line: ' TIME INCREMENT COMPLETED  5.000E-02,  FRACTION OF STEP COMPLETED  5.000E-02'
        s=T{incrementLines(qIncrement)+3};
        d=cellfun(@str2double,regexp(s,'([\+-]?((\d*\.\d+)|\d+))(E[\+-]?\d+)?','match'));
        abaqusData.STEP(qStep).INCREMENT(qIncrement).TIME_INCREMENT_COMPLETED=d(1);
        abaqusData.STEP(qStep).INCREMENT(qIncrement).FRACTION_OF_STEP_COMPLETED=d(2);

        %Scans this type of line: '  STEP TIME COMPLETED       5.000E-02,  TOTAL TIME COMPLETED        5.000E-02'
        s=T{incrementLines(qIncrement)+4};
        d=cellfun(@str2double,regexp(s,'([\+-]?((\d*\.\d+)|\d+))(E[\+-]?\d+)?','match'));
        abaqusData.STEP(qStep).INCREMENT(qIncrement).STEP_TIME_COMPLETED=d(1);
        abaqusData.STEP(qStep).INCREMENT(qIncrement).TOTAL_TIME_COMPLETED=d(2);

        %Get current data containing part
        T_sub=T(incrementLines(qIncrement):incrementLines(qIncrement+1)-1);

        %Get data
        stopTarget=optionStruct.stopTarget;
        if optionStruct.getNodeData==1 %Get node data
            dataTypeStr=optionStruct.nodeData.dataTypeStr;
            startTarget=optionStruct.nodeData.startTarget;
            try
                abaqusData.STEP(qStep).INCREMENT(qIncrement).nodeOutput=getDataField(T_sub,dataTypeStr,startTarget,stopTarget);
            catch
                abaqusData.STEP(qStep).INCREMENT(qIncrement).nodeOutput=NaN;
            end
        end
        if optionStruct.getElementData==1 %Get element data
            dataTypeStr=optionStruct.elementData.dataTypeStr;
            startTarget=optionStruct.elementData.startTarget;
            try
                abaqusData.STEP(qStep).INCREMENT(qIncrement).elementOutput=getDataField(T_sub,dataTypeStr,startTarget,stopTarget);
            catch
                abaqusData.STEP(qStep).INCREMENT(qIncrement).elementOutput=NaN;
            end
        end
    end
end
end

%%

function S=getDataField(T_sub,dataTypeStr,startTarget,stopTarget)

dataOutput = any(~cellfun(@isempty,regexp(T_sub,dataTypeStr)));
if dataOutput==1 %Element data found
    dataOutputStartLines = find(~cellfun(@isempty,regexp(T_sub,startTarget))); %Start lines
    dataOutputEndLines = find(~cellfun(@isempty,regexp(T_sub,stopTarget)))-2; %End lines
    dataOutputEndLines=dataOutputEndLines(dataOutputEndLines>dataOutputStartLines(1));
    ignoreSet={'FOOT-'};
    for q_dataOutput=1:1:numel(dataOutputStartLines)
        s=T_sub{dataOutputStartLines(q_dataOutput)}; %Header line
        s = regexprep(s,'\s+',' '); %Remove excessive spaces
        s=s(2:end-1); %Remove initial and final space
        headerSet=strsplit(s,' '); %Header set in cell
        %         numHeadersAll=numel(headerSet);
        logicRemove=gcontains(headerSet,ignoreSet);
        headerSet=headerSet(~logicRemove);
        numHeadersKeep=numel(headerSet);
        s=T_sub(dataOutputStartLines(q_dataOutput)+3:dataOutputEndLines(q_dataOutput)); %Data text field
        try
            D=cell2mat(cellfun(@(x) sscanf(x,'%f')',s,'uniformOutput',0)); %Data array
        catch %Alternative which ignores character components
            D=cell2mat(cellfun(@str2double,regexp(s,'([\+-]?((\d*\.\d+)|\d+))(E[\+-]?\d+)?','match'),'UniformOutput',0)); %Data array
        end

        %Create output structure with subfields names after headers
        for q=1:1:numHeadersKeep
            fieldNameNow=headerSet{q} ;
            S(q_dataOutput).data.(fieldNameNow)=D(:,q);
        end

    end
end
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON