text:
function [minCostTuples,gain]=assign3DBB(C,maximize,boundType,initMethod,maxIter,epsVal)
%%ASSIGN3DBB Solve the data fusion axial 3D assignment problem using a
%         branch-and-bound algorithm. Such problems are NP-hard and thus
%         cannot be solved in polynomial time, so this function can only
%         solve problems of a moderate size. The optimization problem being
%         solved is minimize (or maximize)
%         \sum_{i=1}^{n1}\sum_{j=1}^{n2}\sum_{k=1}^{n3}C_{i,j,k}*\rho_{i,j,k}
%         subject to
%         \sum_{i=1}^{n1}\sum_{j=1}^{n2}\rho_{i,j,k}<=1 for all k
%         \sum_{i=1}^{n1}\sum_{k=1}^{n3}\rho_{i,j,k}<=1 for all j
%         \sum_{j=1}^{n2}\sum_{k=1}^{n3}\rho_{i,j,k} =1 for all i
%         \rho_{i,j,k} = 0 or 1
%         assuming that n1<=n2<=n3, and C is and n1Xn2Xn3 cost matrix.
%         This is equivalent to the optimization problem
%         min (or max) sum_{i} C(i,phi_2(i),phi_3(i))
%         where phi_2 and phi_3 are length n1 arrangements of n2 and n3
%         items over which the minimization is performed. In general, it is
%         not required that n1<=n2<=n3 (for any other ordering, the
%         equality constraint is on the dimension with the fewest
%         elements), so the solution is given as tuples solving the problem
%         min (or max) sum_{i} C(tuple(1,i),tuple(2,i),tuple(3,i))
%         where the tuples satisfy all of the above constraints.
%
%INPUTS: C An n1Xn2Xn3 cost hypermatrix. C cannot contain any NaNs and the
%          largest finite element minus the smallest element is a finite
%          quantity (does not overflow) when performing minimization and
%          where the smallest finite element minus the largest element is
%          finite when performing maximization.  Forbidden assignments can
%          be given costs of +Inf for minimization and -Inf for
%          maximization. During minimization, there should be no elements
%          with -Inf cost and no elements with +Inf cost during
%          maximization.
% maximize If true, the minimization problem is transformed into a
%          maximization problem. The default if this parameter is omitted
%          or an empty matrix is passed is false.
% boundType This selects the bound to use for the branch and bound method.
%          These correspond to the method input of the assign3DLB function.
%          The default if this parameter is omitted or an empty matrix is
%          passed is 2.
% initMethod A parameter indicating how the branch-and-bound algorithm
%          should be initialized. Possible values are:
%          0 (The default if omitted or an empty matrix is passed) Do not
%            use an initial estimate.
%          1 Use an initial solution (and lower bound) via algorithm 0 of
%            the assign3D function.
%  maxIter If initMethod=1, then this is the maximum number of iterations
%          to allow the initialization routine in assign3D to perform. If
%          initMethod~=1, then this parameter is ignored. The default if
%          omitted or an empty matrix is passed is 10.
%   epsVal If initMethod=1, then this is both the AbsTold and RelTol inputs
%          to assign3D to determine whether the function has converged in
%          terms of the relative duality gap. If initMethod~=1, then this
%          parameter is ignored. The default if omitted or an empty matrix
%          is passed is eps(1).
%
%OUTPUTS: tuples A 3Xn1 matrix where tuples(:,i) is the ith assigned tuple
%                in C as described above.
%           gain The cost value of the optimal assignment found. This is
%                the sum of the values of the assigned elements in C.
%
%The algorithm is similar to the branch-and bound procedure describes in
%[1], which uses the "classical" branching of [2]. There is a total of
%min([n1,n2,n3]) elements to choose out of the entire C matrix. Suppose
%that n1=min([n1,n2,n3]). Then, the tuples chosen will have the first index
%from 1 to n1 and the indices of the remaining two tuple components will
%vary. The branching here is on the tuple chosen for i=1 of C(i,j,k), then
%the next level is for i=2 until branching is done for i=n1. At each level
%of branching there are n2*n3 possibilities, as is illustrated in Fig. 1 of
%[1]. Each time a branch is taken, all tuples of (j,k) sharing any common
%tuples with the branch taken (and any previous branches taken) must be
%removed. At each level, lower bounds for all n2*n3 possible branches are
%computed. Those branches with lower bounds that are less than the best
%found solution are discarded.
%
%REFERENCES:
%[1] W. P. Pierskalla, "The multidimensional assignment problem,"
%    Operations Research, vol. 16, no. 2, pp. 422-431, Mar-Apr. 1968.
%[2] R. E. Burkard and R. Rudolf, "Computational investigations on 3-
%    dimensional axial assignment problems," Belgian Journal of Operations
%    Research, Statistics and Computer Science, vol. 32, no. 1-3, pp.
%    85-98, 1993.
%
%September 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<6||isempty(epsVal))
	epsVal=eps(1);
end

if(nargin<5||isempty(maxIter))
    maxIter=10;
end

if(nargin<4||isempty(initMethod))
	initMethod=0;
end

if(nargin<3||isempty(boundType))
    boundType=2;
end

if(nargin<2||isempty(maximize))
	maximize=false;
end

%For some types of lower bounds, the cost matrix must have all non-negative
%elements for the assignment algorithm to work. This forces all of the
%elements to be positive. The delta is added back in when computing the
%gain in the end.
if(maximize==true)
    CDelta=max(C(:));
    C=-C+CDelta;
else
    CDelta=min(C(:));
    C=C-CDelta;
end

%The algorthm was implemented assuming that n1<=n2<=n3 for the
%dimensionality of C. If a matrix having a different arrangement of indices
%is passed, one could permute the indices to make the assumptions below
%hold.
nVals=size(C);
%Deal with trailing singleton dimensions.
if(length(nVals)==2)
    nVals=[nVals,1];
elseif(length(nVals)==1)
    nVals=[nVals,1,1];
end

if(length(nVals)<3||~(nVals(1)<=nVals(2)&&nVals(2)<=nVals(3)))
    error('It is required that size(C,1)<=size(C,2)<=size(C,3)')
end

n1=nVals(1);
n2=nVals(2);
n3=nVals(3);

%Get an initial estimate.
switch(initMethod)
    case 0%Do not use an initial estimate.
        q=-Inf;
        gain=Inf;
        minCostTuples=[];
    case 1
        subgradMethod=0;
        maximize=false;
        [minCostTuples,gain,q,exitCode]=assign3D(C,false,subgradMethod,maximize,maxIter,epsVal,epsVal);
        if(exitCode==-3||exitCode==-1)
            %If no valid assignment was found, then we just start the
            %algorithm without initialization. This most commonly occurs
            %when the problem is infeasible and no solution will be found
            %in the end anyway.
            q=-Inf;
            gain=Inf;
            minCostTuples=[];
        end
    otherwise
        error('Invalid initMethod specified.');
end

