text:
% Defines a transcription unit
%
% Author: Jonathan Karr
% Affilitation: Covert Lab, Department of Bioengineering, Stanford University
% Last updated: 5/7/2009
classdef TranscriptionUnit < edu.stanford.covert.cell.kb.ssRNA
    properties
        genome           = edu.stanford.covert.cell.kb.Genome.empty(0,0);

        genes            = edu.stanford.covert.cell.kb.Gene.empty(0,0);
        geneCompartments = edu.stanford.covert.cell.kb.Compartment.empty(0,0);

        compartment      = edu.stanford.covert.cell.kb.Compartment.empty(0,0);

        transcriptionFactorProteinMonomers            = edu.stanford.covert.cell.kb.ProteinMonomer.empty(0,0);
        transcriptionFactorProteinMonomerAffinitys    = []; %nM
        transcriptionFactorProteinMonomerActivitys    = [];
        transcriptionFactorProteinMonomerConditions   = [];
        transcriptionFactorProteinMonomerBindingSiteStartCoordinates = [];
        transcriptionFactorProteinMonomerBindingSiteLengths          = [];
        transcriptionFactorProteinMonomerBindingSiteDirections       = [];
        transcriptionFactorProteinMonomerCompartments = edu.stanford.covert.cell.kb.Compartment.empty(0,0);

        transcriptionFactorProteinComplexs            = edu.stanford.covert.cell.kb.ProteinComplex.empty(0,0);
        transcriptionFactorProteinComplexAffinitys    = []; %nM
        transcriptionFactorProteinComplexActivitys    = [];
        transcriptionFactorProteinComplexConditions   = [];
        transcriptionFactorProteinComplexBindingSiteStartCoordinates = [];
        transcriptionFactorProteinComplexBindingSiteLengths          = [];
        transcriptionFactorProteinComplexBindingSiteDirections       = [];
        transcriptionFactorProteinComplexCompartments = edu.stanford.covert.cell.kb.Compartment.empty(0,0);
    end

    properties %(SetAccess = protected)
        halfLife
        sequence
        promoter35Coordinate
        promoter35Length
        promoter10Coordinate
        promoter10Length
        tssCoordinate
    end

    %computed properties
    properties %(SetAccess = protected)
        type
        dnaSequence
        startCoordinate
        endCoordinate
        fivePrimeCoordinate
        threePrimeCoordinate
        direction
        synthesisRate
        expression

        promoter1Sequence
        promoter10Sequence
        promoter35Sequence

        intergenicSequences
        intergenicSequenceLengths
        intergenicSequenceBaseCounts
        intergenicSequenceDecayReactions
        intergenicSequenceMolecularWeights
        intergenicSequenceHalfLives

        numGenes
        features
    end

    methods
        function this = TranscriptionUnit(knowledgeBase, wid,wholeCellModelID,name,...
                promoter35Coordinate, promoter35Length, ...
                promoter10Coordinate, promoter10Length, ...
                tssCoordinate, ...
                comments,crossReferences)

            if nargin == 0; return; end;

            this = edu.stanford.covert.cell.kb.TranscriptionUnit.empty(size(wid,1),0);
            this(size(wid,1),1) = edu.stanford.covert.cell.kb.TranscriptionUnit;
            for i = 1:size(wid,1)
                this(i,1).idx = i;
                this(i,1).knowledgeBase = knowledgeBase;
                this(i,1).wid = wid(i);
                this(i,1).wholeCellModelID = wholeCellModelID{i};
                this(i,1).name = name{i};
                if exist('comments','var') && ~isempty(comments); this(i,1).comments = comments{i}; end;
                if exist('crossReferences','var')
                    if size(crossReferences,1)>1
                        this(i,1).crossReferences = crossReferences(i);
                    else
                        this(i,1).crossReferences = struct;
                        fields = fieldnames(crossReferences);
                        for j = 1:size(fields,1)
                            values = crossReferences.(fields{j});
                            this(i,1).crossReferences.(fields{j}) = values(i);
                        end
                    end
                end

                this(i,1).promoter35Coordinate = promoter35Coordinate(i);
                this(i,1).promoter35Length     = promoter35Length(i);
                this(i,1).promoter10Coordinate = promoter10Coordinate(i);
                this(i,1).promoter10Length     = promoter10Length(i);
                this(i,1).tssCoordinate        = tssCoordinate(i);
            end
        end

        function serializeLinks(this)
            for i = 1:numel(this)
                this(i).genome                                        = this.serializeLinksHelper(this(i).genome);

                this(i).genes                                         = this.serializeLinksHelper(this(i).genes);
                this(i).geneCompartments                              = this.serializeLinksHelper(this(i).geneCompartments);

                this(i).compartment                                   = this.serializeLinksHelper(this(i).compartment);

                this(i).transcriptionFactorProteinMonomers            = this.serializeLinksHelper(this(i).transcriptionFactorProteinMonomers);
                this(i).transcriptionFactorProteinMonomerAffinitys    = this.serializeLinksHelper(this(i).transcriptionFactorProteinMonomerAffinitys);
                this(i).transcriptionFactorProteinMonomerActivitys    = this.serializeLinksHelper(this(i).transcriptionFactorProteinMonomerActivitys);
                this(i).transcriptionFactorProteinMonomerConditions   = this.serializeLinksHelper(this(i).transcriptionFactorProteinMonomerConditions);
                this(i).transcriptionFactorProteinMonomerCompartments = this.serializeLinksHelper(this(i).transcriptionFactorProteinMonomerCompartments);

                this(i).transcriptionFactorProteinComplexs            = this.serializeLinksHelper(this(i).transcriptionFactorProteinComplexs);
                this(i).transcriptionFactorProteinComplexAffinitys    = this.serializeLinksHelper(this(i).transcriptionFactorProteinComplexAffinitys);
                this(i).transcriptionFactorProteinComplexActivitys    = this.serializeLinksHelper(this(i).transcriptionFactorProteinComplexActivitys);
                this(i).transcriptionFactorProteinComplexConditions   = this.serializeLinksHelper(this(i).transcriptionFactorProteinComplexConditions);
                this(i).transcriptionFactorProteinComplexCompartments = this.serializeLinksHelper(this(i).transcriptionFactorProteinComplexCompartments);

                serializeLinks@edu.stanford.covert.cell.kb.ssRNA(this(i));
            end
        end

        function deserializeLinks(this, kb)
            for i = 1:numel(this)
                this(i).genome                                        = this.deserializeLinksHelper(this(i).genome, kb.genome);

                this(i).genes                                         = this.deserializeLinksHelper(this(i).genes, kb.genes);
                this(i).geneCompartments                              = this.deserializeLinksHelper(this(i).geneCompartments, kb.compartments);

                this(i).compartment                                   = this.deserializeLinksHelper(this(i).compartment, kb.compartments);

                this(i).transcriptionFactorProteinMonomers            = this.deserializeLinksHelper(this(i).transcriptionFactorProteinMonomers, kb.proteinMonomers);
                this(i).transcriptionFactorProteinMonomerAffinitys    = this.deserializeLinksHelper(this(i).transcriptionFactorProteinMonomerAffinitys);
                this(i).transcriptionFactorProteinMonomerActivitys    = this.deserializeLinksHelper(this(i).transcriptionFactorProteinMonomerActivitys);
                this(i).transcriptionFactorProteinMonomerConditions   = this.deserializeLinksHelper(this(i).transcriptionFactorProteinMonomerConditions);
                this(i).transcriptionFactorProteinMonomerCompartments = this.deserializeLinksHelper(this(i).transcriptionFactorProteinMonomerCompartments, kb.compartments);

                this(i).transcriptionFactorProteinComplexs            = this.deserializeLinksHelper(this(i).transcriptionFactorProteinComplexs, kb.proteinComplexs);
                this(i).transcriptionFactorProteinComplexAffinitys    = this.deserializeLinksHelper(this(i).transcriptionFactorProteinComplexAffinitys);
                this(i).transcriptionFactorProteinComplexActivitys    = this.deserializeLinksHelper(this(i).transcriptionFactorProteinComplexActivitys);
                this(i).transcriptionFactorProteinComplexConditions   = this.deserializeLinksHelper(this(i).transcriptionFactorProteinComplexConditions);
                this(i).transcriptionFactorProteinComplexCompartments = this.deserializeLinksHelper(this(i).transcriptionFactorProteinComplexCompartments, kb.compartments);

                deserializeLinks@edu.stanford.covert.cell.kb.ssRNA(this(i), kb);
            end
        end

        function deleteLinks(this)
            for i = 1:numel(this)
                this(i).genome                                        = [];

                this(i).genes                                         = [];
                this(i).geneCompartments                              = [];

                this(i).compartment                                   = [];

                this(i).transcriptionFactorProteinMonomers            = [];
                this(i).transcriptionFactorProteinMonomerAffinitys    = [];
                this(i).transcriptionFactorProteinMonomerActivitys    = [];
                this(i).transcriptionFactorProteinMonomerConditions   = [];
                this(i).transcriptionFactorProteinMonomerCompartments = [];

                this(i).transcriptionFactorProteinComplexs            = [];
                this(i).transcriptionFactorProteinComplexAffinitys    = [];
                this(i).transcriptionFactorProteinComplexActivitys    = [];
                this(i).transcriptionFactorProteinComplexConditions   = [];
                this(i).transcriptionFactorProteinComplexCompartments = [];

                deleteLinks@edu.stanford.covert.cell.kb.ssRNA(this(i));
            end
        end

        function value = get.halfLife(this)
            %retrieve
            if ~isempty(this.halfLife)
                value = this.halfLife;
                return;
            end

            %compute
            value = mean([this.genes.halfLife]);

            %store
            this.halfLife = value;
        end

        function value = get.type(this)
            %retrieve
            if ~isempty(this.type)
                value = this.type;
                return;
            end

            %compute
            if ~isempty(this.genes)
                value = this.genes(1).type;
            else
                value = [];
            end

            %store
            this.type = value;
        end

        function value = get.dnaSequence(this)
            %retrieve
            if ~isempty(this.dnaSequence)
                value = this.dnaSequence;
                return;
            end

            %compute
            if this.direction
                value = this.genome.sequence(this.startCoordinate:this.endCoordinate);
            else
                value = seqrcomplement(this.genome.sequence(this.startCoordinate:this.endCoordinate));
            end

            %store
            this.dnaSequence = value;
        end

        function value = get.sequence(this)
            %retrieve
            if ~isempty(this.sequence)
                value = this.sequence;
                return;
            end

            %compute
            value = this.dnaSequence;
            value(this.dnaSequence == 'A') = 'U';
            value(this.dnaSequence == 'C') = 'G';
            value(this.dnaSequence == 'G') = 'C';
            value(this.dnaSequence == 'T') = 'A';

            %store
            this.sequence = value;
        end

        function value = get.startCoordinate(this)
            %retrieve
            if ~isempty(this.startCoordinate)
                value = this.startCoordinate;
                return;
            end

            %compute
            value = min([this.genes.startCoordinate]);

            %store
            this.startCoordinate = value;
        end

        function value = get.endCoordinate(this)
            %retrieve
            if ~isempty(this.endCoordinate)
                value = this.endCoordinate;
                return;
            end

            %compute
            value = max([this.genes.endCoordinate]);

            %store
            this.endCoordinate = value;
        end

        function value = get.fivePrimeCoordinate(this)
            %retrieve
            if ~isempty(this.fivePrimeCoordinate)
                value = this.fivePrimeCoordinate;
                return;
            end

            %compute
            if this.direction
                value = this.startCoordinate;
            else
                value = this.endCoordinate;
            end

            %store
            this.fivePrimeCoordinate = value;
        end

        function value = get.threePrimeCoordinate(this)
            %retrieve
            if ~isempty(this.threePrimeCoordinate)
                value = this.threePrimeCoordinate;
                return;
            end

            %compute
            if this.direction
                value = this.endCoordinate;
            else
                value = this.startCoordinate;
            end

            %store
            this.threePrimeCoordinate = value;
        end

        function value = get.direction(this)
            %retrieve
            if ~isempty(this.direction)
                value = this.direction;
                return;
            end

            %compute
            value = this.genes(1).direction;

            %store
            this.direction = value;
        end

        function value = get.synthesisRate(this)
            %retrieve
            if ~isempty(this.synthesisRate)
                value = this.synthesisRate;
                return;
            end

            %compute
            value = mean(reshape([this.genes.synthesisRate],[],this.numGenes),2)';

            %store
            this.synthesisRate = value;
        end

        function value = get.expression(this)
            %retrieve
            if ~isempty(this.expression)
                value = this.expression;
                return;
            end

            %compute
            value = this.synthesisRate*this.halfLife/log(2);

            %store
            this.expression = value;
        end

        function value = get.promoter1Sequence(this)
            %retrieve
            if ~isempty(this.promoter1Sequence)
                value = this.promoter1Sequence;
                return;
            end

            %compute
            if this.direction
                value = this.genome.sequence(this.fivePrimeCoordinate + this.tssCoordinate + (-1:-1) + 1);
            else
                value = seqcomplement(this.genome.sequence(this.fivePrimeCoordinate - this.tssCoordinate + (1:-1:1) - 1));
            end

            %store
            this.promoter1Sequence = value;
        end

        function value = get.promoter10Sequence(this)
            %retrieve
            if ~isempty(this.promoter10Sequence)
                value = this.promoter10Sequence;
                return;
            end

            %compute
            if this.direction
                value = this.genome.sequence(this.fivePrimeCoordinate + this.promoter10Coordinate + (-this.promoter10Length:-1) + 1);
            else
                value = seqcomplement(this.genome.sequence(this.fivePrimeCoordinate - this.promoter10Coordinate + (this.promoter10Length:-1:1) - 1));
            end

            %store
            this.promoter10Sequence = value;
        end

        function value = get.promoter35Sequence(this)
            %retrieve
            if ~isempty(this.promoter35Sequence)
                value = this.promoter35Sequence;
                return;
            end

            %compute
            if this.direction
                value = this.genome.sequence(this.fivePrimeCoordinate + this.promoter35Coordinate + (-this.promoter35Length:-1) + 1);
            else
                value = seqcomplement(this.genome.sequence(this.fivePrimeCoordinate - this.promoter35Coordinate + (this.promoter35Length:-1:1) - 1));
            end

            %store
            this.promoter35Sequence = value;
        end

        function value = get.intergenicSequences(this)
            %retrieve
            if ~isempty(this.intergenicSequences)
                value = this.intergenicSequences;
                return;
            end

            %compute
            value = cell(0,1);
            for i = 1:numel(this.genes)-1
                if strcmp(this.type,'mRNA')
                    continue;
                end
                if this.genes(i).endCoordinate+1 > this.genes(i+1).startCoordinate-1
                    continue;
                end

                if this.direction
                    seq = this.genome.sequence(this.genes(i).endCoordinate+1:this.genes(i+1).startCoordinate-1);
                else
                    seq = seqrcomplement(this.genome.sequence(this.genes(i).endCoordinate+1:this.genes(i+1).startCoordinate-1));
                end

                value{i,1} = seq;
                value{i,1}(seq == 'A') = 'U';
                value{i,1}(seq == 'C') = 'G';
                value{i,1}(seq == 'G') = 'C';
                value{i,1}(seq == 'T') = 'A';
            end

            %store
            this.intergenicSequences = value;
        end

        function value = get.intergenicSequenceLengths(this)
            %retrieve
            if ~isempty(this.intergenicSequenceLengths)
                value = this.intergenicSequenceLengths;
                return;
            end

            %compute
            value = zeros(0,1);
            for i = 1:numel(this.genes)-1
                if strcmp(this.type,'mRNA')
                    continue;
                end
                if this.genes(i).endCoordinate+1 > this.genes(i+1).startCoordinate-1
                    continue;
                end

                value(i,1) = (this.genes(i+1).startCoordinate-1) - (this.genes(i).endCoordinate+1) + 1;
            end

            %store
            this.intergenicSequenceLengths = value;
        end

        function value = get.intergenicSequenceBaseCounts(this)
            %retrieve
            if ~isempty(this.intergenicSequenceBaseCounts)
                value = this.intergenicSequenceBaseCounts;
                return;
            end

            %compute
            value = zeros(numel(this.intergenicSequences), this.knowledgeBase.numMetabolites);

            for i = 1:numel(this.intergenicSequences)
                value(i, this.knowledgeBase.nmpIndexs) = [...
                    sum('A' == this.intergenicSequences{i}) ...
                    sum('C' == this.intergenicSequences{i}) ...
                    sum('G' == this.intergenicSequences{i}) ...
                    sum('U' == this.intergenicSequences{i})];
            end

            %store
            this.intergenicSequenceBaseCounts = value;
        end

        function value = get.intergenicSequenceDecayReactions(this)
            %retrieve
            if ~isempty(this.intergenicSequenceDecayReactions)
                value = this.intergenicSequenceDecayReactions;
                return;
            end

            %compute
            value = this.intergenicSequenceBaseCounts;
            for i = 1:numel(this.intergenicSequences)
                value(i, this.knowledgeBase.waterIndexs)    = value(i, this.knowledgeBase.waterIndexs)    - max(0, length(this.intergenicSequences{i}) - 1);
                value(i, this.knowledgeBase.hydrogenIndexs) = value(i, this.knowledgeBase.hydrogenIndexs) + max(0, length(this.intergenicSequences{i}) - 1);
            end

            %store
            this.intergenicSequenceDecayReactions = value;
        end

        function value = get.intergenicSequenceMolecularWeights(this)
            %retrieve
            if ~isempty(this.intergenicSequenceMolecularWeights)
                value = this.intergenicSequenceMolecularWeights;
                return;
            end

            %compute
            value = zeros(size(this.intergenicSequences));

            for i = 1:numel(this.intergenicSequences)
                value(i) = this.calculateMolecularWeight(...
                    this.intergenicSequences{i}, length(this.intergenicSequences{i}), 'linear', ...
                    this.knowledgeBase.metaboliteMolecularWeights(this.knowledgeBase.nmpIndexs));
            end

            %store
            this.intergenicSequenceMolecularWeights = value;
        end

        function value = get.intergenicSequenceHalfLives(this)
            %retrieve
            if ~isempty(this.intergenicSequenceHalfLives)
                value = this.intergenicSequenceHalfLives;
                return;
            end

            %compute
            value = zeros(size(this.intergenicSequences));

            %store
            this.intergenicSequenceHalfLives = value;
        end

        function value = get.numGenes(this)
            %retrieve
            if ~isempty(this.numGenes)
                value = this.numGenes;
                return;
            end

            %compute
            value = length(this.genes);

            %store
            this.numGenes = value;
        end

        function value = get.features(this)
            %retrieve
            if ~isempty(this.features)
                value = this.features;
                return;
            end

            %compute
            featureStartCoordinates = [this.genome.features.startCoordinate];
            featureEndCoordinates = [this.genome.features.endCoordinate];
            value = this.genome.features(...
                (this.startCoordinate > featureStartCoordinates & this.startCoordinate < featureEndCoordinates) | ...
                (this.endCoordinate   > featureStartCoordinates & this.endCoordinate   < featureEndCoordinates) | ...
                (this.startCoordinate < featureStartCoordinates & this.endCoordinate   > featureEndCoordinates));

            %store
            this.features = value;
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function c = compute_children(bnet)
% COMPUTE_CHILDREN
% precomputes the children of nodes in a bnet
%
% The return value is a cell array for now

