text:
function theta = polygon3dNormalAngle(points, ind)
%POLYGON3DNORMALANGLE Normal angle at a vertex of the 3D polygon
%
%   THETA = polygon3DNormalAngle(POLYGON, IND)
%   where POLYGON is a set of points, and IND is index of a point in
%   polygon. The function compute the angle of the normal cone localized at
%   this vertex.
%   If IND is a vector of indices, normal angle is computed for each vertex
%   specified by IND.
%
%   Example
%   % create an equilateral triangle in space
%   poly3d = [1 1 0;-1 0 1;0 -1 -1];
%   % compute each normal angle
%   theta = polygon3dNormalAngle(poly3d, 1:size(poly3d, 1));
%   % sum of normal angles must be equal to 2*PI for simple polygons
%   sum(theta)
%
%   IMPORTANT NOTE: works only for convex angles ! ! ! !
%
%   See also
%   polygons3d, faceNormalAngle
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2005-11-30
% Copyright 2005 INRA - CEPIA Nantes - MIAJ (Jouy-en-Josas).

% number of points
np = size(points, 1);

% number of angles to compute
nv = length(ind);

theta = zeros(nv, 1);

for i=1:nv
    p0 = points(ind(i), :);

    if ind(i)==1
        p1 = points(np, :);
    else
        p1 = points(ind(i)-1, :);
    end

    if ind(i)==np
        p2 = points(1, :);
    else
        p2 = points(ind(i)+1, :);
    end

    theta(i) = pi - anglePoints3d(p1, p0, p2);
end

author: rpng


text:
function val = file2mat(a,varargin)
% Function for reading from file_array objects.
% FORMAT val = file2mat(a,ind1,ind2,ind3,...)
% a      - file_array object
% indx   - indices for dimension x (int32)
% val    - the read values
%
% This function is normally called by file_array/subsref
% _______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: file2mat.m 1530 2008-04-30 19:28:22Z guillaume $