%If the initial approximation returned a suboptimal solution.
if(abs(gain-q)>=epsVal*gain)
    numLevels=n1;%The number of things to assign.

    %There are n1*n2*n3 total tuples. However, we define a level by the
    %value of the first index of the tuple. In the initial level, there are
    %n2*n3 possible tuples (fix the value of the first index of C and the
    %other two are free). In subsequent levels, there are
    %(n2-curLevel+1)*(n3-curLevel+1) possibilities.
    maxTuplesPerLevel=n2*n3;

    %Generate all possible tuples of the second and third coordinates that
    %can be assigned. The assignment of the first coordinate is the level
    %of the recursion below.
    origFreeTuples=genAllTuples([n2-1;n3-1],false)+1;

    %The list of tuples that are candidates in the current level.
    freeTuples=zeros(2,maxTuplesPerLevel,numLevels);
    %Save the tuples that could possibly be assigned in the first level.
    freeTuples(:,:,1)=origFreeTuples;

    %These values link the modified tuples at each level to the original
    %tuples. These are important, because the ordering of the tuples in the
    %lower levels will be changed.
    freeTupleOrigIdx=zeros(maxTuplesPerLevel,numLevels);
    freeTupleOrigIdx(:,1)=1:maxTuplesPerLevel;

    %This keeps track of how many tuples are being considered at each
    %level after eliminating those that conflict with prior assignments and
    %eliminating those whose lower bounds are too high.
    numFreeTuples=zeros(numLevels,1);
    %In the top level, all possibilities are initially considered.
    numFreeTuples(1)=maxTuplesPerLevel;

    %The list of tuples that were candidates when entering the current
    %level. In the current level, some tuples might be removed due to a low
    %cost at that level. However, when going to a higher level, all tuples
    %that do not conflict with assignments at lower levels must be
    %considered. Thus, this stores the value of the tuples in the current
    %level before pruning.
    freeTuplesEntering=zeros(2,maxTuplesPerLevel,numLevels);
    freeTuplesEntering(:,:,1)=origFreeTuples;

    %Like freeTupleOrigIdx, these link the values in freeTuplesEntering to
    %the original set of tuples.
    freeTuplesEnteringOrigIdx=zeros(maxTuplesPerLevel,numLevels);
    freeTuplesEnteringOrigIdx(:,1)=1:maxTuplesPerLevel;

    %This keeps track of how many tuples are being considered at each
    %level after eliminating those that conflict with prior assignments but
    %not eliminating anything where the lower bounds are too high.
    numFreeTuplesEntering=zeros(numLevels,1);
    %In the top level, all possibilities are considered.
    numFreeTuplesEntering(1)=maxTuplesPerLevel;

    %The tuple that is assigned at each level, after the indexation has
    %been changed to relate it to the current cost submatrix at that level.
    %We don't need to save the first index here, because it is just the
    %level number. We don't need to save the assigned tuple in the maximum
    %level, because we use it right away.
    assignedTuples=zeros(2,numLevels-1);

    %This links the assigned tuple to the index of the tuple in
    %origFreeTuples. It is simpler to use this than to try to reconstruct
    %the tuple from assignedTuples, where each tuple has been modified
    %based on assignments at lower levels. Unlike assignedTuples, this
    %needs to be numLevels in size.
    assignedTupleOrigIdx=zeros(numLevels,1);

    %The cumulative assigned cost at each level. The final cost is computed
    %at the top level and does not need to be stored in this.
    cumAssignedCost=zeros(numLevels-1,1);

    %Allocate space for the bounded total costs for assigning something at
    %each level. Each level is given by the value in the first index that
    %is assigned. For each first index, there are n2*n3 values. There are
    %no bounds for the deepest level, because the cost matrix can be
    %directly used.
    boundVals=zeros(maxTuplesPerLevel,numLevels-1);

    %This holds what was the minimum cost value before going to a higher
    %level of recursion. Thus, when returning to a particular level, if the
    %minimum cost value changed, one should check the previous bounds and
    %see if any of the candidate tuples should be thrown out and never
    %visited.
    prevMinCost=zeros(numLevels-1,1);

    %Allocate space for the marginal costs at each level (the cost
    %submatrices). We could just allocate this like
    %CostMatCur=zeros(n1,n2,n3,numLevels); and then easily address each
    %level with the final index. However, that wastes a lot of space,
    %because each level removes one element from every dimension. The first
    %matrix is n1Xn2Xn3 in size, the next is (n1-1)*(n2-1)*(n3-1) in size
    %and so on until the final level which is 1X(n2-n1+1)X(n3-n1+1) in
    %size. To save space, we shall allocate the minimum amount of space
    %required, as well as an offset array, so that we know where the
    %elements of the matrices start and end so that we can later consider
    %each matrix individually.
    costMatCurStartIdx=zeros(numLevels+1,1);
    costMatCurStartIdx(1)=1;
    for curLevel=2:(numLevels+1)
        offset=curLevel-2;
        costMatCurStartIdx(curLevel)=costMatCurStartIdx(curLevel-1)+(n1-offset)*(n2-offset)*(n3-offset);
    end
    %The final value in costMatCurStartIdx is one past the end of the
    %matrix.
    CostMatCur=zeros(costMatCurStartIdx(numLevels+1)-1,1);
    %Assign the marginal cost matrix of the first level.
    CostMatCur(costMatCurStartIdx(1):(costMatCurStartIdx(2)-1))=C(:);

    %Enter the recursion (unrolled).
    curLevel=1;
    increasingLevelNumber=true;
    while(curLevel>0)
        if(curLevel==numLevels)
            %If we have reached the top level, the final assignment is just
            %choosing the most likely free tuple in the current cost
            %matrix.

            minVal=Inf;
            minIdx=1;
            for curTuple=1:numFreeTuples(curLevel)
                i=1;
                j=freeTuples(1,curTuple,curLevel);
                k=freeTuples(2,curTuple,curLevel);
                %The dimensionality of the cost matrix at the current
                %level.
                dims=[n1-curLevel+1;n2-curLevel+1;n3-curLevel+1];
                %The index within the cost matrix at the current level.
                idx=sub2ind(dims,i,j,k);
                curCost=CostMatCur(costMatCurStartIdx(curLevel)+idx-1);

                if(curCost<minVal)
                    minVal=curCost;
                    minIdx=curTuple;
                end
            end
            if(numLevels>1)
                minVal=cumAssignedCost(curLevel-1)+minVal;
            end

            if(minVal<gain)%If a new global minimum has been found.
                gain=minVal;

                assignedTupleOrigIdx(curLevel)=freeTupleOrigIdx(minIdx,curLevel);

                minCostTuples=zeros(3,numLevels);
                for i=1:numLevels
                    minCostTuples(:,i)=[i;origFreeTuples(:,assignedTupleOrigIdx(i))];
                end
            end

            increasingLevelNumber=false;
            curLevel=curLevel-1;
            continue;
        elseif(increasingLevelNumber==true)
            %We have just entered this level, so we have to first compute
            %the lower bounds for the current level.

            curTuple=1;
            while(curTuple<=numFreeTuples(curLevel))
                i=1;
                j=freeTuples(1,curTuple,curLevel);
                k=freeTuples(2,curTuple,curLevel);

                %The i,j,k indices of the tuple are with respect to indices
                %in the shrunken cost matrix CostMatCur(i,j,k,curLevel),
                %not to the set of tuples with respect to the global cost
                %matrix C. Since we are picking off rows one at a time, i
                %will always be 1.

                %The dimensionality of the cost matrix at the current
                %level.
                dims=[n1-curLevel+1;n2-curLevel+1;n3-curLevel+1];
                %The index within the cost matrix at the current level.
                idx=sub2ind(dims,i,j,k);
                curCost=CostMatCur(costMatCurStartIdx(curLevel)+idx-1);

                %If the assignment is not allowed.
                if(~isfinite(curCost))
                    curTuple=curTuple+1;
                    continue;
                end

                %Create the cost matrix at the next level up by removing
                %the appropriate entry from each dimension of the current
                %cost matrix. First, we extract the matrix for the current
                %level.
                CCur=reshape(CostMatCur(costMatCurStartIdx(curLevel):(costMatCurStartIdx(curLevel+1)-1)),[n1-curLevel+1,n2-curLevel+1,n3-curLevel+1]);
                %The indices of the cost matrix one level up.
                spanNext=costMatCurStartIdx(curLevel+1):(costMatCurStartIdx(curLevel+2)-1);
                CostMatCur(spanNext)=CCur(2:(n1-curLevel+1),[1:(j-1),(j+1):(n2-curLevel+1)],[1:(k-1),(k+1):(n3-curLevel+1)]);

                if(curLevel>1)
                    boundVals(curTuple,curLevel)=cumAssignedCost(curLevel-1)+curCost+assign3DLB(reshape(CostMatCur(spanNext),[n1-curLevel,n2-curLevel,n3-curLevel]),boundType);
                else
                    boundVals(curTuple,1)=curCost+assign3DLB(reshape(CostMatCur(spanNext),[n1-curLevel,n2-curLevel,n3-curLevel]),boundType);
                end

                %If the bound is no better than the current best solution,
                %then remove that tuple as a contender in this level. We
                %just overwrite it with the last tuple, decriment the
                %number of tuples present and then revisit the tuple with
                %the same index (which will now be the last tuple).
                if(boundVals(curTuple,curLevel)>=gain)
                    CostMatCur(costMatCurStartIdx(curLevel)+idx-1)=Inf;

                    freeTuples(:,curTuple,curLevel)=freeTuples(:,numFreeTuples(curLevel),curLevel);
                    freeTupleOrigIdx(curTuple,curLevel)=freeTupleOrigIdx(numFreeTuples(curLevel),curLevel);
                    numFreeTuples(curLevel)=numFreeTuples(curLevel)-1;
                else
                    curTuple=curTuple+1;
                end
            end

            %If there are not enough tuples left to make a full assignment.
            %The assumption here is that numLevels>1, so if we are here,
            %there must be one tuple left to assign after assigning the one
            %in this level.
            if(numFreeTuples(curLevel)==0)
                increasingLevelNumber=false;
                curLevel=curLevel-1;
                continue;
            end

            %Sort the tuples by cost.
            [boundVals(1:numFreeTuples(curLevel),curLevel),idx]=sort(boundVals(1:numFreeTuples(curLevel),curLevel),'descend');
            freeTuples(:,1:numFreeTuples(curLevel),curLevel)=freeTuples(:,idx,curLevel);
            freeTupleOrigIdx(1:numFreeTuples(curLevel),curLevel)=freeTupleOrigIdx(idx,curLevel);

            %Record the current minimum cost value.
            prevMinCost(curLevel)=gain;
        end

        %If we are here, then either increasingLevelNumber=false and we are
        %backtracking, or increasingLevelNumber=true, and we have just
        %entered this level.
        if(increasingLevelNumber==false)
            %If we are backtracking, then we remove the last tuple visited
            %in this level from consideration.
            i=1;
            j=assignedTuples(1,curLevel);
            k=assignedTuples(2,curLevel);

            %The dimensionality of the cost matrix at the current
            %level.
            dims=[n1-curLevel+1;n2-curLevel+1;n3-curLevel+1];
            %The index within the cost matrix at the current level.
            idx=sub2ind(dims,i,j,k);
            %freeTuples(:,curLevel,numFreeTuples(curLevel),curLevel)
            %contains the last assigned tuple for this level. Here, we
            %eliminate it from further consideration by setting
            %CostMatCur(i,j,k,curLevel) to Inf and reducing the count of
            %the tuple by one.
            CostMatCur(costMatCurStartIdx(curLevel)+idx-1)=Inf;
            numFreeTuples(curLevel)=numFreeTuples(curLevel)-1;

            %If the last tuple visited leaves too few tuples left for a
            %full assignment, then go up another level.
            if(numFreeTuples(curLevel)==0)
                increasingLevelNumber=false;
                curLevel=curLevel-1;
                continue;
            end

            %Also, if the minimum bound changed, remove all tuples whose
            %bounds are larger than the new minimum.
            if(prevMinCost(curLevel)~=gain)
                prevMinCost(curLevel)=gain;

                %Find all of the tuples where the bound is no better than
                %the current best solution and remove them as contenders in
                %this level.

                %Remember that when entering this level, the largest bounds
                %are at the beginning of the boundVals array, because it
                %has been sorted in decreasing order. Thus, we have to find
                %the first index where the bound is <gain.

                idx=find(boundVals(1:numFreeTuples(curLevel),curLevel)<gain,1);

                %If this gets rid of all tuples.
                if(isempty(idx))
                    curLevel=curLevel-1;
                    continue;
                elseif(idx>1)
                    %If here, remove the tuples with bounds that are too
                    %large.
                    sel1=1:(numFreeTuples(curLevel)-idx+1);
                    sel2=idx:numFreeTuples(curLevel);
                    numFreeTuples(curLevel)=numFreeTuples(curLevel)-idx+1;

                    boundVals(sel1,curLevel)=boundVals(sel2,curLevel);

                    %Set the entries in the cost matrix associated with the
                    %tuples that are being removed to Inf, so they won't be
                    %assigned.
                    for curTuple=1:(idx-1)
                        i=1;
                        j=freeTuples(1,curTuple,curLevel);
                        k=freeTuples(2,curTuple,curLevel);
                        %The dimensionality of the cost matrix at the
                        %current level.
                        dims=[n1-curLevel+1;n2-curLevel+1;n3-curLevel+1];
                        %The index within the cost matrix at the current
                        %level.
                        idx=sub2ind(dims,i,j,k);
                        CostMatCur(costMatCurStartIdx(curLevel)+idx-1)=Inf;
                    end

                    freeTuples(:,sel1,curLevel)=freeTuples(:,sel2,curLevel);
                    freeTupleOrigIdx(sel1,curLevel)=freeTupleOrigIdx(sel2,curLevel);
                end
            end

            increasingLevelNumber=true;
        end

        %The tuple having the lowest cost is visited first.
        %Assign the tuple in this current level.
        assignedTuples(:,curLevel)=freeTuples(:,numFreeTuples(curLevel),curLevel);
        assignedTupleOrigIdx(curLevel)=freeTupleOrigIdx(numFreeTuples(curLevel),curLevel);

        %Add the cost of the assigned tuple to the cumulative assigned
        %cost.
        i=1;
        j=assignedTuples(1,curLevel);
        k=assignedTuples(2,curLevel);
        %The dimensionality of the cost matrix at the current
        %level.
        dims=[n1-curLevel+1;n2-curLevel+1;n3-curLevel+1];
        %The index within the cost matrix at the current level.
        idx=sub2ind(dims,i,j,k);
        curCost=CostMatCur(costMatCurStartIdx(curLevel)+idx-1);

        if(curLevel>1)
            cumAssignedCost(curLevel)=cumAssignedCost(curLevel-1)+curCost;
        else
            cumAssignedCost(curLevel)=curCost;
        end

        %Construct the cost matrix for the next level. This means removing
        %the entire row, column, etc. in each dimension that contains the
        %assigned tuple.
        CCur=reshape(CostMatCur(costMatCurStartIdx(curLevel):(costMatCurStartIdx(curLevel+1)-1)),[n1-curLevel+1,n2-curLevel+1,n3-curLevel+1]);
        %The indices of the cost matrix one level up.
        spanNext=costMatCurStartIdx(curLevel+1):(costMatCurStartIdx(curLevel+2)-1);
        CostMatCur(spanNext)=CCur(2:(n1-curLevel+1),[1:(j-1),(j+1):(n2-curLevel+1)],[1:(k-1),(k+1):(n3-curLevel+1)]);

        %Copy the tuples that can be considered for the next and subsequent
        %levels into freeTuples(:,:,curLevel+1) and
        %freeTuplesEntering(:,:,curLevel+1).
        %This draws tuples for the next level down from
        %freeTuplesEntering(:,:,curLevel), because freeTuples(:,:,curLevel)
        %may have had tuples removed only for curLevel due to bounding. At
        %the same time as the tuples are added, the coordinates of all
        %valid tuples are adjusted to deal with the removed elements of
        %CostMatCur.

        numFreeTuples(curLevel+1)=0;
        for curTuple=1:numFreeTuplesEntering(curLevel)
            if(any(freeTuplesEntering(:,curTuple,curLevel)==assignedTuples(:,curLevel)))
                continue;
            else
                numFreeTuples(curLevel+1)=numFreeTuples(curLevel+1)+1;

                %Free tuples are indexed with respect to the shrunken
                %matrix of costs in CostMatCur(:,:,:,curLevel+1). This
                %means that for index values > the assigned index, one must
                %be subtracted.
                freeTuples(:,numFreeTuples(curLevel+1),curLevel+1)=freeTuplesEntering(:,curTuple,curLevel)-(freeTuplesEntering(:,curTuple,curLevel)>assignedTuples(:,curLevel));
                freeTupleOrigIdx(numFreeTuples(curLevel+1),curLevel+1)=freeTuplesEnteringOrigIdx(curTuple,curLevel);
            end
        end
        numFreeTuplesEntering(curLevel+1)=numFreeTuples(curLevel+1);
        freeTuplesEntering(:,1:numFreeTuples(curLevel+1),curLevel+1)=freeTuples(:,1:numFreeTuples(curLevel+1),curLevel+1);
        freeTuplesEnteringOrigIdx(1:numFreeTuples(curLevel+1),curLevel+1)=freeTupleOrigIdx(1:numFreeTuples(curLevel+1),curLevel+1);

        %Go to the next level.
        curLevel=curLevel+1;
    end