ss = size(bnet.dag, 1);
c = cell(ss, 1);
for i = 1:ss
  c{i} = children(bnet.dag, i);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function y = quadcc(fun,a,b,tol)
%QUADCC   Numerical integration using Clenshaw-Curtis quadrature.
%   Y = QUADCC(FUN,A,B) estimates the definite integral of the
%   function FUN from A to B, using an adaptive Clenshaw-Curtis
%   quadrature scheme.  FUN is either a MATLAB expression written
%   as a string or inline function, or a function M-file.  A and
%   B are the lower and upper limits of integration, respectively.
%
%   QUADCC computes the integral by recursively dividing the interval
%   (A,B) into finer subintervals, defined by the zero-points of the
%   Type 1 Chebyshev polynomials.  The recursion continues until
%   either (1) the difference between computed estimates of the
%   integral in successive recursion steps falls below a tolerance
%   (default TOL = 1e-6), or (2) the integral domain has been
%   divided into 1024 subintervals.  A warning is returned if the
%   maximum number of subintervals is reached before the tolerance
%   criterion is met.
%   Y = QUADCC(FUN,A,B,TOL) uses the supplied value of TOL instead
%   of the default.
%
%   QUADCC provides limited handling of certain types of improper
%   integral.  If the integrand has a pole at one (or both) of the
%   integration limits, QUADCC returns a warning and attempts to
%   evaluate the integral by shifting the integration limit a
%   distance EPS inside the integration interval.  If QUADCC
%   encounters a pole inside the integration interval, an error is
%   produced.
%
%   Examples:
%
%       1) String expression
%            Y = quadcc('1./(x.^2-5)',-1,2)
%
%       2) Inline function
%            foo = inline('x.*sin(x)')
%            Y = quadcc(foo,0,2*pi)
%
%       3) Function M-file
%            Y = quadcc(@foo,0,1)
%            with foo.m a function M-file:
%
%            function y = foo(x)
%            y = tan(x);
%
%   See also  QUAD, QUADL, DBLQUAD, TRIPLEQUAD, INLINE

%   References
%   ----------
%   "Numerical Recipes in C. The Art of Scientific Computing, 2nd edition".
%   W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery.
%   Cambridge Unversity Press, 2002.

% Paul Fricker  12/09/2002

if a==b
    y = 0;
    return
end

if ~exist('tol')
    tol = 1e-6;
end

% Make sure 'fun' is an inline function
f = fcnchk(fun);

% Initialize by breaking the integration interval into four regions.
N = 4;
p = (b+a)/2;
q = (b-a)/2;
% Write out [ cos((0:4)/4*pi) ] explicitly to avoid cos(pi/2) roundoff
x = p-q*[1 0.7071067811865475 0 -0.7071067811865475 -1];
F = feval(f,x);

% Simple attempt to handle improper integral
if isinf(F(1))
    F(1) = feval(f,x(1)+eps);
    warning('QUADCC:improperLower', ...
            ['Improper integral: Integrand has pole ' ...
             'at lower integration limit.'])
end

if isinf(F(5))
    F(5) = feval(f,x(5)-eps);
    warning('QUADCC:improperUpper', ...
            ['Improper integral: Integrand has pole ' ...
             'at upper integration limit.'])
end

% Modify the first and last terms of 'F'.
F(1) = F(1)/2;
F(end) = F(end)/2;

% Call the integrator function recursively
[y,warn] = intcc(f,x,F,N,p,q,tol,Inf);

if warn==1
    warning('QUADCC:MaxSubDivide', ...
            ['Computation terminated before tolerance criterion (TOL = ' num2str(tol) ') was met.'])
