text:
function [vert,conn,tria,tnum] = refine2_om(varargin)
%REFINE2_OM (Frontal)-Delaunay-refinement for two-dimensional,
% polygonal geometries.
% Has been slightly edited for OceanMesh2D to make sure boundary edges are
% absolutely preserved and to get rid of non-unique points
%   [VERT,EDGE,TRIA,TNUM] = REFINE2(NODE,EDGE) returns a co-
%   nstrained Delaunay triangulation of the polygonal region
%   {NODE,EDGE}. NODE is an N-by-2 array of polygonal verti-
%   ces and EDGE is an E-by-2 array of edge indexing. Each
%   row in EDGE represents an edge of the polygon, such that
%   NODE(EDGE(JJ,1),:) and NODE(EDGE(JJ,2),:) are the coord-
%   inates of the endpoints of the JJ-TH edge. If the argum-
%   ent EDGE is omitted it assumed that the vertices in NODE
%   are connected in ascending order.
%
%   [...] = REFINE2(NODE,EDGE,PART) computes a triangulation
%   for a multiply-connected geometry. PART is a cell-array
%   of polygonal "parts", where each element PART{KK} is an
%   array of edge indices defining a given polygonal region.
%   EDGE(PART{KK}, :) is the set of edges in the KK-TH part.
%
%   VERT is a V-by-2 array of XY coordinates in the triangu-
%   lation, EDGE is an array of constrained edges, TRIA is a
%   T-by-3 array of triangles, and TNUM is a T-by-1 array of
%   part indices. Each row of TRIA and EDGE define an eleme-
%   nt. VERT(TRIA(II,1),:), VERT(TRIA(II,2),:) and VERT(TRIA
%   (II,3),:) are the coordinates of the II-TH triangle. The
%   edges in EDGE are defined in a similar manner. NUM is an
%   array of part indexing, such that TNUM(II) is the index
%   of the part in which the II-TH triangle resides.
%
%   [...] = REFINE2(..., OPTS) passes an additional options
%   structure OPTS, containing various user-defined paramet-
%   ers, including:
%
% - OPTS.KIND = {'DELFRONT'}, 'DELAUNAY' -- the type of ref-
%   inement employed. The 'DELFRONT' algorithm is typically
%   slower, but produces higher quality output.
%
% - OPTS.RHO2 = {1.025} -- the maximum allowable radius-edge
%   ratio. Refinement proceeds until all interior triangles
%   satisfy the radius-edge threshold. Smaller radius-edge
%   ratios lead to improved triangle shape, with RHO2=1 req-
%   uiring that all angles exceed 30 degrees. Setting RHO2<1
%   may lead to non-convergence.
%
% - OPTS.REF1 = {'REFINE'}, 'PRESERVE' -- refinement 'flag'
%   for 1-dimensional faces (i.e. edges). The 'PRESERVE' op-
%   tion results in minimal refinement, attempting to retain
%   the initial edges without further subdivision. Edges are
%   split only to satisfy basic geomertical conformance.
%
% - OPTS.REF2 = {'REFINE'}, 'PRESERVE' -- refinement 'flag'
%   for 2-dimensional faces (i.e. trias). The 'PRESERVE' op-
%   tion results in minimal refinement, attempting to retain
%   the initial trias without further subdivision. Trias are
%   split only to satisfy basic geomertical conformance.
%
% - OPTS.SIZ1 = {1.333} -- the normalised rel.-length th-
%   reshold for edge-elements. Each exterior edge is refined
%   until LL/HH<SIZ1, where LL is the edge-length, HH is the
%   edge-centred mesh-size value.
%
% - OPTS.SIZ2 = {1.300} -- the normalised rel.-length th-
%   reshold for tria-elements. Each interior tria is refined
%   until RE/HH<SIZ2, where RE is an effective tria length,
%   based on the circumradius, HH is the tria-centred mesh-
%   size value.
%
% - OPTS.DISP = { +10 } -- refinement verbosity. Set to INF
%   for quiet execution.
%
%   [...] = REFINE2(..., HFUN,ARGS) also passes an optional
%   mesh-size function argument. Setting HFUN = HMAX, where
%   HMAX is a scalar value, imposes a constant size constra-
%   int over the full domain. HFUN can also be defined as a
%   general function handle [HH] = HFUN(PP), where PP is an
%   N-by-2 array of XY coordinates and HH is the associated
%   vector of mesh-size values. User-defined HFUN must be
%   fully vectorised. Additional arguments {A1,A2,...AN} for
%   HFUN can be passed as trailing parameters to REFINE2. In
%   such cases, HFUN must adopt a signature [HH] = HFUN(PP,
%   A1,A2,...,AN). HFUN must return positive values.
%
%   See also SMOOTH2, TRIDIV2, TRICOST, TRIDEMO

%   This routine implements a "multi-refinement" variant of
%   Delaunay-refinement type mesh-generation. Both standard
%   Delaunay-refinement and Frontal-Delaunay type algorithms
%   are available. The Frontal-Delaunay approach is a simpl-
%   ified version of the JIGSAW algorithm, described in:
%
% * D. Engwirda, (2014): "Locally-optimal Delaunay-refineme-
%   nt and optimisation-based mesh generation", Ph.D. Thesis
%   School of Mathematics and Statistics, Univ. of Sydney.
%   http://hdl.handle.net/2123/13148
%
% * D. Engwirda & D. Ivers, (2016): "Off-centre Steiner poi-
%   nts for Delaunay-refinement on curved surfaces", Comput-
%   er-Aided Design, (72), 157--171.
%   http://dx.doi.org/10.1016/j.cad.2015.10.007

%   This work is an extension of the "off-centre" type tech-
%   niques introduced in:
%
% * H. Erten & A. Ungor, (2009): "Quality triangulation with
%   locally optimal Steiner points", SIAM Journal on Scient-
%   ific Comp. 31(3), 2103--2130.
%   http://doi.org/10.1137/080716748
%
% * S. Rebay, (1993): "Efficient Unstructured Mesh Generati-
%   on by Means of Delaunay Triangulation and Bowyer-Watson
%   Algorithm, J. Comp. Physics 106(1), 125--138.
%   http://dx.doi.org/10.1006/jcph.1993.1097

%   Generally speaking, the Delaunay-refinement method impl-
%   emented here is a variantion of the "classical" algorit-
%   hm introduced in:
%
% * J. Ruppert, (1995): "A Delaunay refinement algorithm for
%   quality 2-dimensional mesh generation." Journal of Algo-
%   rithms 18(3), 548--585.
%   http://dx.doi.org/10.1006/jagm.1995.1021
%
%   See also: S. Cheng, T. Dey & J. Shewchuk, (2012): "Dela-
%   unay mesh generation", CRC Press, for comprehensive cov-
%   erage of Delaunay-based meshing techniques.

%   A much more advanced, and fully three-dimensional imple-
%   mentation is available in the JIGSAW library. For addit-
%   ional information, see:
%   https://github.com/dengwirda/jigsaw-matlab

%-----------------------------------------------------------
%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 09/07/2018
%-----------------------------------------------------------

    node = []; PSLG = []; part = {}; opts = [] ;
    hfun = []; harg = {};

%---------------------------------------------- extract args
    if (nargin>=+1), node = varargin{1}; end
    if (nargin>=+2), PSLG = varargin{2}; end
    if (nargin>=+3), part = varargin{3}; end
    if (nargin>=+4), opts = varargin{4}; end
    if (nargin>=+5), hfun = varargin{5}; end
    if (nargin>=+6), harg = varargin(6:end); end

   [opts] = makeopt(opts) ;

%---------------------------------------------- default EDGE
    nnod = size(node,1) ;

    if (isempty(PSLG))
        PSLG = [(1:nnod-1)',(2:nnod)'; nnod,1] ;
    end

%---------------------------------------------- default PART
    ncon = size(PSLG,1) ;

    if (isempty(part)), part{1} = (1:ncon)'; end

%---------------------------------------------- basic checks
    if (~isnumeric(node) || ~isnumeric(PSLG) || ...
        ~iscell   (part) || ~isstruct (opts) )
        error('refine2:incorrectInputClass' , ...
            'Incorrect input class.') ;
    end