end

%Adjust the gain for the initial offset of the cost matrix.
if(maximize==true)
    gain=-gain+CDelta*n1;
else
    gain=gain+CDelta*n1;
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function L = noiseLikelihood(noise, mu, varsigma, y);

% NOISELIKELIHOOD Return the likelihood for each point under the noise model.
% FORMAT
% DESC returns the likelihoods for data points under the given noise model.
% ARG noise : the noise structure for which the likelihood is required.
% ARG mu : input mean locations for the likelihood.
% ARG varSigma : input variance locations for the likelihood.
% ARG y : target locations for the likelihood.
%
% SEEALSO : noiseParamInit, noiseLogLikelihood
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% NOISE

fhandle = str2func([noise.type 'NoiseLikelihood']);
L = fhandle(noise, mu, varsigma, y);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
% demos for ch07

%% sparse signal recovery demo
clear; close all;

d = 512; % signal length
k = 20;  % number of spikes
n = 100; % number of measurements
%
% random +/- 1 signal
x = zeros(d,1);
q = randperm(d);
x(q(1:k)) = sign(randn(k,1));

% projection matrix
A = unitize(randn(d,n),1);
% noisy observations
sigma = 0.005;
e = sigma*randn(1,n);
y = x'*A + e;

[model,llh] = rvmRegVb(A,y);
plot(llh);
m = model.w;

h = max(abs(x))+0.2;
x_range = [1,d];
y_range = [-h,+h];
figure;
subplot(2,1,1);plot(x); axis([x_range,y_range]); title('Original Signal');
subplot(2,1,2);plot(m); axis([x_range,y_range]); title('Recovery Signal');

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function afthelp() % autogenerated function wrapper
%In order to analyse the parameters in the "Modified Omori Law" for aftershock sequences, especially the p-value, please take into account the following:
 % turned into function by Celso G Reyes 2017

ZG=ZmapGlobal.Data; % used by get_zmap_globals

%1. The routines are available from the "Cumulative Number" window
%   or from "Seismicity Map Window" - for mapping of parameters.
%   Thus, the file timeplot.m and
%   update(mainmap()).m have options to call the relevant routines. Besides,
%   the cross-sectional mapping is available from "Cross-section Window".
%2. Before starting analysing the data, please execute the following
%   command in "Matlab Command Window": set(0, 'recursionlimit', 510).
%   It can be any other
%   number greater or equal than 502 instead of 510. The reason for this
%   command is that the routine for determining the parameters in the
%   Omori formula, ploop_c_and_p_calcs.m, works by recursion and the default of
%   Matlab has a maximum of 500 recursive calls.
%3. The file opened for analysis should contain the mainshock. Please
%   be sure to open the "Cumulative Number Window" at least ones before
%   starting any analysis.
%   This is because the mainshock information is contained in the matrix
%   "ZG.maepi", which is not defined until you open the "Cumulative Number
%   Window". After this,
%   please eliminate the mainshock, by filtering the catalog in magnitude.
%4. There are several files which permit the analysis of the parameters in
%   Omori formula. You can modify them on your convenience:
%   - src/declus/pvalcat.m - starts computation of p, c and k and their
%     standard deviations for the catalog displayed in the "Cumulative
%     Number Window".
%     You can choose the threshold magnitude, the starting time (in days)
%     after the mainshock and if the value of c is fixed or not. If you
%     want a fix c, you can
%     choose its value.
%   - src/declus/pvalcat2.m - determines p, c and K for different
%     magnitude thresholds and different starting times and presents them
%     as a map. You have to choose
%     the minimum and maximum threshold magnitude, the minimum and maximum
%     starting time and a step for each of them.
%   - src/bpvalgrid.m - computes the spatial (map) distribution of p, b
%     and Mc values for an aftershock sequence. Interactively you can choose
%     several parameters.
%     You can save your grid of values and use it at a later time.
%   - src/view_bpva.m - plots map of b, p, Mc values, computed before with
%     bpvalgrid.m.
%   - src/pcrossnew - compute the spatial (cross-section) distribution of p,
%     b and Mc values. Interactively you can choose several parameters.
%   - src/view_bpvs.m - plots map of b, p, Mc values, computed before with
%     pcrossnew.m.
%   - src/declus/mypval2m.m - interface between the programmes pvalcat.m,
%     pvalcat2.m or bpvalgrid.m on one hand and ploop_c_and_p_calcs.m
%     other hand.
%   - src/declus/ploop2.m - this is the routine which computes the p,c and
%     K values in the "Modified Omori Law".
%   - src/declus/ploop3.m - the same as above, but c is fixed (0 or any
%     other value can be chosen). However, if c = 0, ts should not be 0,
%     because when the time
%     approaches 0 the frequency tends to infinity. Please check the routine
%     mypval2m.m to verify ts, the minimum time in the aftershock sequence.
%   - src/adju2.m - filters the results displayed on the map (for spatial
%     variation of parameters) in function of different parameters, including
%     the standard error
%     in p. It is called from src/view_bpva.m.
%   - src/timabs.m - computes the number of minutes from a reference date.
%   - src/afthlp - opens this file (afthelp.m).
%   - afthelp.m - file with explanations regarding the analysis of aftershock
%     sequences (this file).
%5.  Please read the comments in the above mentioned files for further instructions.

% Bogdan Enescu, 5/2001.

%   References:

%  Ogata Y., Estimation of the parameters in the Modified Omori Formula for aftershock sequencies by the maximum likelihood procedure, J. Phys. Earth, 31, 115-
%            -124, 1983.
%  Wiemer, S. and Katsumata, K., Spatial variability of seismicity parameters in aftershock zones, J. Geophys. Res., 104, B6, 13135-13151, 1999.
%  Utsu, T, Ogata, Y., and Matsu'ura, R.S., The centenary of the Omori Formula for a Decay Law of Aftershock Activity., J. Phys. Earth, 43, 1-33, 1995.
%  Reasenberg, P.A., Lucile M. Jones, Earthquake hazard after a mainshock in California, Science, 243, 1173-1176, 1989.
%  Reasenberg, P.A., Lucile M. Jones, California Aftershock Hazard Forecasts, Science, 247, 345-346, 1990.
%  Reasenberg, P.A., Lucile M. Jones, Earthquake Aftershocks: Update, Science, 265, 1251-1252, 1994.
%  Ogata, Y., Seismicity analysis through point-process modelling: a review, Pure Appl. Geophys., 155, 471-507, 1999.

end

author: CelsoReyes


text:
function svgLuminanceToAlpha(s, source, result)
% Adds a feColorMatrix SVG filter that maps luminance to alpha values
% PRELIMINARY IMPLEMENTATION (Parameters may change)
%
% svgLuminanceToAlpha(s, source, result)
% Parameters:
%   s : Array of plot object handles
%   source : Any previous defined filter result string, 'SourceGraphic',
%            or 'SourceAlpha'.
%   result : String that identifies the filter result for following filter
%            stages.
for i = 1:length(s)
    userdata = get(s(i),'UserData');
    if isfield(userdata, 'svg') && isfield(userdata.svg, 'Filter')
        next = length(userdata.svg.Filter) + 1;
    else
        next = 1;
    end
    userdata.svg.Filter(next).Subfilter.Type = 'feColorMatrix';
    userdata.svg.Filter(next).Subfilter.Source = source;
    userdata.svg.Filter(next).Subfilter.Result = result;
    userdata.svg.Filter(next).Subfilter.ColorType = 'luminanceToAlpha';
    set(s(i),'UserData', userdata);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function D = spm_eeg_average_TF(S)
% Average each channel over trials or trial types, for time-frequency data
% FORMAT D = spm_eeg_average_TF(S)
%
% S         - optional input struct
%    fields of S:
% S.D           - MEEG object or filename of M/EEG mat-file with epoched TF data
% S.circularise - flag that indicates whether average is straight (0) or
%                 vector (1) of phase angles.
% S.robust      - (optional) - use robust averaging (only for power)
%                 .savew  - save the weights in an additional dataset
%                 .bycondition - compute the weights by condition (1,
%                                default) or from all trials (0)
%                 .ks     - offset of the weighting function (default: 3)
%
% Output:
% D         - MEEG object (also written to disk).
%__________________________________________________________________________
% Copyright (C) 2008-2017 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel
% $Id: spm_eeg_average_TF.m 7125 2017-06-23 09:49:29Z guillaume $

SVNrev = '$Rev: 7125 $';

%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename, SVNrev);
spm('FigName','M/EEG TF averaging'); spm('Pointer','Watch');

if ~isfield(S, 'circularise'),       S.circularise = 0;          end

%-Configure robust averaging
%--------------------------------------------------------------------------
if isstruct(S.robust)
    if ~isfield(S.robust, 'savew'),        S.robust.savew =  0;        end
    if ~isfield(S.robust, 'bycondition'),  S.robust.bycondition = 0;   end
    if ~isfield(S.robust, 'ks'),           S.robust.ks =  3;           end

    robust      = 1;
    savew       = S.robust.savew;
    bycondition = S.robust.bycondition;
    ks          = S.robust.ks;
    removebad   = S.robust.removebad;
else
    robust = 0;
end

%-Get MEEG object
%--------------------------------------------------------------------------
D = spm_eeg_load(S.D);

%-Check data type
%--------------------------------------------------------------------------
if ~strcmp(D.type, 'single')
    error('This function can only be applied to single trial data');
end

if ~strncmp(D.transformtype, 'TF', 2) % TF and TFphase
    error('This function can only be applied to time-frequency data.');
end

if strcmp(D.transformtype, 'TFphase')
    if ~isequal(S.robust, 0)
        warning('Robust averaging is not applicable to phase data and will not be used.');
        S.robust = 0;
    end
    plv = S.circularise;
end

%-Generate new MEEG object with new files
%--------------------------------------------------------------------------
Dnew = clone(D, [S.prefix fname(D)], [D.nchannels D.nfrequencies D.nsamples D.nconditions]);

if robust && savew
    Dw = clone(D, ['W' fname(D)]);
end

cl   = D.condlist;

ni = zeros(1,D.nconditions);
for i = 1:D.nconditions
    w = indtrial(D, deblank(cl{i}), 'GOOD');
    ni(i) = length(w);
    if ni(i) == 0
        warning('%s: No trials for trial type %s', D.fname, cl{i});
    end
end

goodtrials  =  indtrial(D, cl, 'GOOD');

if robust && removebad
    bad     = badsamples(D, ':', ':', ':');
end

spm_progress_bar('Init', D.nsamples, 'Samples completed');
if D.nsamples > 100, Ibar = floor(linspace(1, D.nsamples, 100));
else Ibar = [1:D.nsamples]; end
for j = 1:D.nsamples
     if robust && ~bycondition
         Y       = D(:, :, j, goodtrials);
         if removebad
             ibad = reshape(bad(:, j, goodtrials), [size(bad, 1), 1, 1, length(goodtrials)]);
             ibad = repmat(ibad, [1, D.nfrequencies, 1, 1]);
             Y(ibad) = NaN;
         end
         [Y, W1] = spm_robust_average(Y, 4, ks);

         if savew
             Dw(:, :, j, goodtrials) = W1;
         end
         W = zeros([D.nchannels D.nfrequencies 1 D.ntrials]);
         W(:, :, 1, goodtrials) = W1;
     end
    for i = 1:D.nconditions

        w = indtrial(D, deblank(cl{i}), 'GOOD');

         if isempty(w)
             continue;
         end

        %-Straight average
        %------------------------------------------------------------------
        if ~strcmp(D.transformtype, 'TFphase')
            if ~robust
                Dnew(:, :, j, i) = mean(D(:, :, j, w), 4);
             else
                 if bycondition
                     Y      = D(:, :, j, w);
                     if removebad
                         ibad = reshape(bad(:, j, w), [size(bad, 1), 1, 1, length(w)]);
                         ibad = repmat(ibad, [1, D.nfrequencies, 1, 1]);
                         Y(ibad) = NaN;
                     end
                     [Y, W] = spm_robust_average(Y, 4, ks);
                     Dnew(:, :, j, i) = Y;
                     if savew
                         Dw(:, :, j, w)   = W;
                     end
                 else
                     X = D(:, :, j, w);
                     X(isnan(X))      = 0;
                     Dnew(:, :, j, i) = ...
                         sum(W(:, :, 1, w).*X, 4)./sum(W(:, :, 1, w), 4);
                 end
            end

            %-Vector average (eg PLV for phase)
            %--------------------------------------------------------------
        else
            tmp = D(:, :, j, w);
            tmp = exp(sqrt(-1)*tmp);
            if plv
                Dnew(:, :, j, i) = abs(mean(tmp,4));
            else
                Dnew(:, :, j, i) = angle(mean(tmp,4));
            end
        end
    end
    if ismember(j, Ibar), spm_progress_bar('Set', j); end