end

% EOF 'quadcc'

function [y,warn] = intcc(f,x,F,N,p,q,tol,yinit)
% INTCC Recursive integrator function for the QUADCC routine.
%
% Input parameters:
% f     : Integrand
% x     : Chebyshev polynomial zero-points between the integration limits (a,b)
% F     : Values of 'f' evaluated at 'x'
% N     : Initial number of zero-points ('discretization')
% p,q   : Endpoint parameters (from integration limits)
% tol   : Integration tolerance (terminates the recursion)
% yinit : Estimated value of the integral from the previous pass

% Terminate recursion if integral has not converged before 1024 points
if size(x,2) > 2^10
    warn = 1;
    y = yinit;
    return
else
    warn = 0;
end

% Double the size of the approximation
N = 2*N;
Nvec = 0:N;

x2 = zeros(1,N+1);
x2(1:2:end) = x;
x2(2:2:end) = p-q*cos((1:2:N)/N*pi);

G = zeros(1,N+1);
G(1:2:end) = F;
G(2:2:end) = feval(f,x2(2:2:end));
G = G(:);

if any(isinf(G))
    error('Improper integral: Integrand has a pole between integration limits.')
end

% Evaluate the Chebyshev polynomial coefficients for
% approximating the input function 'f'
M = cos(Nvec(1:2:end-1)'*Nvec/N*pi);
M(abs(M)<eps) = 0;
C = 4/N*M*G;

% Compute the quadrature coefficients
coeffs = -1./Nvec(2:2:end)./(Nvec(2:2:end)-2);
coeffs(1) = 0.5;

% Compute the integral estimate
y = q*coeffs*C;

% Continue recursion if the difference between the current value
% of the integral and the previous value is still above 'tol'.
if abs(y-yinit) > tol
    [y,warn] = intcc(f,x2,G,N,p,q,tol,y);
end

% EOF 'intcc'

author: Sable

repo: mcbench-benchmarks


text:
function [traj, infStates] = tapas_hgf_binary_pu(r, p, varargin)
% Calculates the trajectories of the agent's representations under the HGF
%
% This function can be called in two ways:
%
% (1) tapas_hgf_binary_pu(r, p)
%
%     where r is the structure generated by tapas_fitModel and p is the parameter vector in native space;
%
% (2) tapas_hgf_binary_pu(r, ptrans, 'trans')
%
%     where r is the structure generated by tapas_fitModel, ptrans is the parameter vector in
%     transformed space, and 'trans' is a flag indicating this.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2015 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Transform paramaters back to their native space if needed
if ~isempty(varargin) && strcmp(varargin{1},'trans');
    p = tapas_hgf_binary_pu_transp(r, p);
end

% Number of levels
try
    l = r.c_prc.n_levels;
catch
    l = (length(p)+1)/5;

    if l ~= floor(l)
        error('tapas:hgf:UndetNumLevels', 'Cannot determine number of levels');
    end
end

% Unpack parameters
mu_0 = p(1:l);
sa_0 = p(l+1:2*l);
rho  = p(2*l+1:3*l);
ka   = p(3*l+1:4*l-1);
om   = p(4*l:5*l-2);
th   = exp(p(5*l-1));
al   = p(5*l);
eta0 = p(5*l+1);
eta1 = p(5*l+2);

% Add dummy "zeroth" trial
u = [0; r.u(:,1)];

% Number of trials (including prior)
n = length(u);

% Assume that if u has more than one column, the last contains t
try
    if r.c_prc.irregular_intervals
        if size(u,2) > 1
            t = [0; r.u(:,end)];
        else
            error('tapas:hgf:InputSingleColumn', 'Input matrix must contain more than one column if irregular_intervals is set to true.');
        end
    else
        t = ones(n,1);
    end
catch
    if size(u,2) > 1
        t = [0; r.u(:,end)];
    else
        t = ones(n,1);
    end
end

% Initialize updated quantities

% Representations
mu = NaN(n,l);
pi = NaN(n,l);

% Other quantities
muhat = NaN(n,l);
pihat = NaN(n,l);
v     = NaN(n,l);
w     = NaN(n,l-1);
da    = NaN(n,l);

% Representation priors
% Note: first entries of the other quantities remain
% NaN because they are undefined and are thrown away
% at the end; their presence simply leads to consistent
% trial indices.
mu(1,1) = tapas_sgm(mu_0(1), 1);
pi(1,1) = Inf;
mu(1,2:end) = mu_0(2:end);
pi(1,2:end) = 1./sa_0(2:end);

% Pass through representation update loop
for k = 2:1:n
    if not(ismember(k-1, r.ign))

        %%%%%%%%%%%%%%%%%%%%%%
        % Effect of input u(k)
        %%%%%%%%%%%%%%%%%%%%%%

        % 2nd level prediction
        muhat(k,2) = mu(k-1,2) +t(k) *rho(2);

        % 1st level
        % ~~~~~~~~~
        % Prediction
        muhat(k,1) = tapas_sgm(ka(1) *muhat(k,2), 1);

        % Precision of prediction
        pihat(k,1) = 1/(muhat(k,1)*(1 -muhat(k,1)));

        % Mean update
        mu(k,1) = u(k);
        und1 = exp(-(u(k) -eta1)^2/(2*al));
        und0 = exp(-(u(k) -eta0)^2/(2*al));
        mu(k,1) = muhat(k,1) *und1 /(muhat(k,1) *und1 +(1 -muhat(k,1)) *und0);

        % Prediction error
        da(k,1) = mu(k,1) -muhat(k,1);

        % 2nd level
        % ~~~~~~~~~
        % Prediction: see above

        % Precision of prediction
        pihat(k,2) = 1/(1/pi(k-1,2) +exp(ka(2) *mu(k-1,3) +om(2)));

        % Updates
        pi(k,2) = pihat(k,2) +ka(1)^2/pihat(k,1);
        mu(k,2) = muhat(k,2) +ka(1)/pi(k,2) *da(k,1);

        % Implied posterior precision at first level
        sgmmu2 = tapas_sgm(ka(1) *mu(k,2), 1);
        pi(k,1) = pi(k,2)/(sgmmu2*(1-sgmmu2));

        % Volatility prediction error
        da(k,2) = (1/pi(k,2) +(mu(k,2) -muhat(k,2))^2) *pihat(k,2) -1;

        if l > 3
            % Pass through higher levels
            % ~~~~~~~~~~~~~~~~~~~~~~~~~~
            for j = 3:l-1
                % Prediction
                muhat(k,j) = mu(k-1,j) +t(k) *rho(j);

                % Precision of prediction
                pihat(k,j) = 1/(1/pi(k-1,j) +t(k) *exp(ka(j) *mu(k-1,j+1) +om(j)));

                % Weighting factor
                v(k,j-1) = t(k) *exp(ka(j-1) *mu(k-1,j) +om(j-1));
                w(k,j-1) = v(k,j-1) *pihat(k,j-1);

                % Updates
                pi(k,j) = pihat(k,j) +1/2 *ka(j-1)^2 *w(k,j-1) *(w(k,j-1) +(2 *w(k,j-1) -1) *da(k,j-1));

                if pi(k,j) <= 0
                    error('tapas:hgf:NegPostPrec', 'Negative posterior precision. Parameters are in a region where model assumptions are violated.');
                end

                mu(k,j) = muhat(k,j) +1/2 *1/pi(k,j) *ka(j-1) *w(k,j-1) *da(k,j-1);

                % Volatility prediction error
                da(k,j) = (1/pi(k,j) +(mu(k,j) -muhat(k,j))^2) *pihat(k,j) -1;
            end
        end

        % Last level
        % ~~~~~~~~~~
        % Prediction
        muhat(k,l) = mu(k-1,l) +t(k) *rho(l);

        % Precision of prediction
        pihat(k,l) = 1/(1/pi(k-1,l) +t(k) *th);

        % Weighting factor
        v(k,l)   = t(k) *th;
        v(k,l-1) = t(k) *exp(ka(l-1) *mu(k-1,l) +om(l-1));
        w(k,l-1) = v(k,l-1) *pihat(k,l-1);

        % Updates
        pi(k,l) = pihat(k,l) +1/2 *ka(l-1)^2 *w(k,l-1) *(w(k,l-1) +(2 *w(k,l-1) -1) *da(k,l-1));

        if pi(k,l) <= 0
            error('tapas:hgf:NegPostPrec', 'Negative posterior precision. Parameters are in a region where model assumptions are violated.');
        end

        mu(k,l) = muhat(k,l) +1/2 *1/pi(k,l) *ka(l-1) *w(k,l-1) *da(k,l-1);

        % Volatility prediction error
        da(k,l) = (1/pi(k,l) +(mu(k,l) -muhat(k,l))^2) *pihat(k,l) -1;
    else

        mu(k,:) = mu(k-1,:);
        pi(k,:) = pi(k-1,:);

        muhat(k,:) = muhat(k-1,:);
        pihat(k,:) = pihat(k-1,:);

        v(k,:)  = v(k-1,:);
        w(k,:)  = w(k-1,:);
        da(k,:) = da(k-1,:);

    end
end

% Implied learning rate at the first level
sgmmu2 = tapas_sgm(ka(1) *mu(:,2), 1);
lr1    = diff(sgmmu2)./da(2:n,1);
lr1(da(2:n,1)==0) = 0;

% Remove representation priors
mu(1,:)  = [];
pi(1,:)  = [];

% Check validity of trajectories
if any(isnan(mu(:))) || any(isnan(pi(:)))
    error('tapas:hgf:VarApproxInvalid', 'Variational approximation invalid. Parameters are in a region where model assumptions are violated.');
else
    % Check for implausible jumps in trajectories
    dmu = diff(mu(:,2:end));
    dpi = diff(pi(:,2:end));
    rmdmu = repmat(sqrt(mean(dmu.^2)),length(dmu),1);
    rmdpi = repmat(sqrt(mean(dpi.^2)),length(dpi),1);

    jumpTol = 16;
    if any(abs(dmu(:)) > jumpTol*rmdmu(:)) || any(abs(dpi(:)) > jumpTol*rmdpi(:))
        error('tapas:hgf:VarApproxInvalid', 'Variational approximation invalid. Parameters are in a region where model assumptions are violated.');
    end
end

% Remove other dummy initial values
muhat(1,:) = [];
pihat(1,:) = [];
v(1,:)     = [];
w(1,:)     = [];
da(1,:)    = [];

% Create result data structure
traj = struct;

traj.mu     = mu;
traj.sa     = 1./pi;

traj.muhat  = muhat;
traj.sahat  = 1./pihat;

traj.v      = v;
traj.w      = w;
traj.da     = da;

% Updates with respect to prediction
traj.ud = mu -muhat;

% Psi (precision weights on prediction errors)
psi        = NaN(n-1,l);
psi(:,2)   = 1./pi(:,2);
psi(:,3:l) = pihat(:,2:l-1)./pi(:,3:l);
traj.psi   = psi;

% Epsilons (precision-weighted prediction errors)
epsi        = NaN(n-1,l);
epsi(:,2:l) = psi(:,2:l) .*da(:,1:l-1);
traj.epsi   = epsi;

% Full learning rate (full weights on prediction errors)
wt        = NaN(n-1,l);
wt(:,1)   = lr1;
wt(:,2)   = psi(:,2);
wt(:,3:l) = 1/2 *(v(:,2:l-1) *diag(ka(2:l-1))) .*psi(:,3:l);
traj.wt   = wt;

% Create matrices for use by the observation model
infStates = NaN(n-1,l,4);
infStates(:,:,1) = traj.muhat;
infStates(:,:,2) = traj.sahat;
infStates(:,:,3) = traj.mu;
infStates(:,:,4) = traj.sa;

return;

author: translationalneuromodeling

repo: tapas


text:
function opt = mpoption_info_ipopt(selector)
%MPOPTION_INFO_IPOPT  Returns MATPOWER option info for IPOPT.
%
%   DEFAULT_OPTS = MPOPTION_INFO_IPOPT('D')
%   VALID_OPTS   = MPOPTION_INFO_IPOPT('V')
%   EXCEPTIONS   = MPOPTION_INFO_IPOPT('E')
%
%   Returns a structure for IPOPT options for MATPOWER containing ...
%   (1) default options,
%   (2) valid options, or
%   (3) NESTED_STRUCT_COPY exceptions for setting options
%   ... depending on the value of the input argument.
%
%   This function is used by MPOPTION to set default options, check validity
%   of option names or modify option setting/copying behavior for this
%   subset of optional MATPOWER options.
%
%   See also MPOPTION.

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

if nargin < 1
    selector = 'D';
end
if have_feature('ipopt')
    switch upper(selector)
        case {'D', 'V'}     %% default and valid options
            opt = struct(...
                'ipopt',        struct(...
                    'opts',         [], ...
                    'opt_fname',    '', ...
                    'opt',          0 ...
                ) ...
            );
        case 'E'            %% exceptions used by nested_struct_copy() for applying
            opt = struct(...
                'name',         { 'ipopt.opts' }, ...
                'check',        0 ...
                );
%                 'copy_mode',    { @ipopt_options } ...
        otherwise
            error('mpoption_info_ipopt: ''%s'' is not a valid input argument', selector);
    end
else
    opt = struct([]);       %% IPOPT is not available
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function data = LoadPerFrameData(obj,fn,n)

islabelfile = regexp(fn,'[lL]abels','once');
isscoresfile = regexp(fn,'[sS]cores','once');
isspacetime = ~isempty(regexp(fn,'^st_(.*)','tokens','once'));
if islabelfile,
  labelfilestr = [fn,'.mat'];
  filename = fullfile(obj.expdirs{n},labelfilestr);
  if ~exist(filename,'file'),
    scoresfn = regexprep(fn,'labels','scores','preservecase','once');
    scoresfilestr = [scoresfn,'.mat'];
    scoresfilename = fullfile(obj.expdirs{n},scoresfilestr);
    if ~exist(scoresfilename,'file'),
      error('Files %s and %s do not exist',filename,scoresfilename);
    end
    [scores,data] = obj.LoadScoresFromFile(scoresfilestr,n);
    obj.units.(scoresfn) = parseunits('unit');
    obj.SetPerFrameData(scoresfn,scores,n);
  else
    data = obj.LoadLabelsFromFile(labelfilestr,n);
  end
  obj.units.(fn) = parseunits('unit');
  obj.SetPerFrameData(fn,data,n);
elseif isscoresfile,
  scoresfilestr = [fn,'.mat'];
  filename = fullfile(obj.expdirs{n},scoresfilestr);
  if ~exist(filename,'file'),
    error('File %s does not exist',filename);
  end
  [data,labels] = obj.LoadScoresFromFile(scoresfilestr,n);
  obj.units.(fn) = parseunits('unit');
  obj.SetPerFrameData(fn,data,n);
  labelfn = regexprep(fn,'scores','labels','preservecase','once');
  obj.SetPerFrameData(labelfn,labels,n);
else

  filename = obj.GetPerFrameFile(fn,n);
  if exist(filename,'file')
    x = load(filename,'data','units');
    if isspacetime
       if ~isequal(obj.stInfo,x.params),
           error('Space Time settings used to generate the perframe features are differnt than current parameters');
       end
    end
  else
    [x.data,x.units] = obj.ComputePerFrameData(fn,n);
  end
  obj.units.(fn) = x.units;
  obj.SetPerFrameData(fn,x.data,n);
  data = x.data;
end

author: kristinbranson

repo: JAABA


text:
classdef bdepth_ratio < ZmapHGridFunction
% BDEPTH_RATIO compare b values at two different depths
% Stefan Wiemer 1/95

    properties
        topzone_ceiling         = 0;    % Top of TOP zone
        topzone_floor           = 5;    % Bottom of TOP zone
        bottomzone_ceiling      = 7;    % Top of BOTTOM zone
        bottomzone_floor        = 15;   % Bottom of BOTTOM zone
        fBinning                = 0.1;  % magnitude bins
        mc_choice   McMethods   = McMethods.MaxCurvature % magnitude of completion method
        useAutoMcomp McAutoEstimate = true
    end

    properties(Constant)
        PlotTag='bdepth'
        ReturnDetails = cell2table({ ... VariableNames, VariableDescriptions, VariableUnits
            ...
            ... % these are returned from the calculation
            'bv_ratio',     'b-value Ratio Map',                            '';... %valueMap [old] #1
            'magco',        'Magnitude of completion (bottom)',             'mag';... #2
            'bv2',          'b-value Ratio map (with autoestimate Mcomp)',  '';... #3
            'av',           'a-value ratios',                   '';... #4 avm
            'Prmap',        'Utsu Probability map',             '';... #5 Prmap
            'top_b',        'Top Zone b-value map',             '';... #6 tob_b
            'bottom_b',     'Bottom Zone b-value map',          '';... #7 bottom_b
            'per_top',      '% of nodal EQs within TOP zone',   '';... #8 per_in_top
            'per_bot',      '% of nodal EQs within BOTTOM zone','';... #9 per_in_bot
            'Number_of_Events_top', 'Number of events in TOP zone',     '';... #10
            'Number_of_Events_bot', 'Number of events in BOTTOM zone',  '';... #11
            ...
            }, 'VariableNames', {'Names','Descriptions','Units'})

        CalcFields = {...
            'bv_ratio',     'magco',    'bv2',      'av',       'Prmap',...
            'top_b',        'bottom_b', 'per_top',  'per_bot',...
            'Number_of_Events_top',     'Number_of_Events_bot'};

        ParameterableProperties = ["topzone_ceiling" "topzone_floor"....
                    "bottomzone_ceiling" "bottomzone_floor"...
                    "fBinning" "NodeMinEventCount" "mc_choice" "useAutoMcomp"];

        References="";
    end

    methods
        function obj = bdepth_ratio(zap, varargin)
            % BDEPTH_RATIO compare b values at two different depths

            obj@ZmapHGridFunction(zap, 'bv_ratio');
            report_this_filefun();

            obj.parseParameters(varargin);
            obj.StartProcess();

        end

        function InteractiveSetup(obj)

            % get two time periods, along with grid and event parameters
            zdlg=ZmapDialog();
            zdlg.AddMcAutoEstimateCheckbox('useAutoMcomp',  obj.useAutoMcomp);
            zdlg.AddMcMethodDropdown('mc_choice',           obj.mc_choice);
            zdlg.AddHeader('Please define two Depth ranges to compare');
            zdlg.AddEdit('topzone_ceiling','TOP zone ceiling [km]',         obj.topzone_ceiling,'');
            zdlg.AddEdit('topzone_floor','TOP zone floor [km]',        obj.topzone_floor,'');
            zdlg.AddEdit('bottomzone_ceiling','BOTTOM zone ceiling [km]',   obj.bottomzone_ceiling,'');
            zdlg.AddEdit('bottomzone_floor','BOTTOM zone floor [km]',  obj.bottomzone_floor,'');
            obj.AddDialogOption(zdlg, 'NodeMinEventCount');
            obj.AddDialogOption(zdlg,'EventSelector');
            zdlg.Create('Name', 'Depth Ratio Grid input parameters','WriteToObj',obj,'OkFcn',@obj.doIt);
        end

        function modifyGlobals(obj)
            obj.ZG.bvg = obj.Result.values;
        end

        function results=Calculate(obj)

            assert(obj.topzone_ceiling < obj.topzone_floor, 'TOP ZONE: zone floor is above zone ceiling');
            assert(obj.bottomzone_ceiling < obj.bottomzone_floor, 'BOTTOM ZONE: zone floor is above zone ceiling');

            %  make grid, calculate start- endtime etc.  ...
            %

            % find row index of ratio midpoint
            l =  obj.topzone_ceiling <= obj.RawCatalog.Depth & obj.RawCatalog.Depth <  obj.topzone_floor;
            top_zone = obj.RawCatalog.subset(l);

            l = obj.bottomzone_ceiling <= obj.RawCatalog.Depth & obj.RawCatalog.Depth <  obj.bottomzone_floor;
            bot_zone = obj.RawCatalog.subset(l);

            % overall b-value
            [top_b_overall] = bvalca3(top_zone.Magnitude, obj.useAutoMcomp);
            [bottom_b_overall] = bvalca3(bot_zone.Magnitude, obj.useAutoMcomp);

            depth_ratio = top_b_overall/bottom_b_overall;
            disp(depth_ratio);

            [~, mcCalculator] = calc_Mc([], obj.mc_choice,obj.fBinning);

            % loop over all points
            obj.gridCalculations(@calculation_function);

            obj.Result.top.depthrange_km=[obj.topzone_ceiling obj.topzone_floor];
            obj.Result.bottom.depthrange_km=[obj.bottomzone_ceiling obj.bottomzone_floor];

            if nargout
                results=obj.Result.values;
            end

            %catsave3('bdepth_ratio');

            % to View the b-value map : view_bdepth

            function out=calculation_function(b)
                topb = b.subset( obj.topzone_ceiling <= b.Depth & b.Depth <  obj.topzone_floor );
                per_in_top = (topb.Count/b.Count)*100.0;

                botb = b.subset( obj.bottomzone_ceiling<= b.Depth & b.Depth <  obj.bottomzone_floor );
                per_in_bot = (botb.Count/b.Count)*100.0;

                if length(topb) < obj.NodeMinEventCount  || length(botb) < obj.NodeMinEventCount

                    [Mc_valueTop] = mcCalculator(topb);
                    [Mc_valueBot] = mcCalculator(botb);
                    [topbv, topbv2, ~, topav, n1]=calc_bval_both_ways(topb,Mc_valueTop, top_b_overall);
                    [botbv, botbv2, magco, botav, n2]=calc_bval_both_ways(botb,Mc_valueBot, bottom_b_overall);
                else
                    [topbv, topbv2, ~, topav,  n1]=deal(nan);
                    [botbv, botbv2, magco, botav, n2]=deal(nan);
                end
                bv = topbv/botbv;
                bv2 = topbv2/botbv2;
                av = topav/botav;

                n = n1+n2;
                obj.ZG.overall_b_value = topbv;
                da = -2*n*log(n) + 2*n1*log(n1+n2 * topbv/botbv) + 2*n2*log(n1 * botbv/topbv + n2) - 2;
                pr = (1  -  exp(-da/2-2))*100;

                ltopb = topb.Count;
                lbotb = botb.Count;

                out = [bv magco ...
                    bv2 av pr ...
                    topbv botbv ...
                    per_in_top per_in_bot ...
                    ltopb lbotb];

                function [bv, bv2, magco, av, n] = calc_bval_both_ways(mycat,magco,b_overall)
                    % where mycat is already the subset value
                    idx = (magco-0.05) <= mycat.Magnitude5;
                    n=sum(idx);
                    if obj.NodeMinEventCount <= sum(idx)
                        [bv, magco, ~, av] =  bvalca3(mycat.Magnitude(idx), McAutoEstimate.manual, b_overall); %not automatic estimate of Mcomp
                        bv2 =  bvalca3(mycat.Magnitude(idx), McAutoEstimate.auto); % automatic estimate of Mcomp
                    else
                        [bv, bv2, magco, av] = deal(nan);
                    end
                end
            end
        end % Calculate

    end %methods

    methods(Static)
        function h = AddMenuItem(parent, zapFcn, varargin)
            % create a menu item
            label = 'b-value depth ratio grid';
            h = uimenu(parent, 'Label', label,...
                'MenuSelectedFcn', @(~,~)XYfun.bdepth_ratio(zapFcn()),...
                varargin{:});
        end

        function obj=my_load()
            unimplemented_error();
            % Load exist b-grid
            [file1,path1] = uigetfile(['*.mat'],'b-value gridfile');
            if length(path1) > 1

                load([path1 file1])
                normlap2=nan(length(tmpgri(:,1)),1)

                normlap2(ll)= bvg(:,1);
                valueMap=reshape(normlap2,length(yvect),length(xvect));

                normlap2(ll)= bvg(:,5);
                r=reshape(normlap2,length(yvect),length(xvect));

                normlap2(ll)= bvg(:,6);
                meg=reshape(normlap2,length(yvect),length(xvect));

                normlap2(ll)= bvg(:,2);
                old1=reshape(normlap2,length(yvect),length(xvect));

                %  normlap2(ll)= bvg(:,7);
                %  pro=reshape(normlap2,length(yvect),length(xvect));

                normlap2(ll)= bvg(:,7);
                avm=reshape(normlap2,length(yvect),length(xvect));

                %  normlap2(ll)= bvg(:,9);
                % stanm=reshape(normlap2,length(yvect),length(xvect));

                normlap2(ll)= bvg(:,8);
                Prmap=reshape(normlap2,length(yvect),length(xvect));

                normlap2(ll)= bvg(:,9);
                top_b=reshape(normlap2,length(yvect),length(xvect));

                normlap2(ll)= bvg(:,10);
                bottom_b=reshape(normlap2,length(yvect),length(xvect));

                normlap2(ll)= bvg(:,11);
                per_top=reshape(normlap2,length(yvect),length(xvect));

                normlap2(ll)= bvg(:,12);
                per_bot=reshape(normlap2,length(yvect),length(xvect));

                %    normlap2(ll)= bvg(:,13);
                %    ltopb=reshape(normlap2,length(yvect),length(xvect));

                %   normlap2(ll)= bvg(:,14);
                %  lbotb=reshape(normlap2,length(yvect),length(xvect));

                old = valueMap;

                view_bdepth
            else
                return
            end
        end
    end
end

author: CelsoReyes


text:
function [ s ] = geosubset_construct(starttime, stoptime, t_stride, zmin, zmax, z_stride, eastmin, ...
  eastmax, x_stride, northmin, northmax, y_stride)
% SUBSET_STRUCT this is a function that, in a matlab way, creates a subset structure for nctoolbox's
% geosubset or timegeosubset functions.
% >> geosubset_construct(1, 1, [], 1, 1, [], -77, -75.5, [],  36, 38, 1)
%
% Note: the starttime and stoptime arugments also accept matlab time formats datevec and datenum.
%

% Bounding Dimensions
s.lat = [northmin northmax]; % both ends must be explicit
s.lon = [eastmin eastmax]; % both ends must be explicit
s.z_index = [zmin zmax]; % both ends must be explicit
s.time = {starttime stoptime}; %both ends must be explicit

% Stride Arugments
if isempty(z_stride) % Z
  s.z_stride = 1;
else
  s.z_stride = z_stride;
end

if isempty(t_stride) % Time
  s.t_stride = 1;
else
  s.t_stride = t_stride;
end

if isempty(x_stride) %lon
  s.xy_stride(1) = 1;
else
  s.xy_stride(1) = x_stride;
end

if isempty(y_stride) %lat
  s.xy_stride(2) = 1;
else
  s.xy_stride(2) = y_stride;
end

end

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function rational_quadratic_test ( )

%*****************************************************************************80
%
%% RATIONAL_QUADRATIC_TEST examines the rational_quadratic correlation.
%
%  Discussion:
%
%    This code is based substantially on a document by Toby Driscoll.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Toby Driscoll,
%    Mercer's theorem and the Karhunen-Loeve expansion,
%    Oxford University Mathematical Institute,
%    http://www2.maths.ox.ac.uk/chebfun/examples/stats/pdf/MercerKarhunenLoeve.pdf
%
  rmpath ( '/usr/local/matlab/toolbox/datafeed/datafeed' )
  addpath ( '../chebfun' )

  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'RATIONAL_QUADRATIC_TEST:\n' );
  fprintf ( 1, '  Demonstrate Mercer''s theorem and the KL expansion\n' );
  fprintf ( 1, '  for the rational_quadratic kernel.\n' );