%-This is merely the help file for the compiled routine
error('file2mat.c not compiled - see Makefile');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [parameters, ll, ht, VCV, scores, diagnostics] = matrix_garch(data,dataAsym,p,o,q,startingvals,options)
% Estimation of symmetric and asymmetric MATRIX multivariate GARCH models
%
% USAGE:
%  [PARAMETERS,LL,HT,VCV,SCORES,DIAGNOSTICS] = matrix_garch(DATA,DATAASYM,P,O,Q,STARTINGVALS,OPTIONS)
%
% INPUTS:
%   DATA         - A T by K matrix of zero mean residuals -OR-
%                    K by K by T array of covariance estimators (e.g. realized covariance)
%   DATAASYM     - [OPTIONAL] K by K by T array of asymmetric covariance
%                    estimators (e.g. RC scaled by indicator functions)
%   P            - Positive, scalar integer representing the number of lags of the innovation process
%   O            - Non-negative scalar integer representing the number of asymmetric lags to include
%   Q            - Non-negative scalar integer representing the number of lags of conditional covariance
%   STARTINGVALS - [OPTIONAL] INCOMPLETE
%   OPTIONS      - [OPTIONAL] Options to use in the optimization (fminunc)
%
% OUTPUTS:
%   PARAMETERS   - K(K+1)/2*(1+P+O+Q) by 1 vector of parameters of the form
%                    [vech(C)' vech(A(1))' ... vech(A(P))' vech(G(1))' ...
%                    vech(G(O))' vech(B(1))' ... vech(B(Q))']'
%   LL           - The log likelihood at the optimum
%   HT           - A [K K T] dimension matrix of conditional covariances
%   VCV          - A numParams^2 square matrix of robust parameter covariances (A^(-1)*B*A^(-1)*t^(-1))
%   SCORES       - A T by numParams matrix of individual scores
%   DIAGNOSTICS  - Structure containing some diagnostic information
%
% COMMENTS:
%    The conditional variance, H(t), of a MATRIX GARCH is modeled as follows:
%
%      H(t) = CC' + A(1)A(1)'.*r_{t-1}'*r_{t-1} + ... + A(P)A(P)'.*r_{t-P}'*r_{t-P}
%                 + G(1)G(1)'.*n_{t-1}'*n_{t-1} + ... + G(O)G(O)'.*n_{t-P}'*n_{t-P}
%                  + B(1)B(1)'.*H(t-1) +...+ B(Q)B(Q)'.*H(t-q)
%
%    where n_{t} = r_{t} .* (r_{t}<0).  If using realized measures, the
%    RM_{t-1} replaces r_{t-1}'*r_{t-1}, and the asymmetric version
%    replaces n_{t-1}'*n_{t-1}
%
% EXAMPLES:
%     % Estimation of a symmetric MATRIX GARCH(1,0,1) model
%     parameters = matrix_garch(data,[],1,0,1);
%     % Estimation of an asymmetric MATRIX GARCH(1,0,1) model
%     parameters = matrix_garch(data,[],1,1,1);
%     % Estimation of a symmetric MATRIX GARCH(1,0,1) model using realized covariacne
%     data = RC % K by K by T 3D of realized covariacne
%     parameters = matrix_garch(RC,[],1,0,1);

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 3    Date: 10/28/2009

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Argument Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch nargin
    case 3
        o=0;
        q=0;
        startingvals=[];
        options = [];
    case 4
        q=0;
        startingvals=[];
        options = [];
    case 5
        startingvals=[];
        options = [];
    case 6
        options = [];
    case 7
    otherwise
        error('Between 2 and 6 arguments required.')
end

%data should be TxK, T>K
if ndims(data)==2
    [t,k]=size(data);
    if ~isempty(dataAsym)
        error('If DATA is a T by K matrix, DATAASYM must be empty.');
    end
    temp = zeros(k,k,t);
    dataAsym = zeros(k,k,t);
    for i=1:t
        temp(:,:,i) = data(i,:)'*data(i,:);
        dataAsym(:,:,i) = (data(i,:).*(data(i,:)<0))'*(data(i,:).*(data(i,:)<0));
    end
    data = temp;
elseif ndims(data)==3
    [k,m,t] = size(data);
    if m~=k
        error('DATA must be K by K by T is a 3D array.');
    end
    if ~isempty(dataAsym)
        if ndims(dataAsym)~=3
            error('DATAASYM must be a 3D array with the same dimensions as DATA');
        end
        [k2,m2,t2]=size(dataAsym);
        if any([k m t]~=[k2 m2 t2])
            error('DATAASYM must be a 3D array with the same dimensions as DATA');
        end
    end
end
k2 = k*(k+1)/2;
if min(t,k)<2 || t<k
    error('DATA must be a T by K matrix or a K by K by T 3D array, T>K>1');
end

%p, o, q much be non-negative scalars
if length(p)>1 || any(p<1) || floor(p)~=p
    error('P must be a positive scalar');
end
if isempty(o)
    o = 0;
end
if length(o)>1 || any(o<0) || floor(o)~=o
    error('O must be a non-negative scalar');
end
if o>0 && isempty(dataAsym)
    error('DATAASYM must be non-empty if O>0.')
end

if isempty(q)
    q = 0;
end
if length(q)>1 || any(q<0) || floor(q)~=q
    error('Q must be a non-negative scalar');
end

% Startingvals must have (k*(k+1)/2)(1+p+o+q) parameters
if  ~isempty(startingvals)
    if size(startingvals,2)>size(startingvals,1)
        startingvals = startingvals';
    end
    if length(startingvals)<(p+o+q)
        error('STARTINGVALS should be a P+O+Q by 1 vector');
    end
    % Only validate if provided
    % FIXME: Fix Validation
    % TODO: Fix Validation
    kappa = 2;
    A=startingvals(1:p);
    G=startingvals(p+1:p+o);
    B=startingvals(p+o+1:p+o+q);
    if (sum(A)+sum(G)/kappa+sum(B))>=.999998
        error('Weighted sum of STATINGVALUES must be less than 1. See Comments.');
    end
    if any(A<0) || any(B<0) || any(G<0)
        error('STARTINGVALS must all be nonnegative.');
    end
end

%Make sure options is a valid option structure
if isempty(options)
    options=optimset('fminunc');
    options.Display='iter';
    options.Diagnostics='on';
    options.LargeScale='off';
    options.MaxFunEvals = 1000*k2*(1+p+o+q);
end
try
    optimset(options);
catch ME
    error('OPTIONS is not a valid options structure');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Argument Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Compute the backCast
backCast = zeros(k);
backCastAsym = zeros(k);
tau = max(ceil(sqrt(t)),k);
weights = .06 * .94.^(0:tau);
weights = weights / sum(weights);
for i=1:tau
    backCast = backCast + weights(i) * data(:,:,i);
    if o>0
        backCastAsym = backCastAsym + weights(i) * dataAsym(:,:,i);
    end
end

%Get starting values from scalar_vt_vech
if isempty(startingvals)
    startingOptions=optimset('fminunc');
    startingOptions.Display='off';
    startingOptions.Diagnostics='off';
    startingOptions.LargeScale='off';
    startingOptions.TolX = 1e-4;
    startingOptions.TolFun = 1e-4;
    [scalarVechStartingvals,~,~,~,~,~,diagnostics] = scalar_vt_vech(data,dataAsym,p,o,q,[],[],startingOptions);
    CpC = diagnostics.intercept;

    startingvals = zeros(k2 * (1+p+o+q),1);
    startingvals(1:k2) = chol2vec(chol(CpC)');
    index = k2;
    for i=1:(p+o+q)
        matrixParameters = scalarVechStartingvals(i)*(.02*eye(k) + .98*ones(k));
        startingvals(index+1:index+k2) = chol2vec(chol(matrixParameters)');
        index = index + k2;
    end
end

warning('off','MATLAB:illConditionedMatrix')
ll0 = matrix_garch_likelihood(startingvals,data,dataAsym,p,o,q,backCast,backCastAsym);
[parameters,ll,exitflag,output]=fminunc('matrix_garch_likelihood',startingvals,options,data,dataAsym,p,o,q,backCast,backCastAsym);

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Estimation Robustification
%%%%%%%%%%%%%%%%%%%%%%%%%%%
if exitflag<=0 &&  ll<ll0
    %Did not converge, but function improved
    options.MaxFunEvals=4*100*(p+q);
    options.MaxIter=2*100*(p+q);
    parameters=fminunc('matrix_garch_likelihood',parameters,options,data,dataAsym,p,o,q,backCast,backCastAsym);
end
warning('on','MATLAB:illConditionedMatrix')
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Estimation Robustification
%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargout>1
    [ll,lls,ht]=matrix_garch_likelihood(parameters,data,dataAsym,p,o,q,backCast,backCastAsym);
    ll=-ll;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the VCV
%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargout>3
    [VCV,A,B,scores]=robustvcv('matrix_garch_likelihood',parameters,0,data,dataAsym,p,o,q,backCast,backCastAsym);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the VCV
%%%%%%%%%%%%%%%%%%%%%%%%%%%
diagnostics = [];
diagnostics.EXITFLAG=exitflag;
diagnostics.ITERATIONS=output.iterations;
diagnostics.FUNCCOUNT=output.funcCount;
diagnostics.MESSAGE=output.message;
parameterMatrices = zeros(k,k,1+p+o+q);
index = 0;
for i=1:(1+p+o+1)
    temp = vec2chol(parameters(index+1:index+k2));
    parameterMatrices(:,:,i) = temp*temp';
    index=index+k2;
end
diagnostics.C = parameterMatrices(:,:,1);
diagnostics.A = parameterMatrices(:,:,2:p+1);
if o>0
    diagnostics.G = parameterMatrices(:,:,p+2:p+o+1);
else
    diagnostics.G = [];
end

if q>0
    diagnostics.B = parameterMatrices(:,:,p+o+2:p+o+q+1);
else
    diagnostics.B = [];
end

author: bashtage

repo: mfe-toolbox


text:
function dist = sphere01_distance_xyz ( xyz1, xyz2 )

%*****************************************************************************80
%
%% SPHERE01_DISTANCE_XYZ computes great circle distances on a sphere.
%
%  Discussion:
%
%    XYZ coordinates are used.
%
%    We assume the points XYZ1 and XYZ2 lie on the same sphere.
%
%    This computation is a special form of the Vincenty formula.
%    It should be less sensitive to errors associated with very small
%    or very large angular separations.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    "Great-circle distance",
%    Wikipedia.
%
%  Parameters:
%
%    Input, real XYZ1(3), the coordinates of the first point.
%
%    Input, real XYZ2(3), the coordinates of the second point.
%
%    Output, real DIST, the great circle distance between the points.
%
  lat1 = r8_asin ( xyz1(3) );
  lon1 = r8_atan ( xyz1(2), xyz1(1) );

  lat2 = r8_asin ( xyz2(3) );
  lon2 = r8_atan ( xyz2(2), xyz2(1) );

  top = ( cos ( lat2 ) * sin ( lon1 - lon2 ) ).^2 ...
      + ( cos ( lat1 ) * sin ( lat2 ) ...
      -   sin ( lat1 ) * cos ( lat2 ) * cos ( lon1 - lon2 ) ).^2;

  top = sqrt ( top );

  bot = sin ( lat1 ) * sin ( lat2 ) ...
      + cos ( lat1 ) * cos ( lat2 ) * cos ( lon1 - lon2 );

  dist = atan2 ( top, bot );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function GaussianHighpass
a=imread('cameraman.tif');
figure(1)
imshow(a)
[m n]=size(a);
f_transform=fft2(a);
f_shift=fftshift(f_transform);
p=m/2;
q=n/2;
d0=70;
for i=1:m
for j=1:n
distance=sqrt((i-p)^2+(j-q)^2);
low_filter(i,j)=1-exp(-(distance)^2/(2*(d0^2)));
end
end
filter_apply=f_shift.*low_filter;
image_orignal=ifftshift(filter_apply);
image_filter_apply=abs(ifft2(image_orignal));
figure(2)
imshow(image_filter_apply,[])

author: Sable

repo: mcbench-benchmarks


text:
%+ easy11 bbbbbbbbbb
% easy12 chistart
% easy14 bound2
%+ easy17 bbbbbbbbbb
%
%
% bbbbbbbbb rinexa: bbbbbbbbbbbbbbbbbbb in the end

author: yandld


text:
function c = sinh(a)
% SINH for adiff objects.

c = adiff( sinh(a.x), rowmult(cosh(a.x), a.dx), a.root);

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   demo script for mesh generation from binarized volumetric image
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% preparation
% user must add the path of iso2mesh to matlab path list
% addpath('../');

% user need to add the full path to .../iso2mesh/bin directory
% to windows/Linux/Unix PATH environment variable

%% load the sample data
load rat_head.mat

% volimage is a volumetric image such as an X-ray or MRI image
% A,b are registration matrix and vector, respectively
%% perform mesh generation

%% use the alternative 'simplify' method: first create voxel-based
% surface mesh, and then resample it to desired density.
% this method does not guarantee to be free of self-intersecting
% element, as 'cgalsurf' promises.

[node,elem,face]=vol2mesh(volimage>0.05,1:size(volimage,1),1:size(volimage,2),...
                          1:size(volimage,3),0.1,2,1,'simplify');

%% visualize the resulting mesh

plotmesh(node,face);
axis equal;


text:
function [ h, hmin, hmax ] = p18_stepsize ( option )

%*****************************************************************************80
%
%% P18_STEPSIZE returns step sizes for problem 18.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer OPTION, the option index.
%
%    Output, real H, HMIN, HMAX, suggested values for the
%    initial step, the minimum step, and the maximum step.
%
  h =    2.500;
  hmin = 0.001;
  hmax = 5.000;

  return
end


text:
function planC = copyStrToScan_noMesh(structNum,scanNum,planC)
%function planC = copyStrToScan_noMesh(structNum,scanNum,planC)
%
%This function derives a new structure from structNum which is associated
%to scanNum. The naming convention for this new structure is
%[structName assoc scanNum]. If structNum is already associated to scanNum,
%the unchanged planC is returned back.
%
%APA, 01/25/08
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

global stateS

if ~exist('planC','var')
    global planC
end

indexS = planC{end};

%return if structNum is already associated to scanNum
assocScanNum = getAssociatedScan(planC{indexS.structures}(structNum).assocScanUID,planC);
if assocScanNum == scanNum
    warning(['Structure Number ',num2str(structNum),' is already assocoated with scan ',num2str(scanNum)])
    return;
end

%Get transformation matrix
if ~isfield(planC{indexS.scan}(scanNum),'transM') || isempty(planC{indexS.scan}(scanNum).transM)
    transMnew = eye(4);
else
    transMnew = planC{indexS.scan}(scanNum).transM;
end
if ~isfield(planC{indexS.scan}(assocScanNum),'transM') || isempty(planC{indexS.scan}(assocScanNum).transM)
    transMold = eye(4);
else
    transMold = planC{indexS.scan}(assocScanNum).transM;
end
transM = inv(transMnew)*transMold;

%Get z-coordinates of scan
[jnk1,jnk2,scanZv] = getScanXYZVals(planC{indexS.scan}(scanNum));
[jnk, relStructNum] = getStructureAssociatedScan(structNum, planC);
for i = 1:length(scanZv)
    coord = scanZv(i);

    %structUID   = planC{indexS.structures}(structNum).strUID;
    %contourS    = calllib('libMeshContour','getContours',structUID,single(pointOnPlane),single(planeNormal),single([0 1 0]),single([1 0 0]));

    [slcC, sliceXVals, sliceYVals] = getStructureSlice(assocScanNum, 3, coord, transM, planC);
    oneStructM = [];
    if relStructNum <= 52
        cellNum = 1;
        if iscell(slcC) && ~isempty(slcC{cellNum})
            oneStructM = bitget(slcC{cellNum}, relStructNum);
        end
    else
        cellNum = ceil((relStructNum-52)/8)+1; %uint8
        if iscell(slcC) && ~isempty(slcC{cellNum})
            oneStructM = bitget(slcC{cellNum}, relStructNum-52-(cellNum-2)*8); %uint8
        end
    end

    if ~isempty(oneStructM)
        C           = contourc(sliceXVals, sliceYVals, double(oneStructM'),[0.5 0.5]);
        indC = getSegIndices(C);
        if ~isempty(indC)
            for seg = 1:length(indC)
                points = [C(:,indC{seg})' coord*ones(length(C(1,indC{seg})),1)];
                contourS(i).segments(seg).points = points;
            end
        else
            contourS(i).segments.points = [];
        end
    else
        contourS(i).segments.points = [];
    end

end

strname = [planC{indexS.structures}(structNum).structureName,' asoc ',num2str(scanNum)];

newstr = newCERRStructure(scanNum,planC);
newstr.contour = contourS;
newstr.structureName = strname;
newstr.associatedScan = scanNum;
newstr.assocScanUID = planC{indexS.scan}(scanNum).scanUID;
newstr.meshRep = 0;
numStructs = length(planC{indexS.structures});

%Append new structure to planC.
planC{indexS.structures} = dissimilarInsert(planC{indexS.structures}, newstr, numStructs+1, []);

%Create Raster Segments
planC = getRasterSegs(planC, numStructs+1);

%Update uniformized data.
planC = updateStructureMatrices(planC, numStructs+1);

if isfield(stateS,'handle') && ishandle(stateS.handle.CERRSliceViewer)
    stateS.structsChanged = 1;
    CERRRefresh
end

return;

function indC = getSegIndices(C)
% function getSegIndices(C)
%
%This function returns the indices for each segment of input contour C.
%C is output from in-built "contourc" function
%
%APA, 12/15/2006

start = 1;
counter = 1;
indC = [];
while start < length(C(2,:))
    numPts = C(2,start);
    indC{counter} = [(start+1):(start+numPts) start+1];
    start = start + numPts + 1;
    counter = counter + 1;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function result = simplex_unit_05_nd ( func, n )

%*****************************************************************************80
%
%% SIMPLEX_UNIT_05_ND approximates an integral inside a unit simplex in ND.
%
%  Integration region:
%
%    The unit simplex in N dimensions,
%      0 <= X(1:N),
%      Sum ( X(1:N) ) <= 1.
%
%  Discussion:
%
%    An N^2 + 3 N + 3 point formula of degree 5 is used.  This is
%    Stroud formula TN:5-1.
%
%    (For N = 2, the number of points is actually only 7, and
%     for N = 3, the number of points is actually only 15.)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 November 2004
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Arthur H Stroud,
%    A Fifth Degree Integration Formula for the N-Simplex,
%    SIAM Journal on Numerical Analysis,
%    Volume 6, Number 1, March 1969.
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied
%    function which evaluates F(X) at the N-dimensional point
%    X, of the form
%      function value = func ( n, x )
%
%    Input, integer N, the dimension of the space.  For this routine,
%    it must be the case that 2 <= N <= 16.
%
%    Output, real RESULT, the approximate integral of the function.
%
  coef1 = [ ...
    0.0E+00,             0.225E+00, ...
    0.118518518519E+00,  0.0631521898883E+00, ...
    0.235714285714E+00,  0.791575476992E+00, ...
    1.85798728021E+00,   3.53666958042E+00, ...
    5.90844340844E+00,   9.03765432098E+00, ...
    12.9758241758E+00,  17.7645108738E+00, ...
    23.4375030259E+00,  30.0224941950E+00, ...
    37.5423613501E+00,  46.0161454949E+00 ];
  coef21 = [ ...
    0.0E+00,              0.12593918054483E+00, ...
    0.0719370837790E+00,  0.0470456145702E+00, ...
    0.0333009774677E+00,  0.0248633014592E+00, ...
    0.0192679696358E+00,  0.0153322153879E+00, ...
    0.0124316229901E+00,  0.0102112988361E+00, ...
    0.00845730697460E+00, 0.00703433430999E+00, ...
    0.00585330520067E+00, 0.00485356735291E+00, ...
    0.00399261092720E+00, 0.00323988713017E+00 ];
  coef22 = [ ...
    0.0E+00,              0.13239415278851E+00, ...
    0.0690682072263E+00,  0.0371530185868E+00, ...
   -0.0719253160920E+00, -0.264323879461E+00, ...
   -0.537926779961E+00,  -0.886895605701E+00, ...
   -1.30409181465E+00,   -1.78227048964E+00, ...
   -2.31462336314E+00,   -2.89499045158E+00, ...
   -3.51790849765E+00,   -4.17858310668E+00, ...
   -4.87282884913E+00,   -5.59699944261E+00 ];
  coef31 = [ ...
    0.0E+00,             0.0E+00, ...
    0.0529100529100E+00, 0.0261368740713E+00, ...
    0.0499020181331E+00, 0.0782233395867E+00, ...
    0.109041040862E+00,  0.140874828568E+00,  ...
    0.172735353396E+00,  0.203992490408E+00,  ...
    0.234263814181E+00,  0.263332763315E+00,  ...
    0.291091849264E+00,  0.317504208212E+00,  ...
    0.342577872069E+00,  0.366348654344E+00 ];
  coef32 = [ ...
    0.0E+00,              0.0E+00, ...
    0.0E+00,              0.0254485903613E+00, ...
    0.0165000982690E+00,  0.0115218303668E+00,  ...
    0.00850478779483E+00, 0.00655297510968E+00, ...
    0.00522372456259E+00, 0.00428017828134E+00, ...
    0.00358722367033E+00, 0.00306362964360E+00, ...
    0.00265836687133E+00, 0.00233816221525E+00, ...
    0.00208061510846E+00, 0.00187022027571E+00 ];

  if ( n < 2 | 16 < n )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SIMPLEX_UNIT_05_ND - Fatal error!\n' );
    fprintf ( 1, '  Input spatial dimension N out of range.\n' );
    fprintf ( 1, '  N = %d\n', n );
    error ( 'SIMPLEX_UNIT_05_ND - Fatal error!' );
  end

  quad = 0.0;
%
%  S1
%
  x(1:n) = 1.0E+00 / ( n + 1 );
  quad = quad + coef1(n) * feval ( func, n, x );
%
%  S21
%
  r1 = ( ( n + 4 ) - sqrt ( 15.0E+00 ) ) / ( n * n + 8 * n + 1 );
  s1 = 1.0E+00 - n * r1;

  x(1:n) = r1;

  for i = 1 : n + 1

    quad = quad + coef21(n) * feval ( func, n, x );

    if ( 1 < i )
      x(i-1) = r1;
    end

    if ( i < n + 1 )
      x(i) = s1;
    end

  end
%
%  S22
%
  r2 = ( n + 4 + sqrt ( 15.0E+00 ) ) / ( n * n + 8 * n + 1 );
  s2 = 1.0E+00 - n * r2;

  x(1:n) = r2;

  for i = 1 : n + 1

    quad = quad + coef22(n) * feval ( func, n, x );

    if ( 1 < i )
      x(i-1) = r2;
    end

    if ( i < n + 1 )
      x(i) = s2;
    end

  end
%
%  S31
%
  u1 = ( n + 7 + 2.0E+00 * sqrt ( 15.0E+00 ) ) / ( n * n + 14 * n - 11 );
  v1 = ( ( 4 * n - 2 ) - ( n - 1 ) * sqrt ( 15.0E+00 ) ) / ( n * n + 14 * n - 11 );

  for i = 1 : n

    x(1:n) = u1;
    x(i) = v1;

    for j = i : n

      if ( i < j - 1 )
        x(j-1) = u1;
      end

      x(j) = v1;

      quad = quad + coef31(n) * feval ( func, n, x );

    end

  end
%
%  S32
%
  u2 = ( n + 7 - 2.0E+00 * sqrt ( 15.0E+00 ) ) / ( n^2 + 14 * n - 11 );
  v2 = ( ( 4 * n - 2 ) + ( n - 1 ) * sqrt ( 15.0E+00 ) ) / ( n^2 + 14 * n - 11 );

  for i = 1 : n

    x(1:n) = u2;
    x(i) = v2;

    for j = i : n

      if ( i < j - 1 )
        x(j-1) = u2;
      end

      x(j) = v2;

      quad = quad + coef32(n) * feval ( func, n, x );

    end

  end

  volume = simplex_unit_volume_nd ( n );

  result = quad * volume;

  return
end


function [prices, stdErrs] = Price_MC_Barrier_Strikes_func(Spath, call, down, H, Kvec, M, mult, rebate, r, T)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% About: Calculates Knock-out Barrier option prices for vector of strikes, given the simulatd paths
%          This version Allows for a rebate.
%          Note: to price knock-in options, use parity
% Returns: prices and standard errors for each of the supplied strikes
% Author: Justin Lars Kirkby
%
% -----------------
% Params
% -----------------
% Spath = paths of underlying, dimension N_sim x M+1, where M = number of time steps (since includes S_0)
% call = 1 for call (else put)
% down = 1 for down and out, else up and out
% H = Barrier
% Kvec = strike vector
% M = number of monitoring points, e.g. 252 for "daily" monitoring
% mult = time partitioning multiplier to reduce bias (e.g. mult = 2 or 5)
% rebate = rebate which is paid upon barrier breach, e.g. 5.0
% r = Interest rate
% q = dividend yield
% T = Time (in years)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

prices = zeros(length(Kvec),1);
stdErrs = zeros(length(Kvec),1);
N_sim = size(Spath,1);  % number of paths

knock_time = zeros(N_sim,1);  % time of knock-out

M_mult = M*mult;  %time partitioning to reduce bias
dt_mult = T/M_mult;

if down == 1  % down and out
    for n = 1:N_sim
        for m=1:mult:M_mult+1
            if Spath(n,m) < H
                knock_time(n) = (m-1)*dt_mult;
                break
            end
        end
    end
else %up and out
    for n = 1:N_sim
        for m=1:mult:M_mult+1
            if Spath(n,m) > H
                knock_time(n) = (m-1)*dt_mult;
                break
            end
        end
    end
end

if rebate > 0
   disc_rebate = rebate*exp(-r*knock_time).*(knock_time>0);   % discounted rebate
else
   disc_rebate = 0;
end

for k = 1:length(Kvec)
    K = Kvec(k);
    if call ==1
    	payoffs  = exp(-r*T)*max(0, Spath(:,M_mult+1) - K).*(knock_time==0) + disc_rebate;

    else
        payoffs  = exp(-r*T)*max(0, K - Spath(:,M_mult+1)).*(knock_time==0) + disc_rebate;
    end

    prices(k) = mean(payoffs);
    stdErrs(k) = std(payoffs) / sqrt(N_sim);
end

end


text:
%% processing for concentrations
if contains(lower(param.method),'conc')
    if ~isfield(param,'maxConc')
        param.maxConc=1e4;
    end
    if ~isfield(param,'externalNetFluxBounds')
        if isfield(model,'dxl') || isfield(model,'dxu')
            param.externalNetFluxBounds='dxReplacement';
        else
            param.externalNetFluxBounds='original';
        end
    end

    nMetabolitesPerRxn = sum(model.S~=0,1)';
    bool = nMetabolitesPerRxn~=1 & ~model.SConsistentRxnBool;
    if any(bool)
        warning('Exchange reactions involving more than one metabolite, check bounds on x - x0')
        disp(model.rxns(bool))
    end

    if any(~model.SConsistentRxnBool)

        switch param.externalNetFluxBounds
            case 'original'
                if param.printLevel>0
                    fprintf('%s\n','Using existing external net flux bounds without modification.')
                end
                if (isfield(model,'dxl') && any(model.dxl~=0)) || (isfield(model,'dxu') && any(model.dxu~=0))
                    error('Option clash between param.externalNetFluxBounds=''original'' and (isfield(model,''dxl'') && any(model.dxl~=0)) || (isfield(model,''dxu'') && any(model.dxu~=0))')
                end
                %
                vel = lb(~model.SConsistentRxnBool);
                veu = ub(~model.SConsistentRxnBool);
                %force initial and final concentration to be equal
                model.dxl = zeros(m,1);
                model.dxu = zeros(m,1);
            case 'dxReplacement'
                %TODO
                error('revise how net to initial and final conc bounds are dealt with')
                if ~isfield(model,'dxl')
                    %close bounds by default
                    model.dxl = zeros(m,1);
                    if exist('B','var')
                        dxlB =  -B*model.lb(~model.SConsistentRxnBool);
                        model.dxl(dxlB~=0)=dxlB(dxlB~=0);
                    end
                end
                if ~isfield(model,'dxu')
                    %close bounds by default
                    model.dxu = zeros(m,1);
                    if exist('B','var')
                        dxuB =  -B*model.ub(~model.SConsistentRxnBool);
                        model.dxu(dxuB~=0)=dxuB(dxuB~=0);
                    end
                end
                %eliminate all exchange reactions
                B = B*0;
                vel = lb(~model.SConsistentRxnBool)*0;
                veu = ub(~model.SConsistentRxnBool)*0;
            otherwise
                error(['param.externalNetFluxBounds = ' param.externalNetFluxBounds ' is an unrecognised input'])
        end
    else
        model.dxl = -inf*ones(m,1);
        model.dxu =  inf*ones(m,1);
    end

    clear lb ub

    if ~isfield(model,'x0l')
        model.x0l = zeros(m,1);
    end
    if ~isfield(model,'x0u')
        model.x0u = param.maxConc*ones(m,1);
    end
    if ~isfield(model,'xl')
        model.xl = zeros(m,1);
    end
    if ~isfield(model,'xu')
        model.xu = param.maxConc*ones(m,1);
    end

    if ~isfield(model,'u0') || isempty(model.u0)
        model.u0='zero';
    end
    if ischar(model.u0)
        switch model.u0
            case 'rand'
                u0=rand(m,1);
            case 'one'
                u0=ones(m,1);
            case 'zero'
                u0=zeros(m,1);
            otherwise
                error('unrecognised option for model.u0')
        end
    else
        if length(model.u0)==size(model.S,1)
            u0 = columnVector(model.u0);
        else
            if length(model.u0)==1
                u0=ones(m,1)*model.u0;
            else
                error('model.u0 is of incorrect dimension')
            end
        end
        if any(~isfinite(u0))
            error('u0 must be finite')
        end
    end

    %assume concentrations are in uMol
    if ~isfield(model,'concUnit')
        concUnit = 10-6;
    end

    % Define constants
    if isfield(model,'gasConstant') && isfield(model,'T')
        if isfield(model,'gasConstant')
            gasConstant = model.gasConstant;
        else
            %8.31446261815324 J K^-1 mol^-1
            gasConstant = 8.3144621e-3; % Gas constant in kJ K^-1 mol^-1
        end
        if isfield(model,'T')
            temperature = model.T;
        else
            temperature = 310.15;
        end
        %dimensionless
        u0 = u0/(gasConstant*temperature);
    end

    if ~isfield(model,'f') || isempty(model.f)
        model.f='one';
    end
    if ischar(model.f)
        switch model.f
            case 'rand'
                f=N'*rand(m,1);
            case 'one'
                f=ones(m,1);
            case 'two'
                f=ones(m,1)*2;
        end
    else
        if length(model.f)==size(model.S,1)
            f = columnVector(model.f);
        else
            if length(model.f)==1
                f=ones(m,1)*model.f;
            end
        end
        if any(~isfinite(f))
            error('f must all be finite')
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function out = echo_enc_bf(signal, text, d0, d1, alpha, L)
%ECHO_ENC Echo Hiding with Backward and Forward Echo Kernels
%
%   INPUT VARIABLES
%       signal : Cover signal
%       text   : Message to hide
%       d0     : Delay rate for bit0
%       d1     : Delay rate for bit1
%       alpha  : Echo amplitude
%       L      : Length of frames
%
%   OUTPUT VARIABLES
%       out    : Stego signal
%
%   Kadir Tekeli (kadir.tekeli@outlook.com)

if nargin < 4
	d0 = 150;     %Delay rate for bit0
	d1 = 200;     %Delay rate for bit1
end

if nargin < 5
	alpha = 0.5;  %Echo amplitude
end

if nargin < 6
	L = 8*1024;   %Length of frames
end

[s.len, s.ch] = size(signal);
bit = getBits(text);
nframe = floor(s.len/L);
N = nframe - mod(nframe,8);      %Number of frames (for 8 bit)

if (length(bit) > N)
	warning('Message is too long, being cropped!');
	bits = bit(1:N);
else
	warning('Message is being zero padded...');
	bits = [bit, num2str(zeros(N-length(bit), 1))'];
end

[echo_zro, echo_one] = bf_echo(signal, d0, d1, alpha);  %Echo signals
mix = mixer(L, bits, 0, 1, 256) * ones(1, s.ch);        %Mixer signal

%%%%%%%%%%%%%%%%%%%%%%% EMBEDDING MESSAGE... %%%%%%%%%%%%%%%%%%%%%%%
out = signal(1:N*L, :) + echo_zro(1:N*L, :) .* abs(mix-1) ...
                       + echo_one(1:N*L, :) .* mix;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
out = [out; signal(N*L+1:s.len, :)];   %Rest of the signal
end

author: ktekeli

repo: audio-steganography-algorithms


text:
function divdif_test20 ( )

%*****************************************************************************80
%
%% DIVDIF_TEST20 tests DIF_DERIVK_TABLE;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 May 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'DIVDIF_TEST20\n' );
  fprintf ( 1, '  For a divided difference polynomial:\n' );
  fprintf ( 1, '  DIF_DERIVK_TABLE computes the K-th derivative;\n' );
%
%  Set the 0 data points.
%
  n0 = 5;
  x0 = linspace ( -2.0, +2.0, 5 );
%
%  Set data for x^4/24+x^3/3+x^2/2+x+1
%
  f0(1:n0) = 1.0;
  for i = 4 : -1 : 1
    f0(1:n0) = f0(1:n0) .* x0(1:n0) / i + 1.0;
  end
%
%  Compute the difference table.
%
  d0 = data_to_dif ( n0, x0, f0 );
  dif_print ( n0, x0, d0, '  The divided difference polynomial P0:' );

  c0 = dif_to_r8poly ( n0, x0, d0 );

  r8poly_print ( n0, c0, '  Using DIF_TO_R8POLY' );
%
%  Compute the difference table for the K=1 derivative.
%
  k = 1;
  n1 = n0 - k;
  [ x1, d1 ] = dif_derivk_table ( n0, x0, d0, k );
%
%  Compute the difference table for the K=2 derivative.
%
  k = 2;
  n2 = n0 - k;
  [ x2, d2 ] = dif_derivk_table ( n0, x0, d0, k );
%
%  Compute the difference table for the K=3 derivative.
%
  k = 3;
  n3 = n0 - k;
  [ x3, d3 ] = dif_derivk_table ( n0, x0, d0, k );
%
%  Compute the difference table for the K=4 derivative.
%
  k = 4;
  n4 = n0 - k;
  [ x4, d4 ] = dif_derivk_table ( n0, x0, d0, k );
%
%  Evaluate all 5 polynomials.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Evaluate difference tables for the function P0\n' );
  fprintf ( 1, '  and its first four derivatives, P1...P4.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X         P0        P1        P2        P3        P4\n' );
  fprintf ( 1, '\n' );

  for i = 0 : 10
    x = i / 5.0;
    y0 = dif_val ( n0, x0, d0, x );
    y1 = dif_val ( n1, x1, d1, x );
    y2 = dif_val ( n2, x2, d2, x );
    y3 = dif_val ( n3, x3, d3, x );
    y4 = dif_val ( n4, x4, d4, x );
    fprintf ( 1, '  %8.4f  %8.4f  %8.4f  %8.4f  %8.4f  %8.4f\n', ...
      x, y0, y1, y2, y3, y4 );
  end

  return
end


text:
%% Transient diffusion equation
%% PDE and boundary conditions
% The transient diffusion equation reads
%
% $$\alpha\frac{\partial c}{\partial t}+\nabla.\left(-D\nabla c\right)=0,$$
%
% where $c$ is the independent variable (concentration, temperature, etc)
% , $D$ is the diffusion coefficient, and $\alpha$ is a constant.
% Written by Ali A. Eftekhari
% Last checked: June 2021
clc;

%% Define the domain and create a mesh structure
L = 10;  % domain length
Nx = 20; % number of cells
Ntetta = 21;
m = createMeshRadial2D(Nx,Ntetta, L,2*pi);
%% Create the boundary condition structure
BC = createBC(m); % all Neumann boundary condition structure
m.cellcenters.x=m.cellcenters.x+1;
m.facecenters.x=m.facecenters.x+1;
BC.left.a(1:3:20) = 0; BC.left.b(1:3:20)=1; BC.left.c(1:3:20)=10; % left boundary
BC.right.a(:) = 0; BC.right.b(:)=1; BC.right.c(:)=0; % right boundary
% BC.right.a(1:floor(Ntetta/2)) = 0; BC.right.b(1:floor(Ntetta/2))=1; BC.right.c(1:floor(Ntetta/2))=3; % right boundary
BC.top.periodic=1;
% BC.top.a(:) = 0; BC.top.b(:)=1; BC.top.c(:)=0; % top boundary
% BC.bottom.a(:) = 0; BC.bottom.b(:)=1; BC.bottom.c(:)=0; % bottom boundary
%% define the transfer coeffs
D_val = 0.01;
D = createCellVariable(m, D_val);
alfa = createCellVariable(m, 1);
u_val=1;
u = createFaceVariable(m, [u_val, 0]);
%% define initial values
c_init = 0.1;
c_old = createCellVariable(m, c_init,BC); % initial values
c = c_old; % assign the old value of the cells to the current values
%% loop
dt = 0.5; % time step
final_t = 7;
Dave = harmonicMean(D);
Mdiff = diffusionTerm(Dave);
[Mbc, RHSbc] = boundaryCondition(BC);
FL = fluxLimiter('Superbee');
Mconv = convectionTerm(u);
for t=dt:dt:final_t
    [M_trans, RHS_trans] = transientTerm(c, dt, alfa);
    M = M_trans+Mconv-Mdiff+Mbc;
    RHS = RHS_trans+RHSbc;
    c = solvePDE(m,M, RHS);
    c_old = c.value;
    figure(1);visualizeCells(c);shading interp; drawnow;
end

author: simulkade

repo: FVTool

save_path: github-repos/MATLAB/simulkade-FVTool


text:
function [wins] = identify_subj_wins_is(ekg_RSlinB_am_ELF, ekg_RSlinB_bw_ELF, ekg_RSlinB_fm_ELF, up)
%UNTITLED4 Summary of this function goes here
%   Detailed explanation goes here

%% Identify start and end times of each respiratory signal
respSigs.ekg_RSlinB_am_ELF = ekg_RSlinB_am_ELF; respSigs.ekg_RSlinB_bw_ELF = ekg_RSlinB_bw_ELF; respSigs.ekg_RSlinB_fm_ELF = ekg_RSlinB_fm_ELF;
respSig_names = fieldnames(respSigs);
[timings.start, timings.end] = deal(nan(length(respSig_names),1));
for sig_no = 1 : length(respSig_names)
    eval(['rel_data = respSigs.' respSig_names{sig_no} '.t;']);
    timings.start(sig_no) = rel_data(1);
    timings.end(sig_no) = rel_data(end);
end

latest_start = max(timings.start);
earliest_end = min(timings.end);

%% Define windows
duration_of_one_win = up.paramSet.winLeng;
no_of_secs_bet_con_wins = up.paramSet.winStep;
gap_between_win_starts = duration_of_one_win - no_of_secs_bet_con_wins;
wins.t_start = latest_start : gap_between_win_starts : (earliest_end - duration_of_one_win);
wins.t_end = wins.t_start + duration_of_one_win;

end

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function class_map=GetClassMap(label,cmap)
img_size=size(label);
class_map=cmap(label(:)+1,:);
class_map=reshape(class_map,[img_size(1), img_size(2),3]);
class_map=uint8(class_map);

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% function mBlock= PolarToRec(mRadialSlice,cubeSide,cPyramid)
%% Generates the XYZ coordiante of Polar Grid in Cartesian Grid
%%Input:   mRadialSlice: 3-d Radial slcie matrix
%%        cubeSide: order of block matrix in 3-d to be assembled
%%        cPyramid: cell containing one of 3 pyramid information
%%Output:
%%        mBlock: block matrix in 3-d to be assembled
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mBlock= PolarToRec(mRadialIdx,mRadial,cubeSide,pyramid,F)
mBlock=zeros(cubeSide,cubeSide,cubeSide);
for j=mRadialIdx(3,1):mRadialIdx(3,2)
  for i=mRadialIdx(1,1):mRadialIdx(1,2)
    for k=mRadialIdx(2,1):mRadialIdx(2,2)
      mBlock(pyramid.X(i,k,j),pyramid.Y(i,k,j),pyramid.Z(i,k,j))=...
       +mRadial(i-mRadialIdx(1,1)+1, k-mRadialIdx(2,1)+1, j-mRadialIdx(3,1)+1);
    end
  end
end

% FToggle=F==0;
% F=F+FToggle;
%   mBlock=mBlock./F;

author: thomaskuestner


text:
function [U,XI,YI] = takeo_asap(V,F,b,bc)
  % TAKEO_ASAP Solve "As-Similar-As-Possible" according to "As-Rigid-As-Possible
  % Shape Manipulation" by Igarashi et al., this is their "First Step". This is
  % equivalent up to factor of 2.5 to lscm
  %
  % U = takeo_arap(V,F,b,bc)
  %
  % Inputs:
  %   V  #V by 2 list of rest domain positions
  %   F  #F by 3 list of triangle indices into V
  %   b  #b list of indices of constraint (boundary) vertices
  %   bc  #b by 2 list of constraint positions for b
  % Outputs:
  %   U  #V by 2 list of new positions
  %   X1  #F*3 by 1 list of coefficents, so that V(I,:) can be described in
  %     terms of V(J,:) and V(K,:), where I,J,K are the indices of each
  %     triangle
  %   Y1  #F*3 by 1 list of coefficents, see above
  %
  % Note about X1,Y1:
  % V(I,:) =
  %   (V(J,:) + [XI XI].* (V(K,:) - V(J,:)) + [YI -YI].* fliplr(V(K,:)-V(J,:)))
  %
  % See also: arap, takeo_arap, lscm
  %

  % We will build the quadratic system matrix per triangle
  % First we need for each triangle to write each vertex in terms of the other
  % two vertices.

  [XI,YI] = relative_coordinates(V,F);
  % reshape coordinate into single tall column
  XI = reshape(XI,size(F,1)*3,1);
  YI = reshape(YI,size(F,1)*3,1);

  % number of vertices
  n = size(V,1);
  % only works in 2D
  assert(size(V,2) == 2)
  % number of triangles
  nt = size(F,1);

  % Indices of each triangle vertex, I, and its corresponding two neighbors, J
  % and K
  I = [F(:,1);F(:,2);F(:,3)];
  J = [F(:,2);F(:,3);F(:,1)];
  K = [F(:,3);F(:,1);F(:,2)];
  % rename indices to so that I1 and I2 and so on index into vertex positions
  % as single column. Namely V(I,1) = V(I1) and V(I,2) = V(I2) etc.
  I1 = I;
  I2 = I+n;
  J1 = J;
  J2 = J+n;
  K1 = K;
  K2 = K+n;

  % Construct quadratic system matrix each triplet from II,JJ,VV contains an
  % set of elements per vertex per triangle in the quadratic energy summation
  II = [ ...
    J1
    J2
    J1
    K1
    J1
    K2
    J2
    K1
    J2
    K2
    K1
    K2
    I1
    J1
    I1
    J2
    I2
    J1
    I2
    J2
    I1
    I2
    I1
    K1
    I1
    K2
    I2
    K1
    I2
    K2
    ];
  JJ = [
    J1
    J2
    K1
    J1
    K2
    J1
    K1
    J2
    K2
    J2
    K1
    K2
    J1
    I1
    J2
    I1
    J1
    I2
    J2
    I2
    I1
    I2
    K1
    I1
    K2
    I1
    K1
    I2
    K2
    I2
    ];
  VV = [ ...
    (1-2*XI+XI.^2+YI.^2)
    (1-2*XI+XI.^2+YI.^2)
    (XI-XI.^2-YI.^2)
    (XI-XI.^2-YI.^2)
    (YI)
    (YI)
    (-YI)
    (-YI)
    (XI-XI.^2-YI.^2)
    (XI-XI.^2-YI.^2)
    (XI.^2+YI.^2)
    (XI.^2+YI.^2)
    (-1+XI)
    (-1+XI)
    (YI)
    (YI)
    (-YI)
    (-YI)
    (-1+XI)
    (-1+XI)
    ones(nt*3,1)
    ones(nt*3,1)
    (-XI)
    (-XI)
    (-YI)
    (-YI)
    (YI)
    (YI)
    (-XI)
    (-XI)
    ];
  % Assembly matrix
  A = sparse(II,JJ,VV,2*n,2*n);

  % solve
  U = min_quad_with_fixed(A,zeros(2*n,1),[b b+n],bc(:));
  % reshape into columns
  U = reshape(U,n,2);

end

%E = ...
%  sum(sum(((U(J,:) + [XI XI].* (U(K,:) - U(J,:)) + [YI -YI].* fliplr(U(K,:)-U(J,:))) - U(I,:)).^2,2));
%  E
%
%
%E = ...
%  sum( ...
%      (U(J,1) + XI.* U(K,1) - XI.*U(J,1) + YI.* U(K,2) - YI.*U(J,2) - U(I,1)).^2 + ...
%      (U(J,2) + XI.* U(K,2) - XI.*U(J,2) - YI.* U(K,1) + YI.*U(J,1) - U(I,2)).^2 ...
%  );
%
%E = ...
%  sum( ...
%      U(J,1).^2 + ...
%      2.*XI.*U(J,1).*U(K,1) - 2.*XI.*U(J,1).^2 + 2.*YI.*U(J,1).*U(K,2) - 2.*YI.*U(J,1).*U(J,2) - 2.*U(J,1).*U(I,1) + ...
%      (XI.* U(K,1) - XI.*U(J,1) + YI.* U(K,2) - YI.*U(J,2) - U(I,1)).^2 + ...
%      (U(J,2) + XI.* U(K,2) - XI.*U(J,2) - YI.* U(K,1) + YI.*U(J,1) - U(I,2)).^2 ...
%  );
%
%E = ...
%  sum( ...
%      U(J,1).^2 + ...
%      2.*XI.*U(J,1).*U(K,1) - 2.*XI.*U(J,1).^2 + 2.*YI.*U(J,1).*U(K,2) - 2.*YI.*U(J,1).*U(J,2) - 2.*U(J,1).*U(I,1) + ...
%      XI.^2.*U(K,1).^2 + ...
%      -2.*XI.^2.*U(J,1).*U(K,1) + 2.*XI.*YI.*U(K,1).*U(K,2) - 2.*XI.*YI.*U(K,1).*U(J,2) - 2.*XI.*U(K,1).*U(I,1) + ...
%      (- XI.*U(J,1) + YI.* U(K,2) - YI.*U(J,2) - U(I,1)).^2 + ...
%      (U(J,2) + XI.* U(K,2) - XI.*U(J,2) - YI.* U(K,1) + YI.*U(J,1) - U(I,2)).^2 ...
%  );
%
%E = ...
%  sum( ...
%      U(J,1).^2 + ...
%      2.*XI.*U(J,1).*U(K,1) - 2.*XI.*U(J,1).^2 + 2.*YI.*U(J,1).*U(K,2) - 2.*YI.*U(J,1).*U(J,2) - 2.*U(J,1).*U(I,1) + ...
%      XI.^2.*U(K,1).^2 + ...
%      -2.*XI.^2.*U(J,1).*U(K,1) + 2.*XI.*YI.*U(K,1).*U(K,2) - 2.*XI.*YI.*U(K,1).*U(J,2) - 2.*XI.*U(K,1).*U(I,1) + ...
%      XI.^2.*U(J,1).^2 + ...
%      -2.*XI.*U(J,1).*YI.* U(K,2) + 2.*XI.*U(J,1).*YI.*U(J,2) + 2.*XI.*U(J,1).*U(I,1) + ...
%      (YI.* U(K,2) - YI.*U(J,2) - U(I,1)).^2 + ...
%      (U(J,2) + XI.* U(K,2) - XI.*U(J,2) - YI.* U(K,1) + YI.*U(J,1) - U(I,2)).^2 ...
%  );
%
%E = ...
%  sum( ...
%      U(J,1).^2 + ...
%      2.*XI.*U(J,1).*U(K,1) - 2.*XI.*U(J,1).^2 + 2.*YI.*U(J,1).*U(K,2) - 2.*YI.*U(J,1).*U(J,2) - 2.*U(J,1).*U(I,1) + ...
%      XI.^2.*U(K,1).^2 + ...
%      -2.*XI.^2.*U(J,1).*U(K,1) + 2.*XI.*YI.*U(K,1).*U(K,2) - 2.*XI.*YI.*U(K,1).*U(J,2) - 2.*XI.*U(K,1).*U(I,1) + ...
%      XI.^2.*U(J,1).^2 + ...
%      -2.*XI.*U(J,1).*YI.* U(K,2) + 2.*XI.*U(J,1).*YI.*U(J,2) + 2.*XI.*U(J,1).*U(I,1) + ...
%      YI.^2.*U(K,2).^2 + ...
%      - 2.*YI.^2.*U(J,2).*U(K,2) - 2.*YI.*U(I,1).*U(K,2) + ...
%      (- YI.*U(J,2) - U(I,1)).^2 + ...
%      (U(J,2) + XI.* U(K,2) - XI.*U(J,2) - YI.* U(K,1) + YI.*U(J,1) - U(I,2)).^2 ...
%  );
%
%E = ...
%  sum( ...
%      U(J,1).^2 + ...
%      2.*XI.*U(J,1).*U(K,1) - 2.*XI.*U(J,1).^2 + 2.*YI.*U(J,1).*U(K,2) - 2.*YI.*U(J,1).*U(J,2) - 2.*U(J,1).*U(I,1) + ...
%      XI.^2.*U(K,1).^2 + ...
%      -2.*XI.^2.*U(J,1).*U(K,1) + 2.*XI.*YI.*U(K,1).*U(K,2) - 2.*XI.*YI.*U(K,1).*U(J,2) - 2.*XI.*U(K,1).*U(I,1) + ...
%      XI.^2.*U(J,1).^2 + ...
%      -2.*XI.*U(J,1).*YI.* U(K,2) + 2.*XI.*U(J,1).*YI.*U(J,2) + 2.*XI.*U(J,1).*U(I,1) + ...
%      YI.^2.*U(K,2).^2 + ...
%      - 2.*YI.^2.*U(J,2).*U(K,2) - 2.*YI.*U(I,1).*U(K,2) + ...
%      YI.^2.*U(J,2).^2 + ...
%      + 2.*YI.*U(I,1).*U(J,2) + ...
%      U(I,1).^2+...
%      (U(J,2) + XI.* U(K,2) - XI.*U(J,2) - YI.* U(K,1) + YI.*U(J,1) - U(I,2)).^2 ...
%  );
%
%E = ...
%  sum( ...
%      U(J,1).^2 + ...
%      - 2.*XI.*U(J,1).^2 + ...
%      XI.^2.*U(J,1).^2 + ...
%      YI.^2.*U(J,1).^2 + ...
%      ...
%      U(J,2).^2 + ...
%      YI.^2.*U(J,2).^2 + ...
%      - 2.*XI.*U(J,2).^2 + ...
%      XI.^2.*U(J,2).^2 + ...
%      ...
%      - 2.*YI.*U(J,1).*U(J,2) +...
%      2.*XI.*YI.*U(J,1).*U(J,2) + ...
%      2.*YI.*U(J,1).*U(J,2) + ...
%      - 2.*XI.*YI.*U(J,1).*U(J,2) + ...
%      ...
%      2.*XI.*U(J,1).*U(K,1) + ...
%      -2.*XI.^2.*U(J,1).*U(K,1) + ...
%      - 2.*YI.^2.*U(J,1).*U(K,1) + ...
%      ...
%      2.*YI.*U(J,1).*U(K,2) + ...
%      -2.*XI.*YI.*U(J,1).* U(K,2) + ...
%      2.*XI.*YI.*U(J,1).*U(K,2) + ...
%      ...
%      - 2.*XI.*YI.*U(J,2).*U(K,1) + ...
%      - 2.*YI.*U(J,2).*U(K,1) + ...
%      2.*XI.*U(J,2).*YI.* U(K,1) + ...
%      ...
%      - 2.*YI.^2.*U(J,2).*U(K,2) + ...
%      2.*XI.*U(J,2).*U(K,2) + ...
%      -2.*XI.^2.*U(J,2).*U(K,2) + ...
%      ...
%      - 2.*U(I,1).*U(J,1) + ...
%      2.*XI.*U(I,1).*U(J,1) + ...
%      - 2.*YI.*U(I,2).*U(J,1) + ...
%      ...
%      2.*YI.*U(I,1).*U(J,2) + ...
%      ...
%      - 2.*U(I,2).*U(J,2) + ...
%      2.*XI.*U(I,2).*U(J,2) + ...
%      ...
%      XI.^2.*U(K,1).^2 + ...
%      YI.^2.*U(K,1).^2 + ...
%      ...
%      YI.^2.*U(K,2).^2 + ...
%      XI.^2.*U(K,2).^2 + ...
%      ...
%      U(I,1).^2+...
%      ...
%      U(I,2).^2 + ...
%      ...
%      2.*XI.*YI.*U(K,1).*U(K,2) + ...
%      - 2.*XI.*YI.*U(K,1).*U(K,2) + ...
%      ...
%      - 2.*XI.*U(I,1).*U(K,1) + ...
%      ...
%      - 2.*YI.*U(I,1).*U(K,2) + ...
%      ...
%      2.*YI.*U(I,2).*U(K,1) + ...
%      ...
%      - 2.*XI.*U(I,2).*U(K,2) + ...
%  0 ...
%  );
%
%
%E = ...
%  sum( ...
%      (1-2*XI+XI.^2+YI.^2).*U(J,1).^2 + ...
%      ...
%      (1-2*XI+XI.^2+YI.^2).*U(J,2).^2 + ...
%      ...
%      (2*XI-2*XI.^2-2*YI.^2).*U(J,1).*U(K,1) + ...
%      ...
%      (2.*YI).*U(J,1).*U(K,2) + ...
%      ...
%      (-2*YI).*U(J,2).*U(K,1) + ...
%      ...
%      (2*XI-2*XI.^2-2*YI.^2).*U(J,2).*U(K,2) + ...
%      ...
%      (-2+2*XI).*U(I,1).*U(J,1) + ...
%      ...
%      (2.*YI).*U(I,1).*U(J,2) + ...
%      ...
%      (-2*YI).*U(I,2).*U(J,1) + ...
%      ...
%      (-2+2*XI).*U(I,2).*U(J,2) + ...
%      ...
%      (XI.^2+YI.^2).*U(K,1).^2 + ...
%      ...
%      (XI.^2+YI.^2).*U(K,2).^2 + ...
%      ...
%      (1).*U(I,1).^2+...
%      ...
%      (1).*U(I,2).^2 + ...
%      ...
%      (-2*XI).*U(I,1).*U(K,1) + ...
%      ...
%      (-2*YI).*U(I,1).*U(K,2) + ...
%      ...
%      (2*YI).*U(I,2).*U(K,1) + ...
%      ...
%      (-2*XI).*U(I,2).*U(K,2) + ...
%  0 ...
%  );
%
%E = ...
%  sum( ...
%      (1-2*XI+XI.^2+YI.^2).*U(J,1).^2 + ...
%      ...
%      (1-2*XI+XI.^2+YI.^2).*U(J,2).^2 + ...
%      ...
%      2*(XI-XI.^2-YI.^2).*U(J,1).*U(K,1) + ...
%      ...
%      2*(YI).*U(J,1).*U(K,2) + ...
%      ...
%      2*(-YI).*U(J,2).*U(K,1) + ...
%      ...
%      (2*XI-2*XI.^2-2*YI.^2).*U(J,2).*U(K,2) + ...
%      ...
%      2*(-1+XI).*U(I,1).*U(J,1) + ...
%      ...
%      2*(YI).*U(I,1).*U(J,2) + ...
%      ...
%      2*(-YI).*U(I,2).*U(J,1) + ...
%      ...
%      2*(-1+XI).*U(I,2).*U(J,2) + ...
%      ...
%      (XI.^2+YI.^2).*U(K,1).^2 + ...
%      ...
%      (XI.^2+YI.^2).*U(K,2).^2 + ...
%      ...
%      (1).*U(I,1).^2+...
%      ...
%      (1).*U(I,2).^2 + ...
%      ...
%      2*(-XI).*U(I,1).*U(K,1) + ...
%      ...
%      2*(-YI).*U(I,1).*U(K,2) + ...
%      ...
%      2*(YI).*U(I,2).*U(K,1) + ...
%      ...
%      2*(-XI).*U(I,2).*U(K,2) + ...
%  0 ...
%  );

author: alecjacobson

repo: gptoolbox


text:
addpath(genpath('../mutils/My/'));
addpath(genpath('../ptv'));
%%
dat = load('data/molli.mat');
imgs = double(reshape(dat.imgs, [size(dat.imgs, 1), size(dat.imgs, 2), 1, 1, size(dat.imgs,3)]));
imgs = imgs(50:160, 50:160, :, :, :);

opts = [];
opts.check_gradients = 155*0;
opts.mean_penalty = 1e-3;
opts.pix_resolution = [1,1];
opts.metric = 'nuclear';

opts.grid_spacing = [10, 10];
opts.isoTV = 5e-3;

opts.spline_order = 1;
opts.max_iters = 100;

opts.display = 'off';
opts.opt_method = 'lbfsg';

opts.border_mask = 6;
opts.k_down = 0.7;

tic
[voldef_pl, Tmin_pl,  Kmin_pl] = ptv_register(imgs, [], opts);
toc

%%
imgs_exp = [];
for i = 1 : size(imgs, 5)
    hold off;
    imagesc(squeeze(imgs(:,:, 1, :, i)));
    colormap gray;
    hold on;
    plot_wiremesh(squeeze(Tmin_pl(:,:, 1, :, i)), 10, 'r-');
    daspect([1,1,1]);
    F = getframe;
    imgs_exp = cat(4, imgs_exp, double(F.cdata));
%     pause(0.1);
end
%%
savegif('results/molli_register.gif', squeeze([imgs, voldef_pl]), 1/5);
savegif('results/molli_track.gif', squeeze([imgs_exp]/250), 1/5);

repo: pTVreg


text:
function centerOfMass = dtiRoiGetCenterOfMassCoord(roi)

% function centerOfMass = dtiRoiGetCenterOfMassCoord(roi)
%
% For a given ROI this simple function will compute the center of mass and
% return the coordinate for that center point.
%
% User can then use that point to create a new ROI that is centered on the
% center of mass - useful for ROIs that are functionally defined but oddly
% shapped.
%
%% EXAMPLE USAGE SCRIPT
% roiFile = 'LMT.mat';
% % Read in the roi that we want to take the center of mass for
% roi = dtiReadRoi(roiFile);
%
% % Set the radius for the new ROI sphere
% radius = 5;
%
% % Get the coords and find the center of mass
% centerOfMass = dtiRoiGetCenterOfMassCoord(roi);
%
% % Create a new roi with a X mm radius from the center coord
% coords = dtiBuildSphereCoords(centerOfMass,radius);
%
% % Set the coords for the new ROI in the old roi struct and change the name
% roi.coords = coords;
% roi.name = [roi.name '_sphere_' num2str(radius) 'mm'];
%
% % Write out the new roi.
% cd(mrvDirup(roiFile));
% dtiWriteRoi(roi,roi.name);
%
%
% HISTORY:
% 2011.04.08 LMP Wrote the thing
% 2011.06.20 LMP Added example usage script to comments
% 2011.07.21 LMP Now accepts roi as a mat file or struct

if ~isstruct(roi)
 if exist(roi,'file')
    roi = dtiReadRoi(roi);
 else
  keyboard
 end
end

centerOfMass = round(mean(roi.coords,1)*10)/10;
fprintf('Center of mass coordinate: %.1f, %.1f, %.1f\n',centerOfMass);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% This make.m is for MATLAB and OCTAVE under Windows, Mac, and Unix

try
	Type = ver;
	% This part is for OCTAVE
	if(strcmp(Type(1).Name, 'Octave') == 1)
		mex libsvmread.c
		mex libsvmwrite.c
		mex train.c linear_model_matlab.c ../linear.cpp ../tron.cpp ../blas/daxpy.c ../blas/ddot.c ../blas/dnrm2.c ../blas/dscal.c
		mex predict.c linear_model_matlab.c ../linear.cpp ../tron.cpp ../blas/daxpy.c ../blas/ddot.c ../blas/dnrm2.c ../blas/dscal.c
	% This part is for MATLAB
	% Add -largeArrayDims on 64-bit machines of MATLAB
	else
		mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims libsvmread.c
		mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims libsvmwrite.c
		mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims train.c linear_model_matlab.c ../linear.cpp ../tron.cpp ../blas/daxpy.c ../blas/ddot.c ../blas/dnrm2.c ../blas/dscal.c -output liblinear_train
		mex CFLAGS="\$CFLAGS -std=c99" -largeArrayDims predict.c linear_model_matlab.c ../linear.cpp ../tron.cpp ../blas/daxpy.c ../blas/ddot.c ../blas/dnrm2.c ../blas/dscal.c -output liblinear_predict
	end
catch
	fprintf('If make.m fails, please check README about detailed instructions.\n');
end

author: suhangpro

repo: mvcnn

save_path: github-repos/MATLAB/suhangpro-mvcnn


text:
function that = spm_swarp(this,def,M)
% Warp surface
% FORMAT that = spm_swarp(this,def)
% this - a gifti object
% def  - a deformation (nifti object or filename)
% that - the warped gifti object
%
% FORMAT that = spm_swarp(this,def,M)
% this - a gifti object
% def  - a deformation field (nx*ny*nz*1*3)
% M    - mapping from voxels to world, for deformation field
% that - the warped gifti object
%
%__________________________________________________________________________
% Copyright (C) 2009-2015 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_swarp.m 6349 2015-02-26 12:15:06Z guillaume $

%-Input arguments
%--------------------------------------------------------------------------
this = gifti(this);

if nargin<2, that = this; return; end

if ischar(def) || isa(def,'nifti')
    if ischar(def)
        def = nifti(def);
    end
    y   = def(1).dat(:,:,:,:,:);
    M   = def(1).mat;
else
    y = def;
    if nargin<3, M = eye(4); end
end

%-Apply deformation to vertices
%--------------------------------------------------------------------------
v   = double(this.vertices);
iM  = inv(M);
v   = iM(1:3,1:4) * this.mat * [v'; ones(1,size(v,1))];
xyz = {v(1,:)',v(2,:)',v(3,:)'};
v   = [spm_bsplins(y(:,:,:,1,1), xyz{:}, [1 1 1 0 0 0]),...
       spm_bsplins(y(:,:,:,1,2), xyz{:}, [1 1 1 0 0 0]),...
       spm_bsplins(y(:,:,:,1,3), xyz{:}, [1 1 1 0 0 0])];

% Much of the surface data is likely to fall outside the FOV of the
% deformation field.  For this reason, the following code attempts to
% extrapolate the surfaces outside this FOV by replacing NaNs in the
% vertice coordinates by some smooth mesh.
%--------------------------------------------------------------------------
if isfield(this, 'faces')
    f = double(this.faces);
    m = size(v,1);

    % Gradients
    b = double([v(:,1); v(:,2); v(:,3)]);
    w = isfinite(b);
    b(~w) = 0;

    % Hessian
    A = spdiags(w,0,m*3,m*3);
    W = sparse(f(:,1),f(:,1),1,m,m) + sparse(f(:,2),f(:,2),1,m,m) + sparse(f(:,3),f(:,3),1,m,m)...
        - sparse(f(:,1),f(:,2),1,m,m) - sparse(f(:,2),f(:,3),1,m,m) - sparse(f(:,3),f(:,1),1,m,m);
    W = (W'*W + 0.25*W)*1e-2;
    Z = sparse([],[],[],m,m);
    A = A + [W Z Z; Z W Z; Z Z W];

    % Solution to simple linear model
    v = full(reshape(A\b,m,3));
end

%-Generate the gifti structure for the warped data
%--------------------------------------------------------------------------
that = this;
that.vertices = v;
that.mat = eye(4); % this.mat already applied to v

author: spm


text:
% tutorial4_solutions
%   This script contains the solutions for Tutorial 4, see Tutorial 4 in
%   "RAVEN tutorials.docx" for more details.
%
%   NOTE: Many of these changes are easier to do in the Excel sheet. They
%   are done here in code just to avoid having several model files.

%Import the Excel model
model=importExcelModel('smallYeastBad.xlsx');

%Close all uptake and maximize for production
model=setParam(model,'eq',{'glcIN', 'o2IN'},[0 0]);
model=setParam(model,'obj',{'acOUT' 'biomassOUT' 'co2OUT' 'ethOUT' 'glyOUT'},[1 1 1 1 1]);
sol=solveLP(model);
printFluxes(model,sol.x,true); %Nothing produced, good

%Add some fake reactions
rxns.rxns={'FREE_ATP';'FREE_NADH';'FREE_NADPH'};
rxns.equations={'ATP <=> ADP + phosphate';'NAD(+) <=> NADH';'NADP(+) <=> NADPH'};
model=addRxns(model,rxns,2,'c');
sol=solveLP(model,1);

%Lots of ethanol produced. Also plot the equations to make the error easier
%to find
printFluxes(model,sol.x,false,[],[],'%rxnID (%rxnName):%flux\n\t%eqn\n');

%See that ADH1 should only produce one unit of ethanol. Change the reaction
%equation
model=changeRxns(model,'ADH1','acetaldehyde[c] + NADH[c] => ethanol[c] + NAD(+)[c]',3);
sol=solveLP(model,1);
printFluxes(model,sol.x,true); %Nothing produced, good

%Add excretion of all metabolites
model.b=[model.b inf(numel(model.b),1)];
sol=solveLP(model,1);
printFluxes(model,sol.x,false,10^-5,[],'%rxnID (%rxnName):\n\t%eqn\n\t%flux\n');

%By looking at the reactions which were unbalanced and that were in the
%flux list one can see that FBP should be changed to result in only one unit
%of F6P
model=changeRxns(model,'FBP','beta-D-fructofuranose 1,6-bisphosphate[c] => beta-D-fructofuranose 6-phosphate[c] + phosphate[c]',3);
sol=solveLP(model,1);
printFluxes(model,sol.x,false,10^-5,[],'%rxnID (%rxnName):\n\t%eqn\n\t%flux\n');

%The same thing again and one should change PFK to only give one unit of
%F16P
model=changeRxns(model,'PFK','ATP[c] + beta-D-fructofuranose 6-phosphate[c] => ADP[c] + beta-D-fructofuranose 1,6-bisphosphate[c]',3);
sol=solveLP(model,1);
printFluxes(model,sol.x,false,10^-5,[],'%rxnID (%rxnName):\n\t%eqn\n\t%flux\n'); %Now it works

%Set all uptakes and production to 0
model=setParam(model,'eq',getExchangeRxns(model),0);

%Since it is checked which metabolites could be consumed without
%production, one can no longer have free production of all metabolites
model.b=model.b(:,1);
I=canConsume(model);
disp(model.mets(I)); %These 12 metabolites can be consumed without any production

%Allow all uptake
model.b=[ones(numel(model.b),1)*-1000 model.b];

%Pick CO2 and force uptake of it
model=setParam(model,'eq',{'co2OUT'},-1); %Negative output means input
sol=solveLP(model);
printFluxes(model,sol.x,false,10^-5,[],'%rxnID (%rxnName):\n\t%eqn\n\t%flux\n'); %Now it works

%See that PDC converts pyruvate (3 carbons) to acetaldehyde (2 carbons)
%without any other products. If one googles, one may realize that CO2 is
%missing. This would be simpler to change in the Excel file (or using
%changeRxns), but one can change it here as an exercise. One therefore
%needs to find the index of the reactions and the index of cytosolic CO2 in
%order to change the reaction
Irxn=ismember(model.rxns,'PDC');
Imet=ismember(model.mets,'CO2_c');
model.S(Imet,Irxn)=1; %The coefficient is 1.0

%Display the new equation just to be sure
constructEquations(model,Irxn)

%The solution is now not feasible, meaning that it is no longer possible to
%force uptake of CO2 without any output
sol=solveLP(model);

%***Second part of tutorial
model=importExcelModel('smallYeastBad2.xlsx',true,false,true); %This has to be loaded with the setting to ignore error or it would find the error
[reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,false,false,false,true);
disp(deletedReactions);
disp(deletedMetabolites);

%It turned out that G15L_c was spelled G15Lc in one reaction. The best
%solution would be just to change it in the reaction list and remove the
%duplicate metabolite, but one can do it here as an exercise. The indexes
%of the two metabolites are needed.
Igood=ismember(model.mets,'G15L_c');
Ibad=ismember(model.mets,'G15Lc');

%Get all reactions and the coefficients in which the wrong one participates
%move them to be for the right one instead
model.S(Igood,:)=model.S(Igood,:)+model.S(Ibad,:);

%Delete the bad one
model=removeMets(model,'G15Lc');
[reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,false,false,false,true);
disp(deletedReactions);
disp(deletedMetabolites);

%The only difference was that there were 20 deleted metabolites instead of
%21. Nothing too spectacular. Check production can tell the user what one
%needs to connect.
[notProducedMets, ~, neededForProductionMat, minToConnect]=checkProduction(model,true,model.comps,false);

%In order to have production of all 54 metabolites one needs to enable
%production of these 12. This small model does not include net synthesis of
%co-factors, so one should concentrate on the other ones. Glycerone
%phosphate allows for connection 18 others, so it seems like a good target.
disp(minToConnect);

%If one googles around a little bit, and knows the metabolism, one would
%find that DHAP (dihydroxyacetone) and GLYP (glycerone phosphate) are
%actually synonymes. Only use DHAP
Igood=ismember(model.mets,'DHAP_c');
Ibad=ismember(model.mets,'GLYP_c');

%Get all reactions and the coefficients in which the wrong one participates
%move them to be for the right one instead
model.S(Igood,:)=model.S(Igood,:)+model.S(Ibad,:);

%Delete the bad one
model=removeMets(model,'GLYP_c');
[reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(model,false,false,false,true);
disp(deletedReactions);
disp(deletedMetabolites);
[notProducedMets, ~, neededForProductionMat, minToConnect]=checkProduction(model,true,model.comps,false);
disp(minToConnect);

%Still quite a lot of gaps and no immediate way to fix it. One could try
%including reactions from a reference network and see if that helps. Use
%the small yeast model from Tutorial 3
refModel=importExcelModel('smallYeast.xlsx');
[newConnected, cannotConnect, addedRxns, newModel]=fillGaps(model,{refModel},false);
disp(addedRxns);
disp(newConnected);

%By including the ALD6 reaction from the reference model it was possible to
%connect 21 reactions
[reducedModel, deletedReactions, deletedMetabolites]=simplifyModel(newModel,false,false,false,true);
disp(deletedMetabolites);
disp(deletedReactions);

%All the model seems to be connected

%All this stuff may be done in a more automated manner as well
model=importExcelModel('smallYeastBad2.xlsx',true,false,true);
gapReport(model,{refModel});

author: SysBioChalmers

repo: RAVEN


text:
function[varargout]=maternimp(varargin)
%MATERNIMP  Impulse response function for the Matern random process.
%
%   [T,G]=MATERNIMP(N,ALPHA,LAMBDA) returns first N points for the impulse
%   response function G for a *unit-variance* complex-valued Matern process
%   with slope parameter ALPHA and damping or range parameter LAMBDA.
%
%   The impulse response function is also known as the Green's function.
%
%   T is an array of times T=[EPSILON 1:1:N] where the first point is set to
%   a very small number, set to EPSILON=1e-6 to avoid the potential for
%   infinite values at T=0 for some parameter settings.
%
%   The sample interval is taken to be unity.  The damping parameter LAMBDA
%   is understood to have units of the inverse of the sample interval.
%
%   Note that for LAMBDA=0, the case of fractional Brownian motion, G is
%   formally defined but is not useful because it does not decay.
%
%   The input parameters ALPHA and LAMBDA may either be scalars or arrays
%   of the same length M.  If the latter, then the output autocovariance
%   function R will be a matrix with N rows and M columns.
%
%   See MATERNSPEC for a more thorough discussion of the Matern process.
%
%   For details on the Matern process and impulse response function, see:
%
%     Lilly, Sykulski, Early, and Olhede, (2017).  Fractional Brownian
%        motion, the Matern process, and stochastic modeling of turbulent
%        dispersion.  Nonlinear Processes in Geophysics, 24: 481--514.
%   __________________________________________________________________
%
%   Oscillatory Matern
%
%   MATERNIMP(N,ALPHA,C), where C is complex, uses LAMBDA=REAL(C) for the
%   damping parameter and also sets a rotation frequency to OMEGA=-IMAG(C).
%
%   The associated oscillatory process undergoes rotations at frequency
%   OMEGA in addition to the damping at timescale 1/LAMBDA.
%   __________________________________________________________________
%
%   Specified times
%
%   G=MATERNIMP(T,ALPHA,LAMBDA) where the first argument is an *array*
%   rather than a scalar, will alternately return the impulse response at
%   the specified times T.  Note T is not output again in this case.
%   __________________________________________________________________
%
%   Use in fast generation
%
%   [T,G]=MATERNIMP([],ALPHA,LAMBDA) where the first argument is the empty
%   set, evaluates the Green's function at the points T = [0:1:N-1]'+1/2.
%   This form is utilized in the fast generation method implemented by
%   MATERNOISE and discussed in Section 5 of Lilly et al. (2017).
%   __________________________________________________________________
%
%   See also MATERNSPEC, MATERNCOV, MATERNOISE, MATERNFIT.
%
%   'maternimp --t' runs some tests.
%
%   Usage:    [t,G]=maternimp(N,alpha,lambda);
%             G=maternimp(t,alpha,lambda);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2013--2017 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--t')
    maternimp_test,return
end

N=[];
tau=[];

if length(varargin)==3
    if length(varargin{1})==1
        N=varargin{1};
    else
        tau=varargin{1};
    end
end

if isempty(tau)
    tau=[0:1:N-1]'+1/2;
end

alpha=varargin{2};
lambda=real(varargin{3});
omegao=-imag(varargin{3});
arrayify(alpha,omegao,lambda);

G=zeros(size(tau,1),length(lambda));
if size(tau,2)~=length(lambda)
    tau=vrep(tau,length(lambda),2);
end
for i=1:length(lambda)
    U=frac(1,2).*(1+sign(tau(:,i)));
    d=sqrt(frac(lambda(i).^(2*alpha(i)-1),materncfun(alpha(i))));
    G(:,i)=d.*U.*(tau(:,i).^(alpha(i)-1)).*exp(-tau(:,i).*lambda(i)).*exp(sqrt(-1)*omegao(i)*tau(:,i)).*frac(1,gamma(alpha(i)));
end

if ~isempty(N)
    varargout{1}=tau;
    varargout{2}=G;
else
    varargout{1}=G;
end

function[]=maternimp_test
alpha=[1 1.5 2 3 4];
h=[.01 .02 .05 .2 1]/10;
[alpha,h]=meshgrid(alpha,h);
vcolon(alpha,h);
[tau,g]=maternimp(15000,alpha,h);
d=sqrt(frac(h.^(2*alpha-1),materncfun(alpha)));
for i=1:size(g,2)
    g(:,i)=g(:,i)./d(i);
end
reporttest('MATERNIMP integral of Green''s function matches analytic form',aresame((sum(g,1)'.*(h.^alpha))'-1,zeros(size(h')),2.5e-2))

t=[0:1:150]';
alpha=1.4;h=1/10;
g=maternimp(t,alpha,h);
R1=conv(g,flipud(g));
[tau,R2]=materncov(1,length(t),1,alpha,h,'full');
eps=sum(squared(R1-R2))./sum(squared(R1));
reporttest('MATERNIMP Green''s function convolves to autocovariance sample rate',eps<1e-3)

N=1000;
alpha=1.4;h=1/10;
[tau,g]=maternimp(N,alpha,h);
R1=conv(g,flipud(g));
[tau,R2]=materncov(1,N,1,alpha,h,'full');
eps=sum(squared(R1-R2))./sum(squared(R1));
reporttest('MATERNIMP Green''s function convolves to autocovariance without time input',eps<1e-3)

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
%% pathLength
% Below is a demonstration of the features of the |pathLength| function

%%
clear; close all; clc;

%% Syntax
% |D=pathLength(V);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*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 r8mat_solve_3d_test ( )

%*****************************************************************************80
%
%% R8MAT_SOLVE_3D_TEST tests R8MAT_SOLVE_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 June 2011
%
%  Author:
%
%    John Burkardt
%
  n = 3;
  test_num = 5;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8MAT_SOLVE_3D_TEST\n' );
  fprintf ( 1, '  R8MAT_SOLVE_3D solves 3D linear systems.\n' );

  seed = 123456789;

  for test = 1 : test_num

    [ a, seed ] = r8mat_uniform_01 ( n, n, seed );
    [ x, seed ] = r8vec_uniform_01 ( n, seed );
    b(1:n,1) = a(1:n,1:n) * x(1:n);

    [ x2, det ] = r8mat_solve_3d ( a, b );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Solution / Computed:\n' );
    fprintf ( 1, '\n' );

    for i = 1 : n
      fprintf ( 1, '  %14f  %14f\n', x(i), x2(i) );
    end

  end

  return
end

author: johannesgerer


text:
function [yi, ypi] = hermint(x, y, yp, xi, c)

% HERMINT 1-D piecewise hermite interpolation
%    HERMINT(X,Y,YP,XI,C) interpolates to find YI, the values of
%    the underlying function Y at the points in the array XI,
%    using piecewise hermite interpolation.  X and Y must be
%    vectors of length N.
%
%    C specifies the number of data points to use in the
%    interpolation.  The default is to use all points.
%
%    [YI,YPI] = HERMINT() also returns the interpolated derivative
%    of the underlying function Y at points XI.

% Joe Henning - Fall 2011

if (nargin < 5)
   c = 0;
end

n = length(x);

if (n < c)
   fprintf('??? Bad c input to hermint ==> c <= length(x)\n');
   yi = [];
   ypi = [];
   return
end

for i = 1:length(xi)
   % Find the right place in the table by means of a bisection.
   klo = 1;
   khi = n;
   while (khi-klo > 1)
      k = fix((khi+klo)/2.0);
      if (x(k) > xi(i))
         khi = k;
      else
         klo = k;
      end
   end

   h = x(khi) - x(klo);
   if (h == 0.0)
      fprintf('??? Bad x input to hermint ==> x values must be distinct\n');
      yi(i) = NaN;
      ypi(i) = NaN;
      continue;
   end

   isiny = 0;
   for k = 1:n
      if (xi(i) == x(k))
         yi(i) = y(k);
         ypi(i) = yp(k);
         isiny = 1;
         break
      end
   end

   if (isiny)
      continue
   end

   % Evaluate hermite polynomial
   yi(i) = 0;
   ypi(i) = 0;
   if (c == 0)
      for k = 1:n
         f = 0;
         g = 0;
         h = 1;

         for m = 1:n
            if (k ~= m)
               f = f + 1/(x(k)-x(m));
               h = h*(xi(i)-x(m))/(x(k)-x(m));
            end
         end

         for m = 1:n
            if (k ~= m)
               g = g + h/(xi(i)-x(m));
            end
         end

         a = h*h;
         b = 2*h*g;
         u = xi(i) - x(k);

         b1 = a*u;
         b2 = b*u + a;
         a1 = a - 2*f*b1;
         a2 = b - 2*f*b2;

         yi(i) = yi(i) + a1*y(k) + b1*yp(k);
         ypi(i) = ypi(i) + a2*y(k) + b2*yp(k);
      end
   else
      if (mod(c,2) == 0)   % even
         c2 = c/2;
         if (klo < c2)
            klo = c2;
         end
         if (klo > n-c2)
            klo = n-c2;
         end
         khi = klo + 1;
         for k = klo-(c2-1):klo+c2
            f = 0;
            g = 0;
            h = 1;

            for m = klo-(c2-1):klo+c2
               if (k ~= m)
                  f = f + 1/(x(k)-x(m));
                  h = h*(xi(i)-x(m))/(x(k)-x(m));
               end
            end

            for m = klo-(c2-1):klo+c2
               if (k ~= m)
                  g = g + h/(xi(i)-x(m));
               end
            end

            a = h*h;
            b = 2*h*g;
            u = xi(i) - x(k);

            b1 = a*u;
            b2 = b*u + a;
            a1 = a - 2*f*b1;
            a2 = b - 2*f*b2;

            yi(i) = yi(i) + a1*y(k) + b1*yp(k);
            ypi(i) = ypi(i) + a2*y(k) + b2*yp(k);
         end
      else   % odd
         c2 = floor(c/2);
         if (klo < c2+1)
            klo = c2+1;
         end
         if (klo > n-c2)
            klo = n-c2;
         end
         khi = klo + 1;
         for k = klo-c2:klo+c2
            f = 0;
            g = 0;
            h = 1;

            for m = klo-c2:klo+c2
               if (k ~= m)
                  f = f + 1/(x(k)-x(m));
                  h = h*(xi(i)-x(m))/(x(k)-x(m));
               end
            end

            for m = klo-c2:klo+c2
               if (k ~= m)
                  g = g + h/(xi(i)-x(m));
               end
            end

            a = h*h;
            b = 2*h*g;
            u = xi(i) - x(k);

            b1 = a*u;
            b2 = b*u + a;
            a1 = a - 2*f*b1;
            a2 = b - 2*f*b2;

            yi(i) = yi(i) + a1*y(k) + b1*yp(k);
            ypi(i) = ypi(i) + a2*y(k) + b2*yp(k);
         end
      end
   end
end

author: Sable


text:
function mn = r8cvv_size ( m, nr )

%*****************************************************************************80
%
%% R8CVV_SIZE determines the size of an R8CVV.
%
%  Discussion:
%
%    An R8CVV is a "vector of vectors" of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 November 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows in the array.
%
%    Input, integer NR(M), the size of each row.
%
%    Output, integer MN, the size of the cell array.
%
  mn = sum ( nr(1:m) );

  return
end


text:
classdef QuadrilateralNodesCalculator < NodesCalculator

    properties (Access = public)
        nvert
        div
        boundNodes
        totalNodes
    end

    methods (Access = public)

        function obj = QuadrilateralNodesCalculator(cParams)
            obj.init(cParams);
            obj.computeBoundaryNodes();
            obj.computeTotalNodes();
        end

    end

    methods (Access = private)

        function computeTotalNodes(obj)
            divA = obj.div(1);
            divB = obj.div(2);
            obj.totalNodes = obj.boundNodes+(divA-1)*(divB-1);
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
maxiter = 20;
restartProb = 0.50;
dim_drug = 100;
dim_prot = 400;

drugNets = {'Sim_mat_drug_drug', 'Sim_mat_drug_disease', 'Sim_mat_drug_se', 'Sim_mat_Drugs'};
proteinNets = {'Sim_mat_protein_protein', 'Sim_mat_protein_disease', 'Sim_mat_Proteins'};

tic
X = DCA(drugNets, dim_drug, restartProb, maxiter);
toc
tic
Y = DCA(proteinNets, dim_prot, restartProb, maxiter);
toc

dlmwrite(['../feature/drug_vector_d', num2str(dim_drug), '.txt'], X, '\t');
dlmwrite(['../feature/protein_vector_d', num2str(dim_prot), '.txt'], Y, '\t');

author: luoyunan

repo: DTINet

save_path: github-repos/MATLAB/luoyunan-DTINet


text:
function [tfr,t,f,wt]=tfrscalo(X,time,wave,fmin,fmax,N,trace);
%TFRSCALO Scalogram, for Morlet or Mexican hat wavelet.
%	[TFR,T,F,WT]=TFRSCALO(X,T,WAVE,FMIN,FMAX,N,TRACE) computes
%	the scalogram (squared magnitude of a continuous wavelet
%	transform).
%
%	X : signal (in time) to be analyzed (Nx=length(X)). Its
%	    analytic version is used (z=hilbert(real(X))).
%	T : time instant(s) on which the TFR is evaluated
%	     					(default : 1:Nx).
%	WAVE : half length of the Morlet analyzing wavelet at coarsest
% 	    scale. If WAVE = 0, the Mexican hat is used. WAVE can also be
%           a vector containing the time samples of any bandpass
%           function, at any scale.        	(default : sqrt(Nx)).
%	FMIN,FMAX : respectively lower and upper frequency bounds of
%	    the analyzed signal. These parameters fix the equivalent
%	    frequency bandwidth (expressed in Hz). When unspecified, you
%	    have to enter them at the command line from the plot of the
%	    spectrum. FMIN and FMAX must be >0 and <=0.5.
%	N : number of analyzed voices.
%	TRACE : if nonzero, the progression of the algorithm is shown
%	                                 	(default : 0).
%	TFR : time-frequency matrix containing the coefficients of the
%	    decomposition (abscissa correspond to uniformly sampled time,
%	    and ordinates correspond to a geometrically sampled
%	    frequency). First row of TFR corresponds to the lowest
%	    frequency. When called without output arguments, TFRSCALO
%	    runs TFRQVIEW.
%	F : vector of normalized frequencies (geometrically sampled
%	    from FMIN to FMAX).
%	WT : Complex matrix containing the corresponding wavelet
%	    transform. The scalogram TFR is the square modulus of WT.
%
%	Example :
%	 sig=altes(64,0.1,0.45); tfrscalo(sig);
%
%	See also all the time-frequency representations listed in
%	the file CONTENTS (TFR*)

%	P. Goncalves, October 1995 - O. Lemoine, June 1996.
%	Copyright (c) 1995 Rice University - CNRS 1996.
%
%  This program is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published by
%  the Free Software Foundation; either version 2 of the License, or
%  (at your option) any later version.
%
%  This program is distributed in the hope that it will be useful,
%  but WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

if (nargin == 0),
 error('At least one parameter required');
end;

[xrow,xcol] = size(X);
if nargin<=6, trace=0; end

if (nargin == 1),
 time=1:xrow; wave=sqrt(xrow);
elseif (nargin == 2),
 wave=sqrt(xrow);
elseif (nargin==4),
 disp('FMIN will not be taken into account. Determine it with FMAX');
 disp('     from the following plot of the spectrum.');
elseif nargin==5,
 N=xrow;
end;

[trow,tcol] = size(time);
if (xcol==0)|(xcol>2),
 error('X must have one or two columns');
elseif (trow~=1),
 error('TIME must only have one row');
elseif wave<0,
 error('WAVE must be positive');
end;

s = (real(X) - mean(real(X)))';
z = hilbert(s) ;

if trace, disp('Scalogram distribution'); end;

if nargin<=4		        % fmin,fmax,N unspecified
 STF = fft(fftshift(z(min(time):max(time)))); Nstf=length(STF);
 sp = (abs(STF(1:round(Nstf/2)))).^2; Maxsp=max(sp);
 f = linspace(0,0.5,round(Nstf/2)+1) ; f = f(1:round(Nstf/2));
 plot(f,sp) ; grid;
 xlabel('Normalized frequency');
 title('Analyzed signal energy spectrum');
 indmin=min(find(sp>Maxsp/100));
 indmax=max(find(sp>Maxsp/100));
 fmindflt=max([0.01 0.05*fix(f(indmin)/0.05)]);
 fmaxdflt=0.05*ceil(f(indmax)/0.05);
 txtmin=['Lower frequency bound [',num2str(fmindflt),'] : '];
 txtmax=['Upper frequency bound [',num2str(fmaxdflt),'] : '];
 fmin = input(txtmin); fmax = input(txtmax);
 if isempty(fmin), fmin=fmindflt; end
 if isempty(fmax), fmax=fmaxdflt; end
 txt=['Number of frequency samples [',num2str(2^nextpow2(xrow)),'] : '];
 N=input(txt);
 if isempty(N), N=2^nextpow2(xrow); end
end

fmin_s=num2str(fmin); fmax_s=num2str(fmax);
N_s=num2str(N);

if fmin >= fmax
 error('FMAX must be greater or equal to FMIN');
elseif fmin<=0.0 | fmin>0.5,
 error('FMIN must be > 0 and <= 0.5');
elseif fmax<=0.0 | fmax>0.5,
 error('FMAX must be > 0 and <= 0.5');
end
if trace,
 disp(['Frequency runs from ',fmin_s,' to ',fmax_s,' with ',N_s,' points']);
end

f = logspace(log10(fmin),log10(fmax),N);
a = logspace(log10(fmax/fmin),log10(1),N);

wt =zeros(N,tcol);
tfr=zeros(N,tcol);

if wave > 0
 if trace, disp(['using a Morlet wavelet']) ; end
 for ptr=1:N,
  if trace, disprog(ptr,N,10); end
  nha = wave*a(ptr);
  tha = -round(nha) : round(nha);
  ha  = exp(-(2*log(10)/nha^2)*tha.^2).*exp(i*2*pi*f(ptr)*tha);
  detail = conv(z,ha)./sqrt(a(ptr));
  detail = detail(round(nha)+1:length(detail)-round(nha)) ;
  wt(ptr,:)  = detail(time) ;
  tfr(ptr,:) = detail(time).*conj(detail(time)) ;
 end
elseif wave == 0
 if trace, disp(['using a Mexican hat wavelet']) ; end
 for ptr = 1:N
  if trace, disprog(ptr,N,10); end
  ha  = mexhat(f(ptr)) ;
  nha = (length(ha)-1)/2 ;
  detail = conv(z,ha)./sqrt(a(ptr));
  detail = detail(round(nha)+1:length(detail)-round(nha)) ;
  wt(ptr,:)  = detail(time);
  tfr(ptr,:) = detail(time).*conj(detail(time)) ;
 end
elseif length(wave) > 1
 [rwav,cwav]=size(wave);
 if cwav>rwav, wave=wave.'; end
 wavef = fft(wave) ;
 nwave = length(wave) ;
 f0 = find(abs(wavef(1:nwave/2)) == max(abs(wavef(1:nwave/2))));
 f0 = mean((f0-1).*(1/nwave));
 if trace, disp(['mother wavelet centered at f0 = ',num2str(f0)]); end
 a = logspace(log10(f0/fmin),log10(f0/fmax),N);
 B = 0.99;
 R = B/((1.001)/2);
 nscale = max(128,round((B*nwave*(1+2/R)*log((1+R/2)/(1-R/2)))/2));
 if trace, disp('Scale computation :'); end
 wts = scale(wave,a,fmin,fmax,nscale,trace);
 for ptr = 1:N,
  clear detail
  if trace, disprog(ptr,N,10); end
  ha = wts(ptr,:);
  nha = length(ha)/2;
  detail = conv(z,ha)./sqrt(a(ptr));
  detail = detail(fix(nha):length(detail)-round(nha));
  wt(ptr,:) = detail(time);
  tfr(ptr,:) = detail(time).*conj(detail(time));
 end
end

t = time;
f = f';

% Normalization
SP = fft(z);
indmin = 1+round(fmin*(xrow-2));
indmax = 1+round(fmax*(xrow-2));
SPana=SP(indmin:indmax);
tfr=tfr*norm(SPana)^2/integ2d(tfr,t,f)/N;

if (nargout==0),
 tfrqview(tfr,hilbert(real(X)),t,'tfrscalo',wave,N,f);
end;

author: HeLiangHIT


text:
function batch = read_batch(model, file_list, translation)
% load the input batch from file(pre-computed mat files)
% translation: 1 for translating the data by small move and 0 for nothing.
% Hopefully, introducing translation on data would train a better model.

batch_size = length(file_list);
batch = zeros([batch_size, model.layers{1}.layerSize], 'int8');
for i = 1 : batch_size
    load(file_list(i).filename);
    if translation
        orient = randi(7);
        move = 2;
        shifted = zeros(size(instance), 'int8');
        switch orient
            case 1, % 'z'
                shifted(:,:,1:end-move) = instance(:,:,move+1:end);
            case 2, % 'z'
                shifted(:,:,move+1:end) = instance(:,:,1:end-move);
            case 3, % 'y'
                shifted(:,1:end-move,:) = instance(:,move+1:end,:);
            case 4, % 'y'
                shifted(:,move+1:end,:) = instance(:,1:end-move,:);
            case 5, % 'x'
                shifted(1:end-move,:,:) = instance(move+1:end,:,:);
            case 6, % 'x'
                shifted(move+1:end,:,:) = instance(1:end-move,:,:);
            case 7, % 'none'
                shifted = instance;
        end
    else
        shifted = instance;
    end
    batch(i,:,:,:) = shifted;
end

batch = single(batch);

author: zhirongw


text:
function make
% detect platform

compstr = computer;
is64bit = strcmp(compstr(end-1:end),'64');

% compilation parameters

compile_params = cell(0);
if (is64bit)
  compile_params{1} = '-largeArrayDims';
end

% Compile files %

sources = {'mexutils.c'};

cd mex
disp('Compiling sampleTau...');
mex('sampleTau.cpp', sources{:},compile_params{:});
cd ..

author: andrewssobral

repo: lrslibrary


text:
function [ALLEEG cfg] = pop_est_validateMVAR(ALLEEG,typeproc,varargin)
%
% Validate a fitted VAR model. With two inputs, this function generates a
% GUI where the validation scheme can be specified. Validation consists of
% statistical tests for "whiteness" of fitted VAR model residuals [1,2],
% consistency of the fitted model [1,3] and stability of fitted model [1-3]
%
% Input:
%
%   ALLEEG:     Array of EEGLAB data structures containing fitted MODEL
%   typeproc:   reserved for future use. Use 0
%
% Optional:
%
%     'whitenessCriteria':    Cell array containing names of residual whiteness
%                             test criteria to evaluate. See [1, 2] for details.
%                             Possible Values: {'Ljung-Box','ACF','Box-Pierce','Li-McLeod'}
%                             Default Value  : all
%                             Data Input Type: cell array
%
%
%     'checkWhiteness':       Whether or not to check whiteness.
%                             Default Value  : true
%                             Data Input Type: boolean
%
%     'checkConsistency':     Whether or not to check consistency. See [1,3]
%                             for details.
%                             Default Value  : true
%                             Data Input Type: boolean
%
%     'checkStability':       Whether or not to check stability. See
%                             [1-3] for details.
%                             Default Value  : true
%                             Data Input Type: boolean
%
%     'alpha':                significance level for determining whiteness
%                             Data Input Range: [0 1]
%                             Default Value   : 0.05
%                             Data Input Type : real number (double)
%
%     'prctWinToSample':      percent of time windows to randomly select
%                             Data Input Range: [0 100]
%                             Default Value   : 100
%                             Data Input Type : real number (double)
%
%     'verb':                 verbosity level (0=no output, 1=text, 2=gui)
%
%
% Output:
%
%     whitestats:             Structure containing whiteness statistics.
%                             See est_checkMVARWhiteness() for details on
%                             structure format
%
%     PC:                     Vector of percent consistency estimates for
%                             each window.
%
%     stability:              Vector of stability estimates for each window
%     cfg:                    parameter configuration object
%
% See Also: est_checkMVARWhiteness(), est_checkMVARStability(),
%           est_checkMVARConsistency, pop_est_fitMVAR()
%
%
% References:
%
% [1] Mullen T (2010) The Source Information Flow Toolbox (SIFT):
%   Theoretical Handbook and User Manual. Chapter 3.6 and 6.
%   Available at: http://www.sccn.ucsd.edu/wiki/Sift
%
% [2] Lutkepohl, H. (2007) New Introduction to Time Series Analysis.
%   Springer.
%
% [3] Ding M, Bressler SL, Yang W, Liang H (2000) Short-window spectral
%   analysis of cortical event-related potentials by adaptive multivariate
%   autoregressive modeling: data preprocessing, model validation, and
%   variability assessment. Biol. Cybern. 83:35-45
%
% Author: Tim Mullen, 2010, SCCN/INC, UCSD.
% Email:  tim@sccn.ucsd.edu

% This function is part of the Source Information Flow Toolbox (SIFT)
%
% 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, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

if nargin<2
    typeproc = 0;
end

fcnName     = strrep(mfilename,'pop_','');
fcnHandle   = str2func(fcnName);

% check the dataset
res = hlp_checkeegset(ALLEEG,{'model'});
if ~isempty(res)
    error(['SIFT:' fcnName],res{1});
end

if isfield(ALLEEG(1).CAT.configs,fcnName)
    % get default configuration (from prior use) and merge with varargin
    varargin = [hlp_struct2varargin(ALLEEG(1).CAT.configs.(fcnName)) varargin];
end

if strcmpi(typeproc,'nogui')
    % get the config from function
    cfg = arg_tovals(arg_report('rich',fcnHandle,[{'EEG',ALLEEG(1)},varargin]),false);
else
    % render the GUI
    [PGh figh] = feval(['gui_' fcnName],ALLEEG(1),varargin{:});

    if isempty(PGh)
        % user chose to cancel
        cfg = [];
        return;
    end

    % get the specification of the PropertyGrid
    ps = PGh.GetPropertySpecification;
    cfg = arg_tovals(ps,false);
end

drawnow;

if strcmpi(typeproc,'cfg_only')
    return;
end

% Apply model validation routines
for cnd=1:length(ALLEEG)

    [ALLEEG(cnd).CAT.VALIDATION.whitestats ...
     ALLEEG(cnd).CAT.VALIDATION.PCstats    ...
     ALLEEG(cnd).CAT.VALIDATION.stabilitystats ...
     ALLEEG(cnd).CAT.VALIDATION.residualstats] ...
        = est_validateMVAR('EEG',ALLEEG(cnd),cfg);

    if ~isempty(cfg)
        % store the configuration structure
        ALLEEG(cnd).CAT.configs.(fcnName) = cfg;
    end
end

repo: REST


text:
c = 340;                    % Sound velocity (m/s)
fs = 16000;                 % Sample frequency (samples/s)
r = [2 1.5 3 ; 1 1.5 2];    % Receiver positions [x_1 y_1 z_1 ; x_2 y_2 z_2] (m)
s = [2 3.5 2];              % Source position [x y z] (m)
L = [5 4 6];                % Room dimensions [x y z] (m)
beta = 0.4;                 % Reverberation time (s)
n = 4096;                   % Number of samples
mtype = 'omnidirectional';  % Type of microphone
order = -1;                 % -1 equals maximum reflection order!
dim = 3;                    % Room dimension
orientation = 0;            % Microphone orientation (rad)
hp_filter = 1;              % Enable high-pass filter

h = rir_generator(c, fs, r, s, L, beta, n, mtype, order, dim, orientation, hp_filter);

author: ZitengWang

repo: MASP

save_path: github-repos/MATLAB/ZitengWang-MASP


text:
function test_bug1998

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_preprocessing ft_read_data read_neuralynx_ncs

% this bug is detailled on http://bugzilla.fieldtriptoolbox.org/show_bug.cgi?id=1998
% and the workaround is explained on http://www.fieldtriptoolbox.org/getting_started/neuralynx?&#discontinuous_recordings

cd(dccnpath('/home/common/matlab/fieldtrip/data/test/bug1998'));

% start with normal preprocessing of a single channel
cfg         = [];
cfg.dataset = 'CSC1.Ncs';
data        = ft_preprocessing(cfg);

% Warning: discontinuous recording, predicted number of timestamps and observed number of timestamps differ by 1693523717.00
%  Please consult the wiki on http://www.fieldtriptoolbox.org/getting_started/neuralynx?&#discontinuous_recordings
% > In fileio/private/read_neuralynx_ncs at 94
%   In ft_read_header at 1196
%   In ft_preprocessing at 394
%
%
% >> disp(data)
%            hdr: [1x1 struct]
%          label: {'CSC1'}
%           time: {[1x12902912 double]}
%          trial: {[1x12902912 double]}
%        fsample: 1893
%     sampleinfo: [1 12902912]
%            cfg: [1x1 struct]

figure
plot(data.time{1})
xlabel('sample number')
ylabel('time (s)')

ts = ft_read_data(cfg.dataset, 'timestamp', 'true'); % raw timestamps
ts = double(ts); % convert from uint64 into double

figure
plot(ts)
xlabel('sample number')
ylabel('timestamp (a.u.)')

% determine the actual timestamps per sample and interpolate the data
dts       = median(diff(ts));
tsinterp  = ts(1):dts:ts(end);
datinterp = interp1(ts, data.trial{1}, tsinterp);

% you can use NaN to replace the data in the gaps
gaps     = find(diff(ts)>2*dts); % skips at least a sample
for igap = 1:length(gaps)
  sel = tsinterp < ts(gaps(igap)+1) & tsinterp > ts(gaps(igap));
  datinterp(:,sel) = NaN;
end

% update the FieldTrip data structure
data.trial{1} = datinterp;
data.time{1}  = (tsinterp - ts(1)) / (data.hdr.Fs .* dts);

author: fieldtrip

repo: fieldtrip


text:
function X = sethackflag(X,flag)
%GETHACKFLAG Internal function to set constraint type

X.typeflag = flag;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
classdef OptimalExponentComputer < handle

    properties (Access = public)
       qOpt
       qMin
       qMax
    end

    properties (Access = private)
        maxStress
        mx
        my
        rho
        txi
        rhoV
        txiV
        mxMax
        myMax
        phi
        phiV
        samplePoints
        fileName
    end

    methods (Access = public)

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

        function compute(obj)
            nphi = length(obj.phiV);
            npoints = length(obj.rhoV);
            obj.qOpt = zeros(nphi,npoints);
            obj.qMin = zeros(nphi,npoints);
            obj.qMax = zeros(nphi,npoints);
            for iphi = 1:nphi
                obj.phi = obj.phiV(iphi);
                disp([num2str(iphi/nphi*100),'%'])
                for ipoint = 1:npoints
                    obj.rho = obj.rhoV(ipoint);
                    obj.txi = obj.txiV(ipoint);
                    [q,qmin,qmax] = obj.computeOptimalExponent(ipoint,iphi);
                    obj.qOpt(iphi,ipoint) = q;
                    obj.qMin(iphi,ipoint) = qmin;
                    obj.qMax(iphi,ipoint) = qmax;
                end
            end
            x.rho = obj.rhoV;
            x.txi = obj.txiV;
            x.phi = obj.phiV;
            x.q = obj.qOpt;
            save(obj.fileName,'x');
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.fileName = cParams.fileName;
            obj.samplePoints = cParams.samplePoints;
            obj.computeRhoTxiAndPsiSamples();
        end

        function computeRhoTxiAndPsiSamples(obj)
            obj.rhoV = obj.samplePoints.rhoV;
            obj.txiV = obj.samplePoints.txiV;
            obj.phiV = obj.samplePoints.phiV;
        end

        function [q,qMin,qMax] = computeOptimalExponent(obj,ipoint,iphi)
            s.fileName = [obj.fileName,'Txi',num2str(ipoint),'Phi',num2str(iphi)];
            s.rho   = obj.rho;
            s.txi   = obj.txi;
            s.phi   = obj.phi;
            s.pNorm = 'max';
            s.hMesh = 0.1;
            c = OneOptimalExponentComputerAndFunctionVariation(s);
            c.computeOptimalExponent();
            %c.printOptimalMicroStructure();

            qMin = c.qOptIter();
            fMin = c.fOptIter();

            [~,ind] = min(fMin);
            q = qMin(ind);
            qMin = c.qMax;
            qMax = c.qMin;
        end

    end

end


text:
function [bp,wf]=grule(n)
%
% [bp,wf]=grule(n)
%  This function computes Gauss base points and weight factors
%  using the algorithm given by Davis and Rabinowitz in 'Methods
%  of Numerical Integration', page 365, Academic Press, 1975.
%
bp=zeros(n,1); wf=bp; iter=2; m=fix((n+1)/2); e1=n*(n+1);
mm=4*m-1; t=(pi/(4*n+2))*(3:4:mm); nn=(1-(1-1/n)/(8*n*n));
xo=nn*cos(t);
for j=1:iter
   pkm1=1; pk=xo;
   for k=2:n
      t1=xo.*pk; pkp1=t1-pkm1-(t1-pkm1)/k+t1;
      pkm1=pk; pk=pkp1;
   end
   den=1.-xo.*xo; d1=n*(pkm1-xo.*pk); dpn=d1./den;
   d2pn=(2.*xo.*dpn-e1.*pk)./den;
   d3pn=(4*xo.*d2pn+(2-e1).*dpn)./den;
   d4pn=(6*xo.*d3pn+(6-e1).*d2pn)./den;
   u=pk./dpn; v=d2pn./dpn;
   h=-u.*(1+(.5*u).*(v+u.*(v.*v-u.*d3pn./(3*dpn))));
   p=pk+h.*(dpn+(.5*h).*(d2pn+(h/3).*(d3pn+.25*h.*d4pn)));
   dp=dpn+h.*(d2pn+(.5*h).*(d3pn+h.*d4pn/3));
   h=h-p./dp; xo=xo+h;
end
bp=-xo-h;
fx=d1-h.*e1.*(pk+(h/2).*(dpn+(h/3).*(d2pn+(h/4).*(d3pn+(.2*h).*d4pn))));
wf=2*(1-bp.^2)./(fx.*fx);
if ( (m+m) > n )
	bp(m)=0;
end
if ( ~ ((m+m) == n) )
	m=m-1;
end
jj=1:m; n1j=(n+1-jj); bp(n1j)=-bp(jj); wf(n1j)=wf(jj);
bp = reshape(bp,length(bp),1);
wf = reshape(wf,length(wf),1);

end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function p = iptancestor(h,type)
%IPTANCESTOR Get object ancestor.
%   P = IPTANCESTOR(H,TYPE) returns the handle of the ancestor of h of type
%   TYPE.
%
%   Note: IPTANCESTOR does no error checking and assumes that the sought after
%   ancestor is alive and well. If this is not guaranteed, use ANCESTOR
%   instead.
%
%   See also ANCESTOR.

%   Copyright 2005 The MathWorks, Inc.
%   $Revision $  $Date: 2005/06/20 03:10:23 $

p = h;
while ~strcmpi(get(p,'type'),type)
  p = get(p,'parent');
end

repo: mcbench-benchmarks


text:
function [ c, seed ] = i4_uniform_ab ( a, b, seed )

%*****************************************************************************80
%
%% I4_UNIFORM_AB returns a scaled pseudorandom I4.
%
%  Discussion:
%
%    The pseudorandom number will be scaled to be uniformly distributed
%    between A and B.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 November 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Linus Schrage,
%    A Guide to Simulation,
%    Springer Verlag, pages 201-202, 1983.
%
%    Pierre L'Ecuyer,
%    Random Number Generation,
%    in Handbook of Simulation,
%    edited by Jerry Banks,
%    Wiley Interscience, page 95, 1998.
%
%    Bennett Fox,
%    Algorithm 647:
%    Implementation and Relative Efficiency of Quasirandom
%    Sequence Generators,
%    ACM Transactions on Mathematical Software,
%    Volume 12, Number 4, pages 362-376, 1986.
%
%    Peter Lewis, Allen Goodman, James Miller
%    A Pseudo-Random Number Generator for the System/360,
%    IBM Systems Journal,
%    Volume 8, pages 136-143, 1969.
%
%  Parameters:
%
%    Input, integer A, B, the minimum and maximum acceptable values.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, integer C, the randomly chosen integer.
%
%    Output, integer SEED, the updated seed.
%
  i4_huge = 2147483647;

  if ( seed == 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4_UNIFORM_AB - Fatal error!\n' );
    fprintf ( 1, '  Input SEED = 0!\n' );
    error ( 'I4_UNIFORM_AB - Fatal error!' );
  end

  seed = floor ( seed );
  a = round ( a );
  b = round ( b );

  seed = mod ( seed, i4_huge );

  if ( seed < 0 )
    seed = seed + i4_huge;
  end

  k = floor ( seed / 127773 );

  seed = 16807 * ( seed - k * 127773 ) - k * 2836;

  if ( seed < 0 )
    seed = seed + i4_huge;
  end

  r = seed * 4.656612875E-10;
%
%  Scale R to lie between A-0.5 and B+0.5.
%
  r = ( 1.0 - r ) * ( min ( a, b ) - 0.5 ) ...
    +         r   * ( max ( a, b ) + 0.5 );
%
%  Use rounding to convert R to an integer between A and B.
%
  value = round ( r );

  value = max ( value, min ( a, b ) );
  value = min ( value, max ( a, b ) );

  c = value;

  return
end

repo: jburkardt-m