end

spm_progress_bar('Clear');

Dnew = type(Dnew, 'evoked');

%-Update some header information
%--------------------------------------------------------------------------
Dnew = conditions(Dnew, ':', cl);
Dnew = repl(Dnew, ':', ni);

%-Display averaging statistics
%--------------------------------------------------------------------------
fprintf('%s: Number of replications per contrast:\n', Dnew.fname);      %-#
for i = 1:D.nconditions
    fprintf('  average %s: %d trials\n', cl{i}, ni(i));                 %-#
end

%-Save new evoked M/EEG dataset
%--------------------------------------------------------------------------
Dnew = Dnew.history('spm_eeg_average', S);
save(Dnew);

if robust && savew
    Dw = Dw.history('spm_eeg_average', S);
    save(Dw);
end

D = Dnew;

%-Cleanup
%--------------------------------------------------------------------------
fprintf('%-40s: %30s\n','Completed',spm('time'));                       %-#
spm('FigName','M/EEG TF averaging: done'); spm('Pointer', 'Arrow');

author: spm


text:
function tyinfo = slgausstype(GS)
%SLGAUSSTYPE Judges the type of a Gaussian model struct
%
% $ Syntax $
%   - tyinfo = slgausstype(GS)
%
% $ Arguments $
%   - GS:       the Gaussian model struct
%   - tyinfo:   the type information structure with following fields
%               - varform:  the form of variance: 'univar'|'diagvar'|'covar'
%               - sharevar: whether the variance(covariance) is shared
%               - hasinv:   whether is invvars or invcovs exists
%               - hasmw:    whether the mixture weights exist
%
% $ Remarks $
%   - The function will check the validity of the struct and will raise
%     an error for invalid models. So it can be used to check validity
%     even you don't need to know the type of the model.
%
% $ History $
%   - Created by Dahua Lin, on Aug 23rd, 2006
%

%% verify basic fields

if ~isstruct(GS)
    gs_argerror('The Gaussian model should be a struct');
end

if ~all(slisfields(GS, {'dim', 'nmodels', 'means'}))
    gs_argerror('The Gaussian model struct should have all of the fields: dim, nmodels and means');
end

d = GS.dim;
k = GS.nmodels;

if ~isequal(size(GS.means), [d k])
    gs_sizerror('The means field should be an array of size d x k');
end

%% verify variance/covariance field

if isfield(GS, 'vars')
    sizvf = size(GS.vars);
    if isequal(sizvf, [1 1])
        varform = 'univar';
        sharevar = true;
    elseif isequal(sizvf, [1 k])
        varform = 'univar';
        sharevar = false;
    elseif isequal(sizvf, [d 1])
        varform = 'diagvar';
        sharevar = true;
    elseif isequal(sizvf, [d k])
        varform = 'diagvar';
        sharevar = false;
    else
        gs_arrerror('The size of vars is illegal');
    end

    hasinv = isfield(GS, 'invvars');

elseif isfield(GS, 'covs')

    sizcvf = size(GS.covs);
    if isequal(sizcvf, [d d])
        varform = 'covar';
        sharevar = true;
    elseif isequal(sizcvf, [d d k])
        varform = 'covar';
        sharevar = false;
    else
        gs_arrerror('The size of covariance is illegal');
    end

    hasinv = isfield(GS, 'invcovs');

else
    gs_arrerror('The Gaussian struct lacks a field for variance/covariance');
end

hasmw = isfield(GS, 'mixweights');

if nargout >= 1
    tyinfo = struct(...
        'varform', varform, ...
        'sharevar', sharevar, ...
        'hasinv', hasinv, ...
        'hasmw', hasmw);
end

function gs_argerror(errmsg)

error('sltoolbox:invalidarg', errmsg);

function gs_sizerror(errmsg)

error('sltoolbox:sizmismatch', errmsg);

author: lmthang


text:
function RtMul = mulRt(RtLeft, RtRight)

RtMul = [RtLeft(:,1:3)*RtRight(:,1:3) RtLeft(:,1:3)*RtRight(:,4)+ RtLeft(:,4)];

author: jianxiongxiao

repo: ProfXkit


text:
%Revolute Robot manipulator Revolute link class
%
% A subclass of the Link class: holds all information related to a robot
% link such as kinematics parameters, rigid-body inertial parameters, motor
% and transmission parameters.
%
% Notes::
% - This is reference class object
% - Link class objects can be used in vectors and arrays
%
% References::
% - Robotics, Vision & Control, Chap 7
%   P. Corke, Springer 2011.
%
% See also Link, Prismatic, SerialLink.

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

classdef Revolute < Link
    methods
        function L = Revolute(varargin)
            L = L@Link(varargin{:});

            if nargin == 0
                L.theta = [];
            end
            L.sigma = 0;
            if isempty(L.d)
                L.d = 0;
            end
            if ~isempty(L.theta)
                error('theta cannot be specified for a revolute link');
            end
        end

    end
end

author: Allopart

repo: rbpf-gmapping


text:
function twod_to_vtu_test ( )

%*****************************************************************************80
%
%% TWOD_TO_VTU_TEST tests the TWOD_TO_VTU library.
%
%  Discussion:
%
%    The data represents a flow solution in a 10 by 3 channel.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 December 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( '\n' );
  fprintf ( 1, 'TWOD_TO_VTU_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the TWOD_TO_VTU library.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read in TWOD finite element data.\n' );
  fprintf ( 1, '  Write it out as a VTU file.\n' );
%
%  Load the data from text files.
%
  nodes = load ( 'channel_nodes.txt' );
  elements = load ( 'channel_elements.txt' );
  uvp = load ( 'channel_values.txt' );

  output_filename = 'channel.vtu';

  title = 'UVP for 10x3 2D channel problem';
%
%  Have the data written to a file.
%
  twod_to_vtu ( nodes, elements, uvp, output_filename, title );