%
%  Set the interval.
%
  a = 0.0;
  b = +10.0;
  fprintf ( 1, '  Using interval [%g,%g]\n', a, b );
  s_num = 21;
  s_vec = linspace ( a, b, s_num );
%
%  FRED is a function from the CHEBFUN library.
%
%  It constructs a "chebop" representing the Fredholm integral operator
%  with kernel K for functions in domain [A,B].
%
  F = fred ( @rational_quadratic_correlation, domain ( [ a, b ] ) );
%
%  EIGS has been extended to be able to compute the eigenvalues Lambda
%  and eigenfunctions Psi of the Fredholm integral operator represented by F.
%
%  The "LM" switch requests that we return the eigenvalues of largest magnitude.
%
%  Each Psi is a CHEBFUN, that is, it takes a real number argument.
%
  eigen_num = 20;
  [ Psi, Lambda ] = eigs ( F, eigen_num, 'lm' );

  eigen_found = length ( Lambda );
  fprintf ( 1, '  Requested %d eigenmodes, computed %d\n', eigen_num, eigen_found );
  eigen_num = min ( eigen_num, eigen_found );
%
%  Print the eigenvalues.
%
  lambda_vec = diag ( Lambda );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   I       Lambda(I)\n' );
  fprintf ( 1, '\n' );
  for i = 1 : eigen_num
    fprintf ( 1, '  %2d  %10g\n', i, lambda_vec(i) );
  end
