text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spherical Harmonic Modeling and Analysis Toolkit (SPHARM-MAT) is a 3D
% shape modeling and analysis toolkit.
% It is a software package developed at Shenlab in Center for Neuroimaging,
% Indiana University (SpharmMat@gmail.com, http://www.iupui.edu/~shenlab/)
% It is available to the scientific community as copyright freeware
% under the terms of the GNU General Public Licence.
%
% Copyright 2009, 2010, ShenLab, Center for Neuroimaging, Indiana University
%
% This file is part of SPHARM-MAT.
%
% SPHARM-MAT 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.
%
% SPHARM-MAT 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 SPHARM-MAT. If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [inObjs1, inObjs2] = LoadObjList(inputFile)

fin = fopen(inputFile, 'rt');
inObjs1={};
inObjs2={};

[numOut, vals] = Get_Vals(fin, 'inputs', 10^6);
if numOut >0
    inObjs1 = vals;
%    set(handles.IOL_tag, 'String', handles.userdata.inObjs);
end

[numOut, vals] = Get_Vals(fin, 'inputs2', 10^6);
if numOut >0
    inObjs2 = vals;
%    set(handles.IOL_tag, 'String', handles.userdata.inObjs);
end

fclose(fin);

return;

function [numOut, vals] = Get_Vals(fInput, Keyword, NumVals)

%global userdata;

numOut = 0;
vals = {};

fseek(fInput, 0, 'bof');
delim = ' ,';

while 1
    tline = fgetl(fInput);
    if ~ischar(tline),   break,   end;
    [key, remain] = strtok(tline, delim);
    if strcmpi(strtrim(key), strtrim(Keyword))
        while numOut < NumVals
            tline = fgetl(fInput);
            if ~ischar(tline),   break,   end;
            numOut = numOut + 1;
            vals(numOut) = {deblank(tline)};
        end
        break;
    end
end

return;

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [ node_num, tetra_num ] = tet_mesh_order10_example_size ( )

%*****************************************************************************80
%
%% TET_MESH_ORDER10_EXAMPLE_SIZE sizes an example quadratic tet mesh.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NODE_NUM, the number of nodes.
%
%    Output, integer TETRA_NUM, the number of tetrahedrons.
%
  node_num = 27;
  tetra_num = 6;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
 function free(ob)
%function free(ob)
% free static memory associated with object

disp('freeing Gtomo3 object static memory')
f3d_mex('free');

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [ps,ix] = dpsimplify(p,tol)

% Recursive Douglas-Peucker Polyline Simplification, Simplify
%
% [ps,ix] = dpsimplify(p,tol)
%
% dpsimplify uses the recursive Douglas-Peucker line simplification
% algorithm to reduce the number of vertices in a piecewise linear curve
% according to a specified tolerance. The algorithm is also know as
% Iterative Endpoint Fit. It works also for polylines and polygons
% in higher dimensions.
%
% In case of nans (missing vertex coordinates) dpsimplify assumes that
% nans separate polylines. As such, dpsimplify treats each line
% separately.
%
% For additional information on the algorithm follow this link
% http://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm
%
% Input arguments
%
%     p     polyline n*d matrix with n vertices in d
%           dimensions.
%     tol   tolerance (maximal euclidean distance allowed
%           between the new line and a vertex)
%
% Output arguments
%
%     ps    simplified line
%     ix    linear index of the vertices retained in p (ps = p(ix))
%
% Examples
%
% 1. Simplify line
%
%     tol    = 1;
%     x      = 1:0.1:8*pi;
%     y      = sin(x) + randn(size(x))*0.1;
%     p      = [x' y'];
%     ps     = dpsimplify(p,tol);
%
%     plot(p(:,1),p(:,2),'k')
%     hold on
%     plot(ps(:,1),ps(:,2),'r','LineWidth',2);
%     legend('original polyline','simplified')
%
% 2. Reduce polyline so that only knickpoints remain by
%    choosing a very low tolerance
%
%     p = [(1:10)' [1 2 3 2 4 6 7 8 5 2]'];
%     p2 = dpsimplify(p,eps);
%     plot(p(:,1),p(:,2),'k+--')
%     hold on
%     plot(p2(:,1),p2(:,2),'ro','MarkerSize',10);
%     legend('original line','knickpoints')
%
% 3. Simplify a 3d-curve
%
%     x = sin(1:0.01:20)';
%     y = cos(1:0.01:20)';
%     z = x.*y.*(1:0.01:20)';
%     ps = dpsimplify([x y z],0.1);
%     plot3(x,y,z);
%     hold on
%     plot3(ps(:,1),ps(:,2),ps(:,3),'k*-');
%
%
%
% Author: Wolfgang Schwanghart, Markus Greim, 13. July, 2010.
% w.schwanghart[at]unibas.ch

if nargin == 0
    help dpsimplify
    return
end

if nargin ~= 2
    error('wrong number of input arguments')
end

% error checking
if ~isscalar(tol)
    error('tol must be a scalar')
end

% nr of dimensions
nrvertices    = size(p,1);
dims    = size(p,2);

% anonymous function for starting point and end point comparision
% compare = @(a,b) (a+eps >= b && a <= b) || ...
%                 (a-eps <= b && a >= b);
% now in simplifyrec()       for GNU/Octave compatibility

% __________________________________
% what happens, when there are NaNs?
% NaNs divide polylines.
Inan      = any(isnan(p),2);
% any NaN at all?
Inanp     = any(Inan);

% if there is only one vertex
if nrvertices == 1 || isempty(p);
    ps = p;
    ix = 1;

% if there are two
elseif nrvertices == 2 && ~Inanp;
    % when the line has no vertices (except end and start point of the line
    % check if the distance between both is less than the tolerance. If so
    % return the center
    if dims == 2;
        d    = hypot(p(1,1)-p(2,1),p(1,2)-p(2,2));
    else
        d    = sqrt(sum((p(1,:)-p(2,:)).^2));
    end

    if d <= tol;
        ps = sum(p,1)/2;
        ix = 1;
    else
        ps = p;
        ix = [1;2];
    end

elseif Inanp;

    % case: there are nans in the p array
    % --> find start and end indices of contiguous non-nan data
    Inan = ~Inan;
    sIX = strfind(Inan',[0 1])' + 1;
    eIX = strfind(Inan',[1 0])';

    if Inan(end)==true;
        eIX = [eIX;nrvertices];
    end

    if Inan(1);
        sIX = [1;sIX];
    end

    % calculate length of non-nan components
    lIX = eIX-sIX+1;
    % put each component into a single cell
    c   = mat2cell(p(Inan,:),lIX,dims);

    % now call dpsimplify again via cellfun.
    if nargout == 2;
        [ps,ix]   = cellfun(@(x) dpsimplify(x,tol),c,'uniformoutput',false);
        ix        = cellfun(@(x,six) x+six-1,ix,num2cell(sIX),'uniformoutput',false);
    else
        ps   = cellfun(@(x) dpsimplify(x,tol),c,'uniformoutput',false);
    end

    % write the data from a cell array to a matrix
    ps = cellfun(@(x) [x;nan(1,dims)],ps,'uniformoutput',false);
    ps = cell2mat(ps);
    ps(end,:) = [];

    % ix wanted? write ix to a matrix, too.
    if nargout == 2;
        ix = cell2mat(ix);
    end

else

% if there are no nans than start the recursive algorithm
ixe     = size(p,1);
ixs     = 1;

% logical vector for the vertices to be retained
I   = true(ixe,1);

% call recursive function
[p,I]   = simplifyrec(p,tol,ixs,ixe, dims, I);
ps  = p(I,:);

% if desired return the index of retained vertices
if nargout == 2;
    ix  = find(I);
end

end
end
% _________________________________________________________
function [p,I]  = simplifyrec(p,tol,ixs,ixe, dims, I);

mycompare = @(a,b) (a+eps >= b && a <= b) || ...
                 (a-eps <= b && a >= b);

    % check if startpoint and endpoint are the same
    % better comparison needed which included a tolerance eps

    c1 = num2cell(p(ixs,:));
    c2 = num2cell(p(ixe,:));

    % same start and endpoint with tolerance
    sameSE = all(cell2mat(cellfun(mycompare,c1(:),c2(:),'UniformOutput',false)));

    if sameSE;
        % calculate the shortest distance of all vertices between ixs and
        % ixe to ixs only
        if dims == 2;
            d    = hypot(p(ixs,1)-p(ixs+1:ixe-1,1),p(ixs,2)-p(ixs+1:ixe-1,2));
        else
            d    = sqrt(sum(bsxfun(@minus,p(ixs,:),p(ixs+1:ixe-1,:)).^2,2));
        end
    else
        % calculate shortest distance of all points to the line from ixs to ixe
        % subtract starting point from other locations
        pt = bsxfun(@minus,p(ixs+1:ixe,:),p(ixs,:));

        % end point
        a = pt(end,:)';

        beta = (a' * pt')./(a'*a);
        b    = pt-bsxfun(@times,beta,a)';
        if dims == 2;
            % if line in 2D use the numerical more robust hypot function
            d    = hypot(b(:,1),b(:,2));
        else
            d    = sqrt(sum(b.^2,2));
        end
    end

    % identify maximum distance and get the linear index of its location
    [dmax,ixc] = max(d);
    ixc  = ixs + ixc;

    % if the maximum distance is smaller than the tolerance remove vertices
    % between ixs and ixe
    if dmax <= tol;
        if ixs ~= ixe-1;
            I(ixs+1:ixe-1) = false;
        end
    % if not, call simplifyrec for the segments between ixs and ixc (ixc
    % and ixe)
    else
        [p,I]   = simplifyrec(p,tol,ixs,ixc, dims, I);
        [p,I]   = simplifyrec(p,tol,ixc,ixe, dims, I);

    end

end
% end

author: Sable

repo: mcbench-benchmarks


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  out.m           1993-12-28  %       FORMATED OUTPUT OF MATRICES
%%    (c)           M. Balda    %          updated on 2005-12-18
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%          updated on 2006-02-10
%
% The function out.m outputs its first argument to screen or a file in a
% default or required format. The argument can be either empty or a string,
% integer, real or complex numbers or matrices. The format can be either
% user defined or from a previous call of the function, or a default one.
%
% Prerequisity:
% ~~~~~~~~~~~~~
%   function inp.m from  http://www.mathworks.com/matlabcentral/
%
% Forms of calls:
% ~~~~~~~~~~~~~~
%
%   out;        %   Clear the format and close output file if any
%
%   out(arg);   %   Display arg in formats '%s or '%e' or '%d'
%               %   or any other used in the last call of out.m
%       arg   = an argument
%
%   out(arg,form);      %   Display arg in the format form
%       form  = required format of the output (see help to fprintf).
%               A format for complex numbers is generated internally using
%               format for real numbers,
%       fid   = file identifier - handle of the output file
%
%   out(arg,fid);       %   Append arg to the open file of fid handle
%
%   fid = out(arg,form,fname);  % Open file, get handle fid, and write arg
%       fname = name of output file
%       fid   = handle of the file if successfully opened,
%              -1 if not opened,
%               0 if closed.
%
%   out(arg,form,fid);  %   Append arg to fid file
%
% Examples:
% ~~~~~~~~
%  out(rand(2,3));              % display matrix in the '%15.6e' format
%  Z = rand(3,2)+i*rand(3,2);   % create complex matrix Z
%  h = out(Z,'','cmplxZ.dat');  % write Z in '%15.6e' format to the file
%  out('','%s\n',h);            % append new line into the file
%  out(Z,'%8.4f',h);            % append Z once more in '%8.4f' format
%  out;                         % close the file 'cmplxZ.dat'
%  lf = char(10);               % new line character
%  out(['This is a string' lf]);% display a string ending with new line
%  out(1:5,'%4i\n');            % display a column of integers
%  out(char('='*ones(1,50)),'%s\n');% display a line of 50 characters '='
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function fid = out(arg,form,fname)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
persistent fid_ frm_

if nargin<1                     %   Close file & return
    fid = 0;
    if ~isempty(fid_) && fid_~=1
        fid  = fclose(fid_);
    end
    frm_ = [];
    return
end

%                               %   Format test & reconstruction
S = isempty(findstr('s',frm_));
if (ischar(arg) && S) || (isnumeric(arg) && ~S)
    frm_ = [];
end

%                               %   OUTPUT DESTINATION
fid_ = 1;                       %   screen
if nargin>2                     %   write in a file?
    if ischar(fname)             %   file name
        if exist(fname,'file')
            if strcmp(inp(['Overwrite ' fname '?'],'no'),'no')
                fid=-1;
                return          %   without overwriting
            end
        end
        fid_ = fopen(fname, 'w+');% open file
        fid  = fid_;
        if fid_<0, return, end
    else
        if isreal(fname)        %   file identifier
            fid_ = fname;
        end
    end
end

%                               %   OUTPUT FORMAT
if nargin>1
    if isempty(form)
        frm_ = [];
    elseif ischar(form)          %   format
        if isnumeric(arg)       %   for numeric arg
            if isempty(findstr(form,'-') + findstr(form,'+'))
                frmR_ = ['% ' form(2:end)]; %   real
                frmI_ = ['%+' form(2:end)]; %   imag
            else
                frmR_ = ['% ' form(3:end)]; %   real
                frmI_ = ['%+' form(3:end)]; %   imag
            end
            if ~isreal(arg)
                frm_ = [frmR_ frmI_ 'i  '];  %  complex
            else
                frm_ = frmR_;
            end
        else
            frm_ = form;        %   for string arg
        end
    else
        fid_ = form;            %   file identifier
    end
end

if isempty(frm_)
    if ischar(arg),        frm_ = '%s  ';
    elseif isinteger(arg), frm_ = '% d  ';
    elseif isreal(arg),    frm_ = '% e  ';
    else                   frm_ = '% e %+ei  ';
    end
end
%   disp({arg, frm_, fid_})

%                               %   OUTPUT arg
if ischar(arg)
    fprintf(fid_,frm_,arg);     %   text argument
else
    [m,n] = size(arg);
    if isreal(arg)              %   integer or real arg
        for k = 1:m
            for j = 1:n
                fprintf(fid_, frm_, arg(k,j));
            end
            fprintf(fid_,'\n');
        end
    else                        %   complex arg
        for k = 1:m
            Rc = real(arg(k,:));
            Ic = imag(arg(k,:));
            for j = 1:n
                fprintf(fid_, frm_, Rc(j), Ic(j));
            end
            fprintf(fid_,'\n');
        end
    end
end
fid = fid_;


text:
function f = vargplvmObjective(params, model)

% VARGPLVMOBJECTIVE Wrapper function for variational GP-LVM objective.
% FORMAT
% DESC provides a wrapper function for the variational GP-LVM, it
% takes the negative of the log likelihood, feeding the parameters
% correctly to the model.
% ARG params : the parameters of the variational GP-LVM model.
% ARG model : the model structure in which the parameters are to be
% placed.
% RETURN f : the negative of the log likelihood of the model.
%
% SEEALSO : vargplvmCreate, vargplvmLogLikelihood, vargplvmExpandParam
%
% COPYRIGHT : Michalis K. Titsias, 2009-2011
% COPYRIGHT : Neil D. Lawrence, 2009-2011

% VARGPLVM

model = modelExpandParam(model, params);
f = - vargplvmLogLikelihood(model);

author: SheffieldML

repo: vargplvm


text:
function y = poly_env( p, x )

%POLY_ENV Evaluate the convex or concave envelope of a polynomial.
%   POLY_ENV( P, X ) uses a semidefinite program to compute the value of the
%   convex or concave envelope of the polynomial represented by the vector
%   P. The format of the vector P is identical to that required by POLYVAL,
%   with two additional restrictions. First, the elements of P must be real
%   and finite. Second, the length of P must 0, 2, or odd (1,3,5,...)
%
%   If the polynomial described by P is convex or concave, then a call to
%   POLY_ENV( P, X ) produces the same result as POLYVAL( P, X ).
%
%   POLY_ENV looks at the first nonzero element of P (presumably, but not
%   necessarily, P(1)) to determine if a convex or concave envelope is to
%   be selected. If P(1)>0, then a convex envelope is produced---i.e., the
%   function that is the tightest convex lower bound for POLYVAL(P,X).
%   Otherwise, a concave envelope is produced---the function that is the
%   tightest concave upper bound for POLYVAL(P,X).
%
%   If the degree N of P is odd---except for the special case N==1---no
%   proper convex/concave envelope exists. Therefore, POLY_ENV returns an
%   error in such cases.
%
%   If X is an array, the result will be computed elementwise.
%
%   Disciplined convex programming information:
%       POLY_ENV(P,X) is convex or concave and nonmotonic in X; therefore,
%       in CVX expressions, X must be affine, unless the polynomial
%       described by P has a degree of 0 or 1. P must be constant.

%
% Check the polynomial
%

sp = size( p );
if isempty( p ),
    p = zeros( 1, 0 );
elseif ~isa( p, 'double' ) || ~isreal( p ) || length( sp ) > 2 || ~any( sp == 1 ),
    cvx_throw( 'First argument must be a non-empty real vector.' );
elseif any( isnan( p ) | isinf( p ) ),
    cvx_throw( 'First argument must not contain Inf or NaN.' );
end
n = prod( sp );
if n > 2 && rem( n, 2 ) == 0,
    cvx_throw( 'The length of the vector p must be odd.' );
end

%
% Check the second argument
%

if n > 1 && ( ~cvx_isaffine( x ) || ~isreal( x ) ),
    cvx_throw( 'The second argument must be real and affine.' );
end

%
% Handle the special cases
%

sx = size( x );
ndxs = find( p );
if isempty( ndxs ),
    y = zeros( sx );
    return
else
    for k = ndxs(:)',
        pt = p( k : end );
        if ~any( isinf( pt ./ pt(1) ) ),
            p = pt;
            break;
        end
    end
end
n = length( p );
switch n,
    case 0,
        y = zeros(sx);
        return
    case 1,
        y = p(1) * ones(sx);
        return
    case 2,
        y = p(1) * x + p(2);
        return
    case 3,
        % Quadratic
        b2a = p(2) ./ ( 2 * p(1) );
        y = p(1) * ( square( x + b2a ) - b2a * b2a ) + p(3);
        return
    otherwise,
        if all( p(2:end) == 0 ),
            y = p(1) * ( x .^ (n-1) ) + p(end);
            return
        end
end

%
% Build and solve the SDP
%

degr  = n - 1;
deg2  = 0.5 * degr + 1;
nv    = prod( sx );
psign = sign(p(1));
p     = psign * reshape( p, 1, n );
cvx_begin sdp
    variable y(sx);
    variable P(deg2,deg2,sx) hankel;
    if cvx_isconstant( x )
        minimize(sum(y)); %#ok
    else
        minimize(y); %#ok
    end
    P >= 0; %#ok
    1 == P(1,1,:); %#ok
    x == reshape( P(2,1,:), sx ); %#ok
    y == reshape( p(end:-1:1) * [ reshape( P(1,:,:), deg2, nv ) ; reshape( P(2:end,end,:), deg2-1, nv ) ], sx ); %#ok
cvx_end

y = cvx_optval * psign;

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

author: yu-jiang

repo: radpbook

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


text:
function value = p21_i4vec ( action, name, dim_num, value )

%*****************************************************************************80
%
%% P21_I4VEC sets or gets I4VEC parameters for problem 21.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string ACTION, the action.
%    'D' sets the object to a default value.
%    'G' means the current value of the object should be returned.
%    'S' means the input values of the object and its dimension
%    should be stored.
%
%    Input, string NAME, the name of the parameter.
%    'E' is the exponent vector.
%
%    Input, integer DIM_NUM, the dimension of the object.
%
%    Input/output, integer VALUE(DIM_NUM), the value of the object.
%
  persistent e
  persistent dim_num_save

  if ( size ( dim_num_save ) == 0 )
    dim_num_save = 0;
  end

  if ( dim_num_save ~= dim_num )
    dim_num_save = 0;
    e = [];
  end

  if ( dim_num_save == 0 )
    dim_num_save = dim_num;
  end

  if ( action == 'D' | action == 'd' )

    if ( name == 'E' | name == 'e' | name == '*' )
      e(1:dim_num) = 2;
    end

  elseif ( action == 'G' | action == 'g' )

    if ( name == 'E' | name == 'e' )
      value(1:dim_num) = e(1:dim_num);
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P21_I4VEC - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P21_I4VEC - Fatal error!' );
    end

  elseif ( action == 'S' | action == 's' )

    if ( name == 'E' | name == 'e' )
      e(1:dim_num) = value(1:dim_num);
    else
      fprintf ( 1, '\n' );
      fprintf ( 1, 'P21_I4VEC - Fatal error!\n' );
      fprintf ( 1, '  Unrecognized name = "%s".\n', name );
      error ( 'P21_I4VEC - Fatal error!' );
    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P21_I4VEC - Fatal error!\n' );
    fprintf ( 1, '  Unrecognized action = "%s".\n', action );
    error ( 'P21_I4VEC - Fatal error!' );

  end

  return
end


text:
% function to implement the measurement model
function obs = f_h2(x)
% Input:
% x  = [xc xcd zc zcd p1 p2 w].'

load avar % r1,r2, L, and T

% unpack
xc = x(1); zc = x(3); p1 = x(5); p2 = x(6);

% observation 2
obs = L*(xc+r2*cos(p2))/(zc+r2*sin(p2));


text:
function [newnode,newelem]=mergemesh(node,elem,varargin)
%
% [newnode,newelem]=mergemesh(node,elem,varargin)
%
% concatenate two or more tetrahedral meshes or triangular surfaces
%
% author: Qianqian Fang <q.fang at neu.edu>
%
% input:
%      node: node coordinates, dimension (nn,3)
%      elem: tetrahedral element or triangle surface (nn,3) to (nn,5)
%
% output:
%      newnode: the node coordinates after merging, dimension (nn,3)
%      newelem: tetrahedral element or surfaces after merging (nn,4) or (nhn,5)
%
% note: you can call meshcheckrepair for the output newnode and
% newelem to remove the duplicated nodes or elements. mergemesh does
% detect self-intersecting elements when merging; to remove self-intersecting
% elements, you need to use mergesurf().
%
% example:
%
%   [node1,face1,elem1]=meshabox([0 0 0],[10 10 10],1,1);
%   [node2,face2,elem2]=meshasphere([5 5 13.1],3,0.3,3);
%   [newnode,newelem]=mergemesh(node1,elem1,node2,elem2);
%   plotmesh(newnode,newelem);
%   figure;
%   [newnode,newface]=mergemesh(node1,face1,node2,face2);
%   plotmesh(newnode,newface,'x>5');
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

len=length(varargin);
newnode=node;
newelem=elem;
if(len>0 && mod(len,2)~=0)
   error('you must give node and element in pairs');
end

X=mesheuler(newelem);

if(size(newelem,2)==4)
   if(X>=0)
      newelem(:,end+1)=1;
   end
end
if(size(newelem,2)==3)
   newelem(:,end+1)=1;
end
for i=1:2:len
   no=varargin{i};
   el=varargin{i+1};
   baseno=size(newnode,1);
   if(size(no,2)~=size(newnode,2))
        error('input node arrays have inconsistent columns');
   end
   if(size(el,2)==5 || size(el,2)==4)
        el(:,1:4)=el(:,1:4)+baseno;
	if(size(el,2)==4 && X>=0)
	   el(:,5)=1+(i+1)/2;
	end
   	newnode=[newnode;no];
	newelem=[newelem;el];
   elseif(size(el,2)==3 && size(newelem,2)==4)
        el(:,1:3)=el(:,1:3)+baseno;
	if(size(el,2)==3)
	   el(:,4)=1+(i+1)/2;
	end
   	newnode=[newnode;no];
	newelem=[newelem;el];
   else
        error('input element arrays have inconsistent columns');
   end
end

author: fangq


text:
%PEAK2  Find peaks in a matrix
%
% ZP = PEAK2(Z, OPTIONS) are the peak values in the 2-dimensional signal Z.
%
% [ZP,IJ] = PEAK2(Z, OPTIONS) as above but also returns the indices of the
% maxima in the matrix Z.  Use SUB2IND to convert these to row and column
% coordinates
%
% Options::
% 'npeaks',N    Number of peaks to return (default all)
% 'scale',S     Only consider as peaks the largest value in the horizontal
%               and vertical range +/- S points.
% 'interp'      Interpolate peak (default no interpolation)
% 'plot'        Display the interpolation polynomial overlaid on the point data
%
% Notes::
% - A maxima is defined as an element that larger than its eight neighbours.
%   Edges elements will never be returned as maxima.
% - To find minima, use PEAK2(-V).
% - The interp options fits points in the neighbourhood about the peak with
%   a paraboloid and its peak position is returned.  In this case IJ will
%   be non-integer.
%
% See also PEAK, SUB2IND.

% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

% Copyright (c) Peter Corke 1/96

function [zp,xypout, aout] = peak2(z, varargin)

    % process input options
    opt.npeaks = 2;
    opt.scale = 1;
    opt.interp = false;

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

    % create a neighbourhood mask for non-local maxima
    % suppression
    h = opt.scale;
    w = 2*h+1;
    M = ones(w,w);
    M(h+1,h+1) = 0;

    % compute the neighbourhood maximum
    znh = iwindow(double(z), M, 'max', 'wrap');

    % find all pixels greater than their neighbourhood
    k = find(z > znh);

    % sort these local maxima into descending order
    [zpk,ks] = sort(z(k), 'descend');

    k = k(ks);

    npks = min(length(k), opt.npeaks);
    k = k(1:npks);

    [y,x] = ind2sub(size(z), k);
    xy = [x y]';

    % interpolate the peaks if required
    if opt.interp

        xyp = [];
        zp = [];
        ap = [];

        % for each previously identified peak x(i), y(i)
        for xyt=xy
            % fit a polynomial to the local neighbourhood
            try

                x = xyt(1); y = xyt(2);

                % now try to interpolate the peak over a 3x3 window

                zc = z(x,   y);
                zn = z(x,   y-1);
                zs = z(x,   y+1);
                ze = z(x+1, y);
                zw = z(x-1, y);

                dx = (ze - zw)/(2*(2*zc - ze - zw));
                dy = (zn - zs)/(2*(zn - 2*zc + zs));

                zest = zc - (ze - zw)^2/(8*(ze - 2*zc + zw)) - (zn - zs)^2/(8*(zn - 2*zc + zs));

                aest = min(abs([ze/2 - zc + zw/2, zn/2 - zc + zs/2]));

            catch
                % handle situation where neighbourhood falls off the data
                % vector
                warning('Peak at %f too close to edge of image, skipping', x(i));
                continue;
            end
            %

            % store x, y for the refined peak
            xyp = [xyp [x+dx; y+dy]];
            zp = [zp zest];
            ap = [ap aest];

        end
    else
        % no interpolation case
        xyp = xy;
        zp = z(k)';
        ap = [];

    end

    % return values
    if nargout > 1
        xypout = xyp;
    end
    if nargout > 2
        aout = ap;
    end

author: Allopart

repo: rbpf-gmapping


text:
function [ W ] = log_gmm_div_w(X,Priors,Mu,Sigma )
% LOG_GMM_DIV_W Computes the responsibility factor of each GMM component ot
%
%
%   input -----------------------------------------------------------------
%
%       o X         : (D x N),      Samples
%
%       o Priors    : (1 x K),      Weights
%
%       o Mu        : (D x K),      Means
%
%       o Sigma     : (D x D x K),  Covariance
%
%
%   output ----------------------------------------------------------------
%
%       o W         : (N x K),  responsibility factor
%

K = size(Priors,2);
D = size(Mu,1);
N = size(X,2);

% likelihood for all data points
% (N x 1)
Lik = gmm_pdf(X,Priors,Mu,Sigma);

% (N x K)
lik_k = zeros(N,K);

if D==1
    for k=1:K
        lik_k(:,k) = gaussPDF(X,Mu(k),Sigma(k))';
    end
else
    for k=1:K
        lik_k(:,k) = gaussPDF(X,Mu(:,k),squeeze(Sigma(:,:,k)))';
    end
end

%(N x K)
W = (lik_k .* repmat(Priors,N,1)) ./ repmat(Lik,1,K);

% nbStates = K;
% for i=1:nbStates
%     %Compute probability p(x|i)
%     Pxi(:,i) = gaussPDF(X, Mu(:,i), squeeze(Sigma(:,:,i)));
% end
% %Compute posterior probability p(i|x)
% Pix_tmp = repmat(Priors,[N 1]).*Pxi;
% Pix = Pix_tmp ./ repmat(sum(Pix_tmp,2),[1 nbStates]); % -> W
% %Compute cumulated posterior probability
% E = sum(Pix);

end

author: epfl-lasa


text:
function v = integral2(F, S)
%INTEGRAL2   Flux integral of a Chebfun3v object through a 2D-surface.
%   INTEGRAL2(F, S) computes the flux integral of the Chebfun3v object F
%   through the parametric surface S defined as a Chebfun2v object (with 3
%   components):
%                           /
%       INTEGRAL2(F, S) =  |   < F, dS >
%                          /S
%                  //
%               = ||   < F(S(x,y)), cross(S_x(x,y), S_y(x,y)) > dxdy.
%                 //D
%
% See also CHEBFUN3V/INTEGRAL, CHEBFUN3/INTEGRAL, CHEBFUN3/INTEGRAL2,
% CHEBFUN3/INTEGRAL3, CHEBFUN3/SUM, CHEBFUN3/SUM2 and CHEBFUN3/SUM3.

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

% Check that F is a chebfun3v object with 3 components:
if ( F.nComponents ~= 3 )
    error('CHEBFUN:CHEBFUN3V:integral2:dimChebfun3v', ...
        'The chebfun3v object must have 3 components.')
end

% Check that S is a chebfun2v object with 3 components:
if ( isa(S, 'chebfun2v') )
    if ( S.nComponents ~= 3 )
        error('CHEBFUN:CHEBFUN3V:integral2:dimChebfun2v', ...
            'The parametrisation must have 3 components.')
    end
else
    error('CHEBFUN:CHEBFUN3V:integral2:notaChebfun2v', ...
        'The parametrisation must be a chebfun2v object.')
end

% Get components of F:
F1 = F.components{1};
F2 = F.components{2};
F3 = F.components{3};

% Get components of the surface S:
S1 = S(1);
S2 = S(2);
S3 = S(3);

% Build the composition F(S):
FS1_handle = @(x,y) feval(F1, feval(S1, x, y), feval(S2, x, y), ...
    feval(S3, x, y));

FS2_handle = @(x,y) feval(F2, feval(S1, x, y), feval(S2, x, y), ...
    feval(S3, x, y));

FS3_handle = @(x,y) feval(F3, feval(S1, x, y), feval(S2, x, y), ...
    feval(S3, x, y));

FS = chebfun2v(FS1_handle, FS2_handle, FS3_handle, S1.domain);

% Normal vector to the surface:
dS = normal(S);

% By definition:
v = sum2(dot(FS, dS));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%GEOM_CLEAN Cleanup geometry objects.
%
%   [ SOUT ] = GEOM_CLEAN( SIN, BREP2BIN ) Cleans geometry (objects)
%   in SIN by removing polygon data only used for visualization (for
%   BIN, BREP, STEP, and IGES geometry object types). If the BREP2BIN
%   flag is set (default true) ascii brep data is converted to binary
%   format.

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function [textons] = unitex(fb,k)
% function [textons] = unitex(fb,k)
%
% Compute universal textons from the training images.

iids = imgList('train');

n = 100000;
nper = round(n/numel(iids));
n = nper * numel(iids);

d = numel(fb);
data = zeros(d,n);

c = 0;
for i = 1:numel(iids),
  iid = iids(i);
  fprintf(2,'Processing image %d/%d (iid=%d)...\n',i,numel(iids),iid);
  im = imgRead(iid,'gray');
  fim = fbRun(fb,im);
  npix = numel(im);
  p = randperm(npix);
  p = p(1:min(npix,nper));
  m = numel(p);
  for j = 1:d,
    data(j,c+1:c+m) = fim{j}(p);
  end
  c = c + m;
end
data = data(:,1:c);

fprintf(2,'Computing %d universal textons from %d samples...\n',k,c);
[unused,textons] = kmeansML(k,data,'maxiter',30,'verbose',1);

author: zouchuhang

repo: LayoutNet


text:
function trainingModel = loadTrainingData(param)
% Generates the structure that contains all the training data needed for
% Component Contribution.
%
% USAGE:
%
%    trainingModel = loadTrainingData(formation_weight)
%
% INPUT:
%    formation_weight:    the relative weight to give the formation energies (Alberty's data)
%                         compared to the reaction measurements (TECRDB)
%
% OUTPUT:
%    trainingModel:       structure with data for Component Contribution
%                         *.S   `m x n` stoichiometric matrix of training data
%                         *.cids: `m x 1` compound ids
%                         *.dG0_prime: `n x 1`
%                         *.T:  `n x 1`
%                         *.I:  `n x 1`
%                         *.pH:  `n x 1`
%                         *.pMg:  `n x 1`
%                         *.weights:  `n x 1`
%                         *.balance:  `n x 1`
%                         *.cids_that_dont_decompose: k x 1 ids of compounds that do not decompose

if ~exist('param','var')
    formation_weight = 1;
    use_cached_kegg_inchis=true;
    use_model_pKas_by_default=true;
else
    if ~isfield(param,'formation_weight')
        formation_weight = 1;
    end
    if ~isfield(params,'use_cached_kegg_inchis')
        use_cached_kegg_inchis = true;
        % use_cached_kegg_inchis = false;
    else
        use_cached_kegg_inchis=params.use_cached_kegg_inchis;
    end
    if ~isfield(params,'use_model_pKas_by_default')
        use_model_pKas_by_default = true;
    else
        use_model_pKas_by_default=params.use_model_pKas_by_default;
    end
end

TECRDB_TSV_FNAME = 'data/TECRDB.tsv';
FORMATION_TSV_FNAME = 'data/formation_energies_transformed.tsv';
REDOX_TSV_FNAME = 'data/redox.tsv';

WEIGHT_TECRDB = 1;
WEIGHT_FORMATION = formation_weight;
WEIGHT_REDOX = formation_weight;

R=8.31451;
%Energies are expressed in kJ mol^-1.*)
R=R/1000; % kJ/mol/K
%Faraday Constant (kJ/mol)
F=96.48; %kJ/mol

if ~exist(TECRDB_TSV_FNAME, 'file')
    error(['file not found: ', TECRDB_TSV_FNAME]);
end

if ~exist(FORMATION_TSV_FNAME, 'file')
    error(['file not found: ', FORMATION_TSV_FNAME]);
end

if ~exist(REDOX_TSV_FNAME, 'file')
    error(['file not found: ', REDOX_TSV_FNAME]);
end

% Read the raw data of TECRDB (NIST)
reactions = {};
rxns={};
cids = [];
cids_that_dont_decompose = [];
thermo_params = []; % columns are: dG'0, T, I, pH, pMg, weight, balance?

fid = fopen(TECRDB_TSV_FNAME, 'r');
% fields are:
% 1. URL
% 2. REF_ID
% 3. METHOD
% 4. EVAL
% 5. EC
% 6. ENZYME NAME
% 7. REACTION IN KEGG IDS
% 8. REACTION IN COMPOUND NAMES
% 9. K
% 10. K'
% 11. T
% 12. I
% 13. pH
% 14. pMg
% 15. ID

res = textscan(fid, '%s%s%s%s%s%s%s%s%f%f%f%f%f%f%f', 'delimiter','\t');
fclose(fid);

inds = find(~isnan(res{10}) .* ~isnan(res{11}) .* ~isnan(res{13}));

dG0_prime = -R * res{11}(inds) .* log(res{10}(inds)); % calculate dG'0
thermo_params = [dG0_prime, res{11}(inds), res{12}(inds), res{13}(inds), ...
                 res{14}(inds), WEIGHT_TECRDB * ones(size(inds)), ...
                 true(size(inds))];

% parse the reactions in each row
for i = 1:length(inds)
    sprs = reaction2sparse(res{7}{inds(i)});
    cids = unique([cids, find(sprs)]);
    reactions = [reactions, {sprs}];
    rxns = [rxns;strtrim(['TECRDB_' int2str(res{15}(inds(i)))])];
    if strcmp('TECRDB_4403',strtrim(['TECRDB_' int2str(res{15}(inds(i)))]))
        pause(0.1);
    end
end
fprintf('Successfully added %d values from TECRDB\n', length(inds));

% Read the Formation Energy data.
fid = fopen(FORMATION_TSV_FNAME, 'r');
fgetl(fid); % skip the first header line
% fields are:
% 1. cid
% 2. name
% 3. dG'0
% 4. pH
% 5. I
% 6. pMg
% 7. T
% 8. decompose?
% 9. compound_ref
% 10. remark

res = textscan(fid, '%f%s%f%f%f%f%f%f%s%s', 'delimiter','\t');
fclose(fid);

inds = find(~isnan(res{3}));
thermo_params = [thermo_params; [res{3}(inds), res{7}(inds), res{5}(inds), ...
                                 res{4}(inds), res{6}(inds), ...
                                 WEIGHT_FORMATION * ones(size(inds)), ...
                                 false(size(inds))]];

eval(['cid = {' regexprep(sprintf('''C%05d''; ',res{1}),'(;\s)$','') '};']);
for i = 1:length(inds)
    sprs = sparse([]);
    sprs(res{1}(inds(i))) = 1;
    reactions = [reactions, {sprs}];
    rxns = [rxns;['FORM_' cid{i}]];
end

cids = union(cids, res{1}');
cids_that_dont_decompose = res{1}(find(res{8} == 0));

fprintf('Successfully added %d formation energies\n', length(res{1}));

% Read the Reduction potential data.
fid = fopen(REDOX_TSV_FNAME, 'r');
fgetl(fid); % skip the first header line
% fields are:
% 1. name
% 2. CID_ox
% 3. nH_ox
% 4. charge_ox
% 5. CID_red
% 6. nH_red,
% 7. charge_red
% 8. E'0
% 9. pH
% 10. I
% 11. pMg
% 12. T
% 13. ref
res = textscan(fid, '%s%f%f%f%f%f%f%f%f%f%f%f%s', 'delimiter', '\t');
fclose(fid);

delta_e = (res{6} - res{3}) - (res{7} - res{4}); % delta_nH - delta_charge
dG0_prime = -F * res{8} .* delta_e;
thermo_params = [thermo_params; [dG0_prime, res{12}, res{10}, res{9}, ...
                                 res{11}, ...
                                 WEIGHT_REDOX * ones(size(dG0_prime)), ...
                                 false(size(dG0_prime))]];

for i = 1:length(res{1})
    sprs = sparse([]);
    sprs(res{2}(i)) = -1;
    sprs(res{5}(i)) = 1;
    cids = unique([cids, res{2}(i), res{5}(i)]);
    reactions = [reactions, {sprs}];
    rxns = [rxns;['REDOX_' strrep(res{1}{i},' ','_')]];
end

fprintf('Successfully added %d redox potentials\n', length(res{1}));

% convert the list of reactions in sparse notation into a full
% stoichiometric matrix, where the rows (compounds) are according to the
% CID list 'cids'.
S = zeros(length(cids), length(reactions));
for i = 1:length(reactions)
    r = reactions{i};
    S(ismember(cids, find(r)), i) = r(r ~= 0);
    %production of h20 in certain redox reactions
    if any(strcmp(rxns{i},{'REDOX_dimethyl_sulfoxide';'REDOX_Trimethylamine'}))
        S(cids==1,i)=1;
    end
end

trainingModel.S = sparse(S);

if ~isfield(trainingModel,'rxns')
    for i=1:size(trainingModel.S,2)
        trainingModel.rxns{i,1}=['rxn' int2str(i)];
    end
end
if ~isfield(trainingModel,'lb')
    trainingModel.lb=ones(size(trainingModel.S,2),1)*-inf;
end
if ~isfield(trainingModel,'ub')
    trainingModel.ub=ones(size(trainingModel.S,2),1)*inf;
end

trainingModel.rxns = rxns;
trainingModel.cids = cids';
trainingModel.dG0_prime = thermo_params(:, 1);
trainingModel.T = thermo_params(:, 2);
trainingModel.I = thermo_params(:, 3);
trainingModel.pH = thermo_params(:, 4);
trainingModel.pMg = thermo_params(:, 5);
trainingModel.weights = thermo_params(:, 6);
trainingModel.balance = thermo_params(:, 7);
trainingModel.cids_that_dont_decompose = cids_that_dont_decompose;

%remove some problematic entries
boolRemove=false(size(trainingModel.rxns,1),1);

%http://xpdb.nist.gov/enzyme_thermodynamics/enzyme_data1.pl?col=1.&T1=66DED_421	C00089 + C06215 = C00031 + C06215	sucrose(aq) + (2,6--D-fructosyl)n(aq) = D-glucose(aq) + (2,6--D-fructosyl)n+1(aq)
%C06215 present on both sides of the reaction
%Need proper ID's for (2,6--D-fructosyl)n(aq) and (2,6--D-fructosyl)n+1(aq)
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_695');

% http://xpdb.nist.gov/enzyme_thermodynamics/enzyme_data1.pl?col=1.&T1=93AND/BUL_1509	93AND/BUL	radioactivity	B	2.3.1.129	UDP-N-acetylglucosamine acyltransferase	C04688 + C00043 = C03688 + C04738	(R)-3-hydroxytetradecanoyl-[acyl-carrier-protein](aq) + UDP-N-acetyl-D-glucosamine(aq) = acyl-carrier-protein(aq) + UDP-3-O-(3-hydroxytetradecanoyl)-N-acetyl-D-glucosamine(aq)		0.007	296.15		7.4		2905
% http://xpdb.nist.gov/enzyme_thermodynamics/enzyme_data1.pl?col=1.&T1=93AND/BUL_1509	93AND/BUL	radioactivity	B	2.3.1.129	UDP-N-acetylglucosamine acyltransferase	C04688 + C00043 = C03688 + C04738	(R)-3-hydroxytetradecanoyl-[acyl-carrier-protein](aq) + UDP-N-acetyl-D-glucosamine(aq) = acyl-carrier-protein(aq) + UDP-3-O-(3-hydroxytetradecanoyl)-N-acetyl-D-glucosamine(aq)		0.004	296.15		8.5		2906
% http://xpdb.nist.gov/enzyme_thermodynamics/enzyme_data1.pl?col=1.&T1=93AND/BUL_1509	93AND/BUL	radioactivity	B	2.3.1.129	UDP-N-acetylglucosamine acyltransferase	C04688 + C00043 = C03688 + C04738	(R)-3-hydroxytetradecanoyl-[acyl-carrier-protein](aq) + UDP-N-acetyl-D-glucosamine(aq) = acyl-carrier-protein(aq) + UDP-3-O-(3-hydroxytetradecanoyl)-N-acetyl-D-glucosamine(aq)		0.002	296.15		9		2907
% http://xpdb.nist.gov/enzyme_thermodynamics/enzyme_data1.pl?col=1.&T1=93AND/BUL_1509	93AND/BUL	radioactivity	B	2.3.1.129	UDP-N-acetylglucosamine acyltransferase	C04688 + C00043 = C03688 + C04738	(R)-3-hydroxytetradecanoyl-[acyl-carrier-protein](aq) + UDP-N-acetyl-D-glucosamine(aq) = acyl-carrier-protein(aq) + UDP-3-O-(3-hydroxytetradecanoyl)-N-acetyl-D-glucosamine(aq)		0.0095	296.15		8		2908
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_2905');
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_2906');
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_2907');
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_2908');

%These were missing water on the rhs, fixed in TECRDB.tsv
%TECRDB_4540 85LIE		A	4.3.-.-	formaldehyde condensation with THF	C00101 + C00067 = C00143 + C00001	THF(aq) + formaldehyde(aq) = 5,10-CH2-THF(aq) + H2O(l)		30000	311.15	0.25	7		4540
%TECRDB_4541 +66KAL/JEN		E	4.3.-.-	formaldehyde condensation with THF	C00101 + C00067 = C00143 + C00001	THF(aq) + formaldehyde(aq) = 5,10-CH2-THF(aq) + H2O(l)		32000	298.15	1	7		4541
%TECRDB_4544 +59BLA		E	4.3.-.-	formaldehyde condensation with THF	C00101 + C00067 = C00143 + C00001	THF(aq) + formaldehyde(aq) = 5,10-CH2-THF(aq)		7700	293.15		7.2		4544

%The record of 85LIE is problematic i.e. TECRDB_870 - TECRDB_884
%e.g. http://xpdb.nist.gov/enzyme_thermodynamics/enzyme_data1.pl?col=1.&T1=85LIE_292	85LIE	spectrophotometry and enzymatic assay	A	1.4.4.2	glycine dehydrogenase (decarboxylating)	C00037 + C00725 = C80069 + C00288	glycine(aq) + lipoate(aq) = S-aminomethyldihydro--lipoate(aq) + carbon dioxide(aq)		0.031	311.15		6.39		870
%is C00725 really C00248 ???
%C80069 is InChI=1S/C9H19NO2S2/c10-7-14-6-5-8(13)3-1-2-4-9(11)12/h8,13H,1-7,10H2,(H,11,12)
%C80069	but is it  InChI=1S/C9H20N2OS2/c10-7-14-6-5-8(13)3-1-2-4-9(11)12/h8,13H,1-7,10H2,(H2,11,12)/p+1  ???
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_870');
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_871');
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_872');
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_873');
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_874');
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_875');
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_876');
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_877');
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_878');
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_879');
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_880');
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_881');
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_882');
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_883');
boolRemove = boolRemove | ismember(trainingModel.rxns,'TECRDB_884');
% The overall reaction reported in the abstract of the dissertation (https://digitalcommons.library.tmc.edu/dissertations/AAI8516325/) is as follows:
%85LIE		A	1.1.-.-	glycine cleavage system (a series of enzymes)	C00037 + C00003 + C00101 = C00143 + C00004 + C00014 + C00011 	glycine(aq) + NAD+(aq) + THF(aq) = 5,10-CH2-THF(aq) + NADH(aq) + NH3(aq) + CO2(aq)		0.00156	311.15	0.25	7		4539

%http://xpdb.nist.gov/enzyme_thermodynamics/enzyme_data1.pl?col=1.&T1=65STR_520	65STR	spectrophotometry	C	2.6.1.13	ornithine-oxo-acid transaminase	C00077 + C00026 = C03912 + C00025	L-ornithine(aq) + 2-oxoglutarate(aq) = DL-D-1-pyrroline-5-carboxylate(aq) + L-glutamate(aq)		71	310.15		7.1		2719
%is http://xpdb.nist.gov/enzyme_thermodynamics/enzyme_data1.pl?col=1.&T1=65STR_520	65STR	spectrophotometry	C	2.6.1.13	ornithine-oxo-acid transaminase	C00077 + C00026 = C01165 + C00025 + C00001	L-ornithine(aq) + 2-oxoglutarate(aq) = DL-D-1-pyrroline-5-carboxylate(aq) + L-glutamate(aq)		71	310.15		7.1		2719
%i.e. C00001 added to rhs, it is balanced, but is this correct?

%85LIE		A	1.1.-.-	glycine cleavage system	C00037 + C00003 + C00101 = C00143 + C00004 + C00014 + C00288 	glycine(aq) + NAD+(aq) + THF(aq) = 5,10-CH2-THF(aq) + NADH(aq) + NH3(aq) + CO2(aq)		0.00156	311.15	0.25	7		4539
%is 85LIE		A	1.1.-.-	glycine cleavage system	C00037 + C00003 + C00101 = C00143 + C00004 + C00014 + C00011 	glycine(aq) + NAD+(aq) + THF(aq) = 5,10-CH2-THF(aq) + NADH(aq) + NH3(aq) + CO2(aq)		0.00156	311.15	0.25	7		4539
%i.e. C00288 replaced by C00011

%http://xpdb.nist.gov/enzyme_thermodynamics/enzyme_data1.pl?col=1.&T1=07LIN/ALG_1584	07LIN/ALG	spectrophotometry	A	1.1.1.87	homoisocitrate dehydrogenase	 C05662 + C00003 = C00322 + C00288 + C00004	(1R,2S)-1-hydroxybutane-1,2,4-tricarboxylate(aq) + NAD(ox) = 2-oxoadipate(aq) + carbon dioxide(aq) + NAD(red)		0.45	298.15		7.5		1
%is http://xpdb.nist.gov/enzyme_thermodynamics/enzyme_data1.pl?col=1.&T1=07LIN/ALG_1584	07LIN/ALG	spectrophotometry	A	1.1.1.87	homoisocitrate dehydrogenase	C00001 + C05662 + C00003 = C00322 + C00288 + C00004	H2O(l) + (1R,2S)-1-hydroxybutane-1,2,4-tricarboxylate(aq) + NAD(ox) = 2-oxoadipate(aq) + carbon dioxide(aq) + NAD(red)		0.45	298.15		7.5		1
%i.e. lhs missing h20

%http://xpdb.nist.gov/enzyme_thermodynamics/enzyme_data1.pl?col=1.&T1=95PEL/MAC_1595	95PEL/MAC	spectrophotometry and enzymatic assay	B	1.5.1.5	methylenetetrahydrofolate dehydrogenase (NADP+)	C00143 + C00006 = C00234 + C00005	5,10-methylenetetrahydrofolate(aq) + NADP(ox)(aq) = 10-formyltetrahydrofolate(aq) + NADP(red)(aq)		16	303.15		7.3		267
%is http://xpdb.nist.gov/enzyme_thermodynamics/enzyme_data1.pl?col=1.&T1=95PEL/MAC_1595	95PEL/MAC	spectrophotometry and enzymatic assay	B	1.5.1.5	methylenetetrahydrofolate dehydrogenase (NADP+)	C00001 + C00143 + C00006 = C00234 + C00005	H2O(l) + 5,10-methylenetetrahydrofolate(aq) + NADP(ox)(aq) = 10-formyltetrahydrofolate(aq) + NADP(red)(aq)		16	303.15		7.3		267
%i.e. lhs missing h20

%remove some of the training reactions that are problematic, e.g. those involving ACP
trainingModel.S = trainingModel.S(:,~boolRemove);
trainingModel.lb = trainingModel.lb(~boolRemove);
trainingModel.ub = trainingModel.ub(~boolRemove);
trainingModel.rxns = trainingModel.rxns(~boolRemove);
trainingModel.dG0_prime = trainingModel.dG0_prime(~boolRemove);
trainingModel.T = trainingModel.T(~boolRemove);
trainingModel.I = trainingModel.I(~boolRemove);
trainingModel.pH = trainingModel.pH(~boolRemove);
trainingModel.pMg = trainingModel.pMg(~boolRemove);
trainingModel.weights = trainingModel.weights(~boolRemove);
trainingModel.balance = trainingModel.balance(~boolRemove);

author: opencobra

repo: cobratoolbox


text:
function hypograph( varargin )

%HYPOGRAPH Declares a hypograph variable.
%   HYPOGRAPH VARIABLE x
%   where x is a valid MATLAB variable name, declares a scalar
%   variable for the current cvx problem, and specifies it as
%   the objective of a maximization. There is no need to declare
%   a separate MAXIMIZE command.
%
%   This keyword should be used only when attempting to
%   create a new function for the CVX library. Suppose you have a
%   function F(X) whose hypograph {(X,Y)}{F(X)>=Y} can be
%   described as a CVX feasibiliy problem. Then declaring Y
%   as an EPIGRAPH turns the model into one that computes the
%   value of the function itself, maximizing over Y.
%
%   Other uses of this keyword are NOT supported and can lead
%   to numerical errors.
%
%   HYPOGRAPH VARIABLE x(n1,n2,...,nk)
%   declares a vector, matrix, or array hypograph variable with
%   dimensions n1, n2, ..., nk, each of which must be positive
%   integers.
%
%   Structure modifiers such as "symmetric", "toeplitz", etc.
%   are permitted with hypograph variables.
%
%   Examples:
%      hypograph variable x
%      hypograph variable x(100)
%
%   See also VARIABLE, EPIGRAPH.

if ~iscellstr( varargin ),
    error( 'HYPOGRAPH must be used in command mode.' );
elseif nargin ~= 2 || ~strcmpi( varargin{1}, 'variable' ),
    error( 'Syntax: hypograph variable <variable>' );
elseif ~isa( evalin( 'caller', 'cvx_problem', '[]' ), 'cvxprob' ),
    error( 'HYPOGRAPH can only be used within a CVX model.' );
else
    evalin( 'caller', sprintf( '%s ', 'variable', varargin{2:end}, ' hypograph_' ) );
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:
classdef VigdergauzParametersFromVolumeAndStrain < handle

    properties (Access = public)
        parameters
    end

    properties (Access = private)
        ax
        ay
        cx
        cy
        strain
        volume
        mu1
        mu0
        k1
        k0
    end

    methods (Access = public)

        function obj = VigdergauzParametersFromVolumeAndStrain(cParams)
            obj.init(cParams);
        end

        function compute(obj)
            obj.computeAxAy();
            obj.computeVigergauzParameters();
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.strain = cParams.strain;
            obj.volume = cParams.volumeMicro;
            obj.cx     = cParams.cx;
            obj.cy     = cParams.cy;
            k = @(E,nu) E/(2*(1-nu));
            mu = @(E,nu) E/(2*(1+nu));
            obj.mu0 = mu(cParams.E0,cParams.nu0);
            obj.k0  = k(cParams.E0,cParams.nu0);
            obj.mu1 = mu(cParams.E1,cParams.nu1);
            obj.k1  = k(cParams.E1,cParams.nu1);
        end

        function computeAxAy(obj)
            s.strain = obj.strain;
            s.volume = obj.volume;
            s.cx     = obj.cx;
            s.cy     = obj.cy;
            s.mu0    = obj.mu0;
            s.mu1    = obj.mu1;
            s.k0    = obj.k0;
            s.k1    = obj.k1;
            axay = AxAyComputerFromVolumeAndStrain(s);
            [obj.ax,obj.ay] = axay.compute();
        end

        function computeVigergauzParameters(obj)
            s.ax = obj.ax;
            s.ay = obj.ay;
            s.cx = obj.cx;
            s.cy = obj.cy;
            s.type = 'AxAndAy';
            obj.parameters = VigdergauzParameters.create(s);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
%%***************************************************************
%% logcheby: logarithmic Chebyshev approximation problem.
%%           SDP formulation.
%%    minimize   t
%%      x,t
%%    such that  1/t <= (x'*B(i,:))/f(i) <= t,  i = 1:p.
%%
%%    B = pxm matrix,  f = p-vector,  p > m.
%%
%%    Ref: Vanderberghe & Boyd.
%%--------------------------------------------------------------
%% [blk,Avec,C,b,X0,y0,Z0,obj,x] = logcheby(B,f,feas,solve);
%%
%% B = pxm matrix  (p > m).
%% f = px1 vector,  [B(:,j)./f must be positive for each j]
%% feas  = 1 if want feasible starting point
%%       = 0 if otherwise.
%% solve = 0 if just want initialization.
%%       = 1 if want to solve the problem.
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%***************************************************************

  function [blk,Avec,C,b,X0,y0,Z0,objval,x] = logcheby(B,f,feas,solve);

     if (nargin < 4); solve = 0; end;
     if (nargin < 3); feas = 1;  end;

     if ~isreal(B); error('B must be real'); end;

     [p,m] = size(B);
     blk{1,1} = 's'; blk{1,2} = 2*ones(1,p);
     blk{2,1} = 'l'; blk{2,2} = p;

     E = zeros(p,m);
     for j = [1:m];  E(:,j) = B(:,j)./f;  end;
     if any(E < 1e-10);
        error(' B(:,j)./f must have all entry positive');
     end;
     for i = [1:p];  beta(i) = sum(E(i,:)); end;
%%
     temp = zeros(2*p+1,1);
     temp(2:2:2*p) = ones(p,1);
     C{1,1} = spdiags(temp,1,2*p,2*p); C{1,1} = C{1,1} + C{1,1}';
     C{2,1} = zeros(p,1);
     b  = [zeros(m,1); 1];

     A = cell(2,m+1);
     temp = zeros(2*p,1);
     for k = 1:m
         temp(1:2:2*p-1) = -E(:,k);
         A{1,k} = spdiags(temp,0,2*p,2*p);
         A{2,k} = E(:,k);
     end;
     temp = zeros(2*p,1);
     temp(2:2:2*p) = ones(p,1);
     A{1,m+1} = spdiags(temp,0,2*p,2*p);
     A{2,m+1} = ones(p,1);
%%
    Avec = svec(blk,A,ones(size(blk,1),1));
    if (feas == 1);
        X0{1,1} = speye(2*p)/(2*p);
        X0{2,1} = ones(p,1)/(2*p);
        y0 = [ones(m,1); -1.1*max(beta)]/min(beta);
        Z0 = ops(C,'-',Atyfun(blk,Avec,[],[],y0));
    elseif (feas == 0);
        [X0,y0,Z0] = infeaspt(blk,Avec,C,b);
    end;
    if (solve)
       [obj,X,y,Z] = sqlp(blk,Avec,C,b,[],X0,y0,Z0);
       objval = -mean(obj);
       x = y(1:m);
    else
       objval = []; x = [];
    end
%%***************************************************************


text:
% to the value of the slider.
   function slider1_callback(hObject,eventdata)
   global fh slider1_data
      % Get 'slider' appdata.
      slider1_data = getappdata(fh,'slider');
      slider1_data.previous_val = slider1_data.val;
      slider1_data.val = get(hObject,'Value');
%      set(eth,'String',num2str(slider_data.val));
      sprintf('You changed the slider value--> %6.3f',...
              slider1_data.val)
      % Save 'slider' appdata before returning.
      setappdata(fh,'slider',slider1_data)
end


text:
function [h, hdata] = hbayes(net, hdata)
%HBAYES	Evaluate Hessian of Bayesian error function for network.
%
%	Description
%	H = HBAYES(NET, HDATA) takes a network data structure NET together
%	the data contribution to the Hessian for a set of inputs and targets.
%	It returns the regularised Hessian using any zero mean Gaussian
%	priors on the weights defined in NET.  In addition, if a MASK is
%	defined in NET, then the entries in H that correspond to weights with
%	a 0 in the mask are removed.
%
%	[H, HDATA] = HBAYES(NET, HDATA) additionally returns the data
%	component of the Hessian.
%
%	See also
%	GBAYES, GLMHESS, MLPHESS, RBFHESS
%

%	Copyright (c) Ian T Nabney (1996-2001)

if (isfield(net, 'mask'))
  % Extract relevant entries in Hessian
  nmask_rows = size(find(net.mask), 1);
  hdata = reshape(hdata(logical(net.mask*(net.mask'))), ...
     nmask_rows, nmask_rows);
  nwts = nmask_rows;
else
  nwts = net.nwts;
end
if isfield(net, 'beta')
  h = net.beta*hdata;
else
  h = hdata;
end

if isfield(net, 'alpha')
  if size(net.alpha) == [1 1]
    h = h + net.alpha*eye(nwts);
  else
    if isfield(net, 'mask')
      nindx_cols = size(net.index, 2);
      index = reshape(net.index(logical(repmat(net.mask, ...
         1, nindx_cols))), nmask_rows, nindx_cols);
    else
      index = net.index;
    end
    h = h + diag(index*net.alpha);
  end
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
clear
patchSize = 8;

% images courtesy of Stefan Roth
I = double((imread('new.jpg')))/255;
mask = double((imread('new_mask.png')))/255;

if size(I,3)>1
    I = rgb2ycbcr(I);
end

% find which patches are occluded for faster performance
noiseI = I;
mask_inds = find(mask>0);
tt = noiseI(:,:,1);
tt(mask_inds)=NaN;
ttt = im2col(tt,[patchSize patchSize]);
excludeList = find(any(isnan(ttt)));
clear ttt;

noiseI = I;
for i=1:size(I,3)
    tt = I(:,:,i);
    tt(mask_inds) = 0;
    noiseI(:,:,i)=tt;
end

% load ICA model and initialize MAP estimator
load ICAModel
W = W*E;
invW = pinv(W);
prior = @(Z,patchSize,noiseSD,imsize) PatchDCTGG(Z,patchSize,noiseSD,imsize,W,invW,excludeList);

%%
tic

% inpainting
lambda = 1000000*ones([size(I,1) size(I,2)]);
lambda(mask_inds)=0;
cleanI = zeros(size(I));
for i=1:size(I,3)
    cleanI(:,:,i) = EPLLhalfQuadraticSplit(noiseI(:,:,i),lambda,patchSize,10*[1 2 16 128 512],20,prior,I(:,:,i));
end

%%
if (size(I,3)>1)
    cleanI = ycbcr2rgb(cleanI);
    I = ycbcr2rgb(I);
    noiseI = ycbcr2rgb(noiseI);
end
toc

%% output result
figure(1);
imshow(I); title('Original');
figure(2);
imshow(noiseI); title('Corrupted Image');
figure(3);
imshow(cleanI); title('Restored Image');
fprintf('PSNR is:%f\n',20*log10(1/std2(cleanI-I)));

author: lbasek

repo: image-denoising-benchmark


text:
function [cell_case_study,xmesh,ymesh,zmesh] = dicomrt_loadRTdose(filename_list,force)
% dicomrt_loadRTdose(filename_list,force)
%
% Read dicom dose for a case study using MATLAB native function dicomread.
%
% filename contains a list of dose matrices to import.
% filename does not contain link to an RTPLAN in the first line
%
% Dose matrices are stored in a single 3D matrix: case_study
% x-y-z coordinates of the center of the dose-pixel are stored in xmesh,ymesh,zmesh
% cell_case_study is a simple cell array with the following structure:
%
%  ------------------------
%  | [ rtplan structure ] |
%  | ----------------------
%  | [ 3D dose matrix   ] |
%  | ---------------------
%  | [ voxel dimension  ] |
%  ------------------------
%
% NOTE: this function and dicomrt_loaddose differs because in this case the
% RTPLAN is not imported and left black
%
%
% See also dicomrt_loadct, dicomrt_getPatientPosition
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

% Retrieve the number of images
nlines=dicomrt_nASCIIlines(filename_list);

% Get dose images and create 3D Volume
fid=fopen(filename_list);
nfiles=0;
skip_header_mf=0;
DoseGridScaling_mf=1;
case_study=[];
CELLStore=0;

% import plan information
rtplan_location=fgetl(fid);

dictFlg = checkDictUse;
if dictFlg
    rtplan=dicominfo(rtplan_location, 'dictionary', 'ES - IPT4.1CompatibleDictionary.mat');
else
    rtplan=dicominfo(rtplan_location);
end

% retrieve number of beams and beam labels
beams=fieldnames(rtplan.BeamSequence);
nbeams=size(beams,1);

% Define cell array to store info
case_study_info=cell(1);
case_study_info{1}=rtplan;

if exist('force')==0
    [PatientPositionCODE]=dicomrt_getPatientPosition(rtplan);
elseif exist('force')~=0
    if force==0;
        [PatientPositionCODE]=dicomrt_getPatientPosition(rtplan);
    end
else
    warning('dicomrt_loaddose: Dicom load was forced. No PatientPosition check is performed!');
end

try
    fractions=rtplan.FractionGroupSequence.Item_1.NumberOfFractionsPlanned;
    if isempty(fractions)==1
        fractions=1;
    end
catch
    disp('This is a single franction plan');
    fractions=1;
end

% Progress bar
h = waitbar(0,['Loading progress:']);
set(h,'Name','dicomrt_loaddose: loading RTDOSE objects');

% loop until the end-of-file is reached and build 3D dose matrix
while (feof(fid)~=1);
    nfiles=nfiles+1; % counting
    nfilescheck=nfiles; % check for eof
    file_location{1,nfiles}=fgetl(fid);
    if isnumeric(file_location{1,nfiles}), nfiles=nfiles-1, break, end % end of line reached

    dictFlg = checkDictUse;
    if dictFlg
        % import dose information
        if isdeployed
            info_image = dicominfo(file_location{1,nfiles}, 'dictionary', ...
                fullfile(getCERRPath,'bin','ES - IPT4.1CompatibleDictionary.mat'));
        else
            info_image = dicominfo(file_location{1,nfiles}, 'dictionary',...
                'ES - IPT4.1CompatibleDictionary.mat');
        end
    else
        info_image=dicominfo(file_location{1,nfiles});
    end

    try
        image=dicomread(file_location{1,nfiles}); % import pixel data
        dimage=double(squeeze(image)); % convert pixel data to double
    catch
        dictFlg = checkDictUse;
        if dictFlg
            info_image_mf=dicominfo(file_location{1,nfiles}, 'dictionary', 'ES - IPT4.1CompatibleDictionary.mat');
        else
            info_image_mf=dicominfo(file_location{1,nfiles});
        end

        DoseGridScaling_mf=info_image_mf.DoseGridScaling;
        skip_header_mf=1;
        % Give option to load individual beam/segments or total dose
        disp('This appear to be a Multi-frame study with dose contribution');
        disp('Do you want to:');
        disp('(1) load the dose matrix as a sum of all beams [default]');
        disp('(2) keep individual dose contribution');
        % with option (2) it is supposed that only beam contribution is exported
        % no support for segment contribution is given for the time being
        % because no TPS is currently capable of doing this
        STOREoption = input('Option: ');
        if STOREoption == 2
            CELLStore = 1;
            case_study = cell(nbeams,2);
        else
            CELLStore = 0;
        end
    end
    if skip_header_mf==1 & nfiles==1
    elseif skip_header_mf==1 & nfiles>1 & CELLStore == 1
        % Multi-frame with dose contribution
        if strcmpi(info_image.DoseSummationType,'FRACTION')==1 | ...
                strcmpi(info_image.DoseSummationType,'BEAM')==1
            dimage=dimage.*info_image.DoseGridScaling.*DoseGridScaling_mf;
        elseif strcmpi(info_image.DoseSummationType,'PLAN')==1
            dimage=dimage./info_image.DoseGridScaling.*fractions.*DoseGridScaling_mf;
        end
        case_study_info{nfiles}=info_image;
        case_study{nfiles-1,1}=getfield(rtplan,'BeamSequence',char(beams(nfiles-1)),'BeamName');
        case_study{nfiles-1,2}=dimage;
        zmesh=info_image.ImagePositionPatient(3)+info_image.GridFrameOffsetVector;
    elseif skip_header_mf==1 & nfiles>1 & CELLStore == 0
        % Multi-frame with dose contribution
        if strcmpi(info_image.DoseSummationType,'FRACTION')==1 | ...
                strcmpi(info_image.DoseSummationType,'BEAM')==1
            dimage=dimage.*info_image.DoseGridScaling.*DoseGridScaling_mf;
        elseif strcmpi(info_image.DoseSummationType,'PLAN')==1
            dimage=dimage./info_image.DoseGridScaling.*fractions.*DoseGridScaling_mf;
        end
        case_study_info{nfiles+1}=info_image;
        if length(case_study)==0
            case_study=dimage;
        else
            case_study=case_study+dimage;
        end
        zmesh=info_image.ImagePositionPatient(3)+info_image.GridFrameOffsetVector;
    elseif  skip_header_mf==0
        if length(size(image))==4
            % this is a Multi-frame study without dose components
            if strcmpi(info_image.DoseSummationType,'FRACTION')==1 | ...
                    strcmpi(info_image.DoseSummationType,'BEAM')==1
                dimage=dimage.*info_image.DoseGridScaling.*DoseGridScaling_mf;
            elseif strcmpi(info_image.DoseSummationType,'PLAN')==1
                dimage=dimage./info_image.DoseGridScaling.*fractions.*DoseGridScaling_mf;
            end
            case_study_info{nfiles+1}=info_image;
            case_study=dimage;
            zmesh=info_image.ImagePositionPatient(3)+info_image.GridFrameOffsetVector;
        else % not a multi-frame !
            if strcmpi(info_image.DoseSummationType,'FRACTION')==1 | ...
                    strcmpi(info_image.DoseSummationType,'BEAM')==1
                dimage=dimage.*info_image.DoseGridScaling.*DoseGridScaling_mf;
            elseif strcmpi(info_image.DoseSummationType,'PLAN')==1
                dimage=dimage./info_image.DoseGridScaling.*fractions.*DoseGridScaling_mf;
            end
            case_study_info{nfiles+1}=info_image;
            case_study(:,:,nfiles)=dimage;
            zmesh(nfiles)=info_image.ImagePositionPatient(3);
        end
    end
    waitbar(nfiles/nlines,h);
end

% Make zmesh a column vector
zmesh=dicomrt_makevertical(zmesh);

fclose(fid);

if nfilescheck~=nfiles;
    warning('dicomrt_loaddose: End of file was prematurely reached. Check the expected dimensions of your data. It may be OK to continue');
end

% PatientOrientation     ImageOrientationPatient
%
% HFS                    (1,0,0) (0,1,0)
% FFS                    (1,0,0) (0,1,0)
% HFP                    (-1,0,0) (0,-1,0)
% FFP                    (-1,0,0) (0,-1,0)
%
%
if PatientPositionCODE == 1 | PatientPositionCODE == 2
    min_x=info_image.ImagePositionPatient(1);
    pixel_spacing_x=info_image.PixelSpacing(1);

    min_y=info_image.ImagePositionPatient(2);
    pixel_spacing_y=info_image.PixelSpacing(2);

    [xmesh] = dicomrt_create1dmesh(min_x,pixel_spacing_x,size(image,2),0);
    [ymesh] = dicomrt_create1dmesh(min_y,pixel_spacing_y,size(image,1),0);

else
    max_x=info_image.ImagePositionPatient(1);
    pixel_spacing_x=info_image.PixelSpacing(1);

    max_y=info_image.ImagePositionPatient(2);
    pixel_spacing_y=info_image.PixelSpacing(2);

    [xmesh] = dicomrt_create1dmesh(max_x,pixel_spacing_x,size(image,2),1);
    [ymesh] = dicomrt_create1dmesh(max_y,pixel_spacing_y,size(image,1),1);

end

zmesh=dicomrt_mmdigit(zmesh*0.1,7);
ymesh=dicomrt_mmdigit(ymesh*0.1,7);
xmesh=dicomrt_mmdigit(xmesh*0.1,7);

% Check support for current Patient Position
if PatientPositionCODE == 1 % supported Patient Position: HFS
    disp('dicomrt_loadRTdose: Patient Position is Head First Supine (HFS).');
elseif PatientPositionCODE == 2 % supported Patient Position: FFS
    disp('dicomrt_loadRTdose: Patient Position is Feet First Supine (FFS).');
elseif PatientPositionCODE == 3 % supported Patient Position: HFP
    disp('dicomrt_loadRTdose: Patient Position is Head First Prone (HFP).');
elseif PatientPositionCODE == 4 % supported Patient Position: FFP
    disp('dicomrt_loadRTdose: Patient Position is Feet First Prone (FFP).');
else
    warning('dicomrt_loaddose: Unable to determine Patient Position');
    PatientPosition=input('dicomrt_loadRTdose: Please specify Patient Position: HFS(default),FFS,HFP,FFP: ','s');
    if isempty(PatientPosition)==1
        PatientPosition='HFS';
    end

    if strcmpi(PatientPosition, 'HFS')
        PatientPositionCODE = 1;
    elseif strcmpi(PatientPosition, 'FFS')
        PatientPositionCODE = 2;
    elseif strcmpi(PatientPosition, 'HFP')
        PatientPositionCODE = 3;
    elseif strcmpi(PatientPosition, 'FFP')
        PatientPositionCODE = 4;
    end

end

% Create a cell array where to store rtplan and dose matrix
cell_case_study=cell(3,1);
cell_case_study{1,1}=case_study_info;
cell_case_study{2,1}=case_study;
cell_case_study{3,1}=[];
% Completed

% Close progress bar
close(h);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function pn = plane_imp_point_near_3d ( a, b, c, d, p )

%*****************************************************************************80
%
%% PLANE_IMP_POINT_NEAR_3D: nearest point on a implicit plane to a point in 3D.
%
%  Discussion:
%
%    The implicit form of a plane in 3D is:
%
%      A * X + B * Y + C * Z + D = 0
%
%    The normal N to the plane is (A,B,C).
%
%    The line defined by (XN-P(1))/A = (YN-P(2))/B = (ZN-P(3))/C = T
%    goes through P and is parallel to N.
%
%    Solving for the point (XN,YN,ZN) we get
%
%      XN = A*T+P(1)
%      YN = B*T+P(2)
%      ZN = C*T+P(3)
%
%    Now place these values in the equation for the plane:
%
%      A*(A*T+P(1)) + B*(B*T+P(2)) + C*(C*T+P(3)) + D = 0
%
%    and solve for T:
%
%      T = (-A*P(1)-B*P(2)-C*P(3)-D) / (A * A + B * B + C * C )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, D, the implicit plane parameters.
%
%    Input, real P(3), the coordinates of the point.
%
%    Output, real PN(3), the nearest point on the plane.
%
  dim_num = 3;

  if ( plane_imp_is_degenerate_3d ( a, b, c ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PLANE_IMP_POINT_NEAR_3D - Fatal error!\n' );
    fprintf ( 1, '  A = B = C = 0.\n' );
    error ( 'PLANE_IMP_POINT_NEAR_3D - Fatal error!' );
  end

  t = - ( a * p(1) + b * p(2) + c * p(3) + d ) / ( a * a + b * b + c * c );

  pn(1) = p(1) + a * t;
  pn(2) = p(2) + b * t;
  pn(3) = p(3) + c * t;

  return
end


text:
%% Facemark LBF training demo
%
% The user should provides the list of training images accompanied by their
% corresponding landmarks location in separate files.
%
% See below for a description of file formats.
%
% Examples of datasets are available at
% <https://ibug.doc.ic.ac.uk/resources/facial-point-annotations/>.
%
% Sources:
%
% * <https://github.com/opencv/opencv_contrib/blob/3.4.0/modules/face/samples/facemark_demo_lbf.cpp>
% * <https://github.com/opencv/opencv_contrib/blob/3.4.0/modules/face/tutorials/facemark_usage/facemark_usage.markdown>
%

%% Preparation
%
% Before you continue with this tutorial, you should download a training
% dataset of facial landmarks detection.
%
% We suggest you to download the IBUG dataset which can be retrieved at
% <https://ibug.doc.ic.ac.uk/download/annotations/ibug.zip>
%
% First thing to do is to make two text files containing the list of image
% files and annotation files respectively. Make sure that the order of images
% and annotations in both files are matched. Furthermore, it is advised to use
% absolute paths instead of relative paths.
%
% Example to make the file list in Linux machine:
%
%  ls /data/ibug/*.jpg > images.txt
%  ls /data/ibug/*.pts > annotations.txt
%
% Example of content in the |images.txt| file:
%
%  /data/ibug/image_003_1.jpg
%  /data/ibug/image_004_1.jpg
%  /data/ibug/image_005_1.jpg
%  ...
%
% Example of content in the |annotations.txt| file:
%
%  /data/ibug/image_003_1.pts
%  /data/ibug/image_004_1.pts
%  /data/ibug/image_005_1.pts
%  ...
%
% where a |.pts| file contains the position of each face landmark.
% Make sure that the annotation format is supported by the API, where the
% contents should look like the following snippet:
%
%  version: 1
%  n_points:  68
%  {
%  212.716603 499.771793
%  230.232816 566.290071
%  ...
%  }
%
% Once trained, we show how to use the model to detect face landmarks in a
% test image.
%
% You can also download a pre-trained model in this link
% <https://raw.githubusercontent.com/kurnianggoro/GSOC2017/master/data/lbfmodel.yaml>
% (that way you can skip training and simply load the model).
%

%% Options

% [INPUT] path of a text file contains the list of paths to all training images
imgList = fullfile(mexopencv.root(),'test','facemark','ibug','images.lst');
assert(exist(imgList, 'file') == 2, 'missing images list file');

% [INPUT] path of a text file contains the list of paths to all annotations files
ptsList = fullfile(mexopencv.root(),'test','facemark','ibug','annotations.lst');
assert(exist(ptsList, 'file') == 2, 'missing annotations list file');

% [OUTPUT] path for saving the trained model
modelFile = fullfile(tempdir(), 'model_lbf.yaml');

% [INPUT] path to the cascade xml file for the face detector
xmlFace = fullfile(mexopencv.root(),'test','lbpcascade_frontalface.xml');
download_classifier_xml(xmlFace);

% name of user-defined face detector function
faceDetectFcn = 'myFaceDetector';
assert(exist([faceDetectFcn '.m'], 'file') == 2, 'missing face detect function');

% path to test image
testImg = fullfile(mexopencv.root(),'test','lena.jpg');

%% Init
% create the facemark instance
obj = cv.Facemark('LBF', 'NLandmarks',68, 'CascadeFace',xmlFace, ...
    'ModelFilename',modelFile, 'SaveModel',true, 'Verbose',true);

%%
% set user-defined face detector
obj.setFaceDetector(faceDetectFcn);

%% Data
% load the dataset, and add training samples one-by-one
disp('Loading data...')
[imgFiles, ptsFiles] = cv.Facemark.loadDatasetList(imgList, ptsList);
for i=1:numel(imgFiles)
    % load image and its corresponding annotation data, then add pair
    img = cv.imread(imgFiles{i});
    pts = cv.Facemark.loadFacePoints(ptsFiles{i});
    obj.addTrainingSample(img, pts);
end

%% Train
% train the algorithm, model will be saved to specified file
disp('Training...')
tic
obj.training();
toc

%% Test
% run on some test image
img = cv.imread(testImg);
faces = obj.getFaces(img);
landmarks = obj.fit(img, faces);
for i=1:numel(faces)
    img = cv.rectangle(img, faces{i}, 'Color',[255 0 255]);
    img = cv.Facemark.drawFacemarks(img, landmarks{i}, 'Color',[0 0 255]);
end
imshow(img)

%% Helper functions

function download_classifier_xml(fname)
    if exist(fname, 'file') ~= 2
        % attempt to download trained Haar/LBP/HOG classifier from Github
        url = 'https://cdn.rawgit.com/opencv/opencv/3.4.0/data/';
        [~, f, ext] = fileparts(fname);
        if strncmpi(f, 'haarcascade_', length('haarcascade_'))
            url = [url, 'haarcascades/'];
        elseif strncmpi(f, 'lbpcascade_', length('lbpcascade_'))
            url = [url, 'lbpcascades/'];
        elseif strncmpi(f, 'hogcascade_', length('hogcascade_'))
            url = [url, 'hogcascades/'];
        else
            error('File not found');
        end
        urlwrite([url f ext], fname);
    end
end

% The facemark API provides the functionality to the user to use their own
% face detector. The code below implements a sample face detector. This
% function must be saved in its own M-function to be used by the facemark API.
function faces = myFaceDetector(img)
    persistent obj
    if isempty(obj)
        obj = cv.CascadeClassifier();
        obj.load(xmlFace);
    end

    if size(img,3) > 1
        gray = cv.cvtColor(img, 'RGB2GRAY');
    else
        gray = img;
    end
    gray = cv.equalizeHist(gray);
    faces = obj.detect(gray, 'ScaleFactor',1.4, 'MinNeighbors',2, ...
        'ScaleImage',true, 'MinSize',[30 30]);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function q=v_qrmult(q1,q2)
%V_QRMULT multiplies together two real quaternions matrices q=[q1,q2]
%
% Inputs:   q1(4m,n)  Two real quaternions arrays. Either array can
%           q2(4n,r)  also be a scalar quaternion.
%
% Outputs:   q(4m,r)  Matrix product of q1 and q2

%      Copyright (C) Mike Brookes 2000-2012
%      Version: $Id: v_qrmult.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
s1=size(q1);
s2=size(q2);
if isequal(s1,[4 1])
    q=v_qrdotmult(repmat(q1,s2(1)/4,s2(2)),q2);
elseif isequal(s2,[4 1])
    q=v_qrdotmult(q1,repmat(q2,s1(1)/4,s1(2)));
else
    q=v_rotqr2mr(q1)*q2;
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function pass = test_basicArithmetic( pref )
% Check that we can do the basic arithmetic on chebop2 objects.
% Alex Townsend, August 2013.

if ( nargin < 1 )
    pref = chebfunpref();
end
tol = 10*pref.cheb2Prefs.chebfun2eps;

N1 = chebop2(@(x,y,u) diff(u,2,2));
N2 = chebop2(@(x,y,u) x.*diff(u,2,1));
N = N1 + N2;

EXACT = chebop2(@(x,y,u) diff(u,2,2) + x.*diff(u,2,1));
ERROR = N - EXACT;
err = 0; C = ERROR.coeffs;
for jj = 1:size(C,1)
    for kk = 1:size(C,2)
        err = err + norm(C{jj,kk});
    end
end
pass(1) = ( abs(err) < tol );

N1 = chebop2(@(u) diff(u,2,2));
N2 = chebop2(@(x,y,u) x.*diff(u,2,1));
N = N1 + N2;

EXACT = chebop2(@(x,y,u) diff(u,2,2) + x.*diff(u,2,1));
ERROR = N - EXACT;
err = 0; C = ERROR.coeffs;
for jj = 1:size(C,1)
    for kk = 1:size(C,2)
        err = err + norm(C{jj,kk});
    end
end
pass(2) = ( abs(err) < tol );

N1 = chebop2(@(x,y,u) y.*diff(u,2,2));
N2 = chebop2(@(u) diff(u,2,1));
N = N1 + N2;

EXACT = chebop2(@(x,y,u) y.*diff(u,2,2) + diff(u,2,1));
ERROR = N - EXACT;
err = 0; C = ERROR.coeffs;
for jj = 1:size(C,1)
    for kk = 1:size(C,2)
        err = err + norm(C{jj,kk});
    end
end
pass(3) = ( abs(err) < tol );

end


text:
%
% Usage:   AAt =mexCalcAAt(A);
%
% Name: mexCalcAAt
%
% Description: Compute efficiently AAt = A*A', when A is sparse
%   and has a lot more columns than rows. In some cases, it is
%   up to 20 times faster than the equivalent Matlab expression
%   AAt=A*A';
%
% Inputs: A:  double sparse m x n matrix
%
% Output: AAt: double m x m matrix
%
% Author: Julien Mairal, 2009

author: tiepvupsu

repo: DICTOL

save_path: github-repos/MATLAB/tiepvupsu-DICTOL


text:
function BW = PQmovBW (X2)
% Bandwidth tests

% P. Kabal $Revision: 1.1 $  $Date: 2003/12/07 13:34:46 $

persistent kx kl FR FT N

if (isempty (kx))
    NF = 2048;
    Fs = 48000;
    fx = 21586;
    kx = round (fx / Fs * NF);    % 921
    fl = 8109;
    kl = round (fl / Fs * NF);    % 346
    FRdB = 10;
    FR = 10^(FRdB / 10);
    FTdB = 5;
    FT = 10^(FTdB / 10);
    N = NF / 2;     % Limit from pseudo-code
end

Xth = X2(2,kx+1);
for (k = kx+1:N-1)
    Xth = max (Xth, X2(2,k+1));
end

% BWRef and BWTest remain negative if the BW of the test signal
% does not exceed FR * Xth for kx-1 <= k <= kl+1
BW.BWRef = -1;
XthR = FR * Xth;
for (k = kx-1:-1:kl+1)
    if (X2(1,k+1) >= XthR)
        BW.BWRef = k + 1;
        break;
    end
end

BW.BWTest = -1;
XthT = FT * Xth;
for (k = BW.BWRef-1:-1:0)
    if (X2(2,k+1) >= XthT)
        BW.BWTest = k + 1;
        break;
    end
end

author: stephencwelch

repo: Perceptual-Coding-In-Python


text:
function line_ncc_rule_test02 ( )

%*****************************************************************************80
%
%% LINE_NCC_RULE_TEST02 estimates the integral of exp(x) from 0 to 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2014
%
%  Author:
%
%    John Burkardt
%
  a =  0.0;
  b = +1.0;
  exact = exp ( b ) - exp ( a );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST02\n' );
  fprintf ( 1, '  Use a sequence of NCC rules to compute an estimate Q\n' );
  fprintf ( 1, '  of the integral:\n' );
  fprintf ( 1, '    I = integral ( 0 <= x <= 1 ) exp(x) dx.\n' );
  fprintf ( 1, '  The exact value is:\n' );
  fprintf ( 1, '    I = %g\n', exact );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   N       Q             |Q-I|\n' );
  fprintf ( 1, '\n' );

  for n = 1 : 22

    [ x, w ] = line_ncc_rule ( n, a, b );

    q = w(1:n)' * exp ( x(1:n) );
    error = abs ( exact - q );
    fprintf ( 1, '  %2d  %14.6g  %14.6g\n', n, q, error );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function a = combin ( alpha, beta, n )

%*****************************************************************************80
%
%% COMBIN returns the combinatorial matrix.
%
%  Formula:
%
%    If ( I = J ) then
%      A(I,J) = ALPHA + BETA
%    else
%      A(I,J) = BETA
%
%  Example:
%
%    N = 5, ALPHA = 2, BETA = 3
%
%    5 3 3 3 3
%    3 5 3 3 3
%    3 3 5 3 3
%    3 3 3 5 3
%    3 3 3 3 5
%
%  Properties:
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).
%
%    det ( A ) = ALPHA**(N-1) * ( ALPHA + N * BETA ).
%
%    LAMBDA(1:N-1) = ALPHA,
%    LAMBDA(N) = ALPHA + N * BETA.
%
%    The eigenvector associated with LAMBDA(N) is (1,1,1,...,1)/sqrt(N).
%
%    The other N-1 eigenvectors are simply any (orthonormal) basis
%    for the space perpendicular to (1,1,1,...,1).
%
%    A is nonsingular if ALPHA /= 0D+00 and ALPHA + N * BETA /= 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Robert Gregory, David Karney,
%    Example 3.25,
%    A Collection of Matrices for Testing Computational Algorithms,
%    Wiley, New York, 1969, page 53,
%    LC: QA263.G68.
%
%    Donald Knuth,
%    The Art of Computer Programming,
%    Volume 1, Fundamental Algorithms, Second Edition,
%    Addison-Wesley, Reading, Massachusetts, 1973, page 36.
%
%  Parameters:
%
%    Input, real ALPHA, BETA, scalars that define A.
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a(1:n,1:n) = beta;

  for i = 1 :n
    a(i,i) = a(i:i) + alpha;
  end

  return
end


text:
%% Dataset Basics
% Load datasets using cosmo_fmri_dataset
%
% This function loads data stored in a nifti file and return a dataset struct
% where the data are stored in the 2-D array in field dataset.samples
%
% For each of the three masks ('brain','ev','vt'),
% print the number of voxels when loading the dataset with that mask.
%
% Hint: the number of voxels is given by the number of columns in
% dataset.samples
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

%% Set the datapath
%
config=cosmo_config();
data_path=fullfile(config.tutorial_data_path,'ak6','s01');

%% Compute number of voxels in each mask
% Hints:
% -In order not to get overwhelmed, solve a simple problem first and then
% generalize from it.
%   -First write some code in which you load data with full brain mask.
% Generalization:
%   -Next, you could just copy and alter that code such that you do the
%   same thing with the other masks
%   -A more elegant solution would be to put the three mask names into a
%   cell array and to write a loop that performs the same set of operations
%   on the members of the cell array (i.e. the three different mask names

% Let's start with the simple approach.
% Set the filename
% >@@>
mask_fn = fullfile(data_path, 'brain_mask.nii');
% <@@<

% Load the dataset and store in struct 'ds'
% >@@>
ds=cosmo_fmri_dataset(mask_fn);
% <@@<

% Compute number of features that are greater than zero
% hint: use ds.samples
nfeatures=sum(ds.samples>0);

fprintf('There are %d voxels in the whole brain mask\n', nfeatures);
% <@@<

% Now do the same with the EV and VT masks.
% >@@>
mask_fn = fullfile(data_path, 'vt_mask.nii');
ds=cosmo_fmri_dataset(mask_fn);
nfeatures=sum(ds.samples>0);

fprintf('There are %d voxels in the ventral-temporal mask\n', nfeatures);

mask_fn = fullfile(data_path, 'ev_mask.nii');
ds=cosmo_fmri_dataset(mask_fn);
nfeatures=sum(ds.samples>0);
fprintf('There are %d voxels in the early-visual brain mask\n', nfeatures);

% <@@<
%
% And here is space for the more elegant solution in which you define a
% list of mask names and apply the operation in a loop for all masks in the
% list
maskNames = {'brain_mask.nii', 'vt_mask.nii', 'ev_mask.nii'};
data_fn = fullfile(data_path, 'glm_T_stats_perrun.nii');

for iMask=1:numel(maskNames)
    % >@@>
    mask_fn=fullfile(data_path, maskNames{iMask});
    ds=cosmo_fmri_dataset(mask_fn);
    nfeatures=sum(ds.samples>0);
    fprintf('There are %6d voxels in the mask ''%s''\n', nfeatures, maskNames{iMask});
    % <@@<
end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function prev_rng = rcnn_seed_rand()
% seed_rand - Set random number generator to a fixed seed.
%   prev_rng = rcnn_seed_rand(seed)
%
%   Strategic use ensures that results are reproducible.
%
%   To restore the previous rng after calling this do:
%   rng(prev_rng);

% AUTORIGHTS
% ---------------------------------------------------------
% Copyright (c) 2014, Ross Girshick
%
% This file is part of the R-CNN code and is available
% under the terms of the Simplified BSD License provided in
% LICENSE. Please retain this notice and LICENSE if you use
% this file (or any portion of it) in your project.
% ---------------------------------------------------------

% This value works best for me.
seed = 3;
% Just kidding, of course ;-).

prev_rng = rng;
rng(seed, 'twister')

author: s-gupta

repo: rcnn-depth

save_path: github-repos/MATLAB/s-gupta-rcnn-depth


text:
function S=smoothHeaviside(varargin)

switch nargin
    case 1
        x=varargin{1};
        k=6;
        r=0;
    case 2
        x=varargin{1};
        k=varargin{2};
        r=0;
    case 3
        x=varargin{1};
        k=varargin{2};
        r=varargin{3};
end
k=k*2;
S=exp(2*k*(x-r))./(1+exp(2*k*(x-r)));

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [K, runtime, Phi] = spkernel(Graphs, features)
% Compute shortest path kernel for a set of node-labeled graphs
% Copyright 2012 Nino Shervashidze
% Input: Graphs - a 1xN array of graphs
% Output: K - NxN kernel matrix K
%         runtime - scalar
%         features - a boolean: 1 if we want to output the feature
%                    vector representation for each graph, 0 otherwise

N=size(Graphs,2);
Ds = cell(1,N); % shortest path distance matrices for each graph
if nargin<3 features=0; end

t=cputime; % for measuring runtime

%%% PREPROCESSING (mainly in order to find out the size of the node
%%% label alphabet, L, and rename labels as 1 ,..., L)
label_lookup=containers.Map();
label_counter=1;
for i=1:N
  for j=1:length(Graphs(i).nl.values)
    str_label=num2str(Graphs(i).nl.values(j));
    % str_label is the node label of the current node of the
    % current graph converted into a string
    if ~isKey(label_lookup, str_label)
      label_lookup(str_label)=label_counter;
      Graphs(i).nl.values(j)=label_counter;
      label_counter=label_counter+1;
    else
      Graphs(i).nl.values(j)=label_lookup(str_label);
    end
  end
end
L=label_counter-1; % L is the size of the node label alphabet
disp(['the preprocessing step took ', num2str(cputime-t), ' sec']);
t=cputime;

% compute Ds and the length of the maximal shortest path over the dataset
maxpath=0;
for i=1:N
  Ds{i}=floydwarshall(Graphs(i).am);
  aux=max(Ds{i}(~isinf(Ds{i})));
  if aux > maxpath
    maxpath=aux;
  end
 % if rem(i,100)==0 disp(i); end
end

sp=sparse((maxpath+1)*L*(L+1)/2,N);
for i=1:N
  labels_aux=repmat(Graphs(i).nl.values,1,length(Graphs(i).nl.values));
  a=min(labels_aux, labels_aux');
  b=max(labels_aux, labels_aux');
  I=triu(~(isinf(Ds{i})));
  Ind=Ds{i}(I)*L*(L+1)/2+(a(I)-1).*(2*L+2-a(I))/2+b(I)-a(I)+1;
  aux=accumarray(Ind,ones(nnz(I),1));
  sp(Ind,i)=aux(Ind);
end
sp=sp(sum(sp,2)~=0,:);
K=full(sp'*sp);
if features Phi=sp; end
runtime=cputime-t;
disp(['kernel computation took ', num2str(cputime-t), ' sec']);
end

author: muhanzhang

repo: DGCNN

save_path: github-repos/MATLAB/muhanzhang-DGCNN


text:
filename='Cantilever_triangle_fine';
%filename='Cantilever_quad_fine';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'compliance','perimeterInterior'};
%cost = {'compliance','perimeter'};
%cost = {'compliance'};
weights = [1, 0.2];
%weights = 1;
constraint = {'volumeConstraint'};
optimizerUnconstrained = 'SLERP';
incrementFactor = 1;
designVariable = 'LevelSet';
filterType = 'P1';
optimizer = 'DualNestedInPrimal';
%optimizer = 'AlternatingPrimalDual';

nsteps = 5;
Vfrac_final = 0.4;
Perimeter_target = 1;
optimality_final = 1e-4;
constr_final =1e-4;

Vfrac_initial = 0.4;
optimality_initial = 1e-4;
constr_initial = 1e-4;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;

plotting = true;
printing = true;
monitoring = true;
monitoring_interval = 15;

maxiter = 300;