%
%  Terminate.
%
  fprintf ( '\n' );
  fprintf ( 1, 'TWOD_TO_VTU_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [y] = crossprod(x,p,flag)

if nargin<3
  % always augment with constant term
  flag = 1;
end
if nargin<2
  % length of smoothing boxcar
  p = 1;
end
flag = double(flag);

% FIXME only works in case of dim=2
n    = size(x);
indx = tril(ones(n(1)))==1;
N    = 0.5*n(1)*(n(1)+1);
y    = zeros(N+flag, n(2));
for k = 1:n(2)
  tmp    = x(:,k)*x(:,k)';
  y(flag+(1:N),k) = tmp(indx);
end

if flag
  y(1,:) = 1;
end

if p>1
  krn = ones(1,p)./p;
  y   = convn(y, krn, 'valid');
end

% function [y] = crossprod(x, ix)
%
% %FIXME works only in case of dim=2
% n = size(x);
% y = zeros(0.5*n(1)*(n(1)+1), n(2));
% for k = 1:size(ix,1)
%   y(k,:) = x(ix(k,1),:).*x(ix(k,2),:);
% end


text:
function [fYr, nMn, nDay, nHr, nMin, nSec]=decyear2mat(fDy)
    % Calculate decimal year format to matrix with columns year, month, day, hour, minute, and second.
    % [fYr, nMn, nDay, nHr, nMin, nSec] = decyear2mat(fDy)
    % ------------------------------------------------------------------------------------------------
    % Calculate decimal year format (from zmap-function decyear, i.e. 1998.734)
    % to matrix with columns year, month, day, hour, minute, and second.
    % This function was programmed because datevec.m does not work with decimal
    % year input format
    %
    % Input parameters:
    %   fDy     Decimal year (like 1998.2515) as vector or float
    %
    % Output parameters:
    %   fYr     decimal year
    %   nMn     month
    %   nDay    day
    %   nHr     hour
    %   nMin    minute
    %   nSec    second
    %
    % Example [fYr, nMn, nDay, nHr, nMin, nSec]=decyear2mat(decyear([1989 12 31 10 35 44.3]))
    %
    % Thomas van Stiphout
    % Mai 9, 2007

    disp('~/zmap/src/decyear2mat.m')

    % save year in decimal format
    fYr=fDy;
    % define leap years
    bLeapYr = rem(fix(fDy),4) == 0 & rem(fix(fDy),100) ~= 0 | rem(fix(fDy),400) == 0 ;

    % loop over each date
    for i=1:size(bLeapYr,1)
        if bLeapYr(i) % for leap years
            mDay=[0,31,60,91,121,152,182,213,244,274,305,335]'; %leapyear
            nMn(i)=sum(mDay<rem(fDy(i),1).*366); % calculate year
            nDay(i)=ceil(rem(fDy(i),1).*366)-mDay(nMn(i)); % calculate days
            nHr(i)=rem(fDy(i),1).*366.*24-(mDay(nMn(i))+nDay(i)-1).*24;
            nMin(i)=rem(nHr(i),1).*60;
            nHr(i)=fix(nHr(i)); % hours
            nSec(i)=rem(nMin(i),1).*60; % seconds
            nMin(i)=fix(nMin(i)); % minutes

        else
            mDay= [0,31,59,90,120,151,181,212,243,273,304,334]';%cumulative days in one year
            nMn(i)=sum(mDay<rem(fDy(i),1)*365);
            nDay(i)=ceil(rem(fDy(i),1)*365)-mDay(nMn(i));
            nHr(i)=rem(fDy(i),1)*365*24-(mDay(nMn(i))+nDay(i)-1)*24;
            nMin(i)=rem(nHr(i),1)*60;
            nHr(i)=fix(nHr(i));
            nSec(i)=rem(nMin(i),1)*60;
            nMin(i)=fix(nMin(i));
        end
    end
end


text:
function [yi, p, pval] = newtint(x, y, xi, c)

% NEWTINT Interpolation of equally-spaced points.
%    NEWTINT(X,Y,XI,C) interpolates to find YI, the value
%    of the underlying function Y at the point XI, using
%    either Newton's forward interpolation formula or
%    Newton's backward interpolation formula.
%
%    C specifies the interpolation method:
%       0 : Newton's forward or backward (default)
%       1 : Newton's forward
%       2 : Newton's backward
%       3 : Stirling's method
%       4 : Everett's method
%
%    [YI,P,PVAL] = NEWTINT() also returns P, the coefficients
%    of the calculated interpolating polynomial, and PVAL,
%    which specifies the maximum degree of the interpolating
%    polynomial.

% Joe Henning - Fall 2011

if (nargin < 4)
   c = 0;
end

pval = length(x)*(0+1) - 1;

n = length(x)-1;
D = zeros(n+1,n+1);
D(:,1) = y(:);

h = x(2)-x(1);

for i = 1:n
   for j = 1:i
      D(i+1,j+1) = D(i+1,j)-D(i,j);
   end
end

if (c == 0)
   if (xi < x(1))
      c = 2;
   elseif (xi > x(ceil((n+1)/2)) && xi <= x(n+1))
      c = 2;
   else
      c = 1;
   end
end

switch c
   case 1
      % forward difference interpolation
      p = diag(D);
      q = (xi-x(1))/h;
      yi = p(1);
      for i = 1:length(p)-1
         term = 1;
         for k = 0:(i-1)
           term = term*(q-k);
         end
         term = term/factorial(i);
         yi = yi + term*p(i+1);
      end
   case 2
      % backward difference interpolation
      p = D(n+1,:).';
      q = (xi-x(n+1))/h;
      yi = p(1);
      for i = 1:length(p)-1
         term = 1;
         for k = 0:(i-1)
           term = term*(q+k);
         end
         term = term/factorial(i);
         yi = yi + term*p(i+1);
      end
   case 3
      % Stirling's method
      xk = 0;
      ik = 0;
      for i = 1:n+1
         if (xi > x(i))
            xk = x(i);
            ik = i;
         end
      end
      yi = 0;
      q = (xi-xk)/h;
      if (q > 0.5)
         ik = ik + 1;
         xk = x(ik);
         q = (xi-xk)/h;
         for i = 1:2:n+1
            fi = floor(i/2);
            term = 0;
            %[term D(ik+i-1-floor(i/2),i)]s
            for k = fi-1:fi
               term = term + dnchoosek(q+k,i-1)/2.0;
            end
            p(i,1) = term;
            if (ik+i-1-floor(i/2) <= n+1)
               yi = yi + term*D(ik+i-1-floor(i/2),i);
            end
         end
         for i = 2:2:n+1
            fi = floor(i/2);
            term = dnchoosek(q+fi-1,i-1);
            p(i,1) = term;
            %[term D(ik+i-1-floor(i/2),i) D(ik+i-floor(i/2),i)]
            if (ik+i-floor(i/2) > n+1)
               yi = yi + term*D(ik+i-1-floor(i/2),i)/2.0;
            else
               yi = yi + term*(D(ik+i-1-floor(i/2),i)+D(ik+i-floor(i/2),i))/2.0;
            end
         end
      else
         for i = 1:2:n+1
            fi = floor(i/2);
            term = 0;
            for k = fi-1:fi
               term = term + dnchoosek(q+k,i-1)/2.0;
            end
            p(i,1) = term;
            %[term D(ik+i-1-floor(i/2),i)]
            if (ik+i-1-floor(i/2) <= n+1)
               yi = yi + term*D(ik+i-1-floor(i/2),i);
            end
         end
         for i = 2:2:n+1
            fi = floor(i/2);
            term = dnchoosek(q+fi-1,i-1);
            p(i,1) = term;
            %[term D(ik+i-1-floor(i/2),i) D(ik+i-floor(i/2),i)]
            if (ik+i-floor(i/2) > n+1)
               yi = yi + term*D(ik+i-1-floor(i/2),i)/2.0;
            else
               yi = yi + term*(D(ik+i-1-floor(i/2),i)+D(ik+i-floor(i/2),i))/2.0;
            end
         end
      end
   case 4
      % Everett's method
      xk = 0;
      ik = 0;
      for i = 1:n+1
         if (xi > x(i))
            xk = x(i);
            ik = i;
         end
      end
      p = [];
      q = (xi-xk)/h;
      r = 1-q;
      yi = 0;
      for i = 1:2:n+1
         fi = floor(i/2);
         e = dnchoosek(r+fi,i);
         f = dnchoosek(q+fi,i);
         p = [p; e f];
         %[e D(ik+i-1-floor(i/2),i) f D(ik+i-floor(i/2),i)]
         if ((ik+i-floor(i/2)) > n+1)
            yi = yi + e*D(ik+i-1-floor(i/2),i) + f*0;
         else
            yi = yi + e*D(ik+i-1-floor(i/2),i) + f*D(ik+i-floor(i/2),i);
         end
      end
end

function [d] = dnchoosek(n, k)

% DNCHOOSEK Binomial coefficient or all combinations.
%    DNCHOOSEK(N,K) where N and K are non-negative numbers
%    returns N!/K!(N-K)!.

d = gamma(n+1)/(gamma(k+1)*gamma(n-k+1));

author: Sable

repo: mcbench-benchmarks


text:
function [EEG]= selectmodel(EEG, model, type)
startpoints = [];
endpoints = [];
% if type ~= 1
%     newEEG = EEG;
% end

for i = 1:EEG.nbchan
    EEG.icaweights(i,:) = EEG.etc.amica.W(i,:,model);
    EEG.icasphere(i,:)  = EEG.etc.amica.S(i,:);
    EEG.icawinv(i,:)    = EEG.etc.amica.A(i,:,model);
end

disp('Calculating...')

if type == 1  % reject the data where the model is not active
    if EEG.trials > 1
        trials2keep = [];
        for i = 1:EEG.trials
            [dummy in] = max(sum(EEG.etc.amica.v_smooth(:,:,i),2));
            if model == in
                trials2keep = [trials2keep i];
            end
        end

        EEG = pop_select(EEG,'trial',trials2keep);

    else
        nepochs = 0; %number of epochs that will be excluded
        rejecting = 0;
        startpoints = [];
        endpoints = [];
        for i = 1:size(EEG.data,2)

            [y in] = max(EEG.etc.amica.v_smooth(:,i));
            if in ~= model
                if ~rejecting
                    nepochs = nepochs + 1;
                    startpoints(nepochs) = i;
                    rejecting = 1;
                end
            else
                if rejecting == 1
                    endpoints(nepochs) = i-1;
                    rejecting = 0;
                end

            end

        end
        if rejecting == 1
            endpoints(end + 1) = size(EEG.data,2);
        end

        if nepochs == 1 && endpoints == i && startpoints == 1
            error('Selected model is never active. Original data will be preserved.');
            startpoints = [];
            endpoints = [];

            return;
        end
        if nepochs == 0
            disp('Model is always active. Original data will be preserved, ICA weights for the model is loaded.');

            return;
        end

        %-------Rejecting the data where model is not active------------------
        temp_start = startpoints;
        temp_end = endpoints;
        EEG = pop_select(EEG,'nopoint',[startpoints' endpoints']);
%         for i = 1:nepochs
%             newEEG = pop_ozselect(newEEG,'nopoint',[temp_start(i) temp_end(i)]);
%             shift = temp_end(i)-temp_start(i)+1;
%             temp_start = temp_start - shift;
%             temp_end = temp_end - shift;
%         end
        %---------------------------------------------------------------------
    end
    EEG.setname = [EEG.setname ' - Model ' num2str(model)];

else
    if type == 2
        %weight the data w.r.t to the probability of the model chosen
        if isequal(EEG.etc.amica.v_smooth(model,:,:),zeros(size(EEG.etc.amica.v_smooth(model,:,:))))
           disp('Selected model has probability of 0 through the data, old dataset is preserved') ;
        else
            if EEG.trials > 1
                for i = 1:EEG.trials
                    for j = 1:size(EEG.data,2)
                        EEG.data(:,j,i) = EEG.data(:,j,i)*EEG.etc.amica.v_smooth(model,j,i);

                    end
                end

            else
                for i = 1:size(EEG.data,2)
                    EEG.data(:,i) = EEG.data(:,i)*EEG.etc.amica.v_smooth(model,i);
                end

            end
            EEG.etc.amica.weighted = 1;
        end
        EEG.setname = [EEG.setname ' - Model ' num2str(model) ' Weighted'];
    else
        if type == 3

        end

    end

end

author: PatternRecognition

repo: OpenBMI


text:
function setfontsizefigure(fontsize)
% setfontsizefigure(fontsize)
% set font size (axes, labels) of the current figure
set(get(gca,'XLabel'),'FontSize',fontsize)
set(get(gca,'YLabel'),'FontSize',fontsize)
set(get(gca,'ZLabel'),'FontSize',fontsize)
set(gca,'FontSize',fontsize)

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


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

%%
clear; close all; clc;

%% Syntax
% |[subIndOut,L_valid]=snapSubInd(subIndIn,siz);|

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function exactness_test015 ( )

%*****************************************************************************80
%
%% EXACTNESS_TEST015 tests Fejer Type 1 rules for Legendre integrals.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 May 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'EXACTNESS_TEST015\n' );
  fprintf ( 1, '  Test Fejer Type 1 rules on Legendre integrals.\n' );
  fprintf ( 1, '  Density function rho(x) = 1.\n' );
  fprintf ( 1, '  Region: -1 <= x <= +1.\n' );
  fprintf ( 1, '  Exactness: N   for N odd,\n' );
  fprintf ( 1, '             N-1 for N even.\n' );

  for n = 1 : 5

    [ x, w ] = fejer1_set ( n );
    if ( mod ( n, 2 ) == 1 )
      p_max = n + 1;
    else
      p_max = n;
    end
    legendre_exactness ( n, x, w, p_max );

  end

  return
end


text:
%% Classify Text Data Using BERT
% This example shows how to classify text data using a pretrained BERT
% model as a feature extractor.
%
% The simplest use of a pretrained BERT model is to use it as a feature
% extractor. In particular, you can use the BERT model to convert documents
% to feature vectors which you can then use as input to train a deep
% learning classification network.
%
% This example shows how to use a pretrained BERT model to classify failure
% events given a data set of factory reports.

%% Load Pretrained BERT Model
% Load a pretrained BERT model using the |bert| function. The model
% consists of a tokenizer that encodes text as sequences of integers, and
% a structure of parameters.
mdl = bert

%%
% View the BERT model tokenizer. The tokenizer encodes text as sequences of
% integers and holds the details of padding, start, separator and mask
% tokens.
tokenizer = mdl.Tokenizer

%% Load Data
% Load the example data. The file |factoryReports.csv| contains factory
% reports, including a text description and categorical labels for each
% event.

filename = "factoryReports.csv";
data = readtable(filename,"TextType","string");
head(data)

%%
% The goal of this example is to classify events by the label in the
% |Category| column. To divide the data into classes, convert these labels
% to categorical.
data.Category = categorical(data.Category);

%%
% View the number of classes.
classes = categories(data.Category);
numClasses = numel(classes)

%%
% View the distribution of the classes in the data using a histogram.
figure
histogram(data.Category);
xlabel("Class")
ylabel("Frequency")
title("Class Distribution")

%%
% Encode the text data using the BERT model tokenizer using the |encode|
% function and add the tokens to the training data table.
data.Tokens = encode(tokenizer, data.Description);

%%
% The next step is to partition it into sets for training and validation.
% Partition the data into a training partition and a held-out partition for
% validation and testing. Specify the holdout percentage to be 20%.
cvp = cvpartition(data.Category,"Holdout",0.2);
dataTrain = data(training(cvp),:);
dataValidation = data(test(cvp),:);

%%
% View the number of training and validation observations.
numObservationsTrain = size(dataTrain,1)
numObservationsValidation = size(dataValidation,1)

%%
% Extract the text data, labels, and encoded BERT tokens from the
% partitioned tables.
textDataTrain = dataTrain.Description;
textDataValidation = dataValidation.Description;

TTrain = dataTrain.Category;
TValidation = dataValidation.Category;

tokensTrain = dataTrain.Tokens;
tokensValidation = dataValidation.Tokens;

%%
% To check that you have imported the data correctly, visualize the
% training text data using a word cloud.

figure
wordcloud(textDataTrain);
title("Training Data")

%%
% View the BERT token codes of the first few training documents.
tokensTrain{1:5}

%% Prepare Data for Training
% Convert the documents to feature vectors using the BERT model as a
% feature extractor.

% To extract the features of the training data by iterating over
% mini-batches, create a |minibatchqueue| object.

% Mini-batch queues require a single datastore that outputs both the
% predictors and responses. Create array datastores containing the training
% BERT tokens and labels and combine them using the |combine| function.
dsXTrain = arrayDatastore(tokensTrain,"OutputType","same");
dsTTrain = arrayDatastore(TTrain);
cdsTrain = combine(dsXTrain,dsTTrain);

% Create a combined datastore for the validation data using the same steps.
dsXValidation = arrayDatastore(tokensValidation,"OutputType","same");
dsTValidation = arrayDatastore(TValidation);
cdsValidation = combine(dsXValidation,dsTValidation);

%%
% Create a mini-batch queue for the training data. Specify a mini-batch
% size of 32 and preprocess the mini-batches using the
% |preprocessPredictors| function, listed at the end of the example.
miniBatchSize = 32;
paddingValue = mdl.Tokenizer.PaddingCode;
maxSequenceLength = mdl.Parameters.Hyperparameters.NumContext;

mbqTrain = minibatchqueue(cdsTrain,1,...
    "MiniBatchSize",miniBatchSize, ...
    "MiniBatchFcn",@(X) preprocessPredictors(X,paddingValue,maxSequenceLength));

%%%
% Create a mini-batch queue for the validation data using the same steps.
mbqValidation = minibatchqueue(cdsValidation,1,...
    "MiniBatchSize",miniBatchSize, ...
    "MiniBatchFcn",@(X) preprocessPredictors(X,paddingValue,maxSequenceLength));