%
%  Plot the eigenvalues.
%
  figure ( 1 )
  clf
  plot ( lambda_vec, 'Linewidth', 2 );
  hold on
  plot ( lambda_vec, 'b.', 'Markersize', 20 );
  title ( 'rational_quadratic: Mercer eigenvalues' );
  xlabel ( '<--- N --->')
  grid on
  print -dpng 'rational_quadratic_figure1.png'
%
%  Plot selected eigenfunctions.
%
  figure ( 2 )
  subplot ( 4, 1, 1 )
  plot ( Psi(:,1), 'Linewidth', 2 );
  title ( 'rational_quadratic: Mercer eigenfunction PSI(1)' )
  grid on
  subplot ( 4, 1, 2 )
  plot ( Psi(:,2), 'Linewidth', 2 );
  title ( 'rational_quadratic: Mercer eigenfunction PSI(2)' )
  grid on
  subplot ( 4, 1, 3 )
  plot ( Psi(:,5), 'Linewidth', 2 );
  title ( 'rational_quadratic: Mercer eigenfunction PSI(5)' )
  grid on
  subplot ( 4, 1, 4 )
  plot ( Psi(:,10), 'Linewidth', 2 );
  title ( 'rational_quadratic: Mercer eigenfunction PSI(10)' )
  grid on
  print -dpng 'rational_quadratic_figure2.png'
%
%  Orthonormality check.
%
  ptp = Psi' * Psi;
  error_frobenius = r8mat_is_identity ( eigen_num, ptp );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Frobenius norm of I - Psi'' * Psi = %g\n', error_frobenius );
%
%  K(S,S) should be exactly 1.
%  Because we are using a truncated representation of K, our estimate of K(S,S)
%  will be smaller than 1.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Truncated estimate of K(s,s) = 1 for S in the interval.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '        S        K(s,s) estimate\n' );
  fprintf ( 1, '\n' );
  for i = 1 : 21
    s = s_vec(i);
    ptlp = Psi(s,:) * Lambda * Psi(s,:)';
    fprintf ( 1, '  %10g  %14g\n', s, ptlp );
  end
%
%  Look at eigenvalue decay.
%
  x = diff ( log ( ( 1:eigen_num ) ) );
  y = diff ( log ( lambda_vec ) )';
  c = y ./ x;
  figure ( 3 );
  clf
  plot ( c, 'Linewidth', 2 );
  grid on
  hold on
  plot ( c, 'b.', 'Markersize', 20 );
  xlabel ( '<-- N -->' )
  title ( 'rational_quadratic: Eigenvalue decay rate' );
  print -dpng 'rational_quadratic_figure3.png'
%
%  Look at eigenvalue sum.
%
%  The trace of K(s,s) over [a,b] should be the integral of 1 over [a,b],
%  that is, b - a.
%
%  We compare the trace to the partial sums of the lambda's to see how much
%  of the variance of the process we have captured.
%
  trace_K = b - a;
  lambda_cum = cumsum ( lambda_vec ) / trace_K;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Index    Cumulative Eigenvalue sum\n' );
  fprintf ( 1, '\n' );
  for i = 1 : eigen_num
    fprintf ( 1, '  %5d  %10g\n', i, lambda_cum(i) );
  end
%
%  We decide to use the first 10 eigenfunctions.
%
  eigen_use = 10;
%
%  Find 400 realizations of the process by selecting, for each realization,
%  10 random parameters Z in the truncated KL expansion.
%
  Z = randn ( eigen_use, 400 );
  X = Psi(:,1:eigen_use) * ( sqrt ( Lambda(1:eigen_use,1:eigen_use) ) * Z );
%
%  Plot 40 of the realizations;
%  Plot their mean, computed from all 400.
%
  figure ( 4 )
  clf
  plot ( X(:,1:40) )
  mu = sum ( X, 2 ) / 400;
  hold on;
  plot ( mu, 'k', 'linewidth', 3 );
  title ( 'rational_quadratic: 40 Random Realizations X(t,omega), and their Mean.' )
  print -dpng 'rational_quadratic_figure4.png'