%---------------------------------------------- basic checks
    if (ndims(node) ~= +2 || ndims(PSLG) ~= +2)
        error('refine2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end
    if (size(node,2) < +2 || size(PSLG,2) < +2)
        error('refine2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

%---------------------------------------------- basic checks
    if (min([PSLG(:)])<+1 || max([PSLG(:)])>nnod)
        error('refine2:invalidInputs', ...
            'Invalid EDGE input array.') ;
    end

    pmin = cellfun(@min,part);
    pmax = cellfun(@max,part);

    if (min([pmin(:)])<+1 || max([pmax(:)])>ncon)
        error('refine2:invalidInputs', ...
            'Invalid PART input array.') ;
    end

%-------------------------------- prune any non-unique topo.
   [~,ivec,jvec] = ...
        unique(sort(PSLG,+2),'rows') ;

    PSLG = PSLG(ivec,:) ;

    for ppos = +1:length(part)

        if ( ~isnumeric(part{ppos}) )
            error (  ...
            'refine2:incorrectInputClass', ...
                'Incorrect input class. ') ;
        end

        part{ppos} = ...
            unique(jvec(part{ppos})) ;

    end

%-------------------------------- check part "manifold-ness"
    for ppos = +1:length(part)

        eloc = PSLG(part{ppos},:) ;
        nadj = ...
            accumarray(eloc(:),1) ;

        if (any(mod(nadj,2) ~= 0) )
        error('refine2:nonmanifoldInputs', ...
            'Non-manifold PART detected.') ;
        end

    end

%---------------------------------------------- output title
    if (~isinf(opts.disp))
        fprintf(1,'\n') ;
        fprintf(1,' Refine triangulation...\n') ;
        fprintf(1,'\n') ;
        fprintf(1,[...
' -------------------------------------------------------\n', ...
'      |ITER.|          |CDT1(X)|          |CDT2(X)|     \n', ...
' -------------------------------------------------------\n', ...
             ] ) ;
    end

%-------------------------------- PASS 0: inflate box bounds
    vert = node; tria = []; tnum = []; iter = 0 ;
    conn = PSLG;

    vmin = min(vert,[],1);      % inflate bbox for stability
    vmax = max(vert,[],1);

    vdel = vmax - 1.*vmin;
    vmin = vmin - .5*vdel;
    vmax = vmax + .5*vdel;

    vbox = [
        vmin(1), vmin(2)
        vmax(1), vmin(2)
        vmax(1), vmax(2)
        vmin(1), vmax(2)
           ] ;
    vert = [vert ; vbox] ;

%-------------------------------- PASS 0: shield sharp feat.
   [vert,conn,tria,tnum,iter] = ...
        cdtbal0(vert,conn,tria,tnum, ...
            node,PSLG,part,opts,hfun,harg,iter);

%-------------------------------- PASS 1: refine 1-simplexes
   [vert,conn,tria,tnum,iter] = ...
        cdtref1(vert,conn,tria,tnum, ...
            node,PSLG,part,opts,hfun,harg,iter);

%-------------------------------- PASS 2: refine 2-simplexes
   [vert,conn,tria,tnum,iter] = ...
        cdtref2(vert,conn,tria,tnum, ...
            node,PSLG,part,opts,hfun,harg,iter);

    if (~isinf(opts.disp)), fprintf(1,'\n'); end

%-------------------------------- trim extra adjacency info.
    tria = tria( :,1:3) ;

%-------------------------------- trim vert. - deflate bbox.
    keep = false(size(vert,1),1);
    keep(tria(:)) = true;
    keep(conn(:)) = true;

    redo = zeros(size(vert,1),1);
    redo(keep) = ...
        (+1:length(find(keep)))';

    conn = redo(conn);
    tria = redo(tria);

    vert = vert(keep,:) ;

end

function [vert,conn,tria,tnum,iter] = ...
            cdtbal0(vert,conn,tria,tnum, ...
                node,PSLG,part,opts,hfun,harg,iter)
%CDTBAL0 constrained Delaunay-refinement for "sharp" 0-dim.
%features at PSLG vertices.
%   [...] = CDTBAL0(...) refines the set of 1-simplex eleme-
%   nts incident to "sharp" features in the PSLG. Specifica-
%   lly, edges that subtend "small" angles are split about a
%   set of new "collar" vertices, equi-distributed about the
%   centre of "sharp" features. Collar size is computed as a
%   min. of the incident edge-len. and local mesh-size cons-
%   traints.

    if (iter <= opts.iter)

    %------------------------------------- build current CDT
       [vert,conn, ...
        tria,tnum] = deltri2(vert,conn, ...
                             node,PSLG, ...
                             part, ...
                             opts.dtri) ;

    %------------------------------------- build current adj
       [edge,tria] = tricon2(tria,conn) ;

       [feat,ftri] = isfeat2(vert, ...
                             edge,tria) ;

        apex = false(size(vert,1), 1) ;
        apex(tria(ftri)) =  true ;

        if strcmpi(opts.ref1,'preserve')
            return;
        end

    %------------------------------------- eval. length-fun.
        if (~isempty(hfun))
            if (isnumeric(hfun))
            vlen = hfun * ...
              ones(size(vert,1),1) ;
            else
            vlen = feval( ...
                hfun,vert,harg{:}) ;
            vlen = vlen(:) ;
            end
        else
            vlen = +inf * ...
              ones(size(vert,1),1) ;
        end

    %------------------------------------- form edge vectors
        evec = vert(conn(:,2),:) ...
             - vert(conn(:,1),:) ;
        elen = sqrt(sum(evec.^2,2));
        evec = evec./[elen,elen] ;

    %------------------------------------- min. adj. lengths
        for epos = +1 : size(conn,1)

            ivrt = conn(epos,1) ;
            jvrt = conn(epos,2) ;

            vlen(ivrt) = min( ...
            vlen(ivrt), .67*elen(epos)) ;
            vlen(jvrt) = min( ...
            vlen(jvrt), .67*elen(epos)) ;

        end

    %------------------------------------- mark feature edge
        iref = apex(conn(:,1)) ...      %- refine at vert. 1
            & ~apex(conn(:,2)) ;
        jref = apex(conn(:,2)) ...      %- refine at vert. 2
            & ~apex(conn(:,1)) ;
        dref = apex(conn(:,1)) ...      %- refine at both!
            &  apex(conn(:,2)) ;

        keep =~apex(conn(:,1)) ...      %- refine at neither
            & ~apex(conn(:,2)) ;

    %------------------------------------- protecting collar
        ilen = vlen(conn(iref,1)) ;
        inew = vert(conn(iref,1),:) ...
        + [ilen,ilen].*evec(iref,:) ;

        jlen = vlen(conn(jref,2)) ;
        jnew = vert(conn(jref,2),:) ...
        - [jlen,jlen].*evec(jref,:) ;

        Ilen = vlen(conn(dref,1)) ;
        Inew = vert(conn(dref,1),:) ...
        + [Ilen,Ilen].*evec(dref,:) ;

        Jlen = vlen(conn(dref,2)) ;
        Jnew = vert(conn(dref,2),:) ...
        - [Jlen,Jlen].*evec(dref,:) ;

        vnew = [inew; jnew; Inew; Jnew] ;

    %------------------------------------- add new vert/edge
        iset = (1:size(inew,1))' ...
                + size(vert,1) ;

        jset = (1:size(jnew,1))' ...
                + size(inew,1) + ...
                + size(vert,1) ;

        Iset = (1:size(Inew,1))' ...
                + size(inew,1) + ...
                + size(jnew,1) + ...
                + size(vert,1) ;

        Jset = (1:size(Jnew,1))' ...
                + size(inew,1) + ...
                + size(jnew,1) + ...
                + size(Inew,1) + ...
                + size(vert,1) ;

        vert = [vert ; vnew] ;

        cnew = [conn(iref,1), iset ;
                conn(iref,2), iset ;
                conn(jref,2), jset ;
                conn(jref,1), jset ;
                conn(dref,1), Iset ;
                conn(dref,2), Jset ;
                Iset, Jset] ;
        conn = [conn(keep,:); cnew ] ;

    end

end

function [vert,conn,tria,tnum,iter] = ...
            cdtref1(vert,conn,tria,tnum, ...
                node,PSLG,part,opts,hfun,harg,iter)
%CDTREF1 constrained Delaunay-refinement for 1-simplex elem-
%nts embedded in R^2.
%   [...] = CDTREF1(...) refines the set of 1-simplex eleme-
%   nts embedded in the triangulation until all constraints
%   are satisfied. Specifically, edges are refined until all
%   local mesh-spacing and encroachment conditions are met.
%   Refinement proceeds according to either a Delaunay-refi-
%   nement or Frontal-Delaunay type approach, depending on
%   user-settings. In either case, new steiner vertices are
%   introduced to split "bad" edges - those that violate the
%   set of prescribed constraints. In the "-DR" type process
%   edges are split about their circumballs (midpoints). In
%   the "-FD" approach, new vertices are positioned such th-
%   at mesh-spacing constraints are satisfied in a "locally-
%   optimal" fashion.

    tcpu.full = +0. ;
    tcpu.ball = +0. ;
    tcpu.hfun = +0. ;
    tcpu.encr = +0. ;
    tcpu.offc = +0. ;

    vidx = (1:size(vert,1))';     %- "new" vert list to test

    tnow =  tic ;

    ntol = +1.55;

    while (strcmpi(opts.ref1,'refine'))

        iter = iter + 1 ;

        if (iter>=opts.iter),break; end

    %------------------------------------- calc. circumballs
        ttic = tic ;

        bal1 = cdtbal1(vert,conn) ;

        tcpu.ball = ...
            tcpu.ball + toc(ttic) ;

    %------------------------------------- eval. length-fun.
        ttic = tic ;

        if (~isempty(hfun))
            if (isnumeric(hfun))
            fun0 = hfun * ...
              ones(size(vert,1),1);
            fun1 = hfun ;
            else
            fun0(vidx) = ...
                feval(hfun, ...
            vert(vidx,:), harg{:});
            fun0 = fun0(:) ;
            fun1 = fun0(conn(:,1))...
                 + fun0(conn(:,2));
            fun1 = fun1 / +2. ;
            end
        else
            fun0 = +inf * ...
              ones(size(vert,1),1);
            fun1 = +inf ;
        end

        siz1 = ...
         +4. * bal1(:,3)./(fun1.*fun1) ;

        tcpu.hfun = ...
            tcpu.hfun + toc(ttic) ;

    %------------------------------------- test encroachment
        ttic = tic ;

        bal1(:,3) = ...
            (1.-eps^.75) * bal1(:,3) ;

       [vp,vi] = ...
           findball(bal1,vert(:,1:2));

    %------------------------------------- near=>[vert,edge]
        next = +0;
        ebad = false(size(conn,1),1) ;
        near = zeros(size(conn,1),1) ;
        for ii = +1 : size(vp,1)
            for ip = vp(ii,1):vp(ii,2)
                jj = vi(ip);
                if (ii ~= conn(jj,1) ...
                &&  ii ~= conn(jj,2) )
                next = next + 1;
                near(next,1) = ii;
                near(next,2) = jj;
                end
            end
        end

        near = near(1:next-0,:);

        if (~isempty(near))
    %-- mark edge "encroached" if there is a vert within its
    %-- dia.-ball that is not joined to either of its vert's
    %-- via an existing edge...
            ivrt = conn(near(:,2),1);
            jvrt = conn(near(:,2),2);

            pair = [near(:,1), ivrt];
            ivec = setset2(pair,conn) ;

            pair = [near(:,1), jvrt];
            jvec = setset2(pair,conn) ;

            okay = ~ivec & ~jvec ;

            ebad(near(okay,2))=true ;

        end

        tcpu.encr = ...
            tcpu.encr + toc(ttic);

    %------------------------------------- refinement queues
        ref1 = false(size(conn,1),1);
        ref1(ebad)           = true ;   %- edge encroachment
        ref1(siz1>opts.siz1* ...        %- bad equiv. length
                  opts.siz1) = true ;

        num1 = find(ref1)  ;

    %------------------------------------- dump-out progess!
        if (mod(iter,opts.disp)==0)
            numc = size(conn,1) ;
            numt = size(tria,1) ;
            fprintf(+1, ...
            '%11i %18i %18i\n', ...
            [iter,numc,numt]) ;
        end

    %------------------------------------- nothing to refine
        if (isempty(num1)), break; end

    %------------------------------------- refine "bad" tria
        switch (lower(opts.kind))
        case 'delaunay'
    %------------------------------------- do circ-ball pt's
        new1 = bal1(ref1, 1:2) ;

        vidx = (1:size(new1,1))' ...
                + size(vert,1) ;

        cnew = [conn( ref1,1), vidx
                conn( ref1,2), vidx];
        conn = [conn(~ref1,:); cnew];

    %------------------------------------- update vertex set
        vert = [vert; new1(:,1:2)];

        case 'delfront'
    %-- symmetric off-centre scheme:- refine edges from both
    %-- ends simultaneously, placing new vertices to satisfy
    %-- the worst of mesh-spacing and local voronoi constra-
    %-- ints.

        ttic = tic ;

        evec = vert(conn(ref1,2),:) ...
             - vert(conn(ref1,1),:) ;
        elen = sqrt(sum(evec.^2,2)) ;
        evec = evec ./ [elen, elen] ;

    %------------------------------------- "voro"-type dist.
        vlen = sqrt(bal1(ref1,3));

    %------------------------------------- "size"-type dist.
        ihfn = fun0(conn(ref1,1));
        jhfn = fun0(conn(ref1,2));

    %------------------------------------- bind "safe" dist.
        ilen = min(vlen,ihfn) ;
        jlen = min(vlen,jhfn) ;

    %------------------------------------- locate offcentres
        inew = vert(conn(ref1,1),:) ...
             + [ilen,ilen].*evec ;
        jnew = vert(conn(ref1,2),:) ...
             - [jlen,jlen].*evec ;

    %------------------------------------- iter. "size"-type
        for ioff = +1 : +3
    %------------------------------------- eval. length-fun.
        if (~isempty(hfun))
            if (isnumeric(hfun))
            iprj = hfun * ...
              ones(size(inew,1),1);
            jprj = hfun * ...
              ones(size(jnew,1),1);
            else
            iprj = feval( ...
                hfun,inew,harg{:});
            jprj = feval( ...
                hfun,jnew,harg{:});
            iprj = iprj(:);
            jprj = jprj(:);
            end
        else
            iprj = +inf * ...
              ones(size(inew,1),1);
            jprj = +inf * ...
              ones(size(jnew,1),1);
        end

        iprj = 0.5*ihfn + 0.5*iprj;
        jprj = 0.5*jhfn + 0.5*jprj;

    %------------------------------------- bind "safe" dist.
        ilen = min(vlen,iprj) ;
        jlen = min(vlen,jprj) ;

    %------------------------------------- locate offcentres
        inew = vert(conn(ref1,1),:) ...
             + [ilen,ilen].*evec ;
        jnew = vert(conn(ref1,2),:) ...
             - [jlen,jlen].*evec ;

        end

    %------------------------------------- merge i,j if near
        near = ...
            ilen+jlen>=vlen*ntol ;

        znew = inew(near,:) * .5 ...
             + jnew(near,:) * .5 ;

        inew = inew(~near,1:2) ;
        jnew = jnew(~near,1:2) ;

    %------------------------------------- split constraints
        zset = (1:size(znew,1))' ...
                + size(vert,1) ;

        iset = (1:size(inew,1))' ...
                + size(znew,1) + ...
                + size(vert,1) ;

        jset = (1:size(jnew,1))' ...
                + size(znew,1) + ...
                + size(inew,1) + ...
                + size(vert,1) ;

        set1 = num1( near);
        set2 = num1(~near);

        cnew = [conn( set1,1), zset
                conn( set1,2), zset
                conn( set2,1), iset
                conn( set2,2), jset
                iset, jset ] ;
        conn = [conn(~ref1,:); cnew];

    %------------------------------------- update vertex set
        vert = [vert; znew(:,1:2)];
        vert = [vert; inew(:,1:2)];
        vert = [vert; jnew(:,1:2)];

        vidx = [zset; iset; jset] ;

        tcpu.offc = ...
            tcpu.offc + toc(ttic) ;

        end % switch(lower(opts.kind))

    end

    tcpu.full = ...
        tcpu.full + toc(tnow) ;

    if (~isinf(opts.disp) )
    %------------------------------------- print final stats
        numc = size(conn,1) ;
        numt = size(tria,1) ;
        fprintf(+1, ...
        '%11i %18i %18i\n', ...
        [iter,numc,numt]) ;
    end

    if (opts.dbug)
    %------------------------------------- print debug timer
        fprintf(1,'\n') ;
        fprintf(1,' 1-simplex REF. timer...\n');
        fprintf(1,'\n') ;
        fprintf(1, ...
        ' FULL: %f \n', tcpu.full);
        fprintf(1, ...
        ' BALL: %f \n', tcpu.ball);
        fprintf(1, ...
        ' HFUN: %f \n', tcpu.hfun);
        fprintf(1, ...
        ' ENCR: %f \n', tcpu.encr);
        fprintf(1, ...
        ' OFFC: %f \n', tcpu.offc);
        fprintf(1,'\n') ;
    end

end

function [vert,conn,tria,tnum,iter] = ...
            cdtref2(vert,conn,tria,tnum, ...
                node,PSLG,part,opts,hfun,harg,iter)
%CDTREF2 constrained Delaunay-refinement for 2-simplex elem-
%nts embedded in R^2.
%   [...] = CDTREF2(...) refines the set of 2-simplex eleme-
%   nts embedded in the triangulation until all constraints
%   are satisfied. Specifically, triangles are refined until
%   all local mesh-spacing and element-shape conditions are
%   met. Refinement proceeds according to either a Delaunay-
%   refinement or Frontal-Delaunay type approach, depending
%   on user-settings. In either case, new steiner points are
%   introduced to split "bad" triangles - those that violate
%   the set of prescribed constraints. In the "-DR" type pr-
%   ocess triangles are split about their circumballs. In
%   the "-FD" approach, new vertices are positioned such th-
%   at mesh-spacing and element-shape constraints are satis-
%   fied in a "locally-optimal" fashion.

    tcpu.full = +0. ;
    tcpu.dtri = +0. ;
    tcpu.tcon = +0. ;
    tcpu.ball = +0. ;
    tcpu.hfun = +0. ;
    tcpu.offc = +0. ;
    tcpu.filt = +0. ;

    vidx = (1:size(vert,1))';     %- "new" vert list to test

    tnow =  tic ;

    near = +.775;

    while (strcmpi(opts.ref2,'refine'))

        iter = iter + 1 ;

    %------------------------------------- build current CDT
        ttic = tic ;

        nold = size(vert,1) ;

       [vert,conn, ...
        tria,tnum]= deltri2(vert,conn, ...
                            node,PSLG, ...
                            part, ....
                            opts.dtri) ;

        nnew = size(vert,1) ;

        vidx = ...
       [vidx; (nold:nnew)'] ;

        % wjp: in case the new mesh has fewer nodes than old one
        vidx(vidx > nnew) = [];

        tcpu.dtri = ...
            tcpu.dtri + toc(ttic) ;

    %------------------------------------- build current adj
        ttic = tic ;

       [edge,tria]= tricon2(tria,conn) ;

        tcpu.tcon = ...
            tcpu.tcon + toc(ttic) ;

        if (iter>=opts.iter),break; end

    %------------------------------------- calc. circumballs
        ttic = tic ;

        bal1 = cdtbal1(vert,conn) ;
        bal2 = cdtbal2(vert, ...
                       edge,tria) ;
        len2 = minlen2(vert,tria) ;

        rho2 = bal2(:,+3) ./ len2 ;

    %------------------------------------- refinement scores
        scr2 = rho2 .* bal2(:,+3) ;

        tcpu.ball = ...
            tcpu.ball + toc(ttic) ;

    %------------------------------------- eval. length-fun.
        ttic = tic ;

        if (~isempty(hfun))
            if (isnumeric(hfun))
            fun0 = hfun * ...
              ones(size(vert,1),1);
            fun2 = hfun ;
            else
            fun0(vidx) = ...
                feval(hfun, ...
            vert(vidx,:), harg{:});
            fun0 = fun0(:) ;
            fun2 = fun0(tria(:,1))...
                 + fun0(tria(:,2))...
                 + fun0(tria(:,3));
            fun2 = fun2 / +3. ;
            end
        else
            fun0 = +inf * ...
              ones(size(vert,1),1);
            fun2 = +inf ;
        end

        siz2 = ...
         +3. * bal2(:,3)./(fun2.*fun2) ;

        tcpu.hfun = ...
            tcpu.hfun + toc(ttic) ;

    %------------------------------------- refinement queues
        ref1 = false(size(conn,1),1);
        ref2 = false(size(tria,1),1);

        stri = isfeat2(vert,edge,tria) ;

        ref2(rho2>opts.rho2* ...        %- bad rad-edge len.
                  opts.rho2) = true ;
        ref2(stri) = false ;
        ref2(siz2>opts.siz2* ...        %- bad equiv. length
                  opts.siz2) = true ;

        num2 = find(ref2);

    %------------------------------------- dump-out progess!
        if (mod(iter,opts.disp)==0)
            numc = size(conn,1) ;
            numt = size(tria,1) ;
            fprintf(+1, ...
            '%11i %18i %18i\n', ...
            [iter,numc,numt]) ;
        end

    %------------------------------------- nothing to refine
        if (isempty(num2)), break; end

       [scr2,idx2] = sort( ...
            scr2(num2),'descend');
        num2 = num2(idx2);

    %------------------------------------- refine "bad" tria
        switch (lower(opts.kind))
        case 'delaunay'
    %------------------------------------- do circ-ball pt's
        new2 = zeros(length(num2),3);
        new2(:,1:2) = bal2(num2,1:2);

        rmin = ...                      %- min. insert radii
            len2(num2)*(1.-eps^.75)^2 ;

        new2(:,  3) = max( ...
            bal2(num2,3)*near^2,rmin) ;

        case 'delfront'
    %-- off-centre scheme -- refine triangles by positioning
    %-- new vertices along a local segment of the voronoi
    %-- diagram, bounded by assoc. circmballs. New points
    %-- are placed to satisfy the worst of local mesh-length
    %-- and element-shape constraints.

        ttic = tic ;

    %------------------------------------- find frontal edge
       [lmin,emin] = ...
            minlen2(vert,tria(num2,:)) ;

        ftri = false(length(num2),1) ;
        epos = zeros(length(num2),1) ;
        tadj = zeros(length(num2),1) ;

        for ii = +1 : length(epos)
            epos(ii) = tria( ...
                num2(ii),emin(ii)+3) ;
        end

    %------------------------------------- find frontal tria
        for enum = +1 : +3

            eidx = tria(num2,enum+3) ;

            ftri = ...
            ftri | edge(eidx,5) > +0 ;

            ione = ...
                num2 ~= edge(eidx,3) ;
            itwo = ~ione ;

            tadj(ione) = ...
                edge(eidx(ione),3);
            tadj(itwo) = ...
                edge(eidx(itwo),4);

            okay = tadj > +0 ;
            tidx = tadj(okay);

            ftri(okay) = ...
            ftri(okay) | ~ref2(tidx) ;

        end

        if (~any(ftri))                 %- can this happen!?
        ftri = true(length(num2),+1) ;
        end

    %------------------------------------- locate offcentres
        emid = vert(edge(epos,+1),:) ...
             + vert(edge(epos,+2),:) ;
        emid = emid * +0.50 ;

        elen = sqrt(lmin(:));

    %------------------------------------- "voro"-type dist.
        vvec = bal2(num2,1:2)-emid ;
        vlen = sqrt(sum(vvec.^2,2));
        vvec = vvec ./ [vlen,vlen] ;

        hmid = fun0(edge(epos,+1),:) ...
             + fun0(edge(epos,+2),:) ;
        hmid = hmid * +0.50 ;

    %------------------------------------- "ball"-type dist.
        rtri = elen * opts.off2 ;
        rfac = elen * +0.50 ;
        dsqr = rtri.^2 - rfac.^2;
        doff = rtri + ...
            sqrt(max(+0.,dsqr)) ;

    %------------------------------------- "size"-type dist.
        dsiz = +sqrt(3.)/2. * hmid ;

    %------------------------------------- bind "safe" dist.
       [dist,ioff] = ...
          min([dsiz,doff,vlen],[],2) ;

    %------------------------------------- locate offcentres
        off2 = ...
        emid + [dist,dist] .* vvec ;

    %------------------------------------- iter. "size"-type
        for isub = +1 : +3
    %------------------------------------- eval. length-fun.
        if (~isempty(hfun))
            if (isnumeric(hfun))
            hprj = hfun * ...
              ones(size(off2,1),1) ;
            else
            hprj = feval( ...
                hfun,off2,harg{:}) ;
            hprj = hprj(:) ;
            end
        else
            hprj = +inf * ...
              ones(size(off2,1),1) ;
        end

    %------------------------------------- "size"-type dist.
        hprj = .33*hmid + .67*hprj ;

        dsiz = +sqrt(3.)/2. * hprj ;

        dsiz(dsiz<elen*.50) = +inf ;    %- edge-ball limiter
        dsiz(dsiz>vlen*.95) = +inf ;    %- circ-ball limiter

    %------------------------------------- bind "safe" dist.
       [dist,ioff] = ...
          min([dsiz,doff,vlen],[],2) ;

    %------------------------------------- locate offcentres
        off2 = ...
        emid + [dist,dist] .* vvec ;

        end

        orad = ...
        sqrt((elen*.5).^2 + dist.^2) ;

    %------------------------------------- do offcentre pt's
        new2 = ...
        zeros(length(find(ftri)),+3) ;
        new2(:,1:2) = off2(ftri,1:2) ;

        rmin = ...                      %- min. insert radii
            lmin(ftri)*(1.-eps^.75)^2 ;

        new2(:,  3) = max( ...
            (orad(ftri)*near).^2,rmin);

        tcpu.offc = ...
            tcpu.offc + toc (ttic) ;

        end % switch(lower(opts.kind))

    %------------------------------------- inter.-ball dist.
        ttic = tic ;

    %------------------------------------- proximity filters
       [vp,vi] = ...
          findball(new2,new2(:,1:2)) ;

        keep = true (size(new2,1),1) ;
        for ii = size(vp,1):-1:+1
            for ip = vp(ii,1) ...
                   : vp(ii,2)
                jj = vi(ip);
                if (keep(jj) && ...
                    keep(ii) && ...
                    jj < ii )

                keep(ii) = false ;
                break;

                end
            end
        end

        new2 = new2(keep,:);

    %------------------------------------- test encroachment
        bal1(:,3) = ...
            (1.-eps^.75) * bal1(:,3);

       [vp,vi] = ...
          findball(bal1,new2(:,1:2));

        keep = true (size(new2,1),1);
        for ii = +1:+1:size(vp,1)
            for ip = vp(ii,1) ...
                   : vp(ii,2)
                jj = vi(ip);
                ref1(jj) =  true ;
                keep(ii) = false ;
            end
        end

    %------------------------------------- leave sharp edges
        ebnd = false(size(edge,1),1);
        ebnd(tria(stri,4:6)) = true ;

        enot = ...
        setset2(conn,edge(ebnd,1:2));

        ref1(enot) = false ;

    %------------------------------------- preserve boundary
        if (strcmpi(opts.ref1,...
            'preserve'))
        ref1(:)    = false ;
        end

    %------------------------------------- refinement points
        new2 = new2(keep,:);
        new1 = bal1(ref1,:);

        tcpu.filt = ...
            tcpu.filt + toc(ttic) ;

    %------------------------------------- split constraints
        idx1 = ...
       (1:size(new1))'+size(vert,1) ;

        idx2 = ...
       (1:size(new2))'+size(new1,1) ...
                      +size(vert,1) ;

        cnew = [conn( ref1,1), idx1
                conn( ref1,2), idx1];
        conn = [conn(~ref1,:); cnew];

        vidx = [idx1; idx2];

    %------------------------------------- update vertex set
        nold = size(vert,1);
        vert = [vert; new1(:,1:2)];
        vert = [vert; new2(:,1:2)];
        % wjp: make sure vertices are unique
        vert = unique(vert,'rows','stable');
        nnew = size(vert,1);

        if (nnew == nold), break; end   %- we *must* be done

    end

    tcpu.full = ...
        tcpu.full + toc(tnow) ;

    if (~isinf(opts.disp) )
    %------------------------------------- print final stats
        numc = size(conn,1) ;
        numt = size(tria,1) ;
        fprintf(+1, ...
        '%11i %18i %18i\n', ...
        [iter,numc,numt]) ;
    end

    if (opts.dbug)
    %------------------------------------- print debug timer
        fprintf(1,'\n') ;
        fprintf(1,' 2-simplex REF. timer...\n');
        fprintf(1,'\n') ;
        fprintf(1, ...
        ' FULL: %f \n', tcpu.full);
        fprintf(1, ...
        ' DTRI: %f \n', tcpu.dtri);
        fprintf(1, ...
        ' TCON: %f \n', tcpu.tcon);
        fprintf(1, ...
        ' BALL: %f \n', tcpu.ball);
        fprintf(1, ...
        ' HFUN: %f \n', tcpu.hfun);
        fprintf(1, ...
        ' OFFC: %f \n', tcpu.offc);
        fprintf(1, ...
        ' FILT: %f \n', tcpu.filt);
        fprintf(1,'\n') ;
    end

end

function [opts] = makeopt(opts)
%MAKEOPT setup the options structure for REFINE2.

    if (~isfield(opts,'dtri'))
        opts.dtri = 'constrained';
    else
    if (~strcmpi(opts.dtri, 'conforming') && ...
        ~strcmpi(opts.dtri,'constrained') )
        error( ...
    'refine2:invalidOption','Invalid constraint DTRI.');
    end
    end

    if (~isfield(opts,'kind'))
        opts.kind = 'delfront';
    else
    if (~strcmpi(opts.kind, 'delfront') && ...
        ~strcmpi(opts.kind, 'delaunay') )
        error( ...
    'refine2:invalidOption','Invalid refinement KIND.');
    end
    end

    if (~isfield(opts,'ref1'))
        opts.ref1 = 'refine';
    else
    if (~strcmpi(opts.ref1,   'refine') && ...
        ~strcmpi(opts.ref1, 'preserve') )
        error( ...
    'refine2:invalidOption','Invalid refinement REF1.');
    end
    end

    if (~isfield(opts,'ref2'))
        opts.ref2 = 'refine';
    else
    if (~strcmpi(opts.ref2,   'refine') && ...
        ~strcmpi(opts.ref2, 'preserve') )
        error( ...
    'refine2:invalidOption','Invalid refinement REF2.');
    end
    end

    if (~isfield(opts,'iter'))
        opts.iter = +inf;
    else
    if (~isnumeric(opts.iter))
        error('refine2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.iter)~= +1)
        error('refine2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.iter <= +0)
        error('refine2:invalidOptionValues', ...
            'Invalid OPT.ITER selection.') ;
    end
    end

    if (~isfield(opts,'disp'))
        opts.disp = +10 ;
    else
    if (~isnumeric(opts.disp))
        error('refine2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.disp)~= +1)
        error('refine2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.disp <= +0)
        error('refine2:invalidOptionValues', ...
            'Invalid OPT.DISP selection.') ;
    end
    end

    if (~isfield(opts,'rho2'))
        opts.rho2 = 1.025;
    else
    if (~isnumeric(opts.rho2))
        error('refine2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.rho2)~= +1)
        error('refine2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.rho2 < +1.)
        error('refine2:invalidOptionValues', ...
            'Invalid OPT.RHO2 selection.') ;
    end
    end

    if (~isfield(opts,'off2'))
        opts.off2 = 0.933;
    else
    if (~isnumeric(opts.off2))
        error('refine2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.off2)~= +1)
        error('refine2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.off2 < +.7)
        error('refine2:invalidOptionValues', ...
            'Invalid OPT.OFF2 selection.') ;
    end
    end

    if (~isfield(opts,'siz1'))
        opts.siz1 = 1.333;
    else
    if (~isnumeric(opts.siz1))
        error('refine2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.siz1)~= +1)
        error('refine2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.siz1 <= 0.)
        error('refine2:invalidOptionValues', ...
            'Invalid OPT.SIZ1 selection.') ;
    end
    end

    if (~isfield(opts,'siz2'))
        opts.siz2 = 1.300;
    else
    if (~isnumeric(opts.siz2))
        error('refine2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.siz2)~= +1)
        error('refine2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.siz2 <= 0.)
        error('refine2:invalidOptionValues', ...
            'Invalid OPT.SIZ2 selection.') ;
    end
    end

    if (~isfield(opts,'dbug'))
        opts.dbug = false;
    else
    if (~islogical(opts.dbug))
        error('refine2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.dbug)~= +1)
        error('refine2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    end

end

author: CHLNDDEV


text:
function [ResultsAllCellLines] = performPPP(ResultsAllCellLines,mets,step_size,samples,step_num,direct)
% This function performs the a phase plane analysis. The analysis starts from zero and proceed `step_num*step_size` in the specified direction and for two exchanges.
% the results of the analysis will be saved into `ResultsAllCellLines`.
%
% USAGE:
%
%    [ResultsAllCellLines] = performPPP(ResultsAllCellLines, mets, step_size, samples, step_num, direct)
%
% INPUTS:
%    ResultsAllCellLines:
%    samples:                 Conditions
%    mets:                    Set of two metabolites to test e.g.  mets ={`EX_glc(e)`, `EX_gln(e)`; `EX_o2(e)`, `EX_o2(e)`}
%    step_size:               Step size for metabolites specified in mets e.g., `step_size` = 100
%    direct:                  Direction for metabolites specified in mets: uptake (-1) or secretion (1), e.g., direct = [-1,-1,-1,-1];
%    step_num:                Number of steps (1000/20 = 50), step_num = [5,5;5,5]
%
% OUTPUT:
%    ResultsAllCellLines:     The matrix of growth rates are added to the `ResultsAllCellLines` structure alond with the values or the bounds using the same of the manipulated exchanges
%
% .. Author: - Maike K. Aurich 08/07/15 (performed as in Jeff Orth et al. 2010, Supplemental tutorial.)

bound_size = max(max(step_num));
for k=1:length(samples)
    k
    model1 = eval(['ResultsAllCellLines.' samples{k} '.modelPruned']);

    for p=1:size(mets,1)
        bounds = zeros(bound_size,2); % allocate space to fill in bounds
        growthRates=zeros(step_num(p,1),step_num(p,2));
        model=model1;

        for i=1:step_num(p,1)
            if direct(p,1)==1
                bounds(i,1) = direct(p,1)*i*step_size(p,1)- step_size(p,1);
                model=changeRxnBounds(model,mets{p,1},bounds(i,1),'b');
            else
                bounds(i,1) = direct(p,1)*i*step_size(p,1)+ step_size(p,1);
                model=changeRxnBounds(model,mets{p,1},bounds(i,1),'b');
            end

            for j=1:step_num(p,2)
               if direct(p,2)==1
                bounds(j,2) = direct(p,2)*j*step_size(p,2)- step_size(p,2);
                model=changeRxnBounds(model,mets{p,2},bounds(j,2),'b');
               else
                   bounds(j,2) = direct(p,2)*j*step_size(p,2)+ step_size(p,2);
                model=changeRxnBounds(model,mets{p,2},bounds(j,2),'b');
               end
                FBAsolution=optimizeCbModel(model,'max');
                growthRates(i,j)=FBAsolution.f;
                clear FBAsolution
            end
        end
        if contains(mets{p,1},'[')
            name  = ['phasePlane_' strtok(mets{p,1}, '[') '_' strtok(mets{p,2}, '[')];
        else
            name  = ['phasePlane_' strtok(mets{p,1}, '(') '_' strtok(mets{p,2}, '(')];
        end

        ResultsAllCellLines.(samples{k}).(name).growthRates= growthRates;
        ResultsAllCellLines.(samples{k}).(name).bounds =  bounds;
        clear bounds model growthRates
    end

    clear model1
    save('PPP120', '-v7.3');
end
end

author: opencobra

repo: cobratoolbox


text:
function [chunk,p] = onl_filtered(p,desired_length,suppress_output,set_online_scope)
% Obtain processed data from a filter pipeline online.
% [Chunk,Pipeline] = onl_filtered(Pipeline,DesiredLength,SuppressOutput,SetOnlineScope)
%
% This function returns a chunk of most recent filtered output from a filter pipeline data structure
% previously created with onl_newpipeline.
%
% A filter pipeline is a recursive data structure (like a tree) whose nodes are the filter stages and
% whose edges represent the output of one stage going into the input of another stage. The leaf
% nodes refer to raw online streams (structs in the workspace) and are queried via onl_peek, all other
% nodes are evaluated by calling the filter function on its input data (recursively).
%
% In:
%   Pipeline : filter pipeline struct to update and to read data from.
%
%   DesiredLength : number of samples to get (or 0 to get all new samples) (default: 0)
%
%   SuppressOutput : suppress console output (default: true)
%
%   SetOnlineScope : set the regular online-processing scope (can be turned off for efficiency if
%                    that scope is already set for some reason) (default: true)
%
% Out:
%   Chunk : EEGLAB dataset struct representing the desired data chunk
%           Can be shorter than desired length if not enough data is available at the moment; if the
%           chunk is epoched, the desired length is ignored.
%
%   Pipeline : updated filter pipeline struct
%
%
% Example:
%   % load calibration set
%   raw = io_loadset('calib.set')
%
%   % apply a series of filter to it (the processed set now has a filter expression and initial state)
%   processed = exp_eval(flt_iir(flt_resample(raw,128),[0.5 1],'highpass'));
%
%   % start streaming some data
%   run_readdataset('mystream','action.set');
%   % and put a pipeline on top of it that replicates the processing applied to processed and continues it on new data
%   pip = onl_newpipeline(processed,{'mystream'});
%
%   while 1
%      % generate a 200-sample view into the processed stream
%      [EEG,pip] = onl_filtered(pip,200);
%   end
%
% See also:
%   onl_newpipeline, onl_newstream, onl_append, onl_peek
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2012-05-13

% check inputs
if nargin < 4
    set_online_scope = true;
    if nargin < 1
        error('You need to pass a pipeline structure to onl_predict.'); end
    if nargin < 2
        desired_length = 0;
    elseif ~(isnumeric(desired_length) && isscalar(desired_length)) || islogical(desired_length)
        error('The given DesiredLength argument must be a numeric scalar, but was: %s',hlp_tostring(desired_length,10000));
    end
    if nargin < 3
        suppress_output = true;
    elseif ~islogical(suppress_output)
        error('The given SuppressOutput argument must be logical (true or false), but was: %s',hlp_tostring(suppress_output,10000));
    end
else
    if ~(isnumeric(desired_length) && isscalar(desired_length)) || islogical(desired_length)
        error('The given DesiredLength argument must be a numeric scalar, but was: %s',hlp_tostring(desired_length,10000)); end
    if ~islogical(suppress_output)
        error('The given SuppressOutput argument must be logical (true or false), but was: %s',hlp_tostring(suppress_output,10000)); end
    if ~islogical(set_online_scope)
        error('The given SetOnlineScope argument must be logical (true or false), but was: %s',hlp_tostring(set_online_scope,10000)); end
end

if ~all(isfield(p,{'head','parts','subnodes'}))
    error('The given Pipeline argument is not a valid filter pipeline (must be a struct with fields .head, .parts and .subnodes).'); end

% run update_pipeline() with appropriate options
if suppress_output
    if set_online_scope
        [console_output,chunk,p] = evalc('hlp_scope({''disable_expressions'',1,''is_online'',1},@update_pipeline,p)'); %#ok<ASGLU>
    else
        [console_output,chunk,p] = evalc('update_pipeline(p)'); %#ok<ASGLU>
    end
else
    if set_online_scope
        [chunk,p] = hlp_scope({'disable_expressions',1,'is_online',1},@update_pipeline,p);
    else
        [chunk,p] = update_pipeline(p);
    end
end

% if a desired length was specified and the chunk is not epoched
if desired_length && ~isstruct(chunk.epoch)
    try
        % concatenate the chunk with previous buffer contents and trim to desired size
        p.buffer = utl_buffer(chunk,p.buffer,desired_length);
    catch e
        % handle missing .buffer field (first-time use)
        if ~isfield(p,'buffer')
            p.buffer = utl_buffer(chunk,struct('data',[],'event',[]),desired_length);
        else
            rethrow(e);
        end
    end
    chunk = p.buffer;
end

function [chunk,p] = update_pipeline(p)
% Update the given filter pipeline and get a chunk of the newly appended output
% [Chunk,Pipeline] = update_pipeline(Pipeline)
%
% A pipeline is a recursive data structure (like a tree) whose nodes are the filter stages and
% whose edges represent the output of one stage going into the input of another stage. The leaf
% nodes refer to raw data streams (structs in the workspace) and are queried via onl_peek, all other
% nodes are evaluated by calling the filter function on its input data (recursively).
%
% At each node we store the filter function (.head) and its arguments (.parts), some of which may be
% input filter pipelines themselves, plus some miscellaneous book-keeping data. These include:
%  * .israw : true for nodes that represent raw data
%  * .pipeline_indices : indices of those input arguments that are pipelines themselves (if not raw)
%  * .stateful : true if the node has state
%  * .state : previous filter state, if stateful
%  * .smax : number of samples seen so far (if israw)
%
% In:
%   Pipeline : previous filter pipeline struct
%
% Out:
%   Chunk : EEGLAB dataset struct representing newly appended data, filtered
%
%   Pipeline : updated filter pipeline struct

inputs = p.parts;
if p.subnodes
    % update input pipelines to the current node and store the results in inputs
    for k=p.subnodes
        [inputs{k},p.parts{k}] = update_pipeline(p.parts{k}); end
    % process the inputs by calling the respective filter function
    if p.stateful
        [chunk,p.state] = p.head(inputs{:},'state',p.state,'arg_direct',true);
    else
        chunk = p.head(inputs{:});
    end
else
    % get the most recent samples since our buffer's smax from a raw stream:
    % inputs holds the cell array {stream_name,channel_range)
    chunk = onl_peek(inputs{1},p.smax,'index',inputs{2});
    p.smax = chunk.smax;
end

repo: REST


text:
function s = s_sort_a ( s )

%*****************************************************************************80
%
%% S_SORT_A sorts a string into ascending order.
%
%  Discussion:
%
%    The string is assumed to be short, and so a simple bubble sort is used.
%
%    ALL the characters are sorted, including blanks and punctuation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be sorted.
%
%    Output, string S, the sorted string.
%
  s_length = length ( s );

  for i = 1 : s_length - 1

    c = s(i);
    j = i;

    for k = i + 1 : s_length
      if ( s(k) < s(j) )
        j = k;
      end
    end

    s(i) = s(j);
    s(j) = c;

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function rts_ProcInformationUpdtFunc(~,~, hRtsMainGUI, hTable)
%rts_ProcInformationUpdtFunc Summary of this function goes here
%   Detailed explanation goes here
    timers = getappdata(hRtsMainGUI,'timers');

	tableData = {};
    try
        validtimers = timer.empty(1,0);
        for(i=1:length(timers))
            if(isvalid(timers(i)))
                validtimers(end+1) = timers(i);
            else
                delete(timers(i));
            end
        end

        for(i=1:length(validtimers))
            task = validtimers(i);

            if(isvalid(task))
                tableData{i,1} = get(task,'Name');
                tableData{i,2} = cap1stLetter(get(task,'Running'));
                tableData{i,3} = cap1stLetter(get(task,'ExecutionMode'));
                tableData{i,4} = num2str(get(task,'Period'));
            else
                delete(task);
            end
        end

        setTableData(hTable, tableData)
    catch ME

    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function TAU = int_tau(Z)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Program:    Tau Integral Computation for Hydrogen Composition Program
%               in Atmospheric Model
%   Author:     Brent Lewis(RocketLion@gmail.com)
%               University of Colorado-Boulder
%   History:    Original-1/10/2007
%   Input:      Z:      Altitude value
%   Output:     TAU:    Integral Value
%   Note:       This program computes the value of Tau directly with the
%               integral done by hand and only the second integration limit
%               needing to be inputed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   Constants
L_K_9 = 12;
T_10 = 360;
T_inf = 1000;
Z_10 = 120;
g_0 = 9.80665;
r_E = 6.356766e3;
R = 8.31432e3;
lambda = L_K_9/(T_inf-T_10);
M_H = 1.00797;

%   Value of Integration limit computed previously
tau_11 = 8.329503912749350e-004;

tau_Z = M_H*g_0*r_E^2/R*...
    log((exp(lambda*(Z-Z_10)*(r_E+Z_10)/(r_E+Z))-1)*T_inf+T_10)/...
    (lambda*T_inf*(r_E+Z_10)^2);

TAU = tau_Z-tau_11;

author: Sable

repo: mcbench-benchmarks


text:
function f = local ( x )

%*****************************************************************************80
%
%% LOCAL computes the local function.
%
%  Discussion:
%
%    This function has a local minimizer:
%
%      X* = ( 0.2858054412..., 0.2793263206...), F(X*) = 5.9225...
%
%    and a global minimizer:
%
%      X* = ( -21.02667179..., -36.75997872...), F(X*) = 0.
%
%    Suggested starting point for local minimizer:
%
%      X = ( 1, 1 ), F(X) = 3.33 * 10^6.
%
%    Suggested starting point for global minimizer:
%
%      X = ( -15, -35), F(X) = 1.49 * 10^8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 January 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    David Himmelblau,
%    Applied Nonlinear Programming,
%    McGraw Hill, 1972,
%    ISBN13: 978-0070289215,
%    LC: T57.8.H55.
%
%  Parameters:
%
%    Input, real X(2), the argument of the function.
%
%    Output, real F, the value of the function at X.
%
  if ( length ( x ) ~= 2 )
    error ( 'Error: function expects a two dimensional input\n' );
  end

  f = ( x(1)^2 + 12 * x(2) - 1 )^2 ...
    + ( 49 * x(1)^2 + 49 * x(2)^2 + 84 * x(1) + 2324 * x(2) - 681 )^2;

  return
end


text:
function cdf = pareto_cdf ( x, a, b )

%*****************************************************************************80
%
%% PARETO_CDF evaluates the Pareto CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the CDF.
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < A,
%    0.0 < B.
%
%    Output, real CDF, the value of the CDF.
%
  if ( x < a )
    cdf = 0.0;
  else
    cdf = 1.0 - ( a / x )^b;
  end

  return
end


text:
function sizes = getVarSizes(obj, inputSizes)
%GETVARSIZES  Get the size of the variables
%   SIZES = GETVARSIZES(OBJ, INPUTSIZES) computes the SIZES of the
%   DagNN variables given the size of the inputs. `inputSizes` is
%   a cell array of the type `{'inputName', inputSize, ...}`
%   Returns a cell array with sizes of all network variables.
%
%   Example, compute the storage needed for a batch size of 256 for an
%   imagenet-like network:
%   ```
%   batch_size = 256; single_num_bytes = 4;
%   input_size = [net.meta.normalization.imageSize, batch_size];
%   var_sizes = net.getVarSizes({'data', input_size});
%   fprintf('Network activations will take %.2fMiB in single.\n', ...
%     sum(prod(cell2mat(var_sizes, 1))) * single_num_bytes ./ 1024^3);
%   ```

% Copyright (C) 2015 Andrea Vedaldi, Karel Lenc.
% 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).

nv = numel(obj.vars) ;
sizes = num2cell(NaN(nv, 4),2)' ;

for i = 1:2:numel(inputSizes)
  v = obj.getVarIndex(inputSizes{i}) ;
  if isnan(v)
    error('Variable `%s` not found in the network.', inputSizes{i});
  end;
  if isempty(inputSizes{i+1})
    sizes{v} = [0 0 0 0] ;
  else
    sizes{v} = [inputSizes{i+1}(:)' ones(1, 4 - numel(inputSizes{i+1}))] ;
  end
end

for layer = obj.layers(obj.executionOrder)
  in = layer.inputIndexes ;
  out = layer.outputIndexes ;
  sizes(out) = layer.block.getOutputSizes(sizes(in)) ;
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function mnt= mnt_setElectrodePositions(clab, varargin)
%MNT_SETELECTRODEPOSITIONS - Electrode positions of standard named channels
%
%Synopsis:
% MNT= mnt_setElectrodePositions(CLAB, <OPT>);
%
%Input:
% CLAB: Label of channels (according to the extended international
%       10-10 system, see mntutil_posExt1010).
% OPT: Struct or property/value list of optional properties:
%  .PositionFcn - FUNC handle of function that specifies eletrode positions,
%                 default @mntutil_posExt1010
%
%Output:
% MNT: Struct for electrode montage
%  .x     - x coordiante of electrode positions
%  .y     - y coordinate of electrode positions
%  .clab  - channel labels
%
%See also mnt_setGrid

props= {'PositionFcn'   @mntutil_posExt1010   'FUNC'
        };
opt= opt_proplistToStruct(varargin{:});
opt= opt_setDefaults(opt, props, 1);

posSystem= opt.PositionFcn();
x= posSystem.x;
y= posSystem.y;
z= posSystem.z;
elab= posSystem.clab;

maz= max(z(:));
miz= min(z(:));
% ALTERNATIVE PROJECTION:
% This function works with an input of mnt which assumes that the human head "is" a ball with radius 1.
% The lower section of this ball is first projected onto the walls of a cylinder (orthogonal, with radius 1);
% then all (new) points will be projected on the 2d-hyperspace with z=maz.
%ur= [0 0 miz-0.8*(maz-miz)];
% ur is the center of projection. This is why the function only works with radius = 1
% and ur(1:2) = [0 0].
ur= [0 0 -1.5];
if 1==0% old projection. uses center of projection.
  la= (maz-ur(3)) ./ (z(:)-ur(3));
  Ur= ones(length(z(:)),1)*ur;
  % Project the lower halfball onto the wall of the cylinder:
  cx = x;
  cy = y;
  cx(z<0) = cx(z<0)./sqrt(cx(z<0).^2+cy(z<0).^2);
  cy(z<0) = cy(z<0)./sqrt(cx(z<0).^2+cy(z<0).^2);
  % Project everything onto the plane {z = maz}:
  pos2d= Ur + (la*ones(1,3)) .* ([cx(:) cy(:) z(:)] - Ur);
  pos2d= pos2d(:, 1:2);
  %pos2d(z<0,:)= NaN;% TODO: don't throw away the values of the lower halfball!
end

% This projection uses the distance on the "head"-surface to determine the 2d-positions of the electrodes w.r.t. Cz.
alpha = asin(sqrt(x.^2 + y.^2));
stretch = 2-2*abs(alpha)/pi;
stretch(z>0) = 2*abs(alpha(z>0))/pi;
norm = sqrt(x.^2 + y.^2);
norm(norm==0) = 1;
cx = x.*stretch./norm;
cy = y.*stretch./norm;
pos2d = [cx(:) cy(:)];

nChans= length(clab);
mnt.x= NaN*ones(nChans, 1);
mnt.y= NaN*ones(nChans, 1);
mnt.pos_3d= NaN*ones(3, nChans);
for ei= 1:nChans,
  ii= util_chanind(elab, clab{ei});
  if ~isempty(ii),
    mnt.x(ei)= pos2d(ii, 1);
    mnt.y(ei)= pos2d(ii, 2);
    mnt.pos_3d(:,ei)= [x(ii) y(ii) z(ii)];
  end
end
radius = 1.3;
%radius= 1.9;
mnt.x= mnt.x/radius;
mnt.y= mnt.y/radius;
mnt.clab= clab;

author: bbci


text:
function [nodes, edges] = bnMsgPassCreate(M, values, CPT)
% BNMSGPASSCREATE helper function for lungbayesdemo

% Reference: Neapolitan R., "Learning Bayesian Networks", Pearson Prentice Hall,
% Upper Saddle River, New Jersey, 2004.

%=== create a dummy structure for nodes
dummy1.id = [];
dummy1.values = [];
dummy1.parents = [];
dummy1.children = [];
dummy1.peye = [];
dummy1.lambda = [];
dummy1.CPT = [];
dummy1.P = [];

%=== create a dummy structure for edges
dummy2.peyeX = [];
dummy2.lambdaX = [];

%=== create nodes
N = size(M,1); % number of nodes
nodes = repmat(dummy1, N, 1);

%=== create edges
edges = repmat(dummy2, size(M));

%=== populate nodes with data
for i = 1:N
    nodes(i).id = i;
    nodes(i).parents = find(M(:,i));
    nodes(i).children = find(M(i,:));
    nodes(i).CPT = CPT{i};
    nodes(i).values = values{i};
end


text:
function madReadings = removeOLMADx5(readings)
%This function removes outliers from the data set using the median
%absolute deviation method with a x5 factor. The formula used is
%MAD = Median { | y(i) - m | / 0.6745 }

med = median(readings); %get median value of readings
x = abs((readings - med)); %get absolute value of each reading minus median value
MADx5 = median(x/.6745)*5; %get value of x array and multiply it by 5, result is outlier factor
j = 1; %variable for non-outlier array
madReadings = zeros(1,length(readings));%allocate array for storing readings with outliers removed

for i = 1:length(readings) %this loop removes outliears
    if readings(i) <= (med + MADx5) && readings(i) >= (med - MADx5)
        madReadings(j) = readings(i);%if not an outlier add to array
        j = j + 1;
    end

end

madReadings((j):length(madReadings)) = [];%delete empty array items


text:
function y = sample_node(CPD, pev)
% SAMPLE_NODE Draw a random sample from P(Xi | x(pi_i), theta_i)  (gaussian)
% y = sample_node(CPD, parent_evidence)
%
% pev{i} is the value of the i'th parent (if there are any parents)
% y is the sampled value (a scalar or vector)

if length(CPD.dps)==0
  i = 1;
else
  dpvals = cat(1, pev{CPD.dps});
  i = subv2ind(CPD.sizes(CPD.dps), dpvals(:)');
end

if length(CPD.cps) == 0
  y = gsamp(CPD.mean(:,i), CPD.cov(:,:,i), 1);
else
  pev = pev(:);
  x = cat(1, pev{CPD.cps});
  y = gsamp(CPD.mean(:,i) + CPD.weights(:,:,i)*x(:), CPD.cov(:,:,i), 1);
end
y = y(:);

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [vect] = normVector(vect)
%normVector Summary of this function goes here
%   Detailed explanation goes here
    vect = vect / norm(vect);
end


text:
function [] = profiler(data,lonlat,point1,point2,bins,R,units)

% By David Bekaert - December 2010

% give as input the data, the lonlat (2 colums), the profile points [lon lat], the number of
% bins and the radius [m] of the tube used to project points
% by default a 15 bins and a radius of 100m is used.
% for debug figures set debugflag to 1

debugflag=0;

if nargin<7
    units = 'rad';
end
if nargin<6
    R = 100;            % default value of R in m
end
if nargin<5
    bins = 15;          % default number of bins
end
if nargin<4
    error('Specify more input arguments \n');
end
fontsize = 15;          % label fontsize

%%% PLOTTING
profile_fig = figure('name','Cross-sectional plot');
axes('OuterPosition',[0  0  0.7  1])
scatter(lonlat(:,1),lonlat(:,2),3,data,'filled')
hold on
plot([point1(1,1) point2(1,1)]',[point1(1,2) point2(1,2)]','k-*','LineWidth',1.5)
cc =colorbar;
xlabel(cc,units,'fontsize',fontsize)
set(gca,'fontsize',fontsize)
xlabel('Longitude','fontsize',fontsize)
ylabel('Latitude','fontsize',fontsize)
%%% PLOTTING

% crossection computation
% transform to a local reference frame, unit becomes [km]
lonlatXY = llh2local(lonlat',point1);	% transformation of the PS pixels
P1 = llh2local(point1',point1);                 % transformation of the line points
P2 = llh2local(point2',point1);

% rotate the local from such that line becomes horizontal
alpha = atan((P2(2,1)-P1(2,1))/(P2(1,1)-P1(1,1)));		% angle [rad]
rot = [cos(alpha)    -sin(alpha)
       sin(alpha)    cos(alpha)];
XY_new = rot\lonlatXY;                  % X is the first row, Y the second row
P1_new = rot\P1;
P2_new = rot\P2;
clear rot alpha  P1 P2 lonlatXY

if debugflag==1
    figure('name','Debug_plot [1]');
    scatter(XY_new(1,:),XY_new(2,:),3,data,'filled')
    hold on
    scatter([P1_new(1) P2_new(1)],[P1_new(2) P2_new(2)],'ro')
end

% Search for all the PS within distance R
ix_data = find((abs(XY_new(2,:))<=R/1000));

% updating vectors
XY_new = XY_new(:,ix_data);
data_new = data(ix_data);

if debugflag==1
    figure('name','Debug_plot [2]');
    scatter(XY_new(1,:),XY_new(2,:),3,data_new,'filled')
    hold on
    scatter([P1_new(1) P2_new(1)],[P1_new(2) P2_new(2)],'ro')
end

% output to the screen
PS_used = size(ix_data,2);
if (isempty(PS_used)==1 || PS_used==0)
    error('No PS are found within the tube crossection. \n')
else
    fprintf([num2str(PS_used),' PS used to compute projection on crossection \n'])
end

% threshold for minimum PS inside a bin
ps_min = floor(0.5*(PS_used/bins));

% binning of the results
binsize = (max(XY_new(1,:))-min(XY_new(1,:)))/bins;
data_cross_binned = zeros([1 bins]);
bins_xy = zeros([1 bins]);
clear ix

for k=1:bins
    bl = (k-1)*binsize+min(XY_new(1,:));                        % lower bound
    bu = (k)*binsize+min(XY_new(1,:));                          % upper bound
    ix = find(bl<=XY_new(1,:)  & XY_new(1,:)<bu);
    % Only showing bins with a minimum of PS contained
    if size(ix,2)>=ps_min
        bins_xy(1,k) = bl+binsize/2;                            % position center of bin
        data_binned_mean(1,k) = mean(data_new(ix,1));     % take mean value

    else
        bins_xy(1,k) = NaN;
        data_binned_mean(1,k) = NaN;

    end
    clear ix bl bu
end
ix = find(isnan(bins_xy)==1);
if isempty(ix)==0
    bins_xy(ix)=[];
    data_binned_mean(ix)=[];
end

%%% PLOTTING
figure(profile_fig)
axes('OuterPosition',[0.65  0  0.3  1])
plot(data_new,XY_new(1,:),'.','color',[0.9 0.9 0.9])
ylabel('Distance along crossection [km]','fontsize',fontsize)
hold on
plot(data_binned_mean,bins_xy,'k.')
xlabel(units,'fontsize',fontsize)
set(gca,'fontsize',fontsize)
%%% PLOTTING

clear all

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function out = arcLength(f, a, b)
%ARCLENGTH   Compute the length of the arc defined by a CHEBFUN.
%   ARCLENGTH(F) returns the arc length of the curve defined by CHEBFUN F in
%   the x-y plane over the interval where it is defined. If F is complex,
%   the output is the arc length of the curve in the complex plane.
%
%   ARCLENGTH(F, A, B) returns the arc length of F over the interval [A, B],
%   where [A, B] is a subinterval of the domain in which F is defined. In the
%   case of complex-valued F, ARCLENGTH(F, A, B) computes the length of the
%   arc whose ends correspond to A and B.
%
%   If F is a quasimatrix, the arc length of each CHEBFUN in F is computed
%   and a vector is returned.
%
% Examples:
%
%   f = chebfun(@(x) sin(x), [0 1]); L = arcLength(f)
%
%   arcLength(chebfun('exp(1i*pi*s)'))/(2*pi)

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

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

% Check the first input for its type:
if ( ~isa(f, 'chebfun') )
    error('CHEBFUN:CHEBFUN:arcLength:Input', ...
        'The first argument must be a chebfun object.')
end

% Get the interval corresponding to the arc:
if ( nargin == 3 )
    % Full arguments:
    dom(1) = a;
    dom(2) = b;
elseif ( nargin == 2 )
    % Two arguments: the second argument is a vector:
    if ( max( size(a) ) ~= 2 )
        error('CHEBFUN:CHEBFUN:arcLength:Input', ...
            'The second argument must be a 1x2 vector.')
    end
    dom = a;
else
    % Single argument:
    dom = [f(1).domain(1) f(1).domain(end)];
end

% Loop over each column or row:
fPrime = f;
for i = 1:numel(f.funs)
    % This makes sure that no delta functions are
    % generated due to jump discontinuities. pointValues of fPrime are not
    % updated since they are not needed.
    fPrime.funs{i} = diff(f.funs{i});
end

if ( isreal(f) )
    g = sqrt(1 + fPrime.^2);
    out = sum(g, dom(1), dom(2));

else
    out = sum(abs(fPrime), dom(1), dom(2));
end

% Reform the output to accommodate the transposedness with the input:
if ( f(1).isTransposed )
    out = out.';
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [nstate] = tapas_sampler_dlinear_cv_gibbs_node(data, model, ...
    inference, state, node)
%% Samples from a linear model with variance using a gibbs step and a
% leave one out approach.
%
% Input
%
% Output
%

% aponteeduardo@gmail.com
% copyright (C) 2016
%

nstate = state;

[np, nc] = size(state.graph{node - 1}.y);
nm = numel(state.graph{node}.y{1}.mu);
% Compute the means

mu = cell(1, nc);

% Initilize at the prior

y = state.graph{node + 1}.y;

% First store the values somewhere
% np + 1 Number of subjects plus prior
values = zeros(size(y{1}.mu, 1), np + 1);
variance = cell(nc, 1);
for i = 1:nc
    values(:, 1) = y{i}.mu;
    for j = 1:np
        values(:, j + 1)  = state.graph{node - 1}.y{j, i};
    end
    y{i}.mu = mean(values, 2);
    % No bessel correction
    variance{i} = var(values, 1, 2);
end

for i = 1:nc
    % Priors
    alpha = state.graph{node + 1}.y{i}.alpha;
    beta = state.graph{node + 1}.y{i}.beta;
    % Sample variance for each of the components
    % The one degree of freedom comes from the prior
    pe = gamrnd(alpha + np/2, 1./(beta + ((np + 1)/2.0) .* variance{i}));
    nstate.graph{node}.y{i}.pe = pe;
    % Update the mean
    nstate.graph{node}.y{i}.mu = y{i}.mu + ...
        (1./sqrt((np + 1) * pe) .* randn(nm, 1));
end

nstate.llh{node} = model.graph{node}.llh(nstate.graph{node}, ...
    nstate.graph{node + 1}, model.graph{node}.htheta);

nstate.llh{node - 1} = model.graph{node - 1}.llh(...
    nstate.graph{node - 1}, nstate.graph{node}, model.graph{node - 1}.htheta);

end

author: translationalneuromodeling

repo: tapas


text:
function [Hx,Hy,Ez,time] = Maxwell2D(Hx, Hy, Ez, FinalTime)

% function [Hx,Hy,Ez] = Maxwell2D(Hx, Hy, Ez, FinalTime)
% Purpose :Integrate TM-mode Maxwell's until FinalTime starting with initial conditions Hx,Hy,Ez

Globals2D;
time = 0;

% Runge-Kutta residual storage
resHx = zeros(Np,K); resHy = zeros(Np,K); resEz = zeros(Np,K);

% compute time step size
rLGL = JacobiGQ(0,0,N); rmin = abs(rLGL(1)-rLGL(2));
dtscale = dtscale2D; dt = min(dtscale)*rmin*2/3

% outer time step loop
while (time<FinalTime)

  if(time+dt>FinalTime), dt = FinalTime-time; end

   for INTRK = 1:5
      % compute right hand side of TM-mode Maxwell's equations
      [rhsHx, rhsHy, rhsEz] = MaxwellRHS2D(Hx,Hy,Ez);

      % initiate and increment Runge-Kutta residuals
      resHx = rk4a(INTRK)*resHx + dt*rhsHx;
      resHy = rk4a(INTRK)*resHy + dt*rhsHy;
      resEz = rk4a(INTRK)*resEz + dt*rhsEz;

      % update fields
      Hx = Hx+rk4b(INTRK)*resHx; Hy = Hy+rk4b(INTRK)*resHy; Ez = Ez+rk4b(INTRK)*resEz;
   end;
   % Increment time
   time = time+dt;
end
return

repo: Aero-matlab


text:
function ExportAedat2(aedat)

%{
This function exports data to a .aedat file in format version 2.
The .aedat file format is documented here:
http://inilabs.com/support/software/fileformat/

Just merges polarity, frame and imu6 at the moment.

If there is a timestamp reset, events before and after it will be mixed
together. Make sure to cut out time before any resets, or otherwise
monotonise the timestamps for each data type separately before applying
this function.

%}

%% validation of parameters

dbstop if error

if ~exist('aedat', 'var')
    error('Missing input')
end

if ~isfield(aedat, 'exportParams') || ~isfield(aedat.exportParams, 'filePath')
    error('Missing parameter exportParams.filePath')
end

% For source, use an override if it has been given. This allows data from
% one sensor to masquerade as data from another sensor.

if isfield(aedat.exportParams, 'source')
    source = aedat.exportParams.source;
else
    source = aedat.info.source;
end

%% General preparation

% Create overall containers
allTimeStamps = uint32([]);
allSamples = uint32([]);

% Declare function for finding specific event types in eventTypes cell array
cellFind = @(string)(@(cellContents)(strcmp(string, cellContents)));

%% Polarity

if isfield(aedat.data, 'polarity') ...
        && (~isfield(aedat.exportParams, 'dataTypes') ...
        || any(cellfun(cellFind('polarity'), aedat.exportParams.dataTypes)))

    disp('Preparing polarity data ...')

    if strcmp(source, 'Dvs128')
        % In the 32-bit address:
        % bit 1 (1-based) is polarity
        % bit 2-8 is x
        % bit 9-15 is y
        % bit 16 is special

        yShiftBits = 8;
        xShiftBits = 1;
        polShiftBits = 0;
    else % Default to DAVIS
        % In the 32-bit address:
        % bit 32 (1-based) being 1 indicates an APS sample
        % bit 11 (1-based) being 1 indicates a special event
        % bits 11 and 32 (1-based) both being zero signals a polarity event

        yShiftBits = 22;
        xShiftBits = 12;
        polShiftBits = 11;
    end

    y = bitshift(uint32(aedat.data.polarity.y), yShiftBits);
    x = bitshift(uint32(aedat.data.polarity.x), xShiftBits);
    pol = bitshift(uint32(aedat.data.polarity.polarity), polShiftBits);
    allSamples = [allSamples; y + x + pol];
    allTimeStamps = [allTimeStamps; uint32(aedat.data.polarity.timeStamp)];
end

%% Frame

if isfield(aedat.data, 'frame') ...
        && (~isfield(aedat.exportParams, 'dataTypes') ...
        || any(cellfun(cellFind('frame'), aedat.exportParams.dataTypes)))
    disp('Preparing frame data ...')
    yShiftBits = 22;
    xShiftBits = 12;
    % frameShiftBits = 0;
    frameFlagShiftBits = 31;
    signalFlagShiftBits = 10;

    frameData = aedat.data.frame;

    numFrames = frameData.numEvents;
    xDim = aedat.info.deviceAddressSpace(1);
    yDim = aedat.info.deviceAddressSpace(2);
    numPixels = xDim * yDim;

    % Allocate horizontal vectors to hold output data.
    % Why are the vectors that big? The factor of 2 is because
    % we insert dummy 'reset' frames prior to each frame.
    samples = uint32(zeros(1, 2 * numFrames * numPixels));
    timeStamps = uint32(zeros(2 * numFrames * numPixels, 1));

    % The output vector is twice as big again because samples and timeStamps
    % will be interspersed in the 'output' vector.
    y = repmat(uint32(yDim - 1 : -1 : 0)', xDim * numFrames * 2, 1); % y ramps down
    x = repmat(uint32(0 : xDim - 1), yDim, numFrames * 2);           % but x ramps up
    x = x(:);
    % in bit 11 (1-based) 1 means signal read and 0 means reset read.
    signalFlag = repmat([zeros(numPixels, 1, 'uint32');  ...
        bitshift(ones(numPixels, 1, 'uint32'), signalFlagShiftBits)], ...
        numFrames, 1);
    % The last event mask is synonymous with the sample from x=0 y=0; data is
    % therefore ordered backwards.
    for frameIndex = 1 : numFrames
        samplesTemp = frameData.samples{frameIndex}(:);
        samplesTemp = samplesTemp(end : -1 : 1);
        samples((frameIndex * 2 - 1) * numPixels + 1 : frameIndex * 2 * numPixels) ...
            = samplesTemp ;
        timeStamps((frameIndex - 1) * 2 * numPixels + 1 : frameIndex * 2 * numPixels) ...
            = frameData.timeStampStart(frameIndex);
    end
    frameFlag = bitshift(1, frameFlagShiftBits, 'uint32');
    y = bitshift(y, yShiftBits, 'uint32');
    x = bitshift(x, xShiftBits, 'uint32');
    % samples should now be in the range 0-1023 (10-bit).
    % subtract samples from 1023. This has the effect of leaving all the reset
    % frame samples at 1023 - the highest value, against which the signal frames
    % will later be subtracted.
    samples = 1023 - samples';

    allSamples = [allSamples; frameFlag + y + x + signalFlag + samples];
    allTimeStamps = [allTimeStamps; timeStamps];
end

%% IMU6

if isfield(aedat.data, 'imu6') ...
        && (~isfield(aedat.exportParams, 'dataTypes') ...
        || any(cellfun(cellFind('imu6'), aedat.exportParams.dataTypes)))
    disp('Preparing imu6 data ...')

    imuFlag = uint32(0);
    imuFlag = bitset(imuFlag, 31+1);
    imuFlag = bitset(imuFlag, 11+1);
    imuFlag = bitset(imuFlag, 10+1);
    accelXFlag      = bitshift(uint32(0), 28);
    accelYFlag      = bitshift(uint32(1), 28);
    accelZFlag      = bitshift(uint32(2), 28);
    temperatureFlag = bitshift(uint32(3), 28);
    gyroXFlag       = bitshift(uint32(4), 28);
    gyroYFlag       = bitshift(uint32(5), 28);
    gyroZFlag       = bitshift(uint32(6), 28);

    accelX = convertAccelUint32(aedat.data.imu6.accelX);
    accelX = accelX + imuFlag + accelXFlag;

    accelY = convertAccelUint32(aedat.data.imu6.accelY);
    accelY = accelY + imuFlag + accelYFlag;

    accelZ = convertAccelUint32(aedat.data.imu6.accelZ);
    accelZ = accelZ + imuFlag + accelZFlag;

    temp = aedat.data.imu6.temperature; % conversion from g to full scale, and shift bits
    temp = int16((temp - 35) * 340); % conversion from K to full scale 16 range
    temp = [temp zeros(aedat.data.imu6.numEvents, 1, 'int16')];
    temp = temp';
    temp = temp(:);
    temp = typecast(temp, 'uint32');
    temp = bitshift(temp, 12); % shift bits
    temp = temp + imuFlag + temperatureFlag;

    gyroX = convertGyroUint32(aedat.data.imu6.gyroX);
    gyroX = gyroX + imuFlag + gyroXFlag;

    gyroY = convertGyroUint32(aedat.data.imu6.gyroY);
    gyroY = gyroY + imuFlag + gyroYFlag;

    gyroZ = convertGyroUint32(aedat.data.imu6.gyroZ);
    gyroZ = gyroZ + imuFlag + gyroZFlag;

    allData = [accelX accelY accelZ temp gyroX gyroY gyroZ];
    allData = allData';
    allData = allData(:);

    timeStamps = uint32(aedat.data.imu6.timeStamp(:));
    timeStamps = repmat(timeStamps', 7 , 1);
    timeStamps = timeStamps(:);

    allSamples = [allSamples; allData];
    allTimeStamps = [allTimeStamps; timeStamps];
end

%% Sort the events by timestamp

disp('Sorting events ...')

[allTimeStamps, sortIndex] = sort(allTimeStamps);
allSamples = allSamples(sortIndex);

output = zeros(1, length(allSamples) * 2, 'uint32');

output(1:2:end) = allSamples;
output(2:2:end) = allTimeStamps; % set even elements to timestamps

%% Write to file

disp('Writing to file ...')

% Create the file
f = fopen(aedat.exportParams.filePath, 'w', 'b');

if ~isfield(aedat.exportParams, 'noHeader') || aedat.exportParams.noHeader == false

    % CRLF \r\n is needed to not break header parsing in jAER
    fprintf(f,'#!AER-DAT2.0\r\n');
    fprintf(f,'# This is a raw AE data file created by an export function in the AedatTools library\r\n');
    fprintf(f,'# Data format is int32 address, int32 timestamp (8 bytes total), repeated for each event\r\n');
    fprintf(f,'# Timestamps tick is 1 us\r\n');
    fprintf(f,['# AEChip: ' source '\r\n']);
    fprintf(f,'# End of ASCII Header\r\n');
end

% write addresses and timestamps
count = fwrite(f, output, 'uint32', 0, 'b') / 2; % write 4 byte data
fclose(f);
fprintf('wrote %d events to %s\n', count, aedat.exportParams.filePath);

function accel = convertAccelUint32(accel)
% Conversion from g to full scale, and shift bits
accel = int16(accel * 8192); % conversion from g to full scale 16 range
accel = [accel zeros(numel(accel), 1, 'int16')];
accel = accel';
accel = accel(:);
accel = typecast(accel, 'uint32');
accel = bitshift(accel, 12); % shift bits

function gyro = convertGyroUint32(gyro)
gyro = int16(gyro * 65.5);
gyro = [gyro zeros(numel(gyro), 1, 'int16')];
gyro = gyro';
gyro = gyro(:);
gyro = typecast(gyro, 'uint32');
gyro = bitshift(gyro, 12); % shift bits

author: panpanfei


text:
function ncopies = checkworkspace(obj);

stack = dbstack;
stack(1:2) = [];
stack      = rmfield(stack, 'line');
ncopies = 0;
if ~isempty(stack)
    % empty stack means the base workspace
    % check if we are in a different workspace
    if ~isequal(stack, obj.workspace)
        % is the current stack a subset of the obj workspace?
        % if yes, it means that the object was created in a subfunction
        % and that it is OK to modify it (because it does not exist in the
        % subfunction any more)
        subFlag = false;
        for index = 1:length(obj.workspace)
            if isequal(obj.workspace(index), stack(1))
                subFlag = true;
                if length(stack) > 1
                    for index2 = 1:length(obj.workspace)-index
                        if length(stack) < index2+1
                            subFlag = false;
                        elseif ~isequal(obj.workspace(index+index2), stack(index2+1))
                            subFlag = false;
                        end;
                    end;
                end;
                if subFlag, return; end;
            end;
        end;

        % if subfunction, must be a copie
        if ~isempty(obj.workspace) && strcmpi(stack(end).file, obj.workspace(end).file) && ...
                ~strcmpi(stack(end).name, obj.workspace(end).name)
            % We are within a subfunction. The MMO must have
            % been passed as an argument (otherwise the current
            % workspace and the workspace variable would be
            % equal).
            ncopies = 2;
        else
            if ~isscript(stack(1).file)
                ncopies = 2;
                % we are within a function. The MMO must have
                % been passed as an argument (otherwise the current
                % workspace and the workspace variable would be
                % equal).
            else
                % we cannot be in a function with 0 argument
                % (otherwise the current workspace and the workspace
                % variable would be equal). We must assume that
                % we are in a script.
                while ~isempty(stack) && ~isequal(stack, obj.workspace)
                    stack(1) = [];
                end;
                if ~isequal(stack, obj.workspace)
                    ncopies = 2;
                end;
            end;
        end;
    end;
end;

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function O = cluster_cutaways(clusters,textprefix,mycolor,whichcuts,varargin)

% groups clusters to image on brains by the first letter in whichcuts
%
% :Usage:
% ::
%
%    O = cluster_cutaways(clusters,textprefix,mycolor,whichcuts,[coords for center of cuts],[revx])
%
% ..
% groups clusters to image on brains by the first letter in whichcuts
% therefore, if you enter 'y' for whichcuts, the program will select
% clusters with similar y values (w/i 15 mm) and image them on the same
% brain. you need files called brain_render_T1.img/hdr on the path
% these should be scalped anatomicals for the brain image.
%
% you would also need 'brain_render_T1.mat' for the
% transparent brain surface, if you changed this script to get
% the transparent brain surface.
%
% if it cuts from the wrong side, try 'revx' as input into tor_3d.m
% or enter anything as 2nd var arg.
%
%
%
% :Inputs:
%
%   **clusters:**
%        cluster object to display on brain cuts
%
%   **textprefix:**
%        prefix for saving images to files
%
%     **mycolor:**
%         cluster colors - 3 el. vector, single letter, or string of letters
%        enter letter or letter string in single quotes. e.g., O.clcol = 'yrgb';
%
%     **whichcuts:**
%         letter string (no quotes) - which axes to cut along - xyz are choices
%
% :Optional Inputs:
%
%     **bestc**
%     best cut, 3-element vector of coordinates, default [0 0 0]
%
%     **revx**
%     text string to reverse x cut direction, enter 1 to do it.
%
%
% :Outputs:
%
%   **O:**
%        output returned from renderCluster_ui.m
%
%
% :Examples:
% ::
%
%    O = cluster_cutaways(clusters,'myoutname','y','yzx',[0 0 0],[revx])
%
% Uses renderCluster_ui.m
%
% ..
%    by Tor Wager, Jan 2003
% ..

mm = cat(1,clusters.mm_center);
d = tril(dist(mm(:,2)') < 15);  %distances
bestc = [0 0 0];
revx = 0;

done = zeros(1,length(clusters));

if length(varargin) > 0, bestc = varargin{1}; else bestc = [0 0 0]; end
if length(varargin) > 1, revx = varargin{2}; else revx = 0; end

clind = 1;
while any(~done)

    whclust = (d(:,clind));
    whclust(find(done)) = 0;
    whclust = find(whclust);
    done(whclust) = 1;

        if ~isempty(whclust)
            figure('Color','w');
            O = struct('dohead','y','dobrain','n','get_from','workspace', ...
                'clusters',clusters, ...
            'which_cl',whclust,'whichc',whichcuts,'bestCoords',bestc,'clcol',mycolor,'addtext','n', ...
            'head','brain_render_T1','revx',revx);

            O = renderCluster_ui(O); set(O.headp(1:2:end),'FaceColor',[.5 .5 .5]);material dull
            h = findobj('Type','Light');delete(h)

            if whichcuts(1) == 'y', view(171,6),[az,el] = view; hh = lightangle(az,el); %hh = lightangle(az,el);
            elseif whichcuts(1) == 'z', view(0,70);[az,el] = view; hh = lightangle(az,el);
            elseif whichcuts(1) == 'x', view(90,10);[az,el] = view; hh = lightangle(az,el);
            end

            %h = findobj('Type','Light');delete(h),h(1) = camlight(90,0); h(2) = camlight(270,0); h(3) = camlight(270,0);

            if clusters(clind).numVox < 20 & clusters(clind).numVox > 1,
                %hold on, plot3(clusters(1).XYZmm(1,:),clusters(1).XYZmm(2,:),clusters(1).XYZmm(3,:),'o','Color',mycolor,'MarkerFaceColor',mycolor,'MarkerSize',10)
            end

            saveas(gcf,[textprefix '_' whichcuts '_group' num2str(clind)],'tif');
            saveas(gcf,[textprefix '_' whichcuts '_group' num2str(clind)],'fig');
            close
        end

        clind = clind+1;
 end

 return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
classdef(Abstract) AbstractSecondOrderIntegrator < AbstractIntegrator
    %AbstractSecondOrderIntegrator Summary of this class goes here
    %   Detailed explanation goes here

    properties

    end

    methods
        [t,y,yp,te,ye,ype,ie] = integrate(obj, odefun, tspan, y0, yp0, evtsFunc, odeOutputFun)
    end
end


text:
function [intracellRxnList] = findIntracellularRxns(model)
% This function finds all reactions that do not contain metabolites in [e] comparment.
%
% USAGE:
%
%    [intracellRxnList] = findIntracellularRxns(model)
%
% INPUT:
%    model:               Model structure
%
% OUTPUT:
%    intracellRxnList:    List of intracellular reactions in model
%
% .. Author: - Ines Thiele, July 2013, http://thielelab.eu.

intracellRxnList = zeros(length(model.rxns),1);

% exclude extracellular reactions
RL = printRxnFormula(model,model.rxns,false);
for i = 1 : length(RL)
    tmp = regexp(RL{i},'\[e\]');
    if isempty(tmp) % intracellular reaction
        intracellRxnList(i,1)=1;
    end
end
intracellRxnList = model.rxns(find(intracellRxnList));


text:
function K = covADD(cov, hyp, x, z, i)

% Additive covariance function using a 1d base covariance function
% cov(x^p,x^q;hyp) with individual hyperparameters hyp.
%
% k(x^p,x^q) = \sum_{r \in R} sf_r \sum_{|I|=r}
%                 \prod_{i \in I} cov(x^p_i,x^q_i;hyp_i)
%
% hyp = [ hyp_1
%         hyp_2
%          ...
%         hyp_D
%         log(sf_R(1))
%          ...
%         log(sf_R(end)) ]
%
% where hyp_d are the parameters of the 1d covariance function which are shared
% over the different values of R(1) to R(end).
%
% Please see the paper Additive Gaussian Processes by Duvenaud, Nickisch and
% Rasmussen, NIPS, 2011 for details.
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch, 2010-09-10.
%
% See also COVFUNCTIONS.M.

R = cov{1};
nh = eval(feval(cov{2}));           % number of hypers per individual covariance
nr = numel(R);                      % number of different degrees of interaction
if nargin<3                                  % report number of hyper parameters
  K = ['D*', int2str(nh), '+', int2str(nr)];
  return
end
if nargin<4, z = []; end                                   % make sure, z exists
xeqz = isempty(z); dg = strcmp(z,'diag');                       % determine mode

[n,D] = size(x);                                                % dimensionality
sf2 = exp( 2*hyp(D*nh+(1:nr)) );        % signal variances of individual degrees

Kd = Kdim(cov{2},hyp,x,z);                % evaluate dimensionwise covariances K
if nargin<5                                                        % covariances
  EE = elsympol(Kd,max(R));               % Rth elementary symmetric polynomials
  K = 0; for ii=1:nr, K = K + sf2(ii)*EE(:,:,R(ii)+1); end    % sf2 weighted sum
else                                                               % derivatives
  if i <= D*nh                       % individual covariance function parameters
    j = fix(1+(i-1)/nh);              % j is the dimension of the hyperparameter
    if dg, zj='diag'; else if xeqz, zj=[]; else zj=z(:,j); end, end
    dKj = feval(cov{2},hyp(nh*(j-1)+(1:nh)),x(:,j),zj,i-(j-1)*nh);  % other dK=0
    % the final derivative is a sum of multilinear terms, so if only one term
    % depends on the hyperparameter under consideration, we can factorise it
    % out and compute the sum with one degree less
    E = elsympol(Kd(:,:,[1:j-1,j+1:D]),max(R)-1);  %  R-1th elementary sym polyn
    K = 0; for ii=1:nr, K = K + sf2(ii)*E(:,:,R(ii)); end     % sf2 weighted sum
    K = dKj.*K;
  elseif i <= D*nh+nr
    EE = elsympol(Kd,max(R));             % Rth elementary symmetric polynomials
    j = i-D*nh;
    K = 2*sf2(j)*EE(:,:,R(j)+1);                  % rest of the sf2 weighted sum
  else
    error('Unknown hyperparameter')
  end
end

% evaluate dimensionwise covariances K
function K = Kdim(cov,hyp,x,z)
  [n,D] = size(x);                                              % dimensionality
  nh = eval(feval(cov));            % number of hypers per individual covariance
  if nargin<4, z = []; end                                 % make sure, z exists
  xeqz = numel(z)==0; dg = strcmp(z,'diag') && numel(z)>0;      % determine mode

  if dg                                                        % allocate memory
    K = zeros(n,1,D);
  else
    if xeqz, K = zeros(n,n,D); else K = zeros(n,size(z,1),D); end
  end

  for d=1:D
    hyp_d = hyp(nh*(d-1)+(1:nh));                 % hyperparamter of dimension d
    if dg
      K(:,:,d) = feval(cov,hyp_d,x(:,d),'diag');
    else
      if xeqz
        K(:,:,d) = feval(cov,hyp_d,x(:,d));
      else
        K(:,:,d) = feval(cov,hyp_d,x(:,d),z(:,d));
      end
    end
  end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function val=PearsonsGammaInc(u,p)
%%PEARSONSGAMMAINC Evaluate Pearson's incomplete gamma function. This is
%          equal to (1/gamma(p+1))*integral_0^(u*sqrt(p+1)) exp(-t)*t^p dt.
%
%INPUTS: u The positive u value (or matrix for multiple values to be
%          computed in parallel) in the upper bound of the integral. u>=0.
%          It is allowable to pass a scalar for this and a matrix for p.
%        p The exponent of the t term in the integral, or a matrix for
%          multiple values to be done in parallel. p>-1. It is allowable to
%          pass a scalar for this and a matrix for u.
%
%OUTPUTS: val The value or values of Person's gamma function at the desired
%             points.
%
%Pearson's incomplete gamma function is defined in terms of an integral in
%Section 6.5.6 of [1]. In Equation 1.7 of [2], Pearson's incomplete gamma
%function is expressed in terms of Tricomi's incomplete gamma function and
%Equation 1.3 expressed tricomi's gamma function in terms of the type of
%agmma function used in gammainc. Thus, this function uses the
%transformations and the gammainc function to express the result.
%
%REFERENCE
%[1] Abramowitz, M. and Stegun, I. A. (Eds.). "Gamma function and related
%    functions." in Ch. 6 in Handbook of Mathematical Functions with
%    Formulas, Graphs, and Mathematical Tables, 9th printing. New York:
%    Dover, 1972.
%[2] W. Gautschi, "A computational procedure for incomplete gamma
%    functions," ACM Transactions on Mathematical Software, vol. 5, no. 4,
%    pp. 466-481, Dec. 1979.
%
%March 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(isscalar(u)&&~isscalar(p))
    u=u*ones(size(p));
elseif(~isscalar(u)&&isscalar(p))
    p=p*ones(size(u));
end

val=gammainc(sqrt(1+p).*u,1+p);

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:
function trinomial_test ( )

%*****************************************************************************80
%
%% TRINOMIAL_TEST tests TRINOMIAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 April 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRINOMIAL_TEST\n' );
  fprintf ( 1, '  TRINOMIAL evaluates the trinomial coefficient:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  T(I,J,K) = (I+J+K)! / I! / J! / K!\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     I     J     K    T(I,J,K)\n' );
  fprintf ( 1, '\n' );

  for k = 0 : 4
    for j = 0 : 4
      for i = 0 : 4
        t = trinomial ( i, j, k );
        fprintf ( 1, '  %4d  %4d  %4d  %8d\n', i, j, k, t );
      end
    end
  end

  return
end


text:
function [aTS, bTS] = cRecon1(scanParams, fName, slice, compFlag, rotFlag, shifts)

% [aTS, bTS] = cRecon1(scanParams, fName, slice, compFlag, shifts);
%
% Load and crop complex recon data.
% Returns pair of single-slice time-series arrays of size:
% ny x nx x nFrames -- aTS is real and bTS is imaginary;
% the latter is returned only if the optional compFlag
% input is set. Time-duration and spatial cropping are performed,
% but no trend removal. Inputs are pFile name [fName], and
% slice index [slice]. Optional input rotFlag performs a 90
% degree rotation to make functional images match inplanes;
% use MakeMovie.m or similar to check if necessary. Optional
% input shifts allows correction of fallback errors or similar.
% Assumes pwd points to valid session
% directory, and uses global mrSESSION.
%
% DBR 12/00

global mrSESSION

if ~exist('compFlag'), compFlag = 0; end
if ~exist('rotFlag'), rotFlag = 0; end

% Do recon, extract complex time series:
[aTS, bTS] = PaulyRecon(scanParams, fName, slice, compFlag, rotFlag);

% Remove junk frames:
f0 = mrSESSION.junkFirstFrames+1;
nFrames = mrSESSION.nFrames;
fEnd = f0 + nFrames - 1;

% Crop:
x0 = mrSESSION.tseriesCrop(1, 1) + shifts(2);
xN = mrSESSION.tseriesCrop(2, 1) + shifts(2);
y0 = mrSESSION.tseriesCrop(1, 2) + shifts(1);
yN = mrSESSION.tseriesCrop(2, 2) + shifts(1);

% Crop in time/space and shuffle to standard t-series shape:
aTS = aTS(f0:fEnd, :, :);
aTS = aTS(:, y0:yN, x0:xN);
aTS = reshape(aTS, nFrames, (yN-y0+1)*(xN-x0+1));
if compFlag
  bTS = bTS(f0:fEnd, :, :);
  bTS = bTS(:, y0:yN, x0:xN);
  bTS = reshape(bTS, nFrames, (yN-y0+1)*(xN-x0+1));
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function calpak_test686 ( )

%*****************************************************************************80
%
%% CALPAK_TEST686 tests YJF_TO_YMDF_JULIAN and YMDF_TO_YJF_JULIAN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST686\n' );
  fprintf ( 1, '  For the Julian calendar,\n' );
  fprintf ( 1, '  YJF_TO_YMDF_JULIAN: YJF => YMDF\n' );
  fprintf ( 1, '  YMDF_TO_YJF_JULIAN: YMDF => YJF\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  YMDF(in)         YJF        YMDF(out)\n' );
  fprintf ( 1, '\n' );

  jed_epoch = epoch_to_jed_julian ( );

  i = 0;

  while ( 1 )

    i = i + 1;
    jed1 = jed_test ( i );

    if ( jed1 < 0.0 )
      break
    end

    if ( jed_epoch <= jed1 )

      [ y1, m1, d1, f1 ] = jed_to_ymdf_julian ( jed1 );

      s1 = ymdf_to_s_julian ( y1, m1, d1, f1 );

      [ y2, j2, f2 ] = ymdf_to_yjf_julian ( y1, m1, d1, f1 );

      s2 = yjf_to_s_julian ( y2, j2, f2 );

      [ y3, m3, d3, f3 ] = yjf_to_ymdf_julian ( y2, j2, f2 );

      s3 = ymdf_to_s_julian ( y3, m3, d3, f3 );

      fprintf ( 1, '  %10s  %10s  %10s\n', s1, s2, s3 );

    end

  end

  return
end