%%
% To speed up feature extraction. Convert the BERT model weights to
% gpuArray if a GPU is available.
if canUseGPU
    mdl.Parameters.Weights = dlupdate(@gpuArray,mdl.Parameters.Weights);
end

%%
% Convert the training sequences of BERT model tokens to a
% |N|-by-|embeddingDimension| array of feature vectors, where |N| is the
% number of training observations and |embeddingDimension| is the dimension
% of the BERT embedding.

featuresTrain = [];
reset(mbqTrain);
while hasdata(mbqTrain)
    X = next(mbqTrain);
    features = bertEmbed(X,mdl.Parameters);
    featuresTrain = [featuresTrain gather(extractdata(features))];
end

%%
% Transpose the training data to have size |N|-by-|embeddingDimension|.
featuresTrain = featuresTrain.';

%%
% Convert the validation data to feature vectors using the same steps.
featuresValidation = [];

reset(mbqValidation);
while hasdata(mbqValidation)
    X = next(mbqValidation);
    features = bertEmbed(X,mdl.Parameters);
    featuresValidation = cat(2,featuresValidation,gather(extractdata(features)));
end
featuresValidation = featuresValidation.';

%% Define Deep Learning Network
% Define a deep learning network that classifies the feature vectors.

numFeatures = mdl.Parameters.Hyperparameters.HiddenSize;
layers = [
    featureInputLayer(numFeatures)
    fullyConnectedLayer(numClasses)
    softmaxLayer
    classificationLayer];

%% Specify Training Options
% Specify the training options using the |trainingOptions| function.
% * Train with a mini-batch size of 64.
% * Shuffle the data every epoch.
% * Validate the network using the validation data.
% * Display the training progress in a plot and suppress the verbose
%   output.
opts = trainingOptions('adam',...
    "MiniBatchSize",64,...
    "ValidationData",{featuresValidation,dataValidation.Category},...
    "Shuffle","every-epoch", ...
    "Plots","training-progress", ...
    "Verbose",0);

%% Train Network
% Train the network using the |trainNetwork| function.
net = trainNetwork(featuresTrain,dataTrain.Category,layers,opts);

%% Test Network
% Make predictions using the validation data and display the results in a
% confusion matrix.
YPredValidation = classify(net,featuresValidation);

figure
confusionchart(TValidation,YPredValidation)

%%
% Calculate the validation accuracy.
accuracy = mean(dataValidation.Category == YPredValidation)

%% Predict Using New Data
% Classify the event type of three new reports. Create a string array
% containing the new reports.
reportsNew = [ ...
    "Coolant is pooling underneath sorter."
    "Sorter blows fuses at start up."
    "There are some very loud rattling sounds coming from the assembler."];

%%
% Tokenize the text data using the same steps as the training documents.
tokensNew = encode(tokenizer,reportsNew);

%%
% Pad the sequences of tokens to the same length using the |padsequences|
% function and pad using the tokenizer padding code.
XNew = padsequences(tokensNew,2,"PaddingValue",tokenizer.PaddingCode);

%%
% Classify the new sequences using the trained model.
featuresNew = bertEmbed(XNew,mdl.Parameters)';
featuresNew = gather(extractdata(featuresNew));
labelsNew = classify(net,featuresNew)

%% Supporting Functions

%%% Predictors Preprocessing Functions
% The |preprocessPredictors| function truncates the mini-batches to have
% the specified maximum sequence length, pads the sequences to have the
% same length. Use this preprocessing function to preprocess the predictors
% only.
function X = preprocessPredictors(X,paddingValue,maxSeqLen)

X = truncateSequences(X,maxSeqLen);
X = padsequences(X,2,"PaddingValue",paddingValue);

end

%%% BERT Embedding Function
% The |bertEmbed| function maps input data to embedding vectors and
% optionally applies dropout using the "DropoutProbability" name-value
% pair.
function Y = bertEmbed(X,parameters,args)

arguments
    X
    parameters
    args.DropoutProbability = 0
end

dropoutProbabilitiy = args.DropoutProbability;

Y = bert.model(X,parameters, ...
    "DropoutProb",dropoutProbabilitiy, ...
    "AttentionDropoutProb",dropoutProbabilitiy);

% To return single feature vectors, return the first element.
Y = Y(:,1,:);
Y = squeeze(Y);

end

author: matlab-deep-learning

repo: transformer-models


text:
classdef DynamicPressureTermCondition < AbstractEventTerminationCondition
    %AltitudeTermCondition Summary of this class goes here
    %   Detailed explanation goes here

    properties
        dynP(1,1) double = 0; %km
        bodyInfo KSPTOT_BodyInfo
    end

    methods
        function obj = DynamicPressureTermCondition(dynP)
            obj.dynP = dynP;
        end

        function evtTermCondFcnHndl = getEventTermCondFuncHandle(obj)
            evtTermCondFcnHndl = @(t,y) obj.eventTermCond(t,y, obj.dynP, obj.bodyInfo);
        end

        function initTermCondition(obj, initialStateLogEntry)
            obj.bodyInfo = initialStateLogEntry.centralBody;
        end

        function name = getName(obj)
            name = sprintf('Dynamic Pressure (%.3f kPa)', obj.dynP);
        end

        function tf = shouldBeReinitOnRestart(obj)
            tf = true;
        end

        function params = getTermCondUiStruct(obj)
            params = struct();

            params.paramName = 'Dynamic Pressure';
            params.paramUnit = 'kPa';
            params.useParam = 'on';
            params.useStages = 'off';
            params.useTanks = 'off';
            params.useEngines = 'off';
            params.useStopwatches = 'off';

            params.value = obj.dynP;
            params.refStage = LaunchVehicleStage.empty(1,0);
            params.refTank = LaunchVehicleEngine.empty(1,0);
            params.refEngine = LaunchVehicleEngine.empty(1,0);
            params.refStopwatch = LaunchVehicleStopwatch.empty(1,0);
        end

        function optVar = getNewOptVar(obj)
            optVar = DynPressOptimizationVariable(obj);
        end

        function optVar = getExistingOptVar(obj)
            optVar = obj.optVar;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end
    end

    methods(Static)
        function termCond = getTermCondForParams(paramValue, stage, tank, engine, stopwatch)
            termCond = DynamicPressureTermCondition(paramValue);
        end
    end

    methods(Static, Access=private)
        function [value,isterminal,direction] = eventTermCond(ut,y, targetDynP, bodyInfo)
            rVectECI = y(1:3);
            vVectECI = y(4:6);

            altitude = norm(rVectECI) - bodyInfo.radius;

            if(altitude <= bodyInfo.atmohgt && altitude >= 0)
                [lat, long, ~, ~, ~, ~, ~, vVectECEF] = getLatLongAltFromInertialVect(ut, rVectECI, bodyInfo, vVectECI);
                density = getAtmoDensityAtAltitude(bodyInfo, altitude, lat, ut, long);
            elseif(altitude <= 0)
                density = 0;
                vVectECEF = [0;0;0];
            else
                density = 0;
                vVectECEF = [0;0;0];
            end

            vVectEcefMag = norm(vVectECEF);
            vVectEcefMagMS = vVectEcefMag*1000;

            dynP = density * (vVectEcefMagMS^2) / 2; %kg/m^3 * m^2 / s^2  = kg/(m*s^2)
            dynP_kPa = dynP/1000;

            value = dynP_kPa - targetDynP;
            isterminal = 1;
            direction = 0;
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
classdef ExtremumValueConstraint < AbstractConstraint
    %ExtremumValueConstraint Summary of this class goes here
    %   Detailed explanation goes here

    properties
        normFact = 1;
        event LaunchVehicleEvent
        eventNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;
        extremum LaunchVehicleExtrema

        lb(1,1) double = 0;
        ub(1,1) double = 0;

        evalType(1,1) ConstraintEvalTypeEnum = ConstraintEvalTypeEnum.FixedBounds;
        stateCompType(1,1) ConstraintStateComparisonTypeEnum = ConstraintStateComparisonTypeEnum.Equals;
        stateCompEvent LaunchVehicleEvent
        stateCompNode(1,1) ConstraintStateComparisonNodeEnum = ConstraintStateComparisonNodeEnum.FinalState;
    end

    methods
        function obj = ExtremumValueConstraint(event, lb, ub)
            obj.event = event;
            obj.lb = lb;
            obj.ub = ub;

            obj.id = rand();
        end

        function [lb, ub] = getBounds(obj)
            lb = obj.lb;
            ub = obj.ub;
        end

        function [c, ceq, value, lwrBnd, uprBnd, type, eventNum, valueStateComp] = evalConstraint(obj, stateLog, celBodyData)
            type = obj.getConstraintType();

            switch obj.eventNode
                case ConstraintStateComparisonNodeEnum.FinalState
                    stateLogEntry = stateLog.getLastStateLogForEvent(obj.event);

                case ConstraintStateComparisonNodeEnum.InitialState
                    stateLogEntry = stateLog.getFirstStateLogForEvent(obj.event);

                otherwise
                    error('Unknown event node.');
            end

            extremaStates = stateLogEntry.getAllExtremaStates();
            extremaState = extremaStates(extremaStates.extrema == obj.extremum);

            value = extremaState.value;

            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                switch obj.stateCompNode
                    case ConstraintStateComparisonNodeEnum.FinalState
                        stateLogEntryStateComp = stateLog.getLastStateLogForEvent(obj.stateCompEvent);

                    case ConstraintStateComparisonNodeEnum.InitialState
                        stateLogEntryStateComp = stateLog.getFirstStateLogForEvent(obj.stateCompEvent);

                    otherwise
                        error('Unknown event node.');
                end

                extremaStates = stateLogEntryStateComp.getAllExtremaStates();
                extremaState = extremaStates(extremaStates.extrema == obj.extremum);

                valueStateComp = extremaState.value;
            else
                valueStateComp = NaN;
            end

            [c, ceq] = obj.computeCAndCeqValues(value, valueStateComp);

            lwrBnd = obj.lb;
            uprBnd = obj.ub;

            eventNum = obj.event.getEventNum();
        end

        function sF = getScaleFactor(obj)
            sF = obj.normFact;
        end

        function setScaleFactor(obj, sF)
            obj.normFact = sF;
        end

        function tf = usesStage(obj, stage)
            tf = false;
        end

        function tf = usesEngine(obj, engine)
            tf = false;
        end

        function tf = usesTank(obj, tank)
            tf = false;
        end

        function tf = usesEngineToTankConn(obj, engineToTank)
            tf = false;
        end

        function tf = usesEvent(obj, event)
            tf = obj.event == event;
            if(obj.evalType == ConstraintEvalTypeEnum.StateComparison)
                tf = tf || obj.stateCompEvent == event;
            end
        end

        function tf = usesStopwatch(obj, stopwatch)
            tf = false;
        end

        function tf = usesExtremum(obj, extremum)
            tf = [obj.extremum] == extremum;
        end

        function tf = canUseSparseOutput(obj)
            tf = false;
        end

        function event = getConstraintEvent(obj)
            event = obj.event;
        end

        function type = getConstraintType(obj)
            type = 'Extremum Value';
        end

%         function name = getName(obj)
%             name = sprintf('%s - Event %i', obj.getConstraintType(), obj.event.getEventNum());
%         end

        function [unit, lbLim, ubLim, usesLbUb, usesCelBody, usesRefSc] = getConstraintStaticDetails(obj)
            unit = obj.extremum.unitStr;
            lbLim = -Inf;
            ubLim = Inf;
            usesLbUb = true;
            usesCelBody = false;
            usesRefSc = false;
        end

        function addConstraintTf = openEditConstraintUI(obj, lvdData)
%             [listBoxStr, extrema] = lvdData.launchVehicle.getExtremaListBoxStr();
%
%             if(isempty(extrema))
%                 addConstraintTf = false;
%
%                 warndlg('Cannot create extrema value object: no extrema have been created.  Create an extremum first.','Extremum Value Constraint','modal');
%             else
%                 [Selection,ok] = listdlg('PromptString',{'Select an extremum:'},...
%                                 'SelectionMode','single',...
%                                 'Name','Extremum',...
%                                 'ListString',listBoxStr);
%
%                 if(ok == 0)
%                     addConstraintTf = false;
%                 else
%                     ex = extrema(Selection);
%                     obj.extremum = ex;
%
%                     addConstraintTf = lvd_EditGenericMAConstraintGUI(obj, lvdData);
%                 end
%             end

            ex = obj.selectConstraintObj(lvdData);

            if(not(isempty(ex)))
                obj.extremum = ex;