%
%  Estimate the covariance from the data.
%
  figure ( 5 )
  clf
  [ S, T ] = meshgrid ( s_vec, s_vec );
  C = cov ( X(s_vec,:)' );
  mesh ( S, T, C );
  hold on;
  D = rational_quadratic_correlation ( S, T );
  plot3 ( S, T, D, 'k.', 'markersize', 10 )
  title ( 'rational_quadratic: Covariance K(S,T) (dots), and Estimate from Realizations (mesh)' )
  print -dpng 'rational_quadratic_figure5.png'
%
%  Using just 10 functions in the expansion,
%  reduce the correlation length,
%  and examine the sum of the lambda's.
%
  if ( 0 )
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Use a fixed number of eigenfunctions = %d\n', eigen_use );
  fprintf ( 1, '  but vary the correlation length RHOBAR.\n' );
  fprintf ( 1, '  (We used RHOBAR = 1 above.)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The sum of the eigenvalues, divided by (B-A),\n' );
  fprintf ( 1, '  discloses the relative amount of the total variation\n' );
  fprintf ( 1, '  that is captured by the truncated expansion.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      RHOBAR      VARSUM\n' );
  fprintf ( 1, '\n' );
  rhobar = 4.0;
  for i = 1 : 10
    K = @ ( s, t ) exp ( - abs ( s - t ) / rhobar );
    F = fred ( K, domain ( [ a, b ] ) );
    lambda_vec = eigs ( F, 20, 'lm' );
    varsum = sum ( lambda_vec(1:eigen_use) ) / ( b - a );
    fprintf ( 1, '  %10g  %10g\n', rhobar, varsum );
    rhobar = rhobar / 2.0;
  end

  end
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'RATIONAL_QUADRATIC_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../chebfun' )

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef LaunchVehicleAeroState < matlab.mixin.SetGet & matlab.mixin.Copyable
    %LaunchVehicleAeroState Summary of this class goes here
    %   Detailed explanation goes here

    properties
        %drag
        dragCoeffModel DragCoeffModel

        %lift
        liftCoeffModel LiftCoeffModel
    end

    %old properties
    properties(Access=private)
        %drag
        Cd(1,1) double = 0.3;

        %drag
        area(1,1) double = 1; %m^2
        CdInterp %should be (1,1) griddedInterpolant
        CdIndepVar (1,1) AeroIndepVar = AeroIndepVar.Altitude;
        CdInterpMethod (1,1) GriddedInterpolantMethodEnum = GriddedInterpolantMethodEnum.Linear;
        CdInterpPts(1,1) GriddedInterpolantPointSet

        %lift
        useLift(1,1) logical = false;
        areaLift(1,1) double = 16.2;
        Cl_0(1,1) double = 0.731;
        bodyLiftVect(3,1) double = [0;0;-1];
    end

    methods
        function obj = LaunchVehicleAeroState()
            obj.dragCoeffModel = DragCoeffModel();
            obj.liftCoeffModel = LiftCoeffModel();
        end

        function CdA = getDragCoeff(obj, ut, rVect, vVect, bodyInfo, mass, altitude, pressure, density, vVectECEFMag, totalAoA, aoa, sideslip)
            arguments
                obj(1,1) LaunchVehicleAeroState
                ut(1,1) double
                rVect(3,1) double
                vVect(3,1) double
                bodyInfo(1,1) KSPTOT_BodyInfo
                mass(1,1) double
                altitude(1,1) double
                pressure(1,1) double %kPa
                density(1,1) double
                vVectECEFMag(1,1) double
                totalAoA(1,1) double = 0;
                aoa(1,1) double = 0;
                sideslip(1,1) double = 0;
            end

            CdA = obj.dragCoeffModel.getDragCoeff(ut, rVect, vVect, bodyInfo, mass, altitude, pressure, density, vVectECEFMag, totalAoA, aoa, sideslip);
        end

        function [ClS, liftUnitVectInertial] = getLiftCoeffAndDir(obj, ut, rVect, vVect, bodyInfo, mass, altitude, pressure, density, vVectECEF, attState)
            arguments
                obj(1,1) LaunchVehicleAeroState
                ut(1,1) double
                rVect(3,1) double
                vVect(3,1) double
                bodyInfo(1,1) KSPTOT_BodyInfo
                mass(1,1) double
                altitude(1,1) double
                pressure(1,1) double
                density(1,1) double
                vVectECEF(3,1) double
                attState(1,1) LaunchVehicleAttitudeState
            end

            [ClS, liftUnitVectInertial] = obj.liftCoeffModel.getLiftCoeffAndDir(ut, rVect, vVect, bodyInfo, mass, altitude, pressure, density, vVectECEF, attState);
        end

        function newAeroState = deepCopy(obj)
            newAeroState = obj.copy();
        end
    end

    methods(Static)
        function obj = loadobj(obj)
            arguments
                obj(1,1) LaunchVehicleAeroState
            end

            if(isempty(obj.CdInterp))
                pointSet = GriddedInterpolantPointSet();
                obj.CdInterpMethod = GriddedInterpolantMethodEnum.Linear;

                pointSet.addPoint(GriddedInterpolantPoint(0,obj.Cd));
                pointSet.addPoint(GriddedInterpolantPoint(1,obj.Cd));
                obj.CdInterp = pointSet.getGriddedInterpFromPoints(obj.CdInterpMethod, GriddedInterpolantMethodEnum.Nearest);
                obj.CdInterpPts = pointSet;
            end

            if(isempty(obj.dragCoeffModel))
                obj.dragCoeffModel = DragCoeffModel.createFromOldDragInputs(obj.area, obj.CdInterp, obj.CdIndepVar, obj.CdInterpMethod, obj.CdInterpPts);
            end

            if(isempty(obj.liftCoeffModel))
                obj.liftCoeffModel = LiftCoeffModel();
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [pc, x3, y3, z3] = getPointCloud(imName)
	paths = getPaths();
	dt = load(fullfile(paths.pcDir, strcat(imName, '.mat')));
	x3 = dt.x3;
	y3 = dt.y3;
	z3 = dt.z3;
	pc(:,:,1) = x3;
	pc(:,:,2) = y3;
	pc(:,:,3) = z3;
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
function c = upc_check_digit ( p, l, r )

%*****************************************************************************80
%
%% UPC_CHECK_DIGIT returns the check digit of a UPC.
%
%  Discussion:
%
%    UPC stands for Universal Price Code.
%
%    A full UPC is a string of 12 digits, in groups of size 1, 5, 5, and 1,
%    of the form P-LLLLL-RRRRR-C, where:
%
%      P is the one-digit product type code.
%      L is the five-digit manufacturer code.
%      R is the five_digit product code
%      C is the check digit.
%
%  Example:
%
%    0-72890-00011-8
%    0-12345-67890-5
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer P, the one-digit product type code.
%
%    Input, integer L, the five-digit manufacturer code.
%
%    Input, integer R, the five-digit product code.
%
%    Output, integer C, the check digit.
%
  if ( p < 0 | 9 < p )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'UPC_CHECK_DIGIT - Fatal error!\n' );
    fprintf ( 1, '  P < 0 or 9 < P!\n' );
    error ( 'UPC_CHECK_DIGIT - Fatal error!' );
  end

  if ( l < 0 | 99999 < l )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'UPC_CHECK_DIGIT - Fatal error!\n' );
    fprintf ( 1, '  L < 0 or 99999 < L!\n' );
    error ( 'UPC_CHECK_DIGIT - Fatal error!' );
  end

  if ( r < 0 | 99999 < r )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'UPC_CHECK_DIGIT - Fatal error!\n' );
    fprintf ( 1, '  R < 0 or 99999 < R!\n' );
    error ( 'UPC_CHECK_DIGIT - Fatal error!' );
  end

  lc = i4_to_digits_decimal ( l, 5 );
  rc = i4_to_digits_decimal ( r, 5 );

  c = ( p + lc(2) + lc(4) + rc(1) + rc(3) + rc(5) ) * 3 ...
          + lc(1) + lc(3) + lc(5) + rc(2) + rc(4);

  c = mod ( c, 10 );

  c = mod ( 10 - c, 10 );

  return
end


text:
function [flag] = iswithin (X, bounds)
% legacy code
s = warning ('on');
warning ('*** The function `iswithin` is now deprecated and has been renamed `VBA_isInRange` (same syntax).')
warning (s);

% fallback
flag = VBA_isInRange (X, bounds);

author: MBB-team

repo: VBA-toolbox


text:
function [opts,imdb] = DeepInstCoseg(imdb, varargin)
opts.GPUID = 1;
opts.DatasetDir = [];
opts.DatasetName = [];
opts.expDir = [];
opts.ClassName = [];
opts.Resolution = 448;
opts.Lambda = [1 1];
opts.learningRate = 0.000001 * ones(1, 20);
opts.InitModelPath = [];
opts.RootDir = [];
opts.SalDir = [];
opts.solver = @adam;
opts.InitNetPath = [];
opts.RemoveLayerIndex = 34;
[opts, varargin] = vl_argparse(opts, varargin) ;
assert(~isempty(opts.SalDir), 'Provide Sal Data!!!')
if isempty(opts.expDir)
    opts.expDir = [opts.RootDir opts.DatasetName '/' opts.ClassName '/'];
end

if isempty(opts.InitModelPath)
    opts.InitModelPath = [opts.RootDir 'models/'];
end
New_mkdir(opts.expDir);
New_mkdir(opts.InitModelPath);
opts.InitModelName = [opts.InitModelPath '/InitModel.mat'];
opts.ImageDir = [opts.DatasetDir '/' opts.DatasetName '/Image/' opts.ClassName];
opts.ClassifierPath = [opts.RootDir '/models/imagenet-vgg-verydeep-16.mat'];

% training options (SGD)
opts.train.numSubBatches = 1;
opts.train.batchSize = 6;
opts.train.continue = false ;
opts.train.gpus = opts.GPUID;
opts.train.expDir = opts.expDir ;
opts.train.learningRate = opts.learningRate;
opts.train.SaveEachEpoch = false;
opts.train.numEpochs = numel(opts.train.learningRate);
opts.train.derOutputs = {'SalLoss', 1, 'CopeakLoss', opts.Lambda(1), 'AffinityLoss', opts.Lambda(2)};
% the balanced weights presented in Eq. (1)

opts.train.solver = opts.solver;
opts = vl_argparse(opts, varargin) ;

IsModelExist = exist([opts.expDir '/' sprintf('net-epoch-%d.mat', numel(opts.train.learningRate))], 'file');

if ~IsModelExist || nargout == 2
    if isempty(imdb)
        imdb = SetupTrainDataset(opts);
    end
    if length(imdb.images.set) == 1
        return
    end
    train = find(imdb.images.set == 1);
    val = train;

    if IsModelExist
        return;
    end
    if ~exist(opts.InitModelName, 'file')
        Net = ModifyNet(opts);
        Net = Net.saveobj();
        save(opts.InitModelName, '-struct', 'Net', '-v7.3');
    end
    Net = load(opts.InitModelName);
    Net = dagnn.DagNN.loadobj(Net);
    AllImageIndex = train;
    AvgImage = imresize(single(Net.meta.normalization.averageImage), [opts.Resolution opts.Resolution]);
    FunHand = @(imdb,batch)getBatch(imdb, batch, AvgImage, AllImageIndex);
    if length(opts.Lambda) == 3
        Sel = cellfun(@(x) isa(x,'dagnn.AffinityLoss'), {Net.layers.block});
        Net.layers(Sel).block.Alpha = opts.Lambda(end);
    end
    save([opts.expDir '/opts.mat'], 'opts');
    cnn_train_dag(Net, imdb, FunHand, ...
        opts.train, ....
        'train', train, ...
        'val', val, ...
        opts.train);
end
end

function net = ModifyNet(opts)
if isempty(opts.InitNetPath)
    [net, FeatOutput] = fcnInitializeSaliencyModel('sourceModelPath', opts.ClassifierPath, 'RemoveLayerIndex', opts.RemoveLayerIndex) ;
    Sel = cellfun(@(x) isa(x,'dagnn.Loss'), {net.layers.block}) | cellfun(@(x) isa(x,'dagnn.DropOut'), {net.layers.block});
    net.removeLayer({net.layers(Sel).name})
    net.renameVar('prediction', 'SalScore');
    net.addLayer('SigmoidNormalization', dagnn.Sigmoid(), 'SalScore', 'SalMap');
    net.addLayer('SalLoss', dagnn.WeightEuclidLoss, {'SalMap', 'SalImg'}, 'SalLoss');
    net.addLayer('SigmoidNormalization_2', dagnn.Sigmoid(), 'SmallScoreMap', 'SmallSalMap');
    net.rebuild();
    InputStr = 'SmallScoreMap';
