text:
function gamsWrite(prob,filename)
%gamsWrite  Write a GAMS Model from MATLAB data
%
% gamsWrite(prob,filename) writes the OPTI problem to a GAMS model (.gms)
% and optionally a GDX file if SCIP is not used.
%
% You may specify a full path to the file, or if you specify a filename
% only, it will be written to the current directory.

%   Copyright (C) 2013 Jonathan Currie (I2C2)

ext = 'gms';

%Quick Checks
if(~ischar(filename))
    error('Filename must be a char array!');
end
if(~isempty(prob.sdcone))
    error('Semidefinite problems cannot be written by this interface');
end
% Make sure we have file extension for GAMS model
if(isempty(strfind(lower(filename),['.' ext])))
    gdxname = filename;
    filename = [filename '.' ext];
else
    % Remove file extension for GDX file
    ind = strfind(lower(filename),['.' ext]);
    gdxname = filename(1:ind-1);
end

% Check for full filepath
if(isempty(strfind(filename,':')))
    filename = [cd filesep filename];
    gdxname = [cd filesep gdxname];
end

%If we have SCIP installed, use that (much easier, and does nonlinear models)
if(checkSolver('scip',0))
    opts = optiset('solver','scip','warnings','none','display','off','solverOpts',scipset('gamsfile',filename));
    if(isstruct(prob.int) && isfield(prob.int,'str'))
        prob.int = prob.int.str;
    end
    if(isempty(prob.x0) || all(isnan(prob.x0))), prob.x0 = randn(prob.sizes.ndec,1); end
    Opt = opti(prob,opts);
    try
        solve(Opt);
    catch ME
        if(strcmpi(ME.identifier,'OPTI:SCIPVAR'))
            throwAsCaller(MException('OPTI:SCIPVAR','AN ERROR OCCURED WHILE PRE-PROCESSING YOUR MATLAB FUNCTION(S) - PLEASE EXAMINE THE SCIP ERROR BELOW:\n\n%s',ME.message));
        else
            throw(ME);
        end
    end
    return;
end

%Check not nonlinear
if(~isempty(prob.fun) && isa(prob.fun,'function_handle'))
    error('The OPTI GAMS Writer only works for problems up to MIQCQP complexity.\n\nFor general nonlinear problems you must have %s installed to write nonlinear GAMS problems','SCIP');
end

%Check GAMS routines are available
if(~exist(['rgdx.' mexext],'file') || ~exist(['wgdx.' mexext],'file'))
    error('Could not find the required GAMS routines on your PC [wgdx,rgdx]!\n\n%s',...
          'Please download a current GAMS distribution from "http://www.gams.com/download/", then add the GAMS install folder to the MATLAB path.');
end

%Check for SOS constraints (not yet implemented)
if(~isempty(prob.sos))
    optiwarn('OPTI:GamsWriteSOS','SOS constraints are not currently supported by this interface, so they will be skipped!');
end

%Determine number of decision variables
if(isfield(prob,'sizes'))
    ndec = prob.sizes.ndec;
elseif(~isempty(prob.f))
    ndec = numel(prob.f);
elseif(~isempty(prob.x0))
    ndec = numel(prob.x0);
elseif(~isempty(prob.H))
    ndec = size(prob.H,1);
else
    error('Could not determine the number of decision variables!');
end

%Convert linear constraints to row format
if(~isempty(prob.b) || ~isempty(prob.beq))
    [prob.A,prob.rl,prob.ru] = gen2row(prob.A,prob.b,prob.Aeq,prob.beq);
end

%Convert integer constraints to string
if(~isempty(prob.int) && ~isstruct(prob.int))
    if(ischar(prob.int))
        str = prob.int;
        prob.int = struct('str',str);
    elseif(isnumeric(prob.int))
        str = repmat('C',1,ndec);
        str(prob.int) = 'I';
        prob.int = struct('str',str);
    else
        error('Unknown integer format');
    end
elseif(isempty(prob.int))
    prob.int.str = [];
end
prob.int.str = upper(prob.int.str);
nint = sum(prob.int.str == 'I');
nbin = sum(prob.int.str == 'B');

%Starting writing GAMS Model File
fid = fopen(filename,'w');
if(fid < 0), error('Error creating GAMS Model: %s',filename); end

name = regexprep(prob.Name,' ','_');

fprintf(fid,'* OPTI Generated GAMS Model (%s)\n',datestr(now));
fprintf(fid,'*   Model Name             : %s\n',name);
if(isfield(prob,'type') && ~isempty(prob.type))
    fprintf(fid,'*   Model Classification   : %s\n',prob.type);
end
fprintf(fid,'*   Variables              : %d (%d Integer, %d Binary)\n',ndec,nint,nbin);
if(~isempty(prob.rl)),  fprintf(fid,'*   Linear Constraints     : %d\n',numel(prob.rl)); end
if(~isempty(prob.qrl)), fprintf(fid,'*   Quadratic Constraints  : %d\n',numel(prob.qrl)); end
if(~isempty(prob.cl)),  fprintf(fid,'*   Nonlinear Constraints  : %d\n',numel(prob.cl)); end
fprintf(fid,'\n');

try
    model = 'lp';
    %Check for integer and binary constraints
    if(~all(prob.int.str == 'C'))
        int = struct('name','int','type','set','dim',1,'uels',{{find(prob.int.str ~= 'C')}},'form','full','ts','integer and binary indices');
        if(any(prob.int.str == 'B'))
            %Force binary bounds
            ind = (prob.int.str == 'B')';
            if(isempty(prob.lb)), prob.lb = -Inf(ndec,1); end
            if(isempty(prob.ub)), prob.ub = Inf(ndec,1); end
            %Set bin indicies which are infinite or bounds > 1 or < 0 to 1/0
            prob.lb((isinf(prob.lb) | (prob.lb < 0)) & ind) = 0;
            prob.ub((isinf(prob.ub) | (prob.ub > 1)) & ind) = 1;
        end
    else
        int = [];
    end
    %Check for bounds
    if(~isempty(prob.lb))
        lb = struct('name','lb','type','parameter','dim',1,'val',full(prob.lb),'form','full','ts','decision variable lower bound');
    else
        lb = [];
    end
    if(~isempty(prob.ub))
        ub = struct('name','ub','type','parameter','dim',1,'val',full(prob.ub),'form','full','ts','decision variable upper bound');
    else
        ub = [];
    end
    %Check for initial guess
    if(~isempty(prob.x0) && all(~isnan(prob.x0)))
        x0 = struct('name','x0','type','parameter','dim',1,'val',full(prob.x0),'form','full','ts','decision variable initial guess');
    else
        x0 = [];
    end
    %Check for linear objective
    if(~isempty(prob.f) && isnumeric(prob.f))
        f = struct('name','f','type','parameter','dim',1,'val',full(prob.f),'form','full','ts','linear objective vector');
    else
        f = [];
    end
    %Check for quadratic objective
    if(~isempty(prob.H) && isnumeric(prob.H))
        %We require all elements of H, see if we have a saved copy before
        %triu/tril
        if(isfield(prob,'save') && isfield(prob.save,'H') && ~isempty(prob.save.H) && isnumeric(prob.save.H))
            H = prob.save.H;
        else
            H = prob.H;
        end
        [i,j,v] = find(sparse(H));
        H = struct('name','H','type','parameter','dim',2,'val',[i j v],'form','sparse','ts','quadratic objective matrix');
        model = 'qp';
    else
        H = [];
    end
    %Check for linear constraints
    if(~isempty(prob.A))
        [i,j,v] = find(sparse(prob.A));
        A = struct('name','A','type','parameter','dim',2,'val',[i j v],'form','sparse','ts','linear constraint matrix');
        rl = struct('name','rl','type','parameter','dim',1,'val',full(prob.rl),'form','full','ts','linear constraint lower bound');
        ru = struct('name','ru','type','parameter','dim',1,'val',full(prob.ru),'form','full','ts','linear constraint upper bound');
    else
        A = []; rl = []; ru = [];
    end
    %Check for quadratic constraints
    if(~isempty(prob.Q))
        %Concatenate each Q matrix into a 3D array
        if(iscell(prob.Q))
            i = []; j = []; k = []; v = [];
            for n = 1:length(prob.Q) %could preallocate...
                [I,J,V] = find(sparse(prob.Q{n}));
                i = [i;I]; %#ok<AGROW>
                j = [j;J]; %#ok<AGROW>
                k = [k;n*ones(size(I))]; %#ok<AGROW>
                v = [v;V]; %#ok<AGROW>
            end
        else
            [i,j,v] = find(sparse(prob.Q));
            k = ones(size(i));
        end
        Q = struct('name','Q','type','parameter','dim',3,'val',[i j k v],'form','sparse','ts','quadratic constraint matrix');
        [i,j,v] = find(sparse(prob.l));
        l = struct('name','l','type','parameter','dim',2,'val',[i j v],'form','sparse','ts','quadratic constraint linear matrix');
        qrl = struct('name','qrl','type','parameter','dim',1,'val',full(prob.qrl),'form','full','ts','quadratic constraint lower bound');
        qru = struct('name','qru','type','parameter','dim',1,'val',full(prob.qru),'form','full','ts','quadratic constraint upper bound');
        model = 'qcqp';
    else
        Q = []; l = []; qrl = []; qru = [];
    end

    %Write GAMS Data
    writeGDX(gdxname,H,f,lb,ub,x0,int,A,rl,ru,Q,l,qrl,qru);

    %Write set definitions
    fprintf(fid,'Sets\n');
    fprintf(fid,' n /1*%d/\n',ndec);
    if(~isempty(int)), fprintf(fid,' ni(n)\n'); end
    if(~isempty(prob.rl)), fprintf(fid,' mlin /1*%d/\n',numel(prob.rl)); end
    if(~isempty(prob.qrl)), fprintf(fid,' mquad /1*%d/\n',numel(prob.qrl)); end
    if(~isempty(prob.cl)), fprintf(fid,' mnlin /1*%d/\n',numel(prob.cl)); end
    fprintf(fid,';\n');

    %Alias for QP, QC
    if(~isempty(H)), fprintf(fid,'Alias (n,n1);\n\n'); end

    %Write parameter definitions
    fprintf(fid,'Parameters\n');
    if(~isempty(H)), fprintf(fid,' H(n,n1)\n'); end
    if(~isempty(f)), fprintf(fid,' f(n)\n'); end
    if(~isempty(lb) && ~isempty(ub))
        fprintf(fid,' lb(n), ub(n)\n');
    elseif(~isempty(lb))
        fprintf(fid,' lb(n)\n');
    elseif(~isempty(ub))
        fprintf(fid,' ub(n)\n');
    end
    if(~isempty(x0)), fprintf(fid,' x0(n)\n'); end
    if(~isempty(A)), fprintf(fid,' A(mlin,n), rl(mlin), ru(mlin)\n'); end
    if(~isempty(Q)), fprintf(fid,' Q(n,n1,mquad), l(n,mquad), qrl(mquad), qru(mquad)\n'); end
    fprintf(fid,';\n');

    %Enter Variables to Import
    fprintf(fid,'* Load Parameters from GDX\n');
    fprintf(fid,'$gdxin %s.gdx\n',gdxname);
    if(~isempty(H)), fprintf(fid,'$load H\n'); end
    if(~isempty(f)), fprintf(fid,'$load f\n'); end
    if(~isempty(lb) && ~isempty(ub))
        fprintf(fid,'$load lb,ub\n');
    elseif(~isempty(lb))
        fprintf(fid,'$load lb\n');
    elseif(~isempty(ub))
        fprintf(fid,'$load ub\n');
    end
    if(~isempty(int)), fprintf(fid,'$load ni = int\n'); end
    if(~isempty(x0)), fprintf(fid,'$load x0\n'); end
    if(~isempty(A)), fprintf(fid,'$load A,rl,ru\n'); end
    if(~isempty(Q)), fprintf(fid,'$load Q,l,qrl,qru\n'); end
    fprintf(fid,'$gdxin\n\n');

    %Complete GAMS file
    %Variable Definition
    fprintf(fid, 'Variables\n x(n)\n');
    fprintf(fid, ' obj;\n\n');

    %Integer and binary def
    if(~isempty(int))
        fprintf(fid,'Integer Variables\n x;\n');
        fprintf(fid,' x.prior(n) = inf;\n');
        fprintf(fid,' x.prior(ni) = 1;\n\n');
    end

    %Equations
    fprintf(fid, 'Equations\n');
    if(~isempty(A))
        fprintf(fid, ' linconLB(mlin)\n');
        fprintf(fid, ' linconUB(mlin)\n');
    end
    if(~isempty(Q))
        fprintf(fid, ' quadconLB(mquad)\n');
        fprintf(fid, ' quadconUB(mquad)\n');
    end
    fprintf(fid, ' cost;\n\n');

    %Objective
    fprintf(fid, '* Cost Function\n');
    switch(model)
        case 'lp'
            fprintf(fid, ' cost.. obj =e= sum(n, f(n)*x(n));\n\n');
        case {'qp','qcqp'}
            fprintf(fid, ' cost.. obj =e= 0.5*sum(n,x(n)*sum(n1,H(n,n1)*x(n1))) + sum(n,f(n)*x(n));\n\n');
        otherwise
            error('Problem type ''%s'' not yet supported',prob.type);
    end

    %Bounds
    if(~isempty(lb) || ~isempty(ub))
        fprintf(fid, '* Decision Variable Bounds\n');
        if(~isempty(lb)), fprintf(fid, ' x.lo(n) = lb(n);\n'); end
        if(~isempty(ub)), fprintf(fid, ' x.up(n) = ub(n);\n'); end
        fprintf(fid,'\n');
    end

    %Initial Guess
    if(~isempty(x0))
        fprintf(fid, '* Initial Solution Guess\n');
        fprintf(fid, ' x.l(n) = x0(n);\n\n');
    end

    %Linear Constraints
    if(~isempty(A))
        fprintf(fid,'* Linear Constraints\n');
        fprintf(fid,' linconLB(mlin).. rl(mlin) =l= sum(n, A(mlin,n)*x(n));\n');
        fprintf(fid,' linconUB(mlin).. sum(n, A(mlin,n)*x(n)) =l= ru(mlin);\n\n');
    end

    %Quadratic Constraints
    if(~isempty(Q))
        fprintf(fid,'* Quadratic Constraints\n');
        fprintf(fid,' quadconLB(mquad).. qrl(mquad) =l= sum(n,x(n)*sum(n1,Q(n,n1,mquad)*x(n1))) + sum(n,l(n,mquad)*x(n));\n');
        fprintf(fid,' quadconUB(mquad).. sum(n,x(n)*sum(n1,Q(n,n1,mquad)*x(n1))) + sum(n,l(n,mquad)*x(n)) =l= qru(mquad);\n\n');
    end

    %Write model declaration
    fprintf(fid,'\nModel %s / all /;\n\n',name);

    %Determine which GAMS solver type to use
    switch(model)
        case 'lp'
            if(isempty(int))
                fprintf(fid,'Solve %s using %s minimizing obj;\n',name,'LP');
            else
                fprintf(fid,'Solve %s using %s minimizing obj;\n',name,'MIP');
            end
        case {'qp','qcqp'}
            if(isempty(int))
                fprintf(fid,'Solve %s using %s minimizing obj;\n',name,'QCP');
            else
                fprintf(fid,'Solve %s using %s minimizing obj;\n',name,'MIQCP');
            end
    end