%                 addConstraintTf = lvd_EditGenericMAConstraintGUI(obj, lvdData);

                output = AppDesignerGUIOutput({false});
                lvd_EditGenericMAConstraintGUI_App(obj, lvdData, output);
                addConstraintTf = output.output{1};
            else
                addConstraintTf = false;
            end
        end

        function ex = selectConstraintObj(obj, lvdData)
            [listBoxStr, extrema] = lvdData.launchVehicle.getExtremaListBoxStr();

            ex = [];
            if(isempty(extrema))
                warndlg('Cannot create extrema value object: no extrema have been created.  Create an extremum first.','Extremum Value Constraint','modal');
            else
                [Selection,ok] = listdlg('PromptString',{'Select an extremum:'},...
                                'SelectionMode','single',...
                                'Name','Extremum',...
                                'ListString',listBoxStr);

                if(ok == 0)
                    ex = [];
                else
                    ex = extrema(Selection);
                end
            end
        end

        function useObjFcn = setupForUseAsObjectiveFcn(obj,lvdData)
            ex = obj.selectConstraintObj(lvdData);

            if(not(isempty(ex)))
                obj.extremum = ex;
                useObjFcn = true;
            else
                useObjFcn = false;
            end
        end
    end

    methods(Static)
        function constraint = getDefaultConstraint(~, ~)
            constraint = ExtremumValueConstraint(LaunchVehicleEvent.empty(1,0),0,0);
        end
    end
end


text:
function y = acot(x)
%ACOT         Implements  acot(x)  for intervals
%
%   y = acot(x)
%
%interval standard function implementation
%

% written  10/16/98     S.M. Rump
% modified 06/24/99     S.M. Rump  complex allowed, sparse input,
%                                  major revision, improved accuracy
% modified 04/04/04     S.M. Rump  set round to nearest for safety
%                                    accelaration for sparse input
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/06/07     S.M. Rump  improved performance
% modified 10/20/08     S.M. Rump  check for zero
% modified 08/26/12     S.M. Rump  global variables removed
% modified 10/13/12     S.M. Rump  INTLAB_INTVAL_STDFCTS
%

  index = ( x==0 );
  y = atan(1./x);

  if ~isempty(find(index))                    % treat zero indices
    INTLAB_STDFCTS_PI = getappdata(0,'INTLAB_STDFCTS_PI');
    PI2 = intval(INTLAB_STDFCTS_PI.PI2INF,INTLAB_STDFCTS_PI.PI2SUP,'infsup');
    %VVVV  y(index) = PI2;
    s.type = '()'; s.subs = {index}; y = subsasgn(y,s,PI2);
    %AAAA  Matlab bug fix
  end

repo: OpenMAS


text:
function publish_scn_session_spike_id(inputimgs, SUBJDATA,outformat)
% This function is a wrapper function to call scn_session_spike_id in
% 'multi-session' mode, using input data across the runs for a single
% subject.  It runs the program, and generates both a yaml-format text file
% for uploading into the CANlab database, and an html file with all the
% results and images for that subject embedded.
%
% :Inputs:
%
%   **inputimgs:**
%        is a cell array of images (4-D) for each run in a separate cell.
%        the cell array outside this function has to be named 'inputimgs'
%
%   **SUBJDATA:**
%        Input fields of SUBJDATA define the experiment name, subject name,
%        and directories for saving both QC images + yaml and HTML
%        - study. a string
%        - subject. another string with subject ID
%        - html_save_dir can be one directory for all subjects if format 'pdf'
%        - subject_dir. a new subdirectory qc_images will be created within
%        subject_dir. images and .yaml files will be written there.
%
%
%   **ouformat**
%       file format for the report as available in 'publish'. can be
%       'html' [default],'pdf', 'latex', etc.
%
% :Examples:
% ::
%
%    SUBJDATA.study = 'NSF';
%    SUBJDATA.subject = subjects{i};
%    SUBJDATA.html_save_dir = fullfile(output_basedir, 'html_output');
%    SUBJDATA.subject_dir = fullfile(output_basedir, 'SubjectData', 'denoised_canlab', SUBJDATA.subject);
%
% ..
%    This has been tested with SPM8 and Matlab2010a only.
%
%    Tor Wager, Aug 2010
% ..

% Programmer's Notes
%
% 8/25/2017 - Stephan Geuter
% added option to use other fileformats like PDF
%

%% check input
% default report format is html. choose pdf or other options from publish.m
if nargin<3, outformat = 'html'; end

%% Initialize yaml file for database integration
if ~exist(SUBJDATA.subject_dir, 'dir'), mkdir(SUBJDATA.subject_dir); end
cd(SUBJDATA.subject_dir)

yamlfilename = fullfile(SUBJDATA.subject_dir, 'qc_results.yaml');

SUBJDATA.unique_id = [SUBJDATA.study '_' SUBJDATA.subject];

struct2yaml(yamlfilename, SUBJDATA, 'new', 'replace');

%% housekeeping: get rid of any old .png qc images
% these will mess it up, as it checks for existing files in order to
% require minimal input to scn_session_spike_id

qcdir = fullfile(SUBJDATA.subject_dir, 'qc_images');
if exist(qcdir, 'dir')
    disp('Removing old qc png images');

    disp(['!rm ' qcdir filesep '*png'])
    eval(['!rm ' qcdir filesep '*png'])
end

%% Create .m file script to publish that will create HTML file

htmlname = ['qc_scn_session_spike_id_' SUBJDATA.subject];

if ~exist(SUBJDATA.html_save_dir, 'dir'), mkdir(SUBJDATA.html_save_dir); end

diaryfilename = fullfile(SUBJDATA.html_save_dir, [htmlname '.m']);
if exist(diaryfilename, 'file'), eval(['!rm ' diaryfilename]); end

% These commands, written in the file called diaryfilename,
% will be evaluated in the m-file, and image and text
% results published in the HTML file called htmlname.
% ------------------------------------------------------------------

diary(diaryfilename);

disp(['%% SCN session spike ID, subject ' SUBJDATA.subject]);
fprintf('%% *Experiment:* %s\n%% \n', SUBJDATA.study); % start new formatted block if two spaces
fprintf('%% *QC results stored in:* %s\n%% \n', SUBJDATA.subject_dir);
fprintf('%% spike ID, quality control step 1.  writes to ''qc_results.yaml''\n');
%%% NEED COMMAND THAT LOADS INPUTIMGS FROM DISK -- OR PASSES IN AS INPUT?
disp('scn_session_spike_id(inputimgs);')
disp('snapnow');

diary off

%%
addpath(SUBJDATA.html_save_dir)
p = struct('useNewFigure', false, 'maxHeight', 800, 'maxWidth', 500, ...
    'outputDir', SUBJDATA.html_save_dir, 'showCode', false,'format',outformat);
publish(diaryfilename, p);

fprintf('\npublished to\n%s\n\n',diaryfilename);

if strcmp(outformat,'html')
    htmlfilename = fullfile(SUBJDATA.html_save_dir, [htmlname '.html']);
    web(htmlfilename);
end