else
    net = load(opts.InitNetPath);
    net = dagnn.DagNN.loadobj(net.net);
end
Sel = cellfun(@(x) isa(x,'dagnn.PlanePeakGen'), {net.layers.block});
if all(~Sel)
    net.addLayer('PlanePeakMapsGen', ...
        dagnn.PlanePeakGen(), InputStr, {'PosPlainPeakValues'});
end
net.addLayer('FeatNorm', dagnn.LPNormalization(), FeatOutput, 'FeatNorm');

net.addLayer('SelfFeatCorrelation', dagnn.CorrelateCross, 'FeatNorm', 'PixelAffinity');
net.addLayer('SelfSalCorrelation', dagnn.CorrelateCross, 'SmallSalMap', 'SalAffinity');
net.addLayer('SelfSalDiff', dagnn.SalDiffCross(), 'SmallSalMap', 'SalDiff');
net.addLayer('WeightPixelAffinityLoss', dagnn.AffinityLoss, {'SalAffinity', 'PixelAffinity', 'SalDiff'}, 'AffinityLoss');

net.addLayer('WeightFeatNorm', dagnn.FeatMultication(), {'FeatNorm', 'SmallSalMap'}, 'WeightFeatNorm');
net.addLayer('FeatCorrelation', dagnn.Correlate, {'WeightFeatNorm'}, 'FourDTensor');
net.addLayer('CoPeakGen', dagnn.CoPeakGen, {'FourDTensor'}, 'CoPeakValues');
net.addLayer('CopeakLoss', dagnn.CopeakLoss, {'CoPeakValues'}, 'CopeakLoss');
net.rebuild();
end

function y = getBatch(imdb, images, AvgImage, AllImageIndex)
if mod(length(images), 2) == 1
    ImageIndex = find(~ismember(AllImageIndex, images));
    RandomIndex= randperm(length(AllImageIndex) - length(images));
    images = [images ImageIndex(RandomIndex(1))];
end
AvgImage = gpuArray(AvgImage);
ImageList = bsxfun(@minus, single(gpuArray(imdb.images.images(:,:,:,images))), AvgImage);
SalImgList = im2single(gpuArray(imdb.images.SalImg(:,:,:,images)));
y = {'input', ImageList, 'SalImg', SalImgList};
end

function imdb = SetupTrainDataset(opts)
ImageList = dir([opts.ImageDir '/*.jpg']);
SalDir = [opts.SalDir '/' opts.ClassName '/'];
NumImages = length(ImageList);
imdb.images.images = zeros([opts.Resolution opts.Resolution 3 NumImages], 'uint8');
imdb.images.SalImg = zeros([opts.Resolution opts.Resolution 1 NumImages], 'uint8');
for i = 1:NumImages
    disp(['Reading Image: ' num2str(i)])
    Image = imread([opts.ImageDir '/' ImageList(i).name]);
    if size(Image, 3) == 1
        Image = repmat(Image, [1 1 3]);
    end
    imdb.images.images(:,:,:,i) = imresize(Image, [opts.Resolution opts.Resolution]);
    [~,ImgName,~] = fileparts(ImageList(i).name);
    SaImg = load([SalDir ImgName '.mat']);
    imdb.images.SalImg(:,:,:,i) = imresize(SaImg.SalMap, [opts.Resolution opts.Resolution]);

end
imdb.images.set = true(1, NumImages);
imdb.images.name = {ImageList.name};
end

function [net, FeatOutput] = fcnInitializeSaliencyModel(varargin)
NumClasses = 1;
opts.sourceModelUrl = 'http://www.vlfeat.org/matconvnet/models/imagenet-vgg-verydeep-16.mat' ;
opts.sourceModelPath = 'data/models/imagenet-vgg-verydeep-16.mat' ;
opts.RemoveLayerIndex = [];
opts = vl_argparse(opts, varargin) ;

if ~exist(opts.sourceModelPath,'file')
  fprintf('%s: downloading %s\n', opts.sourceModelUrl) ;
  mkdir(fileparts(opts.sourceModelPath)) ;
  urlwrite('http://www.vlfeat.org/matconvnet/models/imagenet-vgg-verydeep-16.mat', opts.sourceModelPath) ;
end
net = vl_simplenn_tidy(load(opts.sourceModelPath)) ;

net.meta.cudnnOpts = {'cudnnworkspacelimit', 512 * 1024^3} ;

if ~isempty(opts.RemoveLayerIndex)
    net.layers = net.layers(1:opts.RemoveLayerIndex-1);
end

net = dagnn.DagNN.fromSimpleNN(net, 'canonicalNames', true) ;

net.layers(5).block.pad = [0 1 0 1] ;
net.layers(10).block.pad = [0 1 0 1] ;
net.layers(17).block.pad = [0 1 0 1] ;
net.layers(24).block.pad = [0 1 0 1] ;
if opts.RemoveLayerIndex > 31
    net.layers(31).block.pad = [0 1 0 1] ;
end
if opts.RemoveLayerIndex > 32
    net.layers(32).block.pad = [3 3 3 3] ;
end

% remove the last layers to reduce the model size
switch opts.RemoveLayerIndex
    case 32
        InputDim = 512;
        FeatOutput = 'x31';
    case 34
        InputDim = 4096;
        FeatOutput = 'x33';
end
for i = 1:numel(net.layers)
  if (isa(net.layers(i).block, 'dagnn.Conv') && net.layers(i).block.hasBias)
    filt = net.getParamIndex(net.layers(i).params{1}) ;
    bias = net.getParamIndex(net.layers(i).params{2}) ;
    net.params(bias).learningRate = 2 * net.params(filt).learningRate ;
  end
end

net.addLayer('Conv1X1', dagnn.Conv('size', [1 1 InputDim 1]), net.layers(end).outputs, 'SmallScoreMap', {'Conv1X1_f', 'Conv1X1_b'});
f = net.getParamIndex('Conv1X1_f') ;
net.params(f).value = zeros(1, 1, InputDim, 1, 'single');
net.params(f).learningRate = 1;
net.params(f).weightDecay = 1 ;
f = net.getParamIndex('Conv1X1_b') ;
net.params(f).value = zeros(1, 1, 1, 1, 'single');
net.params(f).learningRate = 2;
net.params(f).weightDecay = 1 ;

upsample = 32;
filters = single(bilinear_u(upsample * 2, NumClasses, NumClasses)) ;
net.addLayer('deconv32', ...
  dagnn.ConvTranspose(...
  'size', size(filters), ...
  'upsample', upsample, ...
  'crop', [upsample/2 upsample/2 upsample/2 upsample/2], ...
  'numGroups', NumClasses, ...
  'hasBias', false, ...
    'opts', net.meta.cudnnOpts), ...
  'SmallScoreMap', 'prediction', 'deconvf') ;

f = net.getParamIndex('deconvf') ;
net.params(f).value = filters ;
net.params(f).learningRate = 0 ;
net.params(f).weightDecay = 1 ;
net.vars(net.getVarIndex('prediction')).precious = 1 ;

net.addLayer('objective', ...
  dagnn.Loss('loss', 'softmaxlog'), ...
  {'prediction', 'label'}, 'objective') ;

net.addLayer('accuracy', ...
   dagnn.Loss(), ...
  {'prediction', 'label'}, 'accuracy') ;
end

author: KuangJuiHsu


text:
function t = subsasgn(t,s,b)
%SUBSASGN Subscripted assignment for tenmat.
%
%   Examples
%   X = tenmat(rand(3,4,2),1);
%   X(1:2,1:2) = ones(2,2); <-- Calls SUBSASGN
%
%   See also TENMAT, TENMAT/SUBSREF.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

switch s.type
    case '()'
        [m n] = size(t.data);
	t.data(s.subs{:}) = b;
        if ~isequal([m n],size(t.data))
            error('Ambiguous change in size')
        end
    otherwise
        error('Invalid assignment for tenmat.')
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function DEM_demo_fMRI_HMM
% Demonstration of Hidden Markov models for fMRI
%__________________________________________________________________________
%  This demonstration routine illustrates the modelling of state
%  transitions generating resting state fMRI timeseries. The hidden states
%  are modelled as a hidden Markov model, where each state corresponds to a
%  particular  point in the parameter space of effective connectivity. This
%  effective connectivity then generates complex cross spectral data
%  features  of the observed timeseries. Model specification requires prior
%  constraints on the probability transition matrix among hidden states,
%  which implicitly specifies the number of hidden states. The user also
%  has to specify the number of windows for epochs to apply to the
%  timeseries, where each epoch  places a lower bound on the duration of
%  each (discrete) state.
%     We first generate synthetic data using regular transitions among
%  three hidden states  (C.F., a discrete version of a heteroclinic
%  cycle  for orbit). The data are then converted by a routine that
%  combines a parametric empirical Bayesian model and a hidden Markov model
%  (as implemented as a special case of a Markov decision process). This
%  inversion is repeated for each model specified in terms of the
%  transition matrices (as prior Dirichlet concentration parameters).
%  Setting a prior transition parameter to 0 precludes that transition. In
%  this way, several different models of transitions and number of hidden
%  states can be  scored in terms of  the variational free energy.
%     Following inversion, the results are plotted in terms of expected
%  state transitions, fluctuations in connections that are allowed to
%  change (specified in the usual way by DCM.b), the deviations in
%  connectivity associated with each hidden state and the expected
%  probability transition matrix.
%     Finally, we consider Bayesian model comparison in terms of group
%  differences (here, simply the difference between the first and second
%  simulated subject).  Bayesian model comparison is simple to do in this
%  context  by comparing the free energy of a hidden Markov model in which
%  both groups share the same state dependent connections and transition
%  probabilities, with two independent models. These can be evaluated
%  efficiently using Bayesian model reduction implicit in PEB. in this
%  example, we did not introduce any differences between the two groups
%  (i.e., subjects) and therefore expected to infer no group effect.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: DEM_demo_fMRI_HMM.m 7679 2019-10-24 15:54:07Z spm $

% (I) Simulate fMRI timeseries
%==========================================================================
rng('default')

%  Assume we have P sessions with N epochs of T scans with a TR of 2 secs:
% -------------------------------------------------------------------------
%  These epochs could be from a single subject or result from the
%  concatenation of multiple sessions (under the assumption that they share
%  the same modes of connectivity).
% -------------------------------------------------------------------------
P  = 2;                               % number of sessions (e.g., subjects)
S  = 3;                               % number of latent (hidden) states
N  = 9;                               % number of epochs or windows
T  = 128;                             % number of observations (per epoch)
TR = 2;                               % repetition time or timing
t  = (1:(T*N))*TR;                    % observation times (seconds)
n  = 3;                               % number of regions or nodes

% setup model for generating timeseries
% -------------------------------------------------------------------------
options.nonlinear  = 0;
options.two_state  = 0;
options.stochastic = 0;
options.induced    = 1;

% get priors to generate simulated data
% -------------------------------------------------------------------------
a   = ones(n,n);
b   = zeros(n,n,0);
c   = zeros(n,n);
d   = zeros(n,n,0);
pP  = spm_dcm_fmri_priors(a,b,c,d,options);

% average parameters - a simple hierarchy of three nodes
% -------------------------------------------------------------------------
pP.A = [  0  -.3    0;
         .4    0  -.1;
          0   .3    0];
pP.C = eye(n,n);
pP.transit = randn(n,1)/64;

% generate spectral density of neuronal fluctuations and observation noise
% -------------------------------------------------------------------------
[Gu,Gn,Hz,dt] = spm_csd_fmri_gu(pP,TR);
Gu    = Gu(:,1,1)*ones(1,n);
Gn    = Gn(:,1,1)*ones(1,n);

% specify and generate Markovian succession of hidden states
%==========================================================================

% connections associated with hidden states: here, intrinsic connectivity
% -------------------------------------------------------------------------
B        = zeros(n,n,S);
B(1,1,1) = 1/2;
B(2,2,2) = 1/2;
B(3,3,3) = 1/2;
pP.B     = any(B,3);                     % state-dependent connections

% generate sequence of hidden states: here, a simple orbit
% -------------------------------------------------------------------------
bb    = spm_speye(S,S,-1); bb(1,S) = 1;
o     = kron(ones(1,N),1:S);
o     = o(1:N);

% State-dependent deviations in connectivity from average
% -------------------------------------------------------------------------
for i = 1:N
    tB(:,:,i) = B(:,:,o(i));
end
for i = 1:S
    for j = 1:S
        tB(i,j,:) = tB(i,j,:) - mean(tB(i,j,:));
    end
end

%% simulate epoch-specific responses to endogenous fluctuations
%==========================================================================
M.x   = sparse(n,5);
M.f   = 'spm_fx_fmri';
X     = cell(N,1);
Y     = cell(N,1);
E     = cell(N,1);
u     = cell(N,1);
for p = 1:P

    % parameters for this epoch, plus a small random effect
    % ---------------------------------------------------------------------
    gu    = spm_rand_power_law(Gu,Hz,dt,N*T);
    ge    = spm_rand_power_law(Gn,Hz,dt,N*T);
    for s = 1:N

        % parameters for this epoch, plus a small random effect
        % -----------------------------------------------------------------
        tP   = pP;
        tP.A = tP.A + tB(:,:,s);
        tP.A = tP.A.*(1 + randn(n,n)/64);
        tP.C = eye(n,n);

        % integrate states with endogenous fluctuations (gu)
        % -----------------------------------------------------------------
        j    = (1:T) + (s - 1)*T;
        M.f  = 'spm_fx_fmri';
        U.u  = gu(j,:);
        U.dt = TR;
        x    = spm_int_J(tP,M,U);
        M.x  = spm_unvec(x(end,:),M.x);

        % haemodynamic observer function to produce BOLD signal
        % -----------------------------------------------------------------
        for i = 1:T
            y(i,:) = spm_gx_fmri(spm_unvec(x(i,:),M.x),[],tP)';
        end

        % response with observation noise (ge)
        % -----------------------------------------------------------------
        e       = ge(j,:);
        X{s}    = x;
        Y{s}    = y + e;
        E{s}    = e;
        u{s}    = U.u;
        TP(s,p) = tP;

    end

    % concatenate epochs into a single timeseries
    %----------------------------------------------------------------------
    xY.dt = TR;
    xY.y  = spm_cat(Y);
    xY.u  = spm_cat(u);
    xY.X  = spm_cat(X);
    xY.E  = spm_cat(E);

    % and create DCM cell array
    %----------------------------------------------------------------------
    DCM{1,p}.options = options;
    DCM{1,p}.a = logical(pP.A);
    DCM{1,p}.b = logical(pP.B);
    DCM{1,p}.c = zeros(n,0);
    DCM{1,p}.d = zeros(n,n,0);
    DCM{1,p}.Y = xY;

end

%% show simulated responses and windows
%--------------------------------------------------------------------------
spm_figure('Getwin','Figure 1'); clf

subplot(3,2,1), plot(t,xY.u)
title('Endogenous fluctuations','FontSize',16)
xlabel('Time (seconds)'), ylabel('Amplitude'), axis square, spm_axis tight

subplot(3,2,2), hold off
plot(t,xY.X(:,(n + 1):end),'c'), hold on
plot(t,xY.X(:,1:n)),             hold off
title('Hidden states','FontSize',16)
xlabel('Time (seconds)'), ylabel('Amplitude'), axis square, spm_axis tight

subplot(3,2,3)
plot(t,xY.y,t,xY.E,':')
title('Hemodynamic response and noise','FontSize',16)
xlabel('Time (seconds)'), ylabel('Amplitude'), axis square, spm_axis tight

%  This completes the simulation of the data. We now turn to inverting the
%  data to see if one can recover the number of hidden states, the form of
%  the state transitions and the connectivity modes associated with each
%  state:
%--------------------------------------------------------------------------

% (II) Inversion under a hidden Markov model
%==========================================================================
%  Specify model space as a cell array of probability transition matrices:
%  here, the model space at the level of the HMM  will allow all
%  transitions among one to 4 hidden states. These models are specified in
%  terms of Dirichlet priors; starting with a small value of allowable
%  transitions (1/16)
%--------------------------------------------------------------------------
for i = 1:4
    b{i} = ones(i,i)/16;
end

% invert hidden Markov model: this is the routine demonstrated
%--------------------------------------------------------------------------
[HMM,CSD] = spm_dcm_HMM(DCM,N,b);

% This completes the inversion. We now just need to look at the results:
%--------------------------------------------------------------------------

%% (III) report analysis
%==========================================================================
spm_figure('Getwin','Figure 1');

%  plot windows
% -------------------------------------------------------------------------
subplot(3,2,3), hold on
for i = 1:N, plot(t,CSD{i,end}.W - 1), end, hold off

% show estimates for a single session
% -------------------------------------------------------------------------
subplot(3,2,4)
spm_plot_ci(CSD{end}.Ep,CSD{end}.Cp), hold on
bar(TP(end).A(:),1/4), hold off, axis square
title('True and MAP connections (Deterministic)','FontSize',16)

% show state-dependent changes in connectivity over sessions
% -------------------------------------------------------------------------
for i = 1:numel(CSD)
    tp(i,:) = spm_vec(TP(i).A);
    qp(i,:) = spm_vec(CSD{i}.Ep.A);
    pp(i,:) = spm_vec(HMM(S).Ep{i}.A);
end
subplot(3,3,7); imagesc(tp)
title('True connections','FontSize',16), axis square
subplot(3,3,8); imagesc(qp)
title('MAP estimates',   'FontSize',16), axis square
subplot(3,3,9); imagesc(pp)
title('PEB estimates',   'FontSize',16), axis square

% report hidden Markov model
%==========================================================================
spm_dcm_HMM_plot(HMM,S)

% And overlay true values, as cyan dots
%==========================================================================

% true state transitions
%--------------------------------------------------------------------------
x     = sparse(o,1:N,1,S,N);
x     = kron(ones(1,P),x);
N     = size(x,2);

% associate true and discovered states - and reorder
%--------------------------------------------------------------------------
r     = x*HMM(S).X';
j     = zeros(S,1);
for i = 1:S
    [d,m]  = max(r(:,i));
    j(i)   = m;
    r(m,:) = 0;
end
[o,i] = find(x(j,:));
B     = B(:,:,j);
bb    = bb(j,j);

% superimpose true values
%--------------------------------------------------------------------------
spm_figure('Getwin','HMM')

%  hidden states
%--------------------------------------------------------------------------
subplot(4,1,1), hold on
for i = 1:N, plot(i,o(i),'.c','MarkerSize',32), end, hold off

% state-dependent parameters - fluctuations
%--------------------------------------------------------------------------
subplot(4,1,2), hold on
for i = 1:N
    [j,k] = max(spm_vec(TP(i).A - pP.A));
    plot(i,k,'.c','MarkerSize',32)
end, hold off

subplot(4,1,3), hold on
for i = 1:N, pA(:,i) = spm_vec(TP(i).A); end
plot(1:N,pA(HMM(S).iP,:),'-.'), hold off, spm_axis tight

% state-dependent parameters - expectations
%--------------------------------------------------------------------------
subplot(4,2,7), hold on
for i = 1:S
    c     = spm_vec(B(:,:,i));
    [j,k] = max(c(HMM(S).iP));
    plot(i,k,'.c','MarkerSize',32)
end, hold off

% expected transition probabilities
%--------------------------------------------------------------------------
subplot(4,2,8), hold on
for i = 1:S, [j,k] = max(bb(:,i)); plot(i,k,'.c','MarkerSize',32), end
hold off

%% Bayesian model comparison in terms of group (i.e., subject) differences
%==========================================================================

% model as a single group or two separate groups
%--------------------------------------------------------------------------
hmm0 = HMM(S);

hmm1 = spm_dcm_HMM(CSD(:,1),b(S));
hmm2 = spm_dcm_HMM(CSD(:,2),b(S));

% compare the free energy of the combined groups with the combined
% free energy:
%--------------------------------------------------------------------------
F    = [hmm0.F; hmm1.F + hmm2.F];
F    = F - min(F);

% report model comparison in terms of free energy (i.e., log evidence)
%--------------------------------------------------------------------------
spm_figure('Getwin','HMM-F')
subplot(2,2,3)
bar(F,'c'),  title('Group difference','FontSize',16)
xlabel('Effect'), ylabel('Log evidence'), axis square
set(gca,'XTickLabel',{'None','Effect'})

% and show the independent maximum a posteriori estimates of state
% dependent connectivity
%--------------------------------------------------------------------------
subplot(4,2,6)
bar(hmm1.qP),  title('Group 1','FontSize',16)
xlabel('Parameter'), ylabel('Connectivity (log)'), axis square
subplot(4,2,8)
bar(hmm2.qP),  title('Group 2','FontSize',16)
xlabel('Parameter'), ylabel('Connectivity (log)'), axis square

return

author: spm


text:
%              rof_mex_demo.m  by Tom Goldstein
% This code tests the method defined by "splitBregmanROF.c".  This code
% must be compiled using the "mex" command before this demo will work.

% Step 1:  Get the test image
exact = double(imread('cameraman.tif'));
dims = size(exact);

noisy = exact+15*randn(dims);

% Step 2: Denoise the Image

clean = splitBregmanROF(noisy,.05,0.001);

% Step 3:  Display Results

close all;
figure;
subplot(2,2,1);
imagesc(exact);
colormap(gray);
title('Original');

subplot(2,2,2);
imagesc(noisy);
colormap(gray);
title('noisy');

subplot(2,2,3);
imagesc(clean);
colormap(gray);
title('denoised');

subplot(2,2,4);
imagesc(noisy-clean);
colormap(gray);
title('difference');

author: BehnoodRasti


text:
function mrk= convert_markers(mrk_old)

mrk= rmfield(mrk_old, {'pos','fs'});

mrk.time= mrk_old.pos*1000/mrk_old.fs;
mrk.event= struct;
if isfield(mrk_old, 'toe'),
  mrk.event.desc= mrk_old.toe(:);
  mrk= rmfield(mrk, {'toe'});
end

if isfield(mrk, 'indexedByEpochs'),
  mrk= rmfield(mrk, setdiff(mrk.indexedByEpochs, 'time'));
  mrk= rmfield(mrk, 'indexedByEpochs');
  nEvents= length(mrk.time);
  for Fld= mrk_old.indexedByEpochs,
    fld= Fld{1};
    fieldvar= mrk_old.(fld);
    sz= size(fieldvar);
    eventdim= find(sz==nEvents);
    if isempty(eventdim),
      error('no event information found in field %s', fld);
    end
    if length(eventdim)>1,
      error('cannot decide event dimension in field %s', fld);
    end
    if eventdim~=1,
      % permute dimensions to make the first one index events
      dimorder= [eventdim setdiff(1:length(sz), eventdim)];
      fieldvar= permute(fieldvar, dimorder);
    end
    mrk.event= setfield(mrk.event, fld, fieldvar);
  end
end

author: bbci