catch ME
    fclose(fid);
    rethrow(ME);
end

fclose(fid);

function writeGDX(name,H,f,lb,ub,x0,int,A,rl,ru,Q,l,qrl,qru)

args = {};
if(~isempty(H)),    args = [args H]; end
if(~isempty(f)),    args = [args f]; end
if(~isempty(lb)),   args = [args lb]; end
if(~isempty(ub)),   args = [args ub]; end
if(~isempty(x0)),   args = [args x0]; end
if(~isempty(int)),  args = [args int]; end
if(~isempty(A)),    args = [args A]; end
if(~isempty(rl)),   args = [args rl]; end
if(~isempty(ru)),   args = [args ru]; end
if(~isempty(Q)),    args = [args Q]; end
if(~isempty(l)),    args = [args l]; end
if(~isempty(qrl)),  args = [args qrl]; end
if(~isempty(qru)),  args = [args qru]; end

if(~isempty(args))
    wgdx(name,args{:});
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function GCO_SetLabeling(Handle,Labeling)
% GCO_SetLabeling     Sets the current labeling
%     GCO_SetLabeling(Handle,Labeling) sets the entire labeling.

GCO_LoadLib();
if (isnumeric(Labeling))
    NumSites = gco_matlab('gco_getnumsites',Handle);
    NumLabels = gco_matlab('gco_getnumlabels',Handle);
    if (length(Labeling) ~= NumSites)
        error('Labeling must be of length NumSites');
    end
    if (~isa(Labeling,'int32'))
        if (any(floor(Labeling) ~= Labeling))
            error('Labeling was not integer valued');
        end
        Labeling = int32(Labeling);
    end
    if (min(Labeling) < 1 || max(Labeling) > NumLabels)
        error('Label must be in range 1..NumLabels');
    end
    gco_matlab('gco_setlabeling',Handle,Labeling);
end

author: andrewssobral

repo: lrslibrary


text:
function [Ixx, Iyy, Ixy] = polygonSecondAreaMoments(poly)
%POLYGONSECONDAREAMOMENTS Compute second-order area moments of a polygon.
%
%   [IXX, IYY, IXY] = polygonSecondAreaMoments(POLY)
%   Compute the second-order inertia moments of a polygon. The polygon is
%   specified by the N-by-2 list of vertex coordinates.
%
%   Example
%   polygonSecondAreaMoments
%
%   References
%   * http://paulbourke.net/geometry/polygonmesh/
%   * https://en.wikipedia.org/wiki/Second_moment_of_area
%
%   See also
%     polygons2d, polygonEquivalentEllipse, polygonArea, polygonCentroid

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2017-09-08, using Matlab 9.1.0.441655 (R2016b)
% Copyright 2017-2022 INRA - Cepia Software Platform

poly = parsePolygon(poly, 'repetition');

px = poly(:,1);
py = poly(:,2);

% vertex indices
N = length(px);
iNext = [2:N 1];

% compute twice signed area of each triangle
common = px .* py(iNext) - px(iNext) .* py;

% compute each term
Ixx = sum( (py.^2 + py .* py(iNext) + py(iNext).^2) .* common) / 12;
Iyy = sum( (px.^2 + px .* px(iNext) + px(iNext).^2) .* common) / 12;
Ixy = sum( ...
    (px .* py(iNext) + 2 * px .* py + 2 * px(iNext) .* py(iNext) ...
    + px(iNext) .* py ) .* common) / 24;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function f=evalContFrac(bFunc,aFunc,tinyVal,maxIter)
%%EVALCONTFRAC Evaluate a continued fraction of the form
%              f=b(1)+a(1)/(b(2)+a(2)/(b(3)+a(3)/(b(4)+a(4)/...
%              Given either matrices for a finite number of a and b terms
%              of function for the a and b values such that the required
%              number of terms can be requested.
%
%INPUTS: bFunc This can either be a function such that bFunc(i) provides
%              b(i) in the continued fraction, or this can be a maxNX1 or
%              1XmaxN vector of explicit b values.
%        aFunc This has to be of the same type as bFunc, either a function
%              handle or a vector to provide the a values. If a vector,
%              then it is 1X(maxN-1) or (maxN-1)X1 as it takes one fewer
%              call to the a terms than the b terms to evaluate the
%              continued fraction.
%      tinyVal An optional parameter that is less than typical values of
%              eps(b), but is not zero. This parameter is used to mitigate
%              divide-by zero errors. The default if this parameter is
%              omitted or an empty matrix is passed if 2e-100.
%      maxIter An optional parameter considering the maximum number of
%              terms to evaluate in the continued fraction. The default
%              value if omitted is 1000. This parameter is only used if
%              aFunc and bFunc are functions, not arrays. If they are
%              arrays, then the maximum number of terms used is determined
%              by the array length.
%
%OUTPUTS: f The value of the continued fraction. The smallest possible
%           value that this can take if b(0)=0 and all a and b are
%           non-negative is eps, due to how the function is initialized to
%           be numerically robust.
%
%The algorithm is the modified Lentz's method. A good description of the
%method, demonstrating its origin is given in Chapter II and Appendix C of
%[1]. The origin of the algorithm is [2]. However, the algorithm is harder
%to understand outside of the context of Bessel functions when consulting
%[2].
%
%An example of a continued fraction is an approximation to pi:
% bFunc=@(i)((i==1)*0+(i~=1)*(2*(i-1)-1));
% aFunc=@(i)(4*(i==1)+(i-1)^2*(i~=1));
% piApprox=evalContFrac(bFunc,aFunc)
%
%REFERENCES:
%[1] C.-Y. Shih, "An investigation of a stable algorithm for the evaluation
%    of fractional order Bessel functions," Master's thesis, Emporia
%    State University, Emporia, KS, Nov. 1993. [Online]. Available:
%    https://esirc.emporia.edu/bitstream/handle/123456789/1773/Shih%201993.pdf
%[2] W. J. Lentz, "Generating Bessel functions in Mie scattering
%    calculations using continued fractions," Applied optics, vol. 15, no.
%    3, pp. 668-671, Mar. 1976.
%
%October 2014 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%The floating point precision.
epsVal=eps();

if(nargin<3)
    tinyVal=2^(-100);
end

if(isa(bFunc,'double'))
    maxN=length(bFunc);
    maxIter=maxN-1;
elseif(nargin<4)
    maxIter=1000;
end

b0=bFunc(1);
if(b0==0)
    f=tinyVal;
else
    f=b0;
end

CPrev=f;
DPrev=0;

j=1;
while(j<maxIter)
    bj=bFunc(j+1);
    aj=aFunc(j);

    DCur=bj+aj*DPrev;
    if(DCur==0)
        DCur=tinyVal;
    end
    DCur=1/DCur;
    DPrev=DCur;

    CCur=bj+aj/CPrev;
    if(CCur==0)
        CCur=tinyVal;
    end
    CPrev=CCur;

    Delta=CCur*DCur;
    f=f*Delta;

    if(all(abs(Delta-1)<epsVal))
        break;
    end
    j=j+1;
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
classdef Grid1d < handle
    % Grid has all the information related to the staggered grid in a single
    % Cartesian axis.  It does not have physical quantities dependent on
    % frequencies, e.g., omega, eps, mu, and PML s-factors.

    properties (SetAccess = immutable)
		axis  % axis of this grid; one of Axis.x, Axis.y, Axis.z
		unit  % instance of PhysUnit
		unitvalue  %  unit value of length
        l  % {locations for primary vertices, locations for dual vertices}
        dl  % {dl at primary vertex, dl at dual vertex} == { diff( l{GT.dual} ), diff( l{GT.prim} ) }
        bc  % instance of BC
        N  % number of grid cells
		L  % length of the grid (domain)
        Npml  % [Npml at (-) end, Npml at (+) end]
		lpml  % [location of PML interface at (-) end, location of PML interface at (+) end]
		Lpml  % [thickness of PML interface at (-) end, thickness of PML interface at (+) end]
		center  % center of grid without PML
	end

	properties (Access = private)
		lghost  % [location of ghost primary vertex (after the last one), location of ghost dual vertex (before the first one)]
		ldual_ext  % location of dual vertex after the last one
	end

	properties (Dependent, SetAccess = immutable)
		lg  % {l_prim with ghost, l_dual with ghost}: l with ghost vertex (beyond boundary)
		lall  % {l_prim with ghost, l_dual with extra vertices}: l with vertices to interpolate fields at corners of simulation domain
		bound  % [lall_prim(1), lall_prim(end)]
	end

	properties (SetAccess = private)
		kBloch
	end

    methods
        function this = Grid1d(axis, unit, lprim_array, Npml_array, bc)
			chkarg(istypesizeof(axis, 'Axis'), '"axis" should be instance of Axis.');
			this.axis = axis;

			chkarg(istypesizeof(unit, 'PhysUnit'), '"unit" should be instance of PhysUnit.');
			this.unit = unit;
			this.unitvalue = unit.value(PhysQ.L);

			chkarg(istypesizeof(Npml_array, 'int', [1, Sign.count]), ...
				'"Npml" should be length-%d row vector with integer elements.', Sign.count);
			this.Npml = Npml_array;

			chkarg(istypesizeof(bc, 'BC'), '"bc" should be instance of BC.');
			this.bc = bc;

			chkarg(istypesizeof(lprim_array, 'real', [1 0]), ...
				'"lprim_array" should be row vector with real elements.');

			% Set N and L.
			lprim = lprim_array;
			this.N = length(lprim)-1;  % # of grid cells in the axis
			this.L = lprim(end) - lprim(1);

			% Set loc and dl.
			ldual = NaN(1, this.N+1);
			ldual(2:(this.N+1)) = (lprim(1:end-1) + lprim(2:end)) / 2;

			if this.bc == BC.p
				ldual(1) = ldual(end) - (lprim(end)-lprim(1));  % lprim(end) - lprim(1) == ldual(end) - ldual(1)
				this.ldual_ext = ldual(2) + (lprim(end)-lprim(1));  % lprim(end) - lprim(1) = ldual_ext - ldual(2)
			else
				ldual(1) = lprim(1) - (ldual(2)-lprim(1));  % lprim(1) - ldual(1) == ldual(2) - lprim(1)
				this.ldual_ext = lprim(end) + (lprim(end) - ldual(end));  % ldual_ext - lprim(end) = lprim(end) - ldual(end)
			end

			this.l = {lprim(1:end-1), ldual(2:end)};
			this.lghost = [lprim(end), ldual(1)];
			this.dl = {diff(ldual), diff(lprim)};  % not {diff(lprim), diff(ldual)}

			% Set lpml, Lpml, and center.
			this.lpml = [lprim(1 + this.Npml(Sign.n)), lprim(end - this.Npml(Sign.p))];
			this.Lpml = [this.lpml(Sign.n) - lprim(1), lprim(end) - this.lpml(Sign.p)];
			this.center = mean(this.lpml);

			% Initialize kBloch
			this.kBloch = 0;
		end

		function lg = get.lg(this)
			lg = {[this.l{GT.prim}, this.lghost(GT.prim)], [this.lghost(GT.dual), this.l{GT.dual}]};
		end

		function lall = get.lall(this)
			lall = {this.lg{GT.prim}, [this.lg{GT.dual}, this.ldual_ext]};
		end

		function bound = get.bound(this)
			bound = this.lall{GT.prim}([1 end]);
		end

		function set_kBloch(this, blochSrc)
			chkarg(istypesizeof(blochSrc, 'WithBloch'), '"blochSrc" should be instance of WithBloch.');
			this.kBloch = blochSrc.kBloch(this.axis);
		end

		function truth = contains(this, l)
			% This function can handle "l" as an array.
			chkarg(istypeof(l, 'real'), '"l" should be array with real elements.');
			truth = (l >= this.lall{GT.prim}(1)) & (l <= this.lall{GT.prim}(end));  % &: elementwise AND operator
		end

		function bound_plot = bound_plot(this, withpml)
			if withpml
				bound_plot = this.bound;
			else
				bound_plot = this.lpml;
			end
		end

		function lplot = lplot(this, g, withinterp, withpml)
			% Return the locations where data are evaluated for plotting.  If
			% the data do not include the boundaries of the simulation domain (or
			% the PML interfaces for "withpml == false"), the boundary points
			% are added to ensure that the plot is drawn from boundary to
			% boundary.
			chkarg(istypesizeof(g, 'GT') , '"g" should be instance of GT');
			chkarg(istypesizeof(withinterp, 'logical'), '"withinterp" should be logical.');
			chkarg(istypesizeof(withpml, 'logical'), '"withpml" should be logical.');

			if g == GT.prim
				lplot = this.lall{g};
			else  % g == GT.dual
				lplot = this.l{g};
			end

			if ~withpml
				lplot = lplot(1+this.Npml(Sign.n):end-this.Npml(Sign.p));
			end

			if g == GT.dual && withinterp
				lbound = this.bound_plot(withpml);
				lplot = [lbound(1), lplot, lbound(end)];
			end
		end

		function lvoxelbound = lvoxelbound(this, g, withpml)
			% Return the locations of boundaries of voxels drawn.  For data
			% at primary grid points, the voxel centers are within the
			% simulation domain including the boundary.  For data at dual
			% grid points, the voxel centers are within the simulation
			% domain excluding the boundary.
			chkarg(istypesizeof(g, 'GT') , '"g" should be instance of GT');
			chkarg(istypesizeof(withpml, 'logical'), '"withpml" should be logical.');

			lvoxelbound = this.lall{alter(g)};
			if ~withpml
				lvoxelbound = lvoxelbound(1+this.Npml(Sign.n):end-this.Npml(Sign.p));
			end

		end
	end
end

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function [m2] = in22m2(in2)
% Convert area from square inches to square meters.
% Chad A. Greene 2012
m2 = in2*0.00064516;

author: Sable

repo: mcbench-benchmarks


text:
function [IND INT]=imsphere_IND_INT(r,v_low,res,size_IM,Xp,Yp,Zp)

% function [I J K]=imsphere_IND_INT(r,v_low,res,size_IM,Xp,Yp,Zp)
% ------------------------------------------------------------------------
%
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 07/08/2008
% ------------------------------------------------------------------------

%% Setting up voxel dimensions and voxel radius dimensions
v_high=v_low/res;
r_v_high=r/v_high(1);
r_v_low=r/v_low(1);

%% Calculating image coordinate of point in high resolution image
[Ip_high Jp_high Kp_high]=cart2im(Xp,Yp,Zp,v_high);

%% Finding the index of the voxel in which point is found or is closest to
Ir_high=round(Ip_high); Jr_high=round(Jp_high); Kr_high=round(Kp_high);

%% Calculating shift from the centre of this voxel
Ic_high=Ir_high-Ip_high; Jc_high=Jr_high-Jp_high; Kc_high=Kr_high-Kp_high;

%% Finding indices of points in a sphere at specified location
[I_sph_high J_sph_high K_sph_high]=sphere_index(r_v_high,Ic_high,Jc_high,Kc_high);

%Calculating image coordinate of point in low resolution image
[Ip_low Jp_low Kp_low]=cart2im(Xp,Yp,Zp,v_low);
%Finding the index of the voxel in which point is found
Ir_low=round(Ip_low); Jr_low=round(Jp_low); Kr_low=round(Kp_low);
%Calculating shift from the centre of this voxel
Ic_low=Ir_low-Ip_low; Jc_low=Jr_low-Jp_low; Kc_low=Kr_low-Kp_low;

%Calculating shift
I_low_high_shift=((Ir_low*res)-(res/2)+0.5)-Ir_high;
J_low_high_shift=((Jr_low*res)-(res/2)+0.5)-Jr_high;
K_low_high_shift=((Kr_low*res)-(res/2)+0.5)-Kr_high;

size_M=(res+(2*round(r_v_low+1)*res))*ones(1,3);
M=ones(size_M);
IJK_middle=round(size_M/2);

IJK_voxel_center(1)=IJK_middle(1)-I_low_high_shift;
IJK_voxel_center(2)=IJK_middle(2)-J_low_high_shift;
IJK_voxel_center(3)=IJK_middle(3)-K_low_high_shift;
%round_IJK_voxel_center=round(IJK_voxel_center);

I_sph_high=round(I_sph_high+IJK_voxel_center(1));
J_sph_high=round(J_sph_high+IJK_voxel_center(2));
K_sph_high=round(K_sph_high+IJK_voxel_center(3));
IND = sub2ind(size_M,I_sph_high,J_sph_high,K_sph_high);
M(IND)=0;

%disp(['High resolution sphere simulated using ', num2str(numel(IND)), ' voxels']);

%Calculating low resolution image
[M]=voxelate(M,res);

IND=find(M<1);
INT=M(IND);
[I,J,K]=ind2sub(size(M),IND);
IJK_middle=round(size(M)/2);
I=I-IJK_middle(1); J=J-IJK_middle(2); K=K-IJK_middle(3);

I=I+Ir_low;
J=J+Jr_low;
K=K+Kr_low;

[IND]=sub2ind(size_IM,I,J,K);

%% END

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function chu_test ( )

%*****************************************************************************80
%
%% CHU_TEST tests R4_CHU and R8_CHU.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 October 2011
%
%  Author:
%
%    John Burkardt
%
  addpath ( '../test_values' );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHU_TEST:\n' );
  fprintf ( 1, '  Test HYPERGEOMETRIC_U_VALUES, R4_CHU, R8_CHU.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '             A               B               X     CHU(A,B,X)\n' );
  fprintf ( 1, '                                                R4_CHU(A,B,X)' );
  fprintf ( 1, '         Diff\n' );
  fprintf ( 1, '                                                R8_CHU(A,B,X)' );
  fprintf ( 1, '         Diff\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, a, b, x, fx1 ] = hypergeometric_u_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = r4_chu ( single ( a ), single ( b ), single ( x ) );
    fx3 = r8_chu ( a, b, x );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %14.4f  %14.4f  %14.4g  %14.6g\n', a, b, x, fx1 );
    fprintf ( 1, '                                                  %14.6g  %14.6g\n', fx2, abs ( fx1 - fx2 ) );
    fprintf ( 1, '                                                  %14.6g  %14.6g\n', fx3, abs ( fx1 - fx3 ) );

  end

  rmpath ( '../test_values' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%DATASETCONV Convert to dataset if needed
%
%		A = DATASETCONV(A)
%
% If A is not a dataset it is converted to a dataset.
%
% SEE ALSO (PRTools Guide)
% DATASETS, DATASET

function a = datasetconv(a)

% This is just programmed like this for speed, as
% a = prdataset(a) will do the same but involves more checking

if ~isdataset(a)
	a = prdataset(a);
end

author: marianux

repo: ecg-kit

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


text:
function [r_corr, dat, sl_size] = searchlight_modcomp(mask1, rdm, varargin)
% Calculate the local pattern similarity between two pattern maps using
% the searchlight approach.
%
% :Usage:
% ::
%
%     [r_corr, dat, sl_size] = searchlight_correlation(mask1, mask2, [additional_inputs])
%
%
% :Inputs:
%
%   **mask1:**
%        pattern or activation maps 1
%
%   **mask2:**
%        pattern or activation mediaps 2
%
% :Optional inputs:
%
%   **'r':**
%        searchlight sphere radius (in voxel) (default: r = 3 voxels)
%
%   **'type':**
%        This calls corr.m, and can take 'type' option.
%
%        'Pearson' (default), 'Kendall', 'Spearman'.
%
% :Outputs:
%
%   **r_corr:**
%        Correlation between weights of two pattern masks
%
%   **dat:**
%        This contains a statistic_image object that contain
%        correlation values between weights of two pattern masks
%        (=r_corr; in .dat) and p values for the correlation values
%        (in .p).
%
%   **sl_size:**
%        The number of voxels within each searchlight. Based on this
%        number, you can detect searchlights on the edge (searchlights
%        with low sl_size should be on the edge of the brain.
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2017  Phil Kragel
%     Copyright (C) 2014  Wani Woo
%
%     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/>.
% ..
%
% :Examples:
% ::
%
%    mask1 = which('weights_NSF_grouppred_cvpcr.img');
%    mask2 = which('nonnoc_v6_109subjmap_mean.nii');
%
%    [r, dat] = searchlight_correlation(mask1, mask2, 'r', 5);

r = 4; % default radius (in voxel)
corr_type = 'Pearson';

% parsing varargin

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}
            % functional commands
            case 'type'
                corr_type = varargin{i+1};
            case 'r' % radius
                r = varargin{i+1};
        end
    end
end

dat1 = fmri_data(mask1); dat1 = remove_empty(dat1);

dat1 = remove_empty(dat1);
n = size(dat1.dat,1);

r_corr = NaN(n,1);
p_corr = NaN(n,1);
sl_size = zeros(n,1);

fprintf('\n Calculating corrleation for voxel                 ');
xyzlist=dat1.volInfo.xyzlist(~dat1.removed_voxels,:);
dat=dat1.dat;
for i = 1:n %(1):vox_to_run(10)
    fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%07d/%07d', i, n);
    searchlight_indx = searchlight_sphere_prep(xyzlist, i, r);
    [r_corr(i), p_corr(i)] = corr(pdist(dat(searchlight_indx,:)','correlation')',rdm', 'type', corr_type);
    sl_size(i) = sum(searchlight_indx);
end

dat = statistic_image;
eval(['dat.volInfo = dat' num2str(1) '.volInfo;']);
eval(['dat.removed_voxels = dat' num2str(1) '.removed_voxels;']);
dat.dat = r_corr;
dat.p = p_corr;

end

% ========== SUBFUNCTION ===========

function indx = searchlight_sphere_prep(xyz, i, r)
seed = xyz(i,:);
indx = sum([xyz(:,1)-seed(1) xyz(:,2)-seed(2) xyz(:,3)-seed(3)].^2, 2) <= r.^2;
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% Init Script
% Remarks:
%   1.  sa
% TODO:
% 	1.  ds
% Release Notes
% - 1.0.000     08/11/2016  Royi Avital
%   *   First release.
%

%% General Parameters
close('all');
clear();
clc();

% set(0, 'DefaultFigureWindowStyle', 'docked');
% defaultLoosInset = get(0, 'DefaultAxesLooseInset');
% set(0, 'DefaultAxesLooseInset', [0.05, 0.05, 0.05, 0.05]);

figPosSmall     = [100, 100, 400, 300];
figPosMedium    = [100, 100, 800, 600];
figPosLarge     = [100, 100, 960, 720];
figPosDefault   = [100, 100, 560, 420];

fontSizeTitle   = 14;
fontSizeAxis    = 12;
fontSizeString  = 12;

lineWidthThin   = 1;
lineWidthNormal = 3;
lineWidthThick  = 4;

markerSizeSmall     = 4;
markerSizeNormal    = 8;
markerSizeLarge     = 10;

randomNumberStream  = RandStream('mlfg6331_64', 'NormalTransform', 'Ziggurat');
subStreamNumber     = round(sum(clock()));
% subStreamNumber     = 162;
subStreamNumber     = 2122;
set(randomNumberStream, 'Substream', subStreamNumber);
RandStream.setGlobalStream(randomNumberStream);

FALSE   = 0;
TRUE    = 1;

OFF     = 0;
ON      = 1;

author: RoyiAvital

repo: StackExchangeCodes


text:
function str=getSummaryStr(obj, format)
% get summary of tests that were run and their output
%
% str=getSummaryStr(obj, format)
%
% Inputs:
%   obj                 MOxUnitTestReport instance
%   format              Output format: 'xml' or 'text'
%
% Output:
%   str                 text or XML string representation of the results of
%                       the tests
%

    format2converter=struct();
    format2converter.text=@report2summary_text;
    format2converter.xml=@report2summary_xml;

    f=format2converter.(format);
    str=f(obj);

function str=report2summary_xml(obj)
% return summary in XML format
    header=sprintf('<testsuite name="%s" time="%.3f" %s>',...
                    moxunit_util_escape_xml(getName(obj)),...
                    getDuration(obj),...
                    statistics_xml(obj));

    summary_text_cell=get_outcome_summary_cell(obj,'xml');
    body=moxunit_util_strjoin(summary_text_cell,'\n');

    footer='</testsuite>';
    str=sprintf('%s\n%s\n%s',header,body,footer);

function str=report2summary_text(obj)
% return summary in text format
    verbosity=obj.verbosity;

    hor_line=repmat('-',1,50);

    body='';
    if verbosity>=2 || ~wasSuccessful(obj)
        summary_text_cell=get_outcome_summary_cell(obj,'text');
        non_empty_msk=~cellfun(@isempty,summary_text_cell);
        if any(non_empty_msk)
            summary_text_cell_nonempty=summary_text_cell(non_empty_msk);
            body_non_pass=moxunit_util_strjoin(...
                                summary_text_cell_nonempty,'\n\n');
            body=sprintf('\n%s\n\n%s\n%s\n\n',hor_line,body_non_pass);
        end
    end

    if verbosity>=1
        footer=sprintf('\n%s\n\n%s',hor_line,statistics_text(obj));
    else
        footer='';
    end

    str=sprintf('%s%s',body,footer);

function [c,n_tests]=get_outcome_summary_cell(obj, format)
% helper function to get the summary from each test
    n_tests=countTestOutcomes(obj);
    c=cellfun(@(x)getSummaryStr(getTestOutcome(obj,x),format),...
                    num2cell(1:n_tests),...
                    'UniformOutput',false);

function str=statistics_xml(obj)
% return xml string indicating number of tests, and how many were failures,
% skips and errors
%
% Example outputs:
%       'tests="9"'
%       'tests="9" failures="2"'

    total_tests=sprintf('tests="%d"',countTestOutcomes(obj));

    labels_singular_plural={'failure','failures';...
                            'skipped','skipped';...
                            'error','errors'};
    n_labels=size(labels_singular_plural,1);

    label2count=getTestOutputStatistics(obj);
    non_success_tests_cell=cell(n_labels,1);
    for k=1:n_labels
        singular=labels_singular_plural{k,1};
        plural=labels_singular_plural{k,2};

        if isfield(label2count,singular)
            count=label2count.(singular);
        else
            count=0;
        end

        non_success_tests_cell{k}=sprintf('%s="%d"',plural,count);
    end

    str=sprintf('%s %s',total_tests,...
                    moxunit_util_strjoin(non_success_tests_cell,' '));

function str=statistics_text(obj)
% return text string indicating whether all tests passed or not, and
% number of failures, skips and errors.
%
% Example outputs:
%       'OK'
%       'FAILED (skipped=2)'

    if wasSuccessful(obj)
        str='OK';
    else
        str='FAILED';
    end

    label_counts=getTestOutputStatistics(obj);
    labels=fieldnames(label_counts);

    n_labels=numel(labels);
    if n_labels>0
        cell_stats=cellfun(@(x) sprintf('%s=%d',x, label_counts.(x)),...
                                    fieldnames(label_counts),...
                                    'UniformOutput',false);
        str=sprintf('%s (%s)',str,moxunit_util_strjoin(cell_stats,', '));
    end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function f = f_function ( x, y1, y2, Md, Nd, L )

%*****************************************************************************80
%
%% F_FUNCTION: RHS function from exact coefficient Q and the exact solution U.
%
%  Discussion:
%
%    The differential equation has the form:
%
%      ( q(x) u_x(x) )_x = f(x)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Parameters:
%
%    x = physical space point of length one
%
%    y1 = stochastic vector of length Md (from the expansion for u)
%
%    y2 = stochastic vector of length Nd (from the expansion for q)
%
%    Md = dimension of the probability space for u
%
%    Nd = dimension of the probability space for q
%
%    L = the correlation length of the random variables
%
  [n_nodes, dim] = size(x);
  f = zeros( n_nodes, 1 );

  for n=1:Nd

  	for m=1:Md

  	 f = f + qn(n) * um( m ) * pn(y2, n) * pm(y1, m) * ...
  	 ( -(pi/L)^2 * m^2 .* sin((m*pi*x)/L) .* cos((n*pi*x)/L) ...
  	 -n*m * (pi/L)^2 * cos((m*pi*x)/L) .* sin((n*pi*x)/L) );

 	end

  end

  return
end
function p_m = pm( y, m )

%*****************************************************************************80
%
%% PM
%
  p_m = y(m);

  return
end
function p_n = pn( y, n )

%*****************************************************************************80
%
%% PN
%
  p_n = y(n);

  return
end
function q_n = qn( n )

%*****************************************************************************80
%
%% QN
%
  q_n = 1;

  return
end
function u_m = um( m )

%*****************************************************************************80
%
%% UM
%
 u_m = 1;

 return
end


text:
function calpak_test345 ( )

%*****************************************************************************80
%
%% CALPAK_TEST345 tests MONTH_LENGTH_GREEK.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 March 2013
%
%  Author:
%
%    John Burkardt
%
  n_test = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST345\n' );
  fprintf ( 1, '  For the Greek calendar,\n' );
  fprintf ( 1, '  MONTH_LENGTH_GREEK returns month lengths.\n' );

  y_test(1) = 3;
  y_test(2) = 4;

  for i_test = 1 : n_test

    y = y_test(i_test);
    sy = y_to_s_greek ( y );
    months = year_length_months_greek ( y );
    days = year_length_greek ( y );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %d\n', y );
    fprintf ( 1, '  %s\n', sy );
    fprintf ( 1, '  Year length in months = %d\n', months );
    fprintf ( 1, '  Year length in days = %d\n', days );
    fprintf ( 1, '\n' );

    for m = 1 : months
      month_name = month_to_month_name_greek ( y, m );
      fprintf ( 1, '  %12s  %2d\n', month_name, month_length_greek ( y, m ) );
    end

  end

  return
end


text:
function nI= upSmpIm(I,new_imgSize,filtS)

if(~exist('filtS','var'))
  filtS=1;
end
if (filtS==2)
  filt=[1,4,6,4,1]/8;
end
if (filtS==1)
  filt=[1,2,1]/2;
end
if (filtS==0)
  filt=[1];
end

id=floor((new_imgSize(1)-size(I,1)*2+1)/2);
iu=ceil((new_imgSize(1)-size(I,1)*2+1)/2);
jd=floor((new_imgSize(2)-size(I,2)*2+1)/2);
ju=ceil((new_imgSize(2)-size(I,2)*2+1)/2);

nI=zeros(new_imgSize(1)+2*filtS,new_imgSize(2)+2*filtS,size(I,3),size(I, ...
                                                  4));
nI(id+filtS+1:2:end-iu-filtS,jd+filtS+1:2:end-ju-filtS,:,:)=I;
nI(id+filtS-1:-2:1,:,:,:)=repmat(nI(id+filtS+1,:,:,:),ceil((id+filtS-1)/2),1);

nI(end-iu-filtS+2:2:end,:,:,:)=repmat(nI(end-iu-filtS,:,:,:),ceil((iu+filtS-1)/2),1);
nI(:,jd+filtS-1:-2:1,:,:)=repmat(nI(:,jd+filtS+1,:,:),1,ceil((jd+filtS-1)/2));
nI(:,end-ju-filtS+2:2:end,:,:)=repmat(nI(:,end-ju-filtS,:,:),1,ceil((ju+filtS-1)/2));

for i=1:size(nI,3)
  for j=1:size(nI,4)
    nI(:,:,i,j)=conv2(filt,filt',nI(:,:,i,j),'same');
  end
end

nI=nI(filtS+1:end-filtS,filtS+1:end-filtS,:,:);

author: luanfujun

repo: deep-photo-styletransfer


text:
function [exp_avg,exp_avg_null,z]=do_homogeneity(x,mskFile,parcelFile,gmFile,MM)
% Input:
% x: fMRI data matrix,dimension: time x number of all gray matter voxels
% mskFile: subcortex atlas in NIFTI (*.nii)
% parcelFile: random parcellation in NIFTI (*.nii)
% MM: number of randomizations

% Output:
% exp_avg: parcellation homogeneity of the empirical parcellation
% exp_avg_null: parcellation homogeneity of random parcellations
% z: internal connectivity matrix between each pair of parcels

warning off

fprintf('Computing homogeneity for empirical data\n')
% Subcortex
[~,sub_msk]=read(mskFile);

% Gray matter
[~,gm_msk]=read(gmFile);
ind_msk=find(~~gm_msk);

expl=zeros(1,max(sub_msk(:)));
T=size(x,1);
x_roi_mean=zeros(T,max(sub_msk(:)));

for i=1:max(sub_msk(:)) % Loop over every parcel in the atlas

    ind_roi=find(sub_msk==i);

    % index of subcortical voxels into the whole gray matter mask
    ind_ind_roi=zeros(1,length(ind_roi));
    for ii=1:length(ind_roi)
        ind_ind_roi(ii)=find(ind_roi(ii)==ind_msk);
    end
    x_roi=x(:,ind_ind_roi); % x has already been demeaned

    % PCA on time series
    [~, ~, ~, ~, explained] = pca(x_roi,'Rows','complete');
    expl(i)=explained(1); % Variance explained by the 1st PC

    % internal matrix
    x_roi_mean(:,i)=mean(x_roi,2);
end

% Average across all parcels
exp_avg=mean(expl);
fprintf('Mean homogeneity=%0.2f\n',exp_avg)

clear expl

% internal connectivity matrix (optional)
x_roi_mean=detrend(x_roi_mean,'constant'); x_roi_mean=x_roi_mean./repmat(std(x_roi_mean),T,1);
c=corr(x_roi_mean);
z=atanh(c); clear c

% Randomization
fprintf('Loading precomputed random parcellations\n')
load(parcelFile,'parcels_random_all')
exp_avg_null=zeros(MM,1);

for m=1:MM

    parcels_random=parcels_random_all(:,:,:,m);

    expl=zeros(1,max(parcels_random(:)));

    for i=1:max(parcels_random(:))
        ind_roi=find(parcels_random==i);

        % index of subcortical voxels into whole gray matter mask
        ind_ind_roi=zeros(1,length(ind_roi));
        for ii=1:length(ind_roi)
            ind_ind_roi(ii)=find(ind_roi(ii)==ind_msk);
        end

        x_roi=x(:,ind_ind_roi); % x has already been demeaned

        % PCA on time series
        [~, ~, ~, ~, explained] = pca(x_roi,'Rows','complete');
        if ~isempty(explained)
            expl(i)=explained(1);
        else
            fprintf('Warning: region %d empty\n',i)
        end
    end

    exp_avg_null(m)=mean(expl);
    fprintf('randomisation %d of %d,mean homogeneity=%0.2f\n',m,MM,mean(expl))
    clear expl
end

author: yetianmed

repo: subcortex

save_path: github-repos/MATLAB/yetianmed-subcortex


text:
% THIS PROGRAM IS FOR IMPLEMENTATION OF DISCRETE TIME PROCESS UNSCENTED KALMAN FILTER
% FOR GAUSSIAN AND LINEAR STOCHASTIC DIFFERENCE EQUATION.
% (19 JULY 2005).
% UNSCENTED KALMAN FILTER (UKF) AT ITS BEST.
%(Under Nonlinear conditions,UNSCENTED KALMAN FILTER
% performs to a much better extent as compared to EXTENDED KALMAN FILTER).

clc;  close all;  clear all;

format long g;
Xo = [1; 0; 0; 0; 0; 0; 0; 0; 0; 0];
nx = length(Xo);
beta = 2;
elfa = 1*(10^-1);
lambda = (((elfa^2)*(nx)) - nx);
mn1 = mean(Xo);                                   %FIRST STEP
CV1 = (Xo-mn1)*(Xo-mn1)';
%CV1 = (Xo)*(Xo)';
PO = CV1;
FX = size(CV1);
VTt = [1 0 0 0 0 0 0 0 0 0]';
ADP1 = randn(1,200);
mn2 = mean(VTt);
NTt = [1 0 0 0 0 0 0 0 0 0]';
ADP2 = randn(1,200);
mn3 = mean(NTt);
Qo = (VTt-mn2)*(VTt-mn2)';
Ro = (NTt-mn3)*(NTt-mn3)';
%Qo = (VTt)*(VTt)';
%Ro = (NTt)*(NTt)';
FX = zeros(10,10);
PAO = [PO FX FX; FX Qo FX;FX FX Ro];
Xao = [mn1 0 0]';
lam = sqrt((elfa^2)*(nx));

Wmo = lambda/(nx+lambda);
Wmi = 1/(2*(nx+lambda));
Wmin = 1/(2*(nx+lambda));

Wco =  (lambda/(nx+lambda)) + (1-(elfa^2)+ beta);
Wci = Wmi;
Wcin = Wmin;

%CALCULATION OF SIGMA POINTS
[Sg11,Sg12,Sg13,Sg21,Sg22,Sg23,Sg31,Sg32,Sg33] = sigmacal(mn1,CV1,mn2,Qo,mn3,Ro,lam);

%TIME UPDATE EQUATIONS.
[Xinew1,Xinew2,Xinew3,Yinew1,Yinew2,Yinew3,Xbark,Ybark,Pbark] = TMUPDT(Wmo,Wmi,Wmin,Wco,Wci,Wcin,Sg11,Sg12,Sg13,Sg21,Sg22,Sg23,Sg31,Sg32,Sg33);

%MEASUREMENT UPDATE EQUATIONS.
[KGain,XNEW,PT1,PT2,PT3] = MSMTUPDT(Xinew1,Xinew2,Xinew3,Yinew1,Yinew2,Yinew3,Xbark,Ybark,Pbark,Wco,Wci,Wcin,Xo,VTt);

for ii = 1:1:100

    VTt = [ADP1(ii+1) 0 0 0 0 0 0 0 0 0]';
    NTt = [ADP2(ii+1) 0 0 0 0 0 0 0 0 0]';

    %CV1 = (XNEW)*(XNEW)';
    mn1 = mean(XNEW);
    mn2 = mean(VTt);
    mn3 = mean(NTt);
    % Qo = (VTt)*(VTt)';
    % Ro = (NTt)*(NTt)';

    CV1 = (XNEW-mn1)*(XNEW-mn1)';
    Qo = (VTt-mn2)*(VTt-mn2)';
    Ro = (NTt-mn3)*(NTt-mn3)';

    lam = sqrt((elfa^2)*(nx));

    [Sg11,Sg12,Sg13,Sg21,Sg22,Sg23,Sg31,Sg32,Sg33] = sigmacal(mn1,CV1,mn2,Qo,mn3,Ro,lam);

    [Xinew1,Xinew2,Xinew3,Yinew1,Yinew2,Yinew3,Xbark,Ybark,Pbark] = TMUPDT(Wmo,Wmi,Wmin,Wco,Wci,Wcin,Sg11,Sg12,Sg13,Sg21,Sg22,Sg23,Sg31,Sg32,Sg33);

    [KGain,XNEW,PT1,PT2,PT3,YNEW] = MSMTUPDT(Xinew1,Xinew2,Xinew3,Yinew1,Yinew2,Yinew3,Xbark,Ybark,Pbark,Wco,Wci,Wcin,XNEW,VTt);

    INP1(ii) = XNEW(1,1);
    INP2(ii) = YNEW(1,1);

end

T = 1:1:100;
figure(1);     subplot(211);    plot(real(INP1));   title('ORIGINAL SIGNAL');
subplot(212);  plot(real(INP2));   title('ESTIMATED SIGNAL (UNDER Nonlinear MODEL)');

figure(2);    plot(T,abs(INP1),T,abs(INP2));  title('Combined plot'); legend('original','estimated');


text:
function r8_to_r8_discrete_test ( )

%*****************************************************************************80
%
%% R8_TO_R8_DISCRETE_TEST tests R8_TO_R8_DISCRETE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  ndx = 19;
  rhi = 10.0;
  rlo = 1.0;
  test_num = 15;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_TO_R8_DISCRETE_TEST\n' );
  fprintf ( 1, '  R8_TO_R8_DISCRETE maps numbers to a discrete set\n' );
  fprintf ( 1, '  of equally spaced numbers in an interval.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of discrete values = %d\n', ndx );
  fprintf ( 1, '  Real interval: [%f, %f]\n', rlo, rhi );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  R   RD\n' );
  fprintf ( 1, '\n' );

  seed = 123456789;

  rlo2 = rlo - 2.0;
  rhi2 = rhi + 2.0;

  for test = 1 : test_num
    [ r, seed ] = r8_uniform_ab ( rlo2, rhi2, seed );
    rd = r8_to_r8_discrete ( r, rlo, rhi, ndx );
    fprintf ( 1, '  %14f  %14f\n', r, rd );
  end

  return
end


text:
function volume = estimateTissueVolume(curvature, thickness, surfaceArea)
%
% volume = estimateTissueVolume(curvature, thickness, surfaceArea)
%
% * curvature is an array of curvature estimates for each mesh triangle
% * thickness is the estimated tissue thickness
% * surfaceArea is the actual surface area of each triangle (on the plane)
%
% To estimate the volume, if the curvature is zero we can just multiply the
% single surface area of the triangle by the thickness of the tissue. But,
% if the local curvature is not zero, then the surface area at the white
% matter boundary differs from the surface area thickness millimeters away.
% So, we can 't just multiply a single area by the thickness.  Instead, we
% have to measure the volume by taking the difference between two spheres
% that define the inner and outer radius of the local curvature at the two
% boundaries of the gray matter.
%
% This routine measures the volume between these two surfaces rather than
% assuming there is a single surface area and multiplying by the thickness.
%
% Estimates the volume of tissue of a given thickness on a surface with the
% specified mean curvature. The algorithm simply estimates the tissue
% volume on a sphere of radius 1/curvature and radius + thickness.  It
% takes the local volume estimate as the surfaceArea/sphereSurfaceArea
% proportion of the total volume between these spheres.
%
% HISTORY: RFD (bob@white.stanford.edu) wrote it.
%

if(~exist('curvature','var') | isempty(curvature))
    help(mfilename);
    return;
end

curvature(curvature==0) = 0.0000001;
r = 1./curvature;

sphereSurfaceArea = 4*pi*r.^2;
innerV = (4*pi/3).*(r.^3);
outerV = (4*pi/3).*(r + thickness).^3;
totalTissueVolume = outerV - innerV;

volume = totalTissueVolume.*(surfaceArea./sphereSurfaceArea);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [modelComp] = createModelNewCompartment(model,OldComp,NewComp,NewCompName,LB,UB,RemoveExch)
% This function converts a two compartment metabolic model into a three compartment metabolic model model
%
% function [modelComp] = createModelNewCompartment(model,OldComp,NewComp,NewCompName,LB,UB,RemoveExch)
%
% INPUT
% model         Model structure
% OldComp       Name of the current (to be replicated) compartment to be
%               replicated and added
% NewCompName   Name of new compartment e.g., 'lumen'
% LB            Value for Lower bound on new exchange compartment (default:-1000)
% UB            Value for Upper bound on new exchange compartment (default:1000)
%
% NewComp       Abrreviation for the new compartment,m e.g., 'lu'
% RemoveExch    Option (if set to 1) to remove the old compartment
%
% OUTPUT
% modelComp     Model with three compartments
%
% Ines Thiele 2016-2017

warning off
if ~exist('LB','var')
    LB = -1000;
end

if ~exist('UB','var')
    UB = 1000;
end

if ~exist('RemoveExch','var')
    RemoveExch = 0;
end
% find exchange reactions
cnt = 1;
for t=1:length(model.rxns)
    if  strfind(model.rxns{t}, 'EX_')
        Ex_RxnsAll(cnt,1) =model.rxns(t); %make exchange reaction list
        cnt=cnt+1;
    elseif  strfind(model.rxns{t}, 'Ex_')
        Ex_RxnsAll(cnt,1) =model.rxns(t); %make exchange reaction list
        cnt=cnt+1;
    end
end

% duplicate all ExtRxns in a new compartment [i] for the New Compartment and
% for new exchange reactions
%a = printRxnFormulaOri(model, Ex_RxnsAll,[],[],[],[],false);
a = printRxnFormula(model,'rxnAbbrList',Ex_RxnsAll,'printFlag',0,'lineChangeFlag',0,'metNameFlag',0,'fid',0,'directionFlag',0);

aOri = a;
modelComp = model;

for i = 1 : length(Ex_RxnsAll)
    b = find(ismember(model.rxns,Ex_RxnsAll(i)));
    if ~isempty(strfind(a{i},strcat('[',OldComp,']'))) && length(find(model.S(:,b))) ==1
        % add new reactions to model
        [metaboliteList,stoichCoeffList] = parseRxnFormula(a{i});
        metaboliteListA = regexprep(metaboliteList,strcat('\[',OldComp,'\]'),strcat('\[',NewComp,'\]'));
        if UB == 0
            UBO = -1*LB;
            LBO = 0;
            revFlag = 0;
            RxnName = strcat(model.rxnNames{b},' (from ',NewCompName,'to ',strcat('[',OldComp,']'));
            modelComp = addReaction(modelComp,{strcat(Ex_RxnsAll{i},'_[',NewComp,']'),RxnName},[metaboliteListA metaboliteList],...
                [-1 1], revFlag,LBO,UBO,0,strcat('Transport'),'',[],[],0,0);
        else
            metaboliteList = [metaboliteList metaboliteListA];
            stoichCoeffList = [-1 1];
            if LB < 0
                revFlag = 1;
            else
                revFlag = 0;
            end
            RxnName = strcat(model.rxnNames{b},' (from ',strcat('[',OldComp,']'),' to ',NewCompName,')');
            %replaced addReactionOri - Ronan
            modelComp = addReaction(modelComp,{strcat(Ex_RxnsAll{i},'_[',NewComp,']'),RxnName},metaboliteList,stoichCoeffList,...
                revFlag,LB,UB,0,strcat('Transport'),'',[],[],0,0);
            %modelComp = addReaction(modelComp,{strcat(Ex_RxnsAll{i},'_[',NewComp,']'),RxnName},'metaboliteList',metaboliteList,'stoichCoeffList',stoichCoeffList,...
             %'lowerBound',LB,'upperBound',UB,'objectiveCoef',0,'subSystem',{'Transport'},'geneRule','','geneNameList',[],'systNameList',[],'checkDuplicate',0,'printLevel',0);
        end
    end
end
a=aOri;
if RemoveExch == 1
    for i = 1 : length(Ex_RxnsAll)
        % remove original reactions
        b = find(ismember(modelComp.rxns,Ex_RxnsAll(i)));
        if ~isempty(strfind(a{i},strcat('[',OldComp,']'))) && length(find(modelComp.S(:,b))) ==1
            if isfield(modelComp,'rxnGeneMat')
            modelComp = rmfield(modelComp,'rxnGeneMat');
            end
            modelComp = removeRxns(modelComp,Ex_RxnsAll{i});
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
u%%
%Read Example Data
file = 'example_data.xlsx';
num = xlsread(file,'narendra4');
%Inputs and outputs have to be matrices where columns=datapoints
%and rows=inputs
P = num(:,2).';
Y = num(:,3).';
Ptest = num(:,4).';
Ytest = num(:,5).';

%%
%Create NN

%create recurrent neural network with 1 input, 2 hidden layers with
%3 neurons each and 1 output
%the NN uses the input data at timestep t-1 and t-2
%The NN has a recurrent connection with delay of 1,2 and 3 timesteps from the output
% to the first layer (and no recurrent connection of the hidden layers)
nn = [1 3 3 1];
dIn = [1,2];
dIntern=[];
dOut=[1,2,3];
net = CreateNN(nn,dIn,dIntern,dOut); %alternative: net = CreateNN([1,3,3,1],[1,2],[],[1,2,3]);

%%
%Train with LM-Algorithm
% Train NN with training data P=input and Y=target
% Set maximum number of iterations k_max to 200
% Set termination condition for Error E_stop to 1e-3
% The Training will stop after 200 iterations or when the Error <=E_stop
netLM = train_LM(P,Y,net,200,1e-3);
%Calculate Output of trained net (LM) for training and Test Data
y_LM = NNOut(P,netLM);
ytest_LM = NNOut(Ptest,netLM);

%%
%Train with BFGS-Algorithm
% Train NN with training data P=input and Y=target
% Set maximum number of iterations k_max to 400
% Set termination condition for Error E_stop to 1e-3
% The Training will stop after 400 iterations or when the Error <=E_stop
% measure time dt
netBFGS = train_BFGS(P,Y,net,400,1e-3);
%Calculate Output of trained net (LM) for training and Test Data
y_BFGS = NNOut(P,netBFGS);
ytest_BFGS = NNOut(Ptest,netBFGS);

%%
%Plot Results
fig = figure();
set(fig, 'Units', 'normalized', 'Position', [0.2, 0.1, 0.6, 0.6])
axis tight

subplot(311)
set(gca,'FontSize',16)
plot(Y,'r:','LineWidth',2)
hold on
grid on
plot(y_LM,'b','LineWidth',2)
plot(y_BFGS,'g','LineWidth',2)
l1 = legend('Train Data','LM output','BFGS output','Location','northwest');
set(l1,'FontSize',14)

subplot(312)
set(gca,'FontSize',16)
plot(Ytest,'r:','LineWidth',2)
hold on
grid on
plot(ytest_LM,'b','LineWidth',2)
plot(ytest_BFGS,'g','LineWidth',2)
l2 = legend('Test Data','LM output','BFGS output','Location','northwest');
set(l2,'FontSize',14)

subplot(313)
set(gca,'FontSize',16)
plot(netLM.ErrorHistory,'b','LineWidth',2)
hold on
grid on
plot(netBFGS.ErrorHistory,'g','LineWidth',2)
ylim([0,5])
l3 = legend('LM Error','BFGS Error','Location','northeast');
set(l3,'FontSize',14)

author: yabata

repo: pyrenn

save_path: github-repos/MATLAB/yabata-pyrenn


text:
filename = 'ArchTriFine';%'CantileverSquareNewFine';%'LshapeTriFine';%'CantileverSquareNew';'ArchTriFine';%'CantileverSquare';%'ArchTriFine';%'CantileverSquare';%'Lshape';%'LshapeTriFine';%'CantileverSquareSmall';%'';%'Lshape';'CantileverSquare';'Lshape';%'CantileverSquare';%'LshapeFine';%'Bridge_quad_coarse';%'Arch_quad_coarse';%'BridgeCool_Quadrilateral_Bilinear_Structured_Coarse';%'Bridge';%'CantileverSquareSmall';
ptype = 'MACRO';
initial_case = 'given';
m1 = 0.0101;
m2 = 0.0101;
%cost = {'compliance'};
cost = {'stressNorm'};
%cost = {'stressNorm','compliance'};
%weights = [0.55,0.45];
weights = 1;
constraint = {'volumeConstraint'};
filterType = 'PDE';
constraint_case = 'EQUALITY';

Vfrac_initial = 0.3;
optimality_initial = 1e-5;
constr_initial = 1e-5;

Vfrac_final = 0.3;
optimality_final = 1e-5;
constr_final = 1e-5;

stressNormExponent_initial = 2;
stressNormExponent_final = 64;

optimizer = 'DualNestedInPrimal';
optimizerUnconstrained = 'PROJECTED GRADIENT';

designVariable = 'MicroParams';
ub = 0.989;
lb = 0.011;
homegenizedVariablesComputer = 'ByVademecum';
%
%vademecumFileName = 'SuperEllipseQMax';
%vademecumFileName = 'SuperEllipseQ2';
vademecumFileName = 'SuperEllipseQOptAnalytic';
%
% designVariable = 'Density';
% homegenizedVariablesComputer = 'ByInterpolation';
% method = 'SIMPALL';
% materialType = 'ISOTROPIC';

kfrac = 2;
nsteps = 50;

plotting = false;
printing = true;
monitoring = true;
monitoring_interval = 1;
maxiter = 1000;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function value = legendre_integral ( d )

%*****************************************************************************80
%
%% LEGENDRE_INTEGRAL returns the Legendre integral of a test function.
%
%  Discussion:
%
%    The same function, integrated over [-1,+1]^D, has an integral that
%    is a factor of 2^D times as large as this result.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 May 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer D, the spatial dimension.
%
%    Output, real VALUE, the value of the integral.
%
  value = ( 0.5 * erf ( 0.5 / sqrt ( 2.0 ) ) ) .^ d;

  return
end


text:
function io_loadworkspace(filename)
% Load workspace from a file.
% io_loadworkspace(filename)
%
% Brings up an appropriate file chooser dialog, if necessary.
%
% In:
%   Filename : workspace file name
%
% See also:
%   io_saveworkspace
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-11-19

% set default
if ~exist('filename','var') || isempty(filename)
    filename = 'bcilab:/resources/workspaces'; end
if ~ischar(filename)
    error('The file name must be given as a string.'); end

% translate to platform-specific path
filename = env_translatepath(filename);

if exist(filename,'dir')
    % load from given directory
    [filename,filepath] = uigetfile({'*.mat','Workspace files (*.mat)'},'Load workspace',filename);
    filename = [filepath filename];
elseif ~exist(filename,'file')
    % open file chooser from current path
    [filename,filepath] = uigetfile({'*.mat','Workspace files (*.mat)'},'Load workspace');
    filename = [filepath filename];
end

if isempty(filename) || ~ischar(filename)
    % user pressed cancel: nothing to do
elseif exist(filename,'file')
    fprintf('Loading workspace...');
    evalin('base',sprintf('clear; io_load(''%s'');',filename));
    fprintf('done.\n');
else
    error(['Could not find ' filename]);
end

repo: REST


text:
function r8_uniform_01_test ( )

%*****************************************************************************80
%
%% R8_UNIFORM_01_TEST tests R8_UNIFORM_01.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_UNIFORM_01_TEST\n' );
  fprintf ( 1, '  R8_UNIFORM_01 produces a sequence of random values.\n' );

  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Using random seed %d\n', seed );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  SEED  R8_UNIFORM_01(SEED)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : 10
    seed_old = seed;
    [ x, seed ] = r8_uniform_01 ( seed );
    fprintf ( 1, '  %12d  %14f\n', seed, x );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Verify that the sequence can be restarted.\n' );
  fprintf ( 1, '  Set the seed back to its original value, and see that\n' );
  fprintf ( 1, '  we generate the same sequence.\n' );

  seed = 123456789;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  SEED  R8_UNIFORM_01(SEED)\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 10
    seed_old = seed;
    [ x, seed ] = r8_uniform_01 ( seed );
    fprintf ( 1, '  %12d  %14f\n', seed, x );
  end

  return
end


text:
function S = mysign(A)
%MYSIGN True sign function with MYSIGN(0) = 1.

%   Called by various matrices in elmat/private.
%
%   Nicholas J. Higham
%   Copyright 1984-2005 The MathWorks, Inc.
%   $Revision: 1.4.4.1 $  $Date: 2005/11/18 14:15:15 $

S = sign(A);
S(find(S==0)) = 1;

author: andrewssobral

repo: lrslibrary


text:
function X = power(X,d)
% POWER (overloaded)

% Vectorize x if d is vector
if numel(X)==1 & (numel(d)>1)
    X = X.*ones(size(d));
end
% Vectorize if x is a vector
if numel(d)==1 & (numel(X)>1)
    d = d.*ones(size(X));
end
s = size(X);
if isa(X,'sdpvar')
    X = sdpvar(X);
else
    X = X(:);
end
d = reshape(d,[],1);
X = power(X,d);
X = reshape(X,s);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% sph2topo() - Convert from a 3-column headplot file in spherical coordinates
%              to 3-column topoplot() locs file in polar (not cylindrical) coords.
%              Used for topoplot() and other 2-D topographic plotting programs.
%              Assumes a spherical coordinate system in which horizontal angles
%              have a range [-180,180] deg,  with zero pointing to the right ear.
%              In the output polar coordinate system, zero points to the nose.
%              See  >> help readlocs
% Usage:
%          >> [chan_num,angle,radius] = sph2topo(input,shrink_factor,method);
%
% Inputs:
%   input         = [channo,az,horiz] = chan_number, azumith (deg), horiz. angle (deg)
%                   When az>0, horiz=0 -> right ear, 90 -> nose
%                   When az<0, horiz=0 -> left ear, -90 -> nose
%   shrink_factor = arc_length shrinking factor>=1 (deprecated).
%                   1 -> plot edge is 90 deg azimuth {default};
%                   1.5 -> plot edge is +/-135 deg azimuth See
%                   >> help topoplot().
%   method        = [1|2], optional. 1 is for Besa compatibility, 2 is for
%                   compatibility with Matlab function cart2sph(). Default is 2
%
% Outputs:
%   channo  = channel number (as in input)
%   angle   = horizontal angle (0 -> nose; 90 -> right ear; -90 -> left ear)
%   radius  = arc_lengrh from vertex (Note: 90 deg az -> 0.5/shrink_factor);
%             By topoplot() convention, radius=0.5 is the nasion-ear_canal plane.
%             Use topoplot() 'plotrad' to plot chans with abs(az) > 90 deg.
%
% Author: Scott Makeig & Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 6/12/98
%
% See also: cart2topo(), topo2sph()

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

% corrected left/right orientation mismatch, Blair Hicks 6/20/98
% changed name sph2pol() -> sph2topo() for compatibility -sm
% 01-25-02 reformated help & license -ad
% 01-25-02 changed computation so that it works with sph2topo -ad

function [channo,angle,radius] = sph2topo(input,factor, method)

chans = size(input,1);
angle = zeros(chans,1);
radius = zeros(chans,1);

if nargin < 1
   help sph2topo
   return
end

if nargin< 2
  factor = 0;
end
if factor==0
  factor = 1;
end
if factor < 1
  help sph2topo
  return
end

if size(input,2) ~= 3
   help sph2topo
   return
end

channo = input(:,1);
az = input(:,2);
horiz = input(:,3);

if exist('method')== 1 & method == 1
  radius = abs(az/180)/factor;
  i = find(az>=0);
  angle(i) = 90-horiz(i);
  i = find(az<0);
  angle(i) = -90-horiz(i);
else
  angle  = -horiz;
  radius = 0.5 - az/180;
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [varargout] = dcopf(varargin)
%DCOPF  Solves a DC optimal power flow.
%   This is a simple wrapper function around OPF that sets the 'model'
%   option to 'DC' before calling OPF.
%   See OPF for the details of input and output arguments.
%
%   See also RUNDCOPF.

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

[mpc, mpopt] = opf_args(varargin{:});
mpopt = mpoption(mpopt, 'model', 'DC');
[varargout{1:nargout}] = opf(mpc, mpopt);

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function prandtl_values_test ( )

%*****************************************************************************80
%
%% PRANDTL_VALUES_TEST demonstrates the use of PRANDTL_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2005
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'PRANDTL_VALUES_TEST:\n' );
  fprintf ( 1, '  PRANDTL_VALUES stores values of\n' );
  fprintf ( 1, '  the Prandtl number of water\n' );
  fprintf ( 1, '  as a function of temperature and pressure.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      T            P            Pr(T,P)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, tc, p, pr ] = prandtl_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %12f  %24.16f\n', tc, p, pr );

  end

  return
end


text:
function [ node_num, triangle_num, hole_num ] = ...
  triangulation_order6_example2_size ( )

%*****************************************************************************80
%
%% TRIANGULATION_ORDER6_EXAMPLE2_SIZE sets sizes for a sample triangulation.
%
%  Diagram:
%
%   21-22-23-24-25
%    |\  6 |\  8 |
%    | \   | \   |
%   16 17 18 19 20
%    |   \ |   \ |
%    | 5  \| 7  \|
%   11-12-13-14-15
%    |\  2 |\  4 |
%    | \   | \   |
%    6  7  8  9 10
%    | 1 \ | 3 \ |
%    |    \|    \|
%    1--2--3--4--5
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer NODE_NUM, the number of points.
%
%    Output, integer TRIANGLE_NUM, the number of triangles.
%
%    Output, integer HOLE_NUM, the number of holes.
%
  node_num = 25;
  triangle_num = 8;
  hole_num = 0;

  return
end


text:
%VL_NNNORMALIZE CNN Local Response Normalization (LRN)
%   Y = VL_NNORMALIZE(X, PARAM) computes the so-called Local Response
%   Normalization (LRN) operator. This operator performs a
%   channel-wise sliding window normalization of each column of the
%   input array X. The normalized output is given by:
%
%     Y(i,j,k) = X(i,j,k) / L(i,j,k)^BETA
%
%   where the normalization factor is given by
%
%     L(i,j,k) = KAPPA + ALPHA * (sum_{q in Q(k)} X(i,j,k)^2,
%
%   PARAM = [N KAPPA ALPHA BETA], and N is the size of the window. The
%   window Q(k) is defined as:
%
%     Q(k) = [max(1, k-FLOOR((N-1)/2)), min(D, k+CEIL((N-1)/2))].
%
%   where D is the number of feature channels in X. Note in particular
%   that, by setting N >= 2D, the function can be used to normalize
%   all the channels as a single group (useful to achieve L2
%   normalization).
%
%   DZDX = VL_NNORMALIZE(X, PARAM, DZDY) computes the derivative of
%   the block projected onto DZDY. DZDX and DZDY have the same
%   dimensions as X and Y respectively.
%
%   **Remark:** Some CNN libraries (e.g. Caffe) use a slightly
%   different convention for the parameters of the LRN. Caffe in
%   particular uses the convention:
%
%     PARAM_CAFFE = [N KAPPA N*ALPHA BETA]
%
%   i.e. the ALPHA paramter is multiplied by N.

% Copyright (C) 2014 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: willard-yuan

repo: cnn-for-image-retrieval


text:
function test_ft_plot_headshape

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

shape.pnt = randn(500,3);
shape.fid.pnt = randn(3,3);
shape.fid.label = {'nas', 'lpa', 'rpa'};

figure
ft_plot_headshape(shape);


text:
function [x,y,z,s,w,flag] = msquadsolve(Q,c,A,b,C),
% MSQUADSOLVE
%
% USAGE:   [x,y,z,s,w,flag] = msquadsolve(Q,c,A,b,C)
%
% PARAMETERS:  Q -> (n,n) symetric matrix (definite positive)
%              c -> (n,1) vector
%              A -> (m,n) matrix
%              b -> (m,1) vector
%              C -> (m,1) vector
%
%            x -> primal variables
%            y -> lagrangian coeff of equality constraints
%            z -> dual variables of x
%            s -> primal auxiliary variable (only if C < Inf)
%            w -> dual variable of s
%            flag -> set to 0 => no problem, set to 1 => problem
%
% DESCRIPTION: Primal-dual method for quadratic programming
%
%            minimize c'*x + 0.5*x'*Q*x
%
%            subject to  A*x=b
%                        0<= x <= C
%            The method used here is a primal dual method with a predictor-corrector
%            approach and a logarithmic barrier. I used the heuristic from two
%            existing methods LOQO and HOPDM. The method is an iterative method. The
%            maximal number of iteration is stored in the variable 'max_iter'.
%
% ERRORS AND BUGS: 1. There is no test about the conditionning of the matrix Q. If the iteration
%                  50 has been reached, then the optimization may not be finished and the output
%                  may be wrong.
%                  2. If C contains infinite values then, the algorithm will consider that all its
%                  components are actually infinite.
%
% NOTES: 50 iterations have always been sufficient to solve all problems.
%        This code should be read with the tech. report:
%                "Regularized Symmetric Indefinite Systems in Interior Point
%                 Methods for Linear and Quadratic Optimization",
%                 A. Altman and J. Gondzio, Logilab Tech. Report 1998.6
%
% Andre Elisseeff, May. 2001

  %% init
  verbose = 0;
  n = size(Q,1);
  m = size(A,1);
  H = zeros(n+m,n+m);
  flag=1;
  maxC = max(C);

  %% Values of the original HOPDM of Gondzio and Altmann
  dinf = 10^(-14);
  smallz = 10^(-14);
  smallt = 2.3*10^(-16);
  opttol = 10^(-6);
  itref = 1;
  mu = 1;
  maxiter = 50;

  %% init values of the primal and dual variables
  x=ones(n,1);
  z=ones(n,1);
  y=ones(m,1);
  if maxC < Inf,
    s=ones(n,1);
    w=ones(n,1);
  else
    s=[];w=[];
  end;

  %% Description of variables:
  %%
  %%    x,s     -> primal variables
  %%    z,w     -> dual variables
  %%    n       -> number of variables in the initial pb (size of x)
  %%    m       -> number of constraints in A
  %%
  %%    dinf     -> smallest value for all variables
  %%    smallz   -> smallest value of z
  %%    smallt   -> smallest value for t in the computation of the matrix theta
  %%    opttol   -> acceptable tolerance for optimality conditions
  %%    itref    -> iteration counter
  %%    maxiter  -> maximum number of iteration

  %% Analyze the constraints...
  disp(sprintf('Analyzing the equality constraints...\n'));
  [QQ,RR]=qr(A',0);
  [mm,nb] = size(QQ); %% number of eq constraints
  ind = 1:1:nb;
  for i=1:nb,
    if abs(RR(i,i)) < 100*eps
      disp(sprintf('Constraints %d removed because of dependence\n',i));
      ind(i)=0;
    end;
  end;
  indice = find(ind >0);

  if (isempty(indice))
    disp(sprintf('No equality constraints... \n'));
    A=[];
    m=0;
  else
    A = A(indice,:); %% new independent eq constraints
    b = b(indice);
    y = y(indice);
    m=length(indice);
  end;
  clear QQ;clear RR;
  u = C;
  %% init values before looping
  cont = 1;
  objQ = 0.5*x'*Q*x;

  %% init values of primal and dual objective functions
  pobjo = abs(c'*x+objQ) + 1;
  if maxC < inf,
    dobjo = abs(b'*y - u'*w - objQ);
  else
    dobjo = abs(b'*y - objQ);
  end;
%%%%%%%%%%%%%%
%% MAIN LOOP
%%%%%%%%%%%%%%
  while (cont)&(itref<=maxiter)
    %% Compute the primal objective function
    objQ = 0.5*x'*Q*x;
    pobj = c'*x+objQ;
    if (maxC<Inf)
      dobj = b'*y - u'*w - objQ;
    else
      dobj = b'*y - objQ;
    end;
    dlgap = pobj - dobj;
    dp = abs(pobj)/(abs(pobjo)+1);
    dd = abs(dobj)/(abs(dobjo)+1);
    dobjo = dobj;
    pobjo = pobj;
    if verbose,
        disp(sprintf('%d - pobj : %f - dobj : %f\n',itref,pobj,dobj));
    end;
    %% Check if the solution are bounded
    if (dp > 10^6)
      disp(sprintf('Solution not bounded in the primal. Exit.\n'));
      return;
    end;
    if (dd > 10^6)
      disp(sprintf('Solution not bounded in the dual. Exit.\n'));
      return;
    end;

    %% test if optimality
    oldgap = dlgap;
    dp = abs(dobj) + 1;
    if ((abs(dlgap)/dp) <= opttol)
      disp(sprintf('Optimal solution found. Exit.\n'));
      cont = 0;
      break;
    end;

    dp = dp + abs(pobj);
    T = abs(dlgap)/dp;

    %% put the variables away from zero (from HOPDM)
    if (itref <= 3)
      ax = 2*10^(-3);
      az = 10^(-3);
    elseif (T >= 0.8)
      ax = 2*10^(-4);
      az = 10^(-4);
    elseif (T >= 0.1)
      ax = 2*10^(-5);
      az = 10^(-5);
    elseif (T >=0.01)
      ax = 2*10^(-6);
      az = 10^(-6);
    elseif (T>=0.001)
      ax = 2*10^(-7);
      az = 10^(-7);
    elseif (T>=0.0001)
      ax = 2*10^(-8);
      az = 10^(-7);
    elseif (T>=0.00001)
      ax = 2*10^(-9);
      az = 10^(-9);
    else
      ax = T*10^(-5);
      az = ax;
    end;

    %% consider only variables that can be changed
    x = x + ax;
    z = z + az;
    if maxC < Inf
      s = s + ax;
      w = w + az;
    end;
    %% Compute the values of xi_b, xi_c and xi_u
    xi_b = -A*x + b;
    xi_c = c - A'*y - z + Q*x;
    xi_z =  - x.*z;
    if maxC < Inf,
      xi_c=xi_c + w;
      xi_u = u - x - s;
      xi_w = - s.*w;
    end;
    %% Compute theta = (z/x + w/s)

    %% for bounded variables
    if maxC<Inf,
      dp = x;
      if (max(abs(dp))<= smallz)
        disp(sprintf('Conditioning problem to invert theta. Abort.\n'));
        return;
      end;
      dpp= s;
      if (max(abs(dpp))<= smallz)
        disp(sprintf('Conditioning problem to invert theta. Abort.\n'));
        return;
      end;
      theta=z./dp + w./dpp;
    end;
    %% for unbounded variables
    if (maxC==Inf),
      dp = x;
      if (max(abs(dp))<= smallz)
        disp(sprintf('Conditioning problem to invert theta. Abort.\n'));
        return;
      end;
      theta = z./dp;
    end;

    %% neglect small elements of theta array
    neglect = find(theta < smallt);
    if ~isempty(neglect)
      theta(neglect)=zeros(size(neglect));
    end;

    %% and control large elements of theta
        neglect = find(theta >= 10^8);
    if ~isempty(neglect)
      theta(neglect)=(10^4)*sqrt(theta(neglect));
    end;
    %% factorize H = [-Q-theta^(-1)   A^T]
    %%               [ A               0 ]

    H = zeros(n+m,n+m);
    H(1:n,1:n) = -Q-diag(theta);
    H(n+1:n+m,1:n) = A;
    H(1:n,n+1:n+m) = A';

    %% Compute the predictor step
    if maxC < Inf,
      f = xi_c-xi_z./x+(xi_w - xi_u.*w)./s ;
      h = xi_b;
    else
      f = xi_c - xi_z./x;
      h = xi_b;
    end;
    delta=H\[f;h];
    dx = delta(1:n);
    dy = delta(n+1:n+m);
    dz = (xi_z-z.*dx)./x;
    if maxC<Inf,
      ds = xi_u - dx;
      dw = (xi_w-w.*ds)./s;
    end;

   %% determine the maximum step size alpha_p (primal) and
   %% alpha_d (dual) to stay in feasible region
   %% (x,s,z,w must be positive and greater than dinf)
      indz = find(dz<0);
      indx = find(dx<0);
      inds=[];mins=1;
      indw=[];minw=1;
      if maxC < Inf,
        inds = find(ds<0);
        indw = find(dw<0);
        if ~isempty(inds)
          mins = min(-(s(inds)-dinf)./ds(inds));
        else
      mins = 1;
        end;
        if ~isempty(indw)
          minw = min(-(w(indw)-dinf)./dw(indw));
        else
      minw = 1;
        end;
      end;
      if ~isempty(indx)
        minx = min(-(x(indx)-dinf)./dx(indx));
      else
        minx = 1;
      end;
      apk = min([minx,mins,1]);
      if ~isempty(indz),
        minz = min(-(z(indz)-dinf)./dz(indz));
      else
        minz = 1;
      end;
      adk = min([minw,minz,1]);

      ax = sum(x.*z);
      as = sum((x+apk*dx).*(z+adk*dz));
      az = sum(dx.^2+dz.^2);
      if maxC < Inf,
        ax = ax + sum(s.*w);
        as = as + sum((s+apk*ds).*(w+adk*dw));
        az = az + sum(ds.^2+dw.^2);
      end;
      %% check if complementary gap is less than opttol
      if (as <= opttol)
        disp(sprintf('Complementary gap is less than %f\n',opttol));
        cont = 0;
        x = x + apk*dx;
        y = y + adk*dy;
        z = z + adk*dz;
        if maxC < Inf,
         s = s + apk*ds;
         w = w + adk*dw;
        end;
        break;
      end;

      %% Set the barrier parameter : LOQO's heuristic
      ap = min(apk,adk);
      mu = (ax/(2*n))*(0.95*(1/ap) -1)^2/(0.95*(1/ap)+10)^2;

      %% Compute the new direction (algo. of Mehrotra) of order 1 (corrector step)
      xi_z =  - x.*z + mu*ones(size(x)) - dx.*dz;
      f = xi_c - xi_z./x;
      if maxC < Inf
        xi_w = - s.*w + mu*ones(size(s)) - ds.*dw;
        f=f+ xi_w./s - (w.*xi_u)./s;
      end;
      h = xi_b;

      delta=H\[f;h];
      dx = delta(1:n);
      dy = delta(n+1:n+m);
      dz = (xi_z-z.*dx)./x;
      if maxC<Inf,
        ds = xi_u - dx;
        dw = (xi_w-w.*ds)./s;
      end;

      %% determine the maximum step size alpha_p (primal) and
      %% alpha_d (dual) to stay in feasible region
      %% (x,s,z,w must be positiv)
      indz = find(dz<0);
      indx = find(dx<0);
      if maxC < Inf,
        inds = find(ds<0);
        indw = find(dw<0);
        if ~isempty(inds)
          mins = min(-s(inds)./ds(inds));
        else
        mins=1;
        end;
        if ~isempty(indw)
          minw = min(-w(indw)./dw(indw));
        else
        minw = 1;
        end;
      end;
      if ~isempty(indx)
        minx = min(-x(indx)./dx(indx));
      else
        minx = 1;
      end;
      alpha_p = min([minx,mins,1]);
      if ~isempty(indz)
        minz = min(-z(indz)./dz(indz));
      else
        minz = 1;
      end;
      alpha_d = min([minw,minz,1]);

      %% Compute step factors
      fp = 0.9*min(alpha_p,alpha_d);
      fd = 0.9*min(alpha_d,alpha_p);
      x = x +fp*dx;
      y = y +fd*dy;
      z = z + fd*dz;
      if maxC < Inf,
        w = w + fd*dw;
        s = s +fp*ds;
      end;
      itref=itref+1;
  end;
%%%%%%%%%%%%%%%%%%%%
%% End of main loop
%%%%%%%%%%%%%%%%%%%%

    if (cont==0)
      disp(sprintf('Optimal Solution found after %d iteration.\n',itref));
      disp(sprintf('Value of the objective : %f.\n',pobj));
      flag=0;
    end;
end; %% function


text:
function [xf,wftilda,wftilda2,wftilda3]=pfilterj(a,b,c,d,g,h,lambda,X0,y)
T=length(y);
n=length(X0);
xf(:,1)=X0;
wftilda(:,1)=1/n*ones(n,1);
for t=2:T
xf(:,t)=gaussmix(mx(a,xf(:,t-1)),b,mx([a,g],xf(:,t-1)),sqrt(1+h)*b,lambda);
%Compute and normalize filter weights
wf(:,t)=normpdf(y(t),my(c,xf(:,t)),d);
%Ensures a nonzero summation in denominator
[n1,val]=zerotest(y(t),my(c,xf(:,t)),wf(:,t));
wf(n1,t)=val;
wftilda(:,t)=wf(:,t)/sum(wf(:,t));
end

for i=1:n
wf3(i,T)=wftilda(i,T)*weight(xf(i,T),xf(:,T),wftilda(i,T),wftilda(:,T),i,a,b,g,h,lambda);
end
wftilda3(:,T)=wf3(:,T)/sum(wf3(:,T));

for t=T-1:-1:1
%Compute joint smoother weight p(x(t),x(t+1)|Y) and p(x(t),x(t)|Y)
for i=1:n
wf2(i,t)=wftilda(i,t)*weight(xf(i,t+1),xf(:,t),wftilda(i,t+1),wftilda(:,t),i,a,b,g,h,lambda);
wf3(i,t)=wftilda(i,t)*weight(xf(i,t),xf(:,t),wftilda(i,t),wftilda(:,t),i,a,b,g,h,lambda);
end

%Ensures a nonzero summation in denominator
if sum(wf2(:,t))==0;
wftilda2(:,t)=1/n*ones(n,1);
else
wftilda2(:,t)=wf2(:,t)/sum(wf2(:,t));
end

if sum(wf3(:,t))==0;
wftilda3(:,t)=1/n*ones(n,1);
else
wftilda3(:,t)=wf3(:,t)/sum(wf3(:,t));
end

end

author: Sable

repo: mcbench-benchmarks


text:
function plotPointClouds(M, B, markerm, markerb)
    if (nargin < 4)
        markerm = 'r.';
    end
    if (nargin < 3)
        markerb = 'b.';
    end
    %figure;
    plot3(M(1,:), M(2,:), M(3,:), markerm, B(1,:), B(2,:), B(3,:), markerb, 'MarkerSize',0.8);
    axis off;

end

author: intellhave

repo: SDRSAC

save_path: github-repos/MATLAB/intellhave-SDRSAC


text:
function [M,Q]=community_louvain(W,gamma,M0,B)
%COMMUNITY_LOUVAIN     Optimal community structure
%
%   M     = community_louvain(W);
%   [M,Q] = community_louvain(W,gamma);
%   [M,Q] = community_louvain(W,gamma,M0);
%   [M,Q] = community_louvain(W,gamma,M0,'potts');
%   [M,Q] = community_louvain(W,gamma,M0,'negative_asym');
%   [M,Q] = community_louvain(W,[],[],B);
%
%   The optimal community structure is a subdivision of the network into
%   nonoverlapping groups of nodes which maximizes the number of within-
%	group edges, and minimizes the number of between-group edges.
%
%   This function is a fast and accurate multi-iterative generalization of
%   the Louvain community detection algorithm. This function subsumes and
%   improves upon,
%		modularity_louvain_und.m, modularity_finetune_und.m,
%		modularity_louvain_dir.m, modularity_finetune_dir.m,
%       modularity_louvain_und_sign.m
%	and additionally allows to optimize other objective functions (includes
%	built-in Potts-model Hamiltonian, allows for custom objective-function
%	matrices).
%
%   Inputs:
%       W,
%           directed/undirected weighted/binary connection matrix with
%           positive and possibly negative weights.
%       gamma,
%           resolution parameter (optional)
%               gamma>1,        detects smaller modules
%               0<=gamma<1,     detects larger modules
%               gamma=1,        classic modularity (default)
%       M0,
%           initial community affiliation vector (optional)
%       B,
%           objective-function type or custom objective matrix (optional)
%           'modularity',       modularity (default)
%           'potts',            Potts-model Hamiltonian (for binary networks)
%           'negative_sym',     symmetric treatment of negative weights
%           'negative_asym',    asymmetric treatment of negative weights
%           B,                  custom objective-function matrix
%
%           Note: see Rubinov and Sporns (2011) for a discussion of
%           symmetric vs. asymmetric treatment of negative weights.
%
%   Outputs:
%       M,
%           community affiliation vector
%       Q,
%           optimized community-structure statistic (modularity by default)
%
%   Example:
%       % Iterative community finetuning.
%       % W is the input connection matrix.
%       n  = size(W,1);             % number of nodes
%       M  = 1:n;                   % initial community affiliations
%       Q0 = -1; Q1 = 0;            % initialize modularity values
%       while Q1-Q0>1e-5;           % while modularity increases
%           Q0 = Q1;                % perform community detection
%           [M, Q1] = community_louvain(W, [], M);
%       end
%
%   References:
%       Blondel et al. (2008)  J. Stat. Mech. P10008.
%       Reichardt and Bornholdt (2006) Phys. Rev. E 74, 016110.
%       Ronhovde and Nussinov (2008) Phys. Rev. E 80, 016109
%       Sun et al. (2008) Europhysics Lett 86, 28004.
%       Rubinov and Sporns (2011) Neuroimage 56:2068-79.
%
%   Mika Rubinov, U Cambridge 2015-2016

%   Modification history
%   2015: Original
%   2016: Included generalization for negative weights.
%         Enforced binary network input for Potts-model Hamiltonian.
%         Streamlined code and expanded documentation.

W=double(W);                                % convert to double format
n=length(W);                                % get number of nodes
s=sum(sum(W));                              % get sum of edges

if ~exist('B','var') || isempty(B)
    type_B = 'modularity';
elseif ischar(B)
    type_B = B;
else
    type_B = 0;
    if exist('gamma','var') && ~isempty(gamma)
        warning('Value of gamma is ignored in generalized mode.')
    end
end
if ~exist('gamma','var') || isempty(gamma)
    gamma = 1;
end

if strcmp(type_B,'negative_sym') || strcmp(type_B,'negative_asym')
    W0 = W.*(W>0);                          %positive weights matrix
    s0 = sum(sum(W0));                      %weight of positive links
    B0 = W0-gamma*(sum(W0,2)*sum(W0,1))/s0; %positive modularity

    W1 =-W.*(W<0);                          %negative weights matrix
    s1 = sum(sum(W1));                      %weight of negative links
    if s1                                   %negative modularity
        B1 = W1-gamma*(sum(W1,2)*sum(W1,1))/s1;
    else
        B1 = 0;
    end
elseif min(min(W))<-1e-10
    err_string = [
        'The input connection matrix contains negative weights.\nSpecify ' ...
        '''negative_sym'' or ''negative_asym'' objective-function types.'];
    error(sprintf(err_string))              %#ok<SPERR>
end
if strcmp(type_B,'potts') && any(any(W ~= logical(W)))
    error('Potts-model Hamiltonian requires a binary W.')
end

if type_B
    switch type_B
        case 'modularity';      B = (W-gamma*(sum(W,2)*sum(W,1))/s)/s;
        case 'potts';           B =  W-gamma*(~W);
        case 'negative_sym';    B = B0/(s0+s1) - B1/(s0+s1);
        case 'negative_asym';   B = B0/s0      - B1/(s0+s1);
        otherwise;              error('Unknown objective function.');
    end
else                            % custom objective function matrix as input
    B = double(B);
    if ~isequal(size(W),size(B))
        error('W and B must have the same size.')
    end
end
if ~exist('M0','var') || isempty(M0)
    M0=1:n;
elseif numel(M0)~=n
    error('M0 must contain n elements.')
end

[~,~,Mb] = unique(M0);
M = Mb;

B = (B+B.')/2;                                          % symmetrize modularity matrix
Hnm=zeros(n,n);                                         % node-to-module degree
for m=1:max(Mb)                                         % loop over modules
    Hnm(:,m)=sum(B(:,Mb==m),2);
end

Q0 = -inf;
Q = sum(B(bsxfun(@eq,M0,M0.')));                        % compute modularity
first_iteration = true;
while Q-Q0>1e-10
    flag = true;                                        % flag for within-hierarchy search
    while flag;
        flag = false;
        for u=randperm(n)                               % loop over all nodes in random order
            ma = Mb(u);                                 % current module of u
            dQ = Hnm(u,:) - Hnm(u,ma) + B(u,u);
            dQ(ma) = 0;                                 % (line above) algorithm condition

            [max_dQ,mb] = max(dQ);                      % maximal increase in modularity and corresponding module
            if max_dQ>1e-10;                            % if maximal increase is positive
                flag = true;
                Mb(u) = mb;                             % reassign module

                Hnm(:,mb) = Hnm(:,mb)+B(:,u);           % change node-to-module strengths
                Hnm(:,ma) = Hnm(:,ma)-B(:,u);
            end
        end
    end
    [~,~,Mb] = unique(Mb);                              % new module assignments

    M0 = M;
    if first_iteration
        M=Mb;
        first_iteration=false;
    else
        for u=1:n                                       % loop through initial module assignments
            M(M0==u)=Mb(u);                             % assign new modules
        end
    end

    n=max(Mb);                                          % new number of modules
    B1=zeros(n);                                        % new weighted matrix
    for u=1:n
        for v=u:n
            bm=sum(sum(B(Mb==u,Mb==v)));                % pool weights of nodes in same module
            B1(u,v)=bm;
            B1(v,u)=bm;
        end
    end
    B=B1;

    Mb=1:n;                                             % initial module assignments
    Hnm=B;                                              % node-to-module strength

    Q0=Q;
    Q=trace(B);                                         % compute modularity
end


text:
load -ascii soybeanD
load -ascii soybeanTD

% soybeanD(2:36,:) = soybeanD(2:36,:)+1;
% soybeanTD(2:36,:) = soybeanTD(2:36,:)+1;
% soybeanD(find(soybeanD==-9998))=-9999;
% soybeanTD(find(soybeanTD==-9998))=-9999;

[N, Napp] = size(soybeanD);
[N, Ntest] = size(soybeanTD);

N, m=Napp+Ntest, Napp, Ntest,
class = 1

%rand('state',0); randn('state',0);
%abalone = abalone(:,randperm(m));

app  = soybeanD;
test = soybeanTD;

unique(app(class,:))
unique(test(class,:))

ns = max(soybeanD')
clear soybeanD soybeanTD

% N, ns(class), Napp, Ntest, mean(ns),

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [lrp] = ft_lateralizedpotential(cfg, avgL, avgR)

% FT_LATERALIZEDPOTENTIAL computes lateralized potentials such as the
% lateralized readiness potential (LRP)
%
% Use as
%   [lrp] = ft_lateralizedpotential(cfg, avgL, avgR)
%
% where the input datasets should come from FT_TIMELOCKANALYSIS
% and the configuration should contain
%   cfg.channelcmb = Nx2 cell-array
%
% An example channelcombination containing the homologous channels
% in the 10-20 standard system is
%    cfg.channelcmb = {'Fp1'   'Fp2'
%                      'F7'    'F8'
%                      'F3'    'F4'
%                      'T7'    'T8'
%                      'C3'    'C4'
%                      'P7'    'P8'
%                      'P3'    'P4'
%                      'O1'    'O2'}
%
% The lateralized potential is computed on combinations of channels and
% not on indivudual channels. However, if you want to make a topographic
% plot with e.g. FT_MULTIPLOTER, you can replace the output lrp.label
% with lrp.plotlabel.
%
% The concept for the LRP was introduced approximately simultaneously in the
% following two papers
% - M. G. H. Coles. Modern mind-brain reading - psychophysiology,
%   physiology, and cognition. Psychophysiology, 26(3):251-269, 1988.
% - R. de Jong, M. Wierda, G. Mulder, and L. J. Mulder. Use of
%   partial stimulus information in response processing. J Exp Psychol
%   Hum Percept Perform, 14:682-692, 1988.
% and it is discussed in detail on a technical level in
% - R. Oostenveld, D.F. Stegeman, P. Praamstra and A. van Oosterom.
%   Brain symmetry and topographic analysis of lateralized event-related
%   potentials. Clin Neurophysiol. 114(7):1194-202, 2003.
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a *.mat
% file on disk and/or the output data will be written to a *.mat file. These mat
% files should contain only a single variable, corresponding with the
% input/output structure.
%
% See also FT_TIMELOCKANALYSIS, FT_MULTIPLOTER

% Copyright (C) 2004, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% these are used by the ft_preamble/ft_postamble function and scripts
ft_revision = '$Id$';
ft_nargin   = nargin;
ft_nargout  = nargout;

% do the general setup of the function
ft_defaults
ft_preamble init
ft_preamble debug
ft_preamble loadvar avgL avgR
ft_preamble provenance avgL avgR

% the ft_abort variable is set to true or false in ft_preamble_init
if ft_abort
  return
end

% check if the input data is valid for this function
avgL = ft_checkdata(avgL, 'datatype', 'timelock');
avgR = ft_checkdata(avgR, 'datatype', 'timelock');

% check if the input cfg is valid for this function
cfg = ft_checkconfig(cfg, 'forbidden',  {'channels'}); % prevent accidental typos, see issue 1729

% set the defaults
if ~isfield(cfg, 'channelcmb')
  cfg.channelcmb = {
    'Fp1'   'Fp2'
    'F7'    'F8'
    'F3'    'F4'
    'T7'    'T8'
    'C3'    'C4'
    'P7'    'P8'
    'P3'    'P4'
    'O1'    'O2'
    };
end

if ~isequal(avgL.time, avgR.time)
  ft_error('the time axes are not the same');
end

% start with an empty output structure
lrp.label     = {};
lrp.plotlabel = {};
lrp.avg       = [];
lrp.time      = avgL.time;

% add timelock signature
if isfield(avgL, 'dimord') && isfield(avgR, 'dimord')
    if ~strcmp(avgL.dimord, avgR.dimord)
        ft_error('The input data are of different dimord types');
    else
        lrp.dimord = avgL.dimord;
    end
else
    ft_error('''dimord'' not found. The function expects timelock data');
end

% compute the lateralized potentials
Nchan = size(cfg.channelcmb,1);
for i=1:Nchan
  % here the channel names "C3" and "C4" are used to clarify the
  % computation of the lateralized potential on all channel pairs
  C3R = strcmp(cfg.channelcmb{i,1}, avgR.label);
  C4R = strcmp(cfg.channelcmb{i,2}, avgR.label);
  C3L = strcmp(cfg.channelcmb{i,1}, avgL.label);
  C4L = strcmp(cfg.channelcmb{i,2}, avgL.label);
  if any(C3R) && any(C4R) && any(C3L) && any(C4L)
    lrp.label{end+1}     = sprintf('%s/%s', cfg.channelcmb{i,1}, cfg.channelcmb{i,2});
    lrp.plotlabel{end+1} = cfg.channelcmb{i,1};
    erpC3L = avgL.avg(C3L,:);
    erpC4L = avgL.avg(C4L,:);
    erpC3R = avgR.avg(C3R,:);
    erpC4R = avgR.avg(C4R,:);
    lrp.avg(end+1,:) = 1/2 * ((erpC3R - erpC4R) + (erpC4L - erpC3L));
  end
end

% do the general cleanup and bookkeeping at the end of the function
ft_postamble debug
ft_postamble previous   avgL avgR
ft_postamble provenance lrp
ft_postamble history    lrp
ft_postamble savevar    lrp


text:
function [lab, err] = iccvTransferLabels(labels, regions, area)
% [lab, err] = transferLabels(labels, regions, area)

nlab = max(labels);

nr = numel(regions);
lab = zeros(nr, 1);

npix = sum(area.*(labels>0));
total = 0;

count = zeros(nr,1);
for k = 1:nr
    count(k) = sum(labels(regions{k})>0);
end
for k = find(count==1)'
    valid = labels(regions{k})>0;
    lab(k) = labels(regions{k}(valid));
    total = total + area(regions{k}(valid));
end

for k = find(count>1)'
    rcount = zeros(nlab, 1);
    origlabs = labels(regions{k});
    valid = find(origlabs>0);
    for k2 = valid(:)'
        rcount(origlabs(k2)) = rcount(origlabs(k2)) + area(regions{k}(k2));
    end
    [tmp, lab(k)] = max(rcount);
    total = total + tmp;
end

err = (npix-total) / npix;

author: Cloud-CV

repo: object-proposals