% clean up the original diary; no longer needed
eval(['!rm ' diaryfilename]);

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [LeakMets,modelClosed,FluxExV] = fastLeakTest(model, testRxns,test)
% This function performs a leak test, i.e., it tests whether the model can produce mass from nothing (i.e., when all lower bound on boundary reactions are set to be 0 corresponding to no uptake).
%
% [LeakMets,modelClosed,FluxExV] = fastLeakTest(model, testRxns,test)
%
% INPUT
% model         Model structure
% testRxns      Reactiosn to be tested
% test          If empty (''), the input model structure will be tested as is.
%               If 'demand' is provided as input or no entry is provided, then demand
%               reactions for each internal metabolite will be added and the function
%               will test whether any of the internal metabolites can be produced from
%               nothing.
%
% OUTPUT
% LeakMets      List of metabolites that can be produced from nothing. If
%               empty the model is leak free.
% modelClosed   Model structure with closed boundary reactions (lower bound
%               = 0.
% FluxExV       List of boundary reactions that produce metabolites from
%               nothing.
%
% Ines Thiele Jan 2015

if nargin<3
    test = 'demands';
end
tol = 1e-06;
modelClosed = model;
% find all reactions that have only one entry in S
clear count ExR
count=false(size(modelClosed.S,2),1);
for i = 1 :size(modelClosed.S,2)
    if length(find(modelClosed.S(:,i))<0)==1
        count(i,1)=1;
        modelClosed.lb(i)=0;
        % modelClosed.ub(i)=0;
    elseif length(find(modelClosed.S(:,i))>0)==1
        count(i,1)=1;
        %   modelClosed.ub(i)=0;
        modelClosed.lb(i)=0;
    end
end
ExR = modelClosed.rxns(count);

modelexchangesAbbr = unique([testRxns;ExR]);
FluxEx = [];
cnt =1;
%% test for all demand reactions is an option
if strcmp(test,'demands')
    % add demand reactions for all metabolites in model to check for those too
    [modelClosed,rxnNames] = addDemandReaction(modelClosed,modelClosed.mets,0);
else
    rxnNames = '';
end
modelexchangesAbbr = unique([modelexchangesAbbr;rxnNames']);
TestRxnNum = length(modelexchangesAbbr)
FluxExV =[];
while cnt == 1
    modelClosed = changeObjective(modelClosed,modelexchangesAbbr);
    FF2=optimizeCbModel(modelClosed,'max');
    ObjValue = FF2.f
    if FF2.f >= tol
        FluxR = modelClosed.rxns(find(abs(FF2.x)>tol));
        FluxEx = [FluxEx;intersect(modelexchangesAbbr,FluxR)];
        FluxExV = [FluxExV;FF2.x(find(ismember( modelClosed.rxns,intersect(modelexchangesAbbr,FluxR))))];
        modelexchangesAbbr = setdiff(modelexchangesAbbr, FluxEx);
        length(unique(FluxEx))
    else
        cnt = 2;
    end
end

LeakMets = FluxEx;

author: opencobra

repo: cobratoolbox


text:
%CAMSHIFT  Finds an object center, size, and orientation
%
%     box = cv.CamShift(probImage, window)
%     [box,window] = cv.CamShift(probImage, window)
%     [...] = cv.CamShift(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __probImage__ Back projection of the object histogram. See
%   cv.calcBackProject for details.
% * __window__ Initial search window `[x,y,w,h]`.
%
% ## Output
% * __box__ Output rectangle with rotation. A scalar structure of the form:
%   `struct('center',[x,y], 'size',[w,h], 'angle',a)`
% * __window__ Converged CAMSHIFT window `[x,y,w,h]`
%
% ## Options
% * __Criteria__ Stop criteria for the underlying cv.meanShift. Accepts a
%   struct with 'type', 'maxCount', and 'epsilon' fields. Default
%   `struct('type','Count+EPS', 'maxCount',100, 'epsilon',1.0)`
%
% The function implements the CAMSHIFT object tracking algorithm
% [Bradski98]. First, it finds an object center using cv.meanShift and then
% adjusts the window size and finds the optimal rotation. The function
% returns the rotated rectangle structure that includes the object
% position, size, and orientation. The next position of the search window
% can be obtained with cv.RotatedRect.boundingRect.
%
% ## References
% [Bradski98]:
% > Gary R Bradski. Computer vision face tracking for use in a perceptual
% > user interface. 1998.
%
% See also: cv.meanShift, cv.calcBackProject, vision.HistogramBasedTracker
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function r = mldivide(a,b)
% o \ v
%
% Syntax
%   h = o \ r
%
% Input
%  o - @orientation
%  r - @vector3d
%
% Output
%  h - @Miller
%

r = inv(a) * b; %#ok<MINV>

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function outstring =  demorse(wavfile);

%demorse a morsed input wav file

vis_on = 0;

threshold = 0.05;

x = wavread(wavfile);

% half-wave rectify x
x2 = abs(x);

% slow-wave filter
y = filter(ones(1,20)/20,1, x2);

% threshold (digitize) y
z = y > threshold;
% z is now effectively our morse signal

if vis_on
   figure(1);
   subplot(3,1,1);
   plot(x, 'r');
   title('original signal');

   subplot(3,1,2);
   plot(y);
   title('HWR + Slow-wave filter -> envelope');
   subplot(3,1,3);
   plot(z, 'o', 'MarkerSize', 2);
   title('Digitized Morse signal')

end

%zero pad z so we always start with an onset
z = [zeros(10,1); z];

% id tones/spaces   -----------------------
% --> find changes between 0/1 and 1/0

b = diff(z);
% figure(3); plot(b, '.');
% 1: change from 1 to 0
% 0: no change
% -1: change from 0 to 1

c = b(b~=0);
c2 = find(b~=0);

tokens = -c .* diff([0; c2]);
% value == length of token
% sign == tone/space

% id shorts/longs  -----------------------

% since short/long should be bi-modal dist, a regular average should give
% us a good cutoff point to distinguish between the two? (assuming equal
% counts of short and long...)
% use mean as simple cutoff point; smarter algorithms can get smarter about
% this classification if they want to.

% 1: short, 2: long, +: tone, -: space
tokens2 = tokens;

% cutoff tones, cutoff spaces;
cut_t = mean(tokens2(tokens2>0));
cut_s = mean(tokens2(tokens2<0));

tokens2(tokens > 0 & tokens < cut_t) = 1;
tokens2(tokens > 0 & tokens > cut_t) = 2;
tokens2(tokens < 0 & tokens > cut_s) = -1;
tokens2(tokens < 0 & tokens < cut_s) = -2;

% now tokens 2 is a string of -1s, -2s, 1s, 2s, can trim first known space;
% put final endstop at end
tokens2 = [tokens2(2:end); -2];

% can drop little spaces, b/c they don't matter when parsing;
tokens2(tokens2 == -1) = [];
tokens3 = tokens2;
tokens4 = {};
ctr = 1;
start_idx = 1;

%parse
toparse = find(tokens3(start_idx:end) == -2);

for j=1:length(toparse)
   a = toparse(j);
   temp = tokens3(start_idx:a-1);
   tokens4{j} = temp;
   % zeropad for easy comparison
   %tokens4{j} = [tokens4{j}; zeros(length(tokens4{j}), 1)];
   start_idx = a+1;

end

% now tokens4 is de-codeable tokens... proceed to setup lookups
% letters
code{1} = [1 2 ];
code{2} = [2 1 1 1];
code{3} = [2 1 2 1];
code{4} = [2 1 1];
code{5} = [1];
code{6} = [1 1 2 1];
code{7} = [2 2 1];
code{8} = [1 1 1 1];
code{9} = [1 1];
code{10} = [1 2 2 2];
code{11} = [2 1 2];
code{12} = [1 2 1 1];
code{13} = [2 2];
code{14} = [2 1];
code{15} = [2 2 2];
code{16} = [1 2 2 1];
code{17} = [1 2 1 2];
code{18} = [1 2 1];
code{19} = [1 1 1];
code{20} = [2];
code{21} = [1 1 2];
code{22} = [1 1 1 2];
code{23} = [1 2 2];
code{24} = [2 1 1 2];
code{25} = [2 1 2 2];
code{26} = [2 2 1 1];

% punct
code{27} = [1 2 1 2 1 2];
code{28} = [2 2 1 1 2 2];
code{29} = [1 1 2 2 1 1];
code{30} = [2 1 1 2 1];

% numbers

code{31} = [1 2 2 2 2];
code{32} = [1 1 2 2 2];
code{33} = [1 1 1 2 2];
code{34} = [1 1 1 1 2];
code{35} = [1 1 1 1 1];
code{36} = [2 1 1 1 1];
code{37} = [2 2 1 1 1];
code{38} = [2 2 2 1 1];
code{39} = [2 2 2 2 1];
code{40} = [2 2 2 2 2];

decode{1} = 'A';
decode{2} = 'B';
decode{3} = 'C';
decode{4} = 'D';
decode{5} = 'E';
decode{6} = 'F';
decode{7} = 'G';
decode{8} = 'H';
decode{9} = 'I';
decode{10} = 'J';
decode{11} = 'K';
decode{12} = 'L';
decode{13} = 'M';
decode{14} = 'N';
decode{15} = 'O';
decode{16} = 'P';
decode{17} = 'Q';
decode{18} = 'R';
decode{19} = 'S';
decode{20} = 'T';
decode{21} = 'U';
decode{22} = 'V';
decode{23} = 'W';
decode{24} = 'X';
decode{25} = 'Y';
decode{26} = 'Z';
decode{27} = '.';
decode{28} = ',';
decode{29} = '?';
decode{30} = '/';
decode{31} = '1';
decode{32} = '2';
decode{33} = '3';
decode{34} = '4';
decode{35} = '5';
decode{36} = '6';
decode{37} = '7';
decode{38} = '8';
decode{39} = '9';
decode{40} = '0';

% compare tokens to tables

out1 = [];

for j = 1:length(tokens4)
    %zero pad temp_tok
    temp_tok = [tokens4{j}; zeros(6 - length(tokens4{j}), 1)];
    for k = 1:length(code)
        if (temp_tok == [code{k}'; zeros(6 - length(code{k}), 1)]);
            out1(j) = char(decode{k});
            %display(decode{k})
        end

    end

    % if didn't find a match
    if isempty(out1(j))
        out1(j) = '_';
    end

end

% semi-prettify
outstring = 32*ones(2*length(out1),1);
outstring(2:2:end) = out1;
outstring = char(outstring');

%display('demorsed message:')
%display(outstring);


text:
function [VM, h, fit, e, param] = Det_Logit(V0,t,tc,Run)
%
% [theta,HH,C,P] = Anneal_Logit(theta0,t,tc,Run)
%
% Estimate inverse logit (IL) HRF model
% Creates fitted curve - 3 logistic functions to be summed together - from parameter estimates
%
% INPUT: V0, t, tc, Run
% Run = stick function
% tc = time course
% t = vector of time points
% V0 = initial value for the parameter vector
%
% By Martin Lindquist and Tor Wager
% Edited 10/01/09
%

% Find optimal values
options = optimset('MaxFunEvals',10000000,'Maxiter',10000000,'TolX',1e-8,'TolFun',1e-8,'Display','off');
VM = fminsearch(@msq_logit,V0,options,Run,t,tc);
VM

% Use optimal values to fit hemodynamic response functions
h = il_hdmf_tw2(t,VM(1:7));

%[param] = get_parameters2(h,t);
[param] = get_parameters_logit(h,t,VM(1:7));

len = length(Run);
fit = conv(Run, h);
fit = fit(1:len);

e = tc-fit;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   SUBFUNCTIONS
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function m=msq_logit(V,Run, t, tc)

HR = il_hdmf_tw2(t,V(1:7));
len = length(Run);
timecourse = conv(Run, HR);
timecourse = timecourse(1:len);

m=sum((tc-timecourse).^2);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [h,base] = il_hdmf_tw2(t,V)
% inverse logit -- creates fitted curve from parameter estimates
%
% t = vector of time points
% V = parameters

% 3 logistic functions to be summed together
base = zeros(length(t),3);
A1 = V(1);
T1 = V(2);
d1 = V(3);
A2 = V(4);
T2 = V(5);
A3 = V(6);
T3 = V(7);
d2 = -d1*(ilogit(A1*(1-T1)) - ilogit(A3*(1-T3)))/(ilogit(A2*(1-T2)) + ilogit(A3*(1-T3)));
d3 = abs(d2)-abs(d1);

base(:,1)= d1*ilogit(A1*(t-T1))';
base(:,2)= d2*ilogit(A2*(t-T2))';
base(:,3)= d3*ilogit(A3*(t-T3))';
h = sum(base,2)';

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [L] = ilogit(t)
L = exp(t)./(1+exp(t));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% function [param] = get_parameters2(hdrf,t)
% % Find model parameters
% %
% % Height - h
% % Time to peak - p (in time units of TR seconds)
% % Width (at half peak) - w
%
% % Calculate Heights and Time to peak:
%
% n = ceil(t(end)*0.8);
% [h,p] = max(abs(hdrf(1:n)));
% h = hdrf(p);
%
% if (h >0)
%     v = (hdrf >= h/2);
% else
%     v = (hdrf <= h/2);
% end;
%
% [a,b] = min(diff(v));
% v(b+1:end) = 0;
% w = sum(v);
%
% cnt = p-1;
% g =hdrf(2:end) - hdrf(1:(end-1));
% while((cnt > 0) & (abs(g(cnt)) <0.001)),
%     h = hdrf(cnt);
%     p = cnt;
%     cnt = cnt-1;
% end;
%
% param = zeros(3,1);
% param(1) = h;
% param(2) = p;
% param(3) = w;
%
% end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


text:
function op = proj_l1( q, d )
%PROJ_L1   Projection onto the scaled 1-norm ball.
%    OP = PROJ_L1( Q ) returns an operator implementing the
%    indicator function for the 1-norm ball of radius q,
%    { X | norm( X, 1 ) <= q }. Q is optional; if omitted,
%    Q=1 is assumed. But if Q is supplied, it must be a positive
%    real scalar.
%
%    OP = PROJ_L1( Q, D ) uses a scaled 1-norm ball of radius q,
%    { X | norm( D.*X, 1 ) <= 1 }. D should be the same size as X
%    and non-negative (some zero entries are OK).
%
% Dual: prox_linf.m
% See also: prox_linf, prox_l1, proj_linf

if nargin == 0,
	q = 1;
elseif ~isnumeric( q ) || ~isreal( q ) || numel( q ) ~= 1 || q <= 0,
	error( 'Argument must be positive.' );
end
if nargin < 2 || isempty(d) || numel(d)==1
    if nargin>=2 && ~isempty(d)
        % d is a scalar, so norm( d*x ) <= q is same as norm(x)<=q/d
        if d==0
            error('If d==0 in proj_l1, the set is just {0}, so use proj_0');
        elseif d < 0
            error('Require d >= 0');
        end
        q = q/d;
    end
    op = @(varargin)proj_l1_q(q, varargin{:} );
else
    if any(d<0)
        error('All entries of d must be non-negative');
    end
    op = @(varargin)proj_l1_q_d(q, d, varargin{:} );
end

function [ v, x ] = proj_l1_q( q, x, t )
v = 0;
switch nargin,
case 2,
	if nargout == 2,
		error( 'This function is not differentiable.'  );
	elseif norm( x(:), 1 ) > q,
		v = Inf;
	end
case 3,
    s      = sort(abs(nonzeros(x)),'descend');
    cs     = cumsum(s);
    % ndx    = find( cs - (1:numel(s))' .* [ s(2:end) ; 0 ] >= q, 1 );
    ndx    = find( cs - (1:numel(s))' .* [ s(2:end) ; 0 ] >= q+2*eps(q), 1 ); % For stability
    if ~isempty( ndx )
        thresh = ( cs(ndx) - q ) / ndx;
        x      = x .* ( 1 - thresh ./ max( abs(x), thresh ) ); % May divide very small numbers
    end
otherwise,
    error( 'Not enough arguments.' );
end

% Allows scaling. Added Feb 21 2014
function [ v, x ] = proj_l1_q_d( q, d,  x, t )
v = 0;
switch nargin,
case 3,
	if nargout == 2,
		error( 'This function is not differentiable.'  );
	elseif norm( d(:).*x(:), 1 ) > q,
		v = Inf;
	end
case 4,
    [lambdas,srt]      = sort(abs(nonzeros(x./d)),'descend');
    s   = abs(x(:).*d(:));  s = s(srt);
    dd  = d(:).^2;          dd= dd(srt);
    cs  = cumsum(s);
    cd  = cumsum(dd);
    ndx    = find( cs - lambdas.*cd >= q+2*eps(q), 1, 'first');
    if ~isempty( ndx )
        ndx     = ndx - 1;
        lambda  = ( cs(ndx) - q )/cd(ndx);
        x       = sign(x).*max( 0, abs(x) - lambda*d );
    end
otherwise,
    error( 'Not enough arguments.' );
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS