text:
function [u,sigma,eqn,info] = PoissonRT0(node,elem,pde,bdFlag,option)
%% POISSONRT0 Poisson equation: lowest order RT element.
%
%  [u,sigma] = PoissonRT0(node,elem,pde,bdFlag) produces an approximation of
%  the Poisson equation
%
%       -div(d*grad(u))=f  in \Omega, with
%       Dirichlet boundary condition u=g_D on \Gamma_D,
%       Neumann boundary condition   d*grad(u)*n=g_N on \Gamma_N
%
%  in the mixed formulation:
%
%  ------------------------------------------------------------------------
%  Find (\sigma , u) in H_{g_N,\Gamma_N}(div,\Omega)\times L^2(\Omega) s.t.
%
%  (d^-1\sigma,\tau) - (div \tau, u)  = <\tau*n,g_D>_{\Gamma_D}
%  \forall \tau in H_{0,\Gamma_N}(div,\Omega)
%    - (div \sigma, v)                =  -(f,v)
%  \forall v in L^2(\Omega)
%
%  where
%  H_{g,\Gamma}(div,\Omega) = {\sigma \in H(div,\Omega); \sigma*n = g
%  on \Gamma \subset \partial\Omega }.
%  ------------------------------------------------------------------------
%
%  The unknown sigma = d*grad(u) is approximated using the lowest order
%  Raviart-Thomas element and u by piecewise constant element (with basis 1).
%
%  [u,sigma] = PoissonRT0(node,elem,pde,bdFlag,option) specifies options
%   - option.solver
%     'direct': the built in direct solver \ (mldivide)
%     'dmg':     multigrid-type solvers mg is used.
%     'uzawapcg': PCG for the Schur complement equation
%     'none': only assemble the matrix equation but not solve
%
%   The default setting is to use the direct solver for small size problems
%   and transforming based multigrid solvers for large size problems.
%
%  Example
%
%    exampleRT0
%
% Created by Ming Wang. Reorganized by Long Chen. Change basis for u.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if ~exist('option','var'), option = []; end

%% Diffusion coefficient
if ~isfield(pde,'d'), pde.d = []; end
if isfield(pde,'d') && ~isempty(pde.d)
   if isnumeric(pde.d)
      K = pde.d;                   % d is an array
   else                            % d is a function
      center = (node(elem(:,1),:) + node(elem(:,2),:) + node(elem(:,3),:))/3;
      K = pde.d(center);  % take inverse sequencil.
   end
else
    K = [];
end

%% Data structure
elemold = elem;
[elem,bdFlag] = sortelem(elem,bdFlag);  % ascend ordering
[elem2edge,edge] = dofedge(elem);
NT = size(elem,1); NE = size(edge,1);
[Dlambda,area,elemSign] = gradbasis(node,elem);

%% Assemble matrix
Nsigma = NE; Nu = NT; Ndof = Nsigma + Nu;

% M. Mass matrix for RT0 element
M = getmassmatvec(elem2edge,area,Dlambda,'RT0',K);

% B. negative divergence operator
B = icdmat(double(elem2edge),elemSign*[1 -1 1]);

% C. zero matrix.
C = sparse(Nu,Nu);

A = [M B';B C];

%% Assemble right hand side.
fu = zeros(Nu,1);
if ~isfield(pde,'f') || (isreal(pde.f) && (pde.f==0))
    pde.f = [];
end
if ~isfield(option,'fquadorder')
    option.fquadorder = 2;   % default order
end
if ~isempty(pde.f)
	[lambda,weight] = quadpts(option.fquadorder);
	nQuad = size(lambda,1);
    for p = 1:nQuad
		% quadrature points in the x-y coordinate
		pxy = lambda(p,1)*node(elem(:,1),:) ...
			+ lambda(p,2)*node(elem(:,2),:) ...
			+ lambda(p,3)*node(elem(:,3),:);
		fp = pde.f(pxy);
		fu = fu - fp*weight(p);
    end
    fu = fu.*area;
end
clear fp area
F((Nsigma+1):(Ndof),1) = fu;

%% Boundary Conditions
if ~exist('bdFlag','var'), bdFlag = []; end
[AD,F,bigu,freeDof,isPureNeumannBC] = getbdRT0(F);
eqn = struct('M',AD(1:NE,1:NE),'B',AD(NE+1:end,1:NE),'C',AD(NE+1:end,NE+1:end),...
             'f',F(1:NE),'g',F(NE+1:end),'freeDof',freeDof);

%% Solve the linear system.
% Set up solver type
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if Ndof <= 2e3  % Direct solver for small size systems
        option.solver = 'direct';
    else         % MGCG  solver for large size systems
        option.solver = 'dmg';
    end
elseif strcmp(option.solver,'mg')
    option.solver = 'dmg';
end
solver = option.solver;
% solve
switch lower(solver);
    case 'direct'
      bigu(freeDof) = AD(freeDof,freeDof)\F(freeDof);
        sigma = bigu(1:NE);
        u = bigu(NE+1:end); info =[];
    case 'none'
        sigma = []; u = []; info =[];
    case 'dmg'
        [sigma,u] = dmg(eqn.M,eqn.B,eqn.C,eqn.f,eqn.g,elemold);
    case 'uzawapcg'
        [sigma,u] = uzawapcg(eqn.M,eqn.B,eqn.C,eqn.f,eqn.g,elemold);
end
info.solverTime = 0; info.assembleTime = 0; info.itStep = 0;
info.stopErr = 0; info.flag = 0;
if isPureNeumannBC == true % post process for u.
    ubar = sum(u.*area)/sum(area);
    u = u - ubar;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function [AD,F,bigu,freeDof,isPureNeumannBC] = getbdRT0(F)
    %% GETBDRT0 Boundary conditions for Poisson equation: RT0 element.
    %
    %  Created by Ming Wang. Improved the check of edgeSign by Long Chen.

    %%
    bigu = zeros(Ndof,1);

    %% Boundary conditions
    if ~isfield(pde,'g_D'), pde.g_D = []; end
    if ~isfield(pde,'g_N'), pde.g_N = []; end

    %% Set up bdFlag
    if isempty(bdFlag) % no bdFlag information
       if ~isempty(pde.g_N) % case: Neumann
           bdFlag = setboundary(node,elem,'Neumann');
       elseif ~isempty(pde.g_D) % case: Dirichlet
           bdFlag = setboundary(node,elem,'Dirichlet');
       end
    end

    %% Find Dirichlet and Neumann dofs
    if ~isempty(bdFlag)
        isDirichlet(elem2edge(bdFlag(:)==1)) = true;
        isNeumann(elem2edge(bdFlag(:)==2)) = true;
        % Direction of boundary edges may not be the outwards normal
        % direction of the domain. edgeSign is introduced to record this
        % inconsistency.
        edgeSign = ones(NE,1);
        idx = (bdFlag(:,1) ~= 0) & (elemSign == -1);% first edge is on boundary
        edgeSign(elem2edge(idx,1)) = -1;
        idx = (bdFlag(:,2) ~= 0) & (elemSign == 1); % second edge is on boundary
        edgeSign(elem2edge(idx,2)) = -1;
        idx = (bdFlag(:,3) ~= 0) & (elemSign == -1);% first edge is on boundary
        edgeSign(elem2edge(idx,3)) = -1;
    end
    Dirichlet = edge(isDirichlet,:);
    Neumann = edge(isNeumann,:);
    isBdDof = false(Ndof,1);
    isBdDof(isNeumann) = true;   % for mixed method, Neumann edges are fixed
    freeDof = find(~isBdDof);

    %% Dirichlet boundary condition (Neumann BC in mixed form)
    %   We need only modify the rhs on dof associated with Dirichlet
    %   boundary. Compute the int_e \Phi\cdot n g_D on the boundary using
    %   quadrature rules.
    if ~isempty(pde.g_D) && isnumeric(pde.g_D) && (pde.g_D==0)
        pde.g_D = [];
    end
    if ~isempty(pde.g_D) && any(isDirichlet)
        if ~isfield(option,'gNquadorder')
            option.gNquadorder = 2;   % default order exact for linear gN
        end
        [lambda,weight] = quadpts1(option.gNquadorder);
        nQuad = size(lambda,1);
        for ip = 1:nQuad
        	pxy = lambda(ip,1)*node(Dirichlet(:,1),:)+...
                  lambda(ip,2)*node(Dirichlet(:,2),:);
            F(isDirichlet) = F(isDirichlet) + weight(ip)*pde.g_D(pxy);
        end
        F(isDirichlet) = F(isDirichlet).*edgeSign(isDirichlet);
        % no edge length since the basis of sigma contains it.
    end

    %% Neumann boundary condition (Dirichlet BC in mixed form)
    if ~isempty(pde.g_N) && any(isNeumann)
        % modify the rhs to include Dirichlet boundary condition
        mid = 1/2*(node(Neumann(:,1),:)+node(Neumann(:,2),:));
        ve = node(Neumann(:,1),:)-node(Neumann(:,2),:);
        edgeLength = sqrt(sum(ve.^2,2));
        if isnumeric(pde.g_N)
            evalg_N = pde.g_N;
        else
            evalg_N = pde.g_N(mid);
        end
        bigu(isNeumann) = edgeLength.*evalg_N;
        if ~isempty(pde.d)
            bigu(isNeumann) = pde.d(mid).*bigu(isNeumann);
        end
        bigu(isNeumann) = bigu(isNeumann).*edgeSign(isNeumann);
        F = F - A*bigu;
        F(isNeumann) = bigu(isNeumann);
    end

    %% Pure Neumann boundary condition
    isPureNeumannBC = false;
    if ~any(isDirichlet) && any(isNeumann)
        freeDof = freeDof(1:end-1);  % eliminate the kernel by enforcing u(NT) = 0;
        isBdDof(end) = true;
        isPureNeumannBC = true;
%         F(end) = 0;
    end

    %% Modify the matrix
    %  Build Neumann boundary condition(Dirichlet BC in mixed form) into the
    %  matrix AD by enforcing  |AD(bdNode,bdNode)=I,
    %  AD(bdNode,FreeNode)=0, AD(FreeNode,bdNode)=0|.
    if any(isBdDof)
       bdidx = zeros(Ndof,1);
       bdidx(isBdDof) = 1;
       Tbd = spdiags(bdidx,0,Ndof,Ndof);
       T = spdiags(1-bdidx,0,Ndof,Ndof);
       AD = T*A*T + Tbd;
    else
       AD = A;
    end
    end
end

repo: Aero-matlab


text:
function val=intCosPow(u,n)
%%INTCOSPOWER Evaluate the integral of cos(u)^n du.  A definite integral
%             can be evaluated, or an indefinite integral (with a
%             particular additive constant).
%
%INPUTS: u A 2XN (for definite integral) or a 1XN (for indefinite
%          integrals) set of N points. For definite integrals, u(1,:) are
%          the real lower bounds and u(2,:) are the real upper bounds. For
%          indefinite integrals, the integral is evaluated at the points in
%          u. The values in u should be real.
%        n The positive integer exponent of the cosine term.
%
%OUTPUTS: val The 1XN set of values of the integral of cos(u)^n.
%
%This function simply implements the recursion that arises from integration
%by parts from basic calculus, as are given in the tables in the back of
%[1].
%
%REFERENCES:
%[1] J. Stewart, Calculus, 7th ed. Belmont, CA: Brooks/Cole, 2012.
%
%October 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numDim=size(u,1);

if(isempty(u))
   val=[];
   return;
end

if(numDim==1)%An indefinite integral
    val=indefIntCosPow(u,n);
else%A definite integral
    val=indefIntCosPow(u(2,:),n)-indefIntCosPow(u(1,:),n);
end

end

function val=indefIntCosPow(u,n)

%Special cases for n=0 and n=1.
if(n==0)
    val=u;
    return;
elseif(n==1)
    val=sin(u);
    return;
end

%If here, n>=1
sinVal=sin(u);
cosVal=cos(u);

val=0;
if(mod(n,2)==0)%If n is even
    endVal=2;
else
    endVal=3;
end

coeffProd=1;
for k=n:-2:endVal
    val=val+coeffProd*(1/k)*cosVal.^(k-1).*sinVal;
    coeffProd=coeffProd*(k-1)/k;
end

if(endVal==2)
    %The final integral is over cos^0
    val=val+coeffProd*u;
else
    %The final integral is over cos^1
    val=val+coeffProd*sinVal;
end

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
classdef BinomialD
%%BINOMIALD Functions to handle the binomial distribution. The binomial
%           distribution is the distribution of the number of successes in
%           a sequence of n independent trials, each with a probability of
%           success p.
%Implemented methods are: mean, var, PMF, CDF, rand
%
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

methods(Static)

function val=mean(n,p)
%%MEAN Obtain the mean of the binomial distribution for the given number of
%      trials and probability of success.
%
%INPUTS: n The number of trials performed.
%        p The probability of success for the underlying Bernoulli trials.
%          0<=p<=1.
%
%OUTPUTS: val The mean number of successes for the binomial distribution.
%
%EXAMPLE:
%We verify the computed mean by comparing it to the sample mean.
% p=0.75;
% n=61;
% meanVal=BinomialD.mean(n,p)
% numSamp=1e5;
% sampMeanVal=mean(BinomialD.rand([1,numSamp],n,p))
%One will find both values are about 45.75.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=n*p;
end

function val=var(n,p)
%%VAR Obtain the variance of the binomial distribution for the given number
%     of trials and probability of success.
%
%INPUTS: n The number of trials performed.
%        p The probability of success for the underlying Bernoulli trials.
%          0<=p<=1.
%
%OUTPUTS: val The variance of the binomial distribution with the given
%             parameters.
%
%EXAMPLE:
%We verify that the computed variance by comparing it to the sample
%variance.
% p=0.75;
% n=61;
% varVal=BinomialD.var(n,p)
% numSamp=1e5;
% sampVarVal=var(BinomialD.rand([1,numSamp],n,p))
%One will find both values are about 11.4375.
%
%October 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=n*p*(1-p);
end

function val=PMF(k,n,p)
%%PMF Evaluate the binomial probability mass function at a given point with
%     a specified number of trials and success probability per trial.
%
%INPUTS: k The nonnegative integer point(s) at which the binomial PMF is to
%          be evaluated (the number of successes in n trials). Note that
%          0<=k<=n.
%        n The number of trials of which k successes are observed.
%        p The probability of success for the underlying Bernoulli trials.
%          0<=p<=1.
%
%OUTPUT: val The value(s) of the binomial PMF with success probability p.
%
%The binomial PMF provides the probability of k successes from n trials,
%each with a success probability of p.
%
%EXAMPLE:
%Here, we validate the PMF by generating random samples and comparing the
%PMF plot with a histogram of the random samples.
% p=0.75;
% n=61;
% numSamples=1e5;
% figure(1)
% clf
% histogram(BinomialD.rand([numSamples,1],n,p),'Normalization','pdf')
% hold on
% x=0:61;
% vals=BinomialD.PMF(x,n,p);
% stem(x,vals,'linewidth',2)
% h1=xlabel('x');
% h2=ylabel('PDF(x)');
% set(gca,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h1,'FontSize',14,'FontWeight','bold','FontName','Times')
% set(h2,'FontSize',14,'FontWeight','bold','FontName','Times')
%One will see that the histogram matches well with the plot.
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    numPoints=length(k);
    val=zeros(size(k));
    for curK=1:numPoints
        if(k(curK)>=0&&k(curK)<=n)
            val(curK)=binomial(n,k(curK))*p^k(curK)*(1-p)^(n-k(curK));
        end
    end
end

function val=CDF(k,n,p)
%%CDF Evaluate the cumulative distribution function (CDF) of the binomial
%     distribution at a desired point.
%
%INPUTS: k The nonnegative integer point(s) at which the binomial CDF is to
%          be evaluated. Note that 0<=k<=n.
%        n The number of trials of which k or fewer successes are observed.
%        p The probability of success for the underlying Bernoulli trials.
%          0<=p<=1.
%
%OUTPUTS: val The CDF of the binomial distribution evaluated at the desired
%             point(s).
%
%Rather than summing over the values of the PMF, the equivalency between
%the binomial distribution and the regularized incomplete beta function,
%described at [1] is used. Thus, the built-in function betainc can be used
%for efficient evaluation of the CDF with large values.
%
%EXAMPLE:
%We validate the CDF value by comparing it to a value computed from random
%samples.
% p=0.75;
% n=61;
% x=45;
% numSamples=1e5;
% prob=BinomialD.CDF(x,n,p)
% probSamp=mean(BinomialD.rand([numSamples,1],n,p)<=x)
%One will find the values ot both be about 0.460.
%
%REFERENCES:
%[1] Weisstein, Eric W. "Binomial Distribution." From MathWorld--A Wolfram
%    Web Resource. http://mathworld.wolfram.com/BinomialDistribution.html
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    val=zeros(size(k));
    sel=k>=0&k<=n;

    val(sel)=betainc(1-p,n-k(sel),k(sel)+1);
end

function k=invCDF(probVals,n,p,choice)
%%INVCDF Evaluate the inverse sumulative distribution function (CDF) of the
%        binomial distirbution for given probabilities.
%
%INPUTS: probVals The CDF probability value(s) at which the argument of the
%          binomial CDF is to be determined. Note that 0<=prob<=1.
%        n The number of trials of which k or fewer successes are observed.
%        p The probability of success for the underlying Bernoulli trials.
%          0<=p<=1.
%   choice An optional parameter indicating what to do when probVals does
%          not exactly equal a CDF value (which will presumably be the case
%          most of the time). Possible values are:
%          0 Return the k that is the closest value.
%          1 (The default if omitted or an empty matrix is passed) return
%            the next lower k if it exists, otherwise return 0.
%          2 Return the next higher k if there is one, otherwise return n.
%
%OUTPUTS: k The set of non-negative integer points that are the inverse CDF
%           values (number of trials) corresponding to the values in
%           probVals.  This has the same dimensions as probVals.
%
%Though one can sum the PMF values until the resulting CDF values exceed
%prob, this can be slow as n gets large and the recursive computation of
%costs can be subject to finite precision limitations. These issues can be
%avoided by expressing the CDF in terms of the regularized incomplete beta
%function, which is given in [1]. The CDF can be evaluated using betainc
%without summing all of the terms. The inverse is found using the binSearch
%function.
%
%REFERENCES:
%[1] Weisstein, Eric W. "Binomial Distribution." From MathWorld--A Wolfram
%    Web Resource. http://mathworld.wolfram.com/BinomialDistribution.html
%
%September 2013 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(nargin<4||isempty(choice))
        choice=1;
    end

   f=@(idx)BinomialD.CDF(idx,n,p);
   k=zeros(size(probVals));
   numPts=numel(probVals);

    for i=1:numPts
        [~,k(i)]=binSearch(f,probVals(i),choice,[0;n]);
    end
end

function vals=rand(N,n,p)
%%RAND Generate binomial random variables with a given number of trials and
%      probability of success.
%
%INPUTS: N If N is a scalar, then rand returns an NXN matrix of random
%          variables. If N=[M,N1] is a two-element row  vector, then rand
%          returns an MXN1 matrix of random variables.
%        n The number of trials performed.
%        p The probability of success for the underlying Bernoulli trials.
%          0<=p<=1.
%
%The algorithm used depends on the value of n*min(p,1-p). If
%n*min(p,1-p)>=10, then the binomial, triangle, parallelogram, exponential
%algorithm of [1] is used. If n*min(p,1-p)<10, then a basic inverse
%transformation method, as described in [1] is used. This assures that the
%random number generation is fast for a wide variety of values of n and p.
%
%REFERENCES:
%[1] V. Kachitvichyanukul and B. W. Schmeiser, "Binomial random variate
%    generation," Communications of the ACM, vol. 31, no. 2, pp. 216-222,
%    Feb. 1988.
%
%August 2015 David F. Crouse, Naval Research Laboratory, Washington D.C.

    if(isscalar(N))
        dims=[N N];
    else
        dims=N;
    end

    vals=zeros(dims);
    numVals=numel(vals);

    %Do not use the BTPE algorithm if n*min(p,1-p) is small. Rather, the
    %inverse transformation algorithm from Section 2 is used. However, it
    %has been modified replacing p with pMin=min(p,1-p) as described in the
    %text to reduce the execution time.
    if(n*min(p,1-p)<10)
        pMin=min(p,1-p);
        q=1-pMin;
        s=pMin/q;
        a=(n+1)*s;

        for curVal=1:numVals
            u=rand(1);
            r=q^n;
            x=0;
            gotoNextVal=false;

            %With n*pMin<10, the mean of the PMF is below 10 and the
            %standard deviation is below sqrt(10). To assure a limited
            %runtime, we thus cap the maximum value of the iterated x at
            %1000, or about 316 standard deviations from the mean.
            while(x<1000)
                if(u<=r)
                    %This part deals with using pMin instead of p.
                    if(p>1/2)
                        x=n-x;
                    end
                    vals(curVal)=x;
                    gotoNextVal=true;
                    break;
                end
                u=u-r;
                x=x+1;
                r=((a/x)-s)*r;
            end

            if(gotoNextVal)
                continue;
            end

            if(p>1/2)
                x=n-x;
            end
            vals(curVal)=x;
        end
        return;
    end

%If we are here, then the BTPE algorithm should be used. Everything below
%is the BTPE algorithm from Section 3 of [1].

%Step 0: Set up constants as functions of n and p.
    r=min(p,1-p);
    q=1-r;
    fM=n*r+r;
    M=floor(fM);
    p1=floor(2.195*sqrt(n*r*q)-4.6*q)+0.5;
    xM=M+0.5;
    xL=xM-p1;
    xR=xM+p1;
    c=0.134+20.5/(15.3+M);
    a=(fM-xL)/(fM-xL*r);
    lambdaL=a*(1+a/2);
    a=(xR-fM)/(xR*q);
    lambdaR=a*(1+a/2);
    p2=p1*(1+2*c);
    p3=p2+c/lambdaL;
    p4=p3+c/lambdaR;

    for curVal=1:numVals
        gotoNextVal=false;
        while(1)
            if(gotoNextVal)
                break;
            end

    %Step 1: Generate a random variable for selecting the region. If region
    %   1 is selected, generate a triangularly distributed random variable.
            u=rand(1)*p4;
            v=rand(1);
            if(u<p1)%Generate a triangular variate.
                y=floor(xM-p1*v+u);
                if(p>1/2)%Step 6
                    y=n-y;
                end
                vals(curVal)=y;
                gotoNextVal=true;
                continue;
            elseif(u<=p2)%Step 2: Parallelograms.
                x=xL+(u-p1)/c;
                v=v*c+1-abs(M-x+0.5)/p1;
                if(v<=0||v>1)
                    continue;
                end

                y=floor(x);
            elseif(u<=p3)%Step 3: Left exponential tail.
                y=floor(xL+log(v)/lambdaL);
                if (y<0)
                    continue;
                end
                v=v*(u-p2)*lambdaL;
            else%Step 4: Right exponential tail.
                y=floor(xR-log(v)/lambdaR);
                if(y>n)
                    continue;
                end
                v=v*(u-p3)*lambdaR;
            end

    %Step 5: Acceptace/ Rejection Comparison.

            k=abs(y-M);
        %5.0 Test for appropriate method of evaluating f(y)
            if(k<=20||k>=n*r*q/2-1)
                %5.1 Evaluate f(y) via recursive relationship
                s=r/q;
                a=s*(n+1);
                F=1;
                if(M<y)
                    i=M;
                    while(1)
                        i=i+1;
                        F=F*(a/i-s);
                        if(i==y)
                            break;
                        end
                    end
                elseif(M>y)
                    i=y;
                    while(1)
                       i=i+1;
                       F=F/(a/i-s);
                       if(i==M)
                           break;
                       end
                    end
                end

                if(v<=F)
                    if(p>1/2)%Step 6
                        y=n-y;
                    end
                    vals(curVal)=y;
                    gotoNextVal=true;
                    continue;
                else
                    continue;
                end
            else
                %5.2 Test for squeezing. Check the value of ln(v) against
                %upper and lower bounds of ln(f(y))
                rho=(k/(n*r*q))*((k*(k/3+0.625)+(1/6))/(n*r*q)+0.5);
                t=-k^2/(2.0*n*r*q);
                A=log(v);
                if(A<t-rho)
                    if(p>1/2)%Step 6
                        y=n-y;
                    end
                    vals(curVal)=y;
                    gotoNextVal=true;
                    continue;
                elseif(A>t+rho)
                    continue;
                end
            %5.3 Final Acceptance/ Rejection Test
                x1=y+1;
                f1=M+1;
                z=n+1-M;
                w=n-y+1;
                x2=x1^2;
                f2=f1^2;
                z2=z^2;
                w2=w^2;

                temp=xM*log(f1/x1)+(n-M+0.5)*log(z/w) ...
                +(y-M)*log(w*r/(x1*q)) ...
                +(13860-(462-(132-(99 ...
                -140/f2)/f2)/f2)/f2)/f1/166320 ...
                +(13860-(462-(132-(99 ...
                -140/z2)/z2)/z2)/z2)/z/166320 ...
                +(13860-(462-(132-(99 ...
                -140/x2)/x2)/x2)/x2)/x1/166320 ...
                +(13860-(462-(132-(99 ...
                -140/w2)/w2)/w2)/w2)/w/166320;
                if(A<=temp)
                    if(p>1/2)%Step 6
                        y=n-y;
                    end
                    vals(curVal)=y;
                    gotoNextVal=true;
                    continue;
                else
                    continue;
                end
            end
        end
    end
end

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.


text:
function [dLdp] = spm_mci_adjoint_sun (Pr,M,U,Y)
% Gradient of log joint from adjoint method (via Sundials)
% FORMAT [dLdp] = spm_mci_adjoint_sun (Pr,M,U,Y)
%
% Pr        Parameters (vectorised and in M.V subspace)
% M         Model structure
% U         Inputs  [Nin x N]
% Y         Data
%
% dLdp      Gradient [Np x 1]
%
% For M.adjlike=1, dLdp is gradient of log likelihood not log joint
% (useful for debugging).
%
% For M.backint=1 (default), compute the integral underlying dLdp
% *during* backwards integration of adjoint. For M.backint=0, this
% integral is computed *after* adjoint (useful for debugging).
%
% B. Sengupta, K. Friston and W. Penny (2014) Efficient Gradient
% Computation for Dynamical Models. Neuroimage,98, 521-527.
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny and Biswa Sengupta
% $Id: spm_mci_adjoint_sun.m 6697 2016-01-27 14:57:28Z spm $

try, backint=M.backint; catch, backint=1; end
try, adjlike=M.adjlike; catch, adjlike=0; end

% Tolerances
tol_scale = 1e-3;
reltol = 1e-3;
abstol = 1e-4;

% Parameters in original space
P = M.V*Pr+M.vpE;

%tDur=[0 M.t(end)];
%tDur=[M.t(1) M.t(end)];

t0=0;
%tf=M.T;
%dt=M.t(2)-M.t(1);
%t0=M.t(1)-dt;
tf=M.T;

if isempty(U)
    U=zeros(1,M.N);
end

data.U=U;
data.P=P;
data.M=M;
data.Y=Y;

options = CVodeSetOptions('UserData', data,...
    'RelTol',reltol, ...
    'AbsTol',abstol, ...
    'LinearSolver','Dense', ...
    'JacobianFn',@djacfn);

CVodeInit(@spm_mci_flow_sun, 'BDF', 'Newton', t0, M.x0, options);

CVodeAdjInit(150, 'Hermite');

optionsB = CVodeSetOptions('UserData',data,...
    'MaxNumSteps',50000, ...
    'RelTol',reltol,...
    'AbsTol',abstol,...
    'LinearSolver','Dense',...
    'JacobianFn',@djacBfn);

lambda_init=zeros(M.n,1);
iB = CVodeInitB(@rhsadjoint, 'BDF', ...
    'Newton', tf, lambda_init, optionsB);

if backint
    % Use CVODE to compute the integral underlying dLdp
    % *during* backwards integration of adjoint (backint=1)

    [status,t,y] = CVode(tf,'Normal');

    optionsQB = CVodeQuadSetOptions('ErrControl',true,...
        'RelTol',reltol,...
        'AbsTol',abstol);

    qB1 = zeros(M.Np,1);
    CVodeQuadInitB(iB, @paramgrad, qB1, optionsQB);

    % Backward integration of the adjoint equation
    % (yB is the adjoint vector)
    [status,t,yB,qB] = CVodeB(M.t(2),'Normal');
    %[status,t,yB,qB] = CVodeB(t0,'Normal');
    if status == -1
        error('Adjoint integration failed');
    end
    qB = -qB';
    dt = M.t(2)-M.t(1);
    dLdp = qB/dt;
else
    % Compute integral underlying dLdp *after* backwards
    % integration of adjoint

    ntout = 1000;
    dt = (tf-t0)/ntout;
    tt = linspace(t0+dt,tf,ntout-1);

    %[status,ttf,x] = CVode(M.t(2:end),'Normal');
    %[status,ttf,x] = CVode(tt,'Normal');
    %[status,ttf,x] = CVode(M.t(2:end),'Normal');
    [status,ttf,x] = CVode(M.t,'Normal');
    interp_x = (interp1q(ttf',x',M.t))';

    %tm(1) = tDur(2);
    %t = tDur(2);
    tm(1) = tf;
    t = tf;
    it = 1;
    % Integrate adjoint equation
    while t > M.t(2)
    %while t > t0+2*dt
        it = it+1;
        % The adjoint vector is yB
        %[status,t,yB] = CVodeB(tDur(1),'OneStep');
        [status,t,yB] = CVodeB(M.t(1),'OneStep');
        %[status,t,yB] = CVodeB(t0+2*dt,'OneStep');
        if status == -1
            error('Adjoint integration failed');
        end
        tm(it) = t;
        dldp(it,:) = feval(@paramgrad,t, interp_x, yB, data);
    end
    dldp_int=interp1(tm,dldp,M.t,'spline');
    dLdp=sum(dldp_int);
end
CVodeFree;

if ~adjlike
    dlogpriordp = spm_mci_gprior_deriv (Pr,M);
    dLdp=dLdp+dlogpriordp;
end

end

% -----------------------------------------------------------
function [qBd, flag, new_data] = paramgrad(t, x, yB, data)
% Np quadratures for parametric gradient
% t     time
% x     state
% yB    adjoint vector, lambda^T
% data  contains P,M,U,Y
%
% yBd   dlambda^T/dt

P = data.P;
M = data.M;
U = data.U;

% If observation function becomes dependent on parameters
% we'll need to update term1
term1 = 0;

% Find nearest time point for which we have pre-computed input
if isempty(U)
    ut=[];
else
    [tmp,ind]=min(abs(t-M.t));
    ut=U(:,ind);
end

% Evaluate parameter Jacobian
if isfield(M,'dfdp')
    dfdp = feval(M.dfdp,x,ut,P,M);
else
    dfdp = spm_diff(M.f,x,ut,P,M,3);
end
term2 = -yB'*dfdp;
qBd = term1 + term2;

flag = 0;
new_data = [];

end

% -----------------------------------------------------------
function [yBd, flag, new_data] = rhsadjoint(t, x, yB, data)
% Adjoint equation
% FORMAT [yBd, flag, new_data] = rhsadjoint(t, x, yB, data)
%
% t     time
% x     state
% yB    adjoint vector, lambda^T
% data  contains P,M,U,Y
%
% yBd   dlambda^T/dt

Y = data.Y;
M = data.M;
U = data.U;
P = data.P;

% Interpolate data to required time point
for j=1:M.l,
    ydata(j)=interp1q(M.t,Y(:,j),t);
end
[y,L]=feval(M.g,x,U(:,1),P,M);
e=ydata-y';

term1 = (djacfn(t,x,[],data))'*yB; % (f_y)^T \lambda

% When computing output sensitivities, assume dydx=L, ie not a
% function of x. Generalise later
dydx=L;

term2 = e*M.iCe*dydx;
yBd = (-term1 + term2')';

flag = 0;
new_data = [];

end

% -----------------------------------------------------------
function [JB, flag, new_data] = djacBfn(t, y, yB, fyB, data)
% Backward problem Jacobian function

J               = djacfn(t,y,[],data);
JB              = -J';

flag            = 0;
new_data        = [];

end

% -----------------------------------------------------------
function [J, flag, new_data] = djacfn(t, y, fy, data)
% State Jacobian at time t

P=data.P;
M=data.M;
U=data.U;

% Find nearest time point for which we have pre-computed input
if isempty(U)
    ut=[];
else
    [tmp,ind]=min(abs(t-M.t));
    ut=U(:,ind);
end

% Evaluate state Jacobian
if isfield(M,'dfdx')
    J = feval(M.dfdx,y,ut,P,M);
else
    J = spm_diff(M.f,y,ut,P,M,1);
end

flag            = 0;
new_data        = [];

end

author: spm


text:
function prob = problem_type(om, recheck)
%PROBLEM_TYPE  Return a string identifying the type of mathematical program
%   PROB_TYPE = OM.PROBLEM_TYPE()
%   PROB_TYPE = OM.PROBLEM_TYPE(RECHECK)
%
%   Returns a string identifying the type of mathematical program
%   represented by the current model, based on the variables, costs,
%   and constraints that have been added to the model. Used to
%   automatically select an appropriate solver.
%
%   Linear and nonlinear equations are models with no costs, no inequality
%   constraints, and an equal number of continuous variables and equality
%   constraints. If the number of variables in a nonlinear equation model
%   is one more than the number of constraints, it is a parameterized
%   nonlinear equation.
%
%   Outputs:
%       PROB_TYPE : problem type, one of the following strings:
%           'LEQ'   - linear equations
%           'NLEQ'  - nonlinear equations
%           'PNE'   - parameterized nonlinear equations
%           'LP'    - linear program
%           'QP'    - quadratic program
%           'NLP'   - nonlinear program
%           'MILP'  - mixed-integer linear program
%           'MIQP'  - mixed-integer quadratic program
%           'MINLP' - mixed-integer nonlinear program
%
%   The output value is cached for future calls, but calling with a true
%   value for the optional RECHECK argument will force it to recheck in
%   case the problem type has changed due to modifying the variables,
%   constraints or costs in the model.
%
%   See also OPT_MODEL

%   MP-Opt-Model
%   Copyright (c) 2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

if isempty(om.prob_type) || nargin > 1 && recheck
    nleN = om.getN('nle');      %% nonlinear equalities
    nliN = om.getN('nli');      %% nonlinear inequalities
    nlcN = om.getN('nlc');      %% general nonlinear costs
    qdcN = om.getN('qdc');      %% quadratic costs
    linN = om.getN('lin');      %% linear constraints
    varN = om.getN('var');      %% variables
    if varN == 0
        prob = '';
    elseif nlcN || qdcN         %% problem has costs
        if nliN || nleN || nlcN %% nonlinear
            prob = 'NLP';           %% nonlinear program
        else                    %% linear constraints, no general nonlinear costs
            %% get quadratic cost coefficients
            H = om.params_quad_cost();
            if isempty(H) || ~any(any(H))
                prob = 'LP';        %% linear program
            else
                prob = 'QP';        %% quadratic program
            end
        end
    else                    %% problem has no costs
        if nliN
            error('@opt_model/problem_type: invalid problem - nonlinear inequality constraints with no costs');
        end
        if nleN + linN == varN || nleN + linN == varN - 1   %% square (or almost) system
            if linN > 0
                %% get lower & upper bounds
                [A, l, u] = om.params_lin_constraint();
                if any(l ~= u)
                    error('@opt_model/problem_type: invalid problem - linear inequality constraints with no costs');
                end
            end
            if nleN + linN == varN  %% square system
                if nleN
                    prob = 'NLEQ';      %% square nonlinear set of equations
                else
                    prob = 'LEQ';       %% square linear set of equations
                end
            elseif nleN + linN + 1 == varN  %% square + 1 extra (parameterization) variable
                if nleN
                    prob = 'PNE';       %% parameterized nonlinear set of equations
                else
                    prob = 'PLEQ';      %% parameterized linear set of equations
                    error('@opt_model/problem_type: invalid problem - PNE not implemented for for linear constraints only');
                end
            else
                error('@opt_model/problem_type: invalid problem - PNE must have num of vars = num of constraints + 1');
            end
        else
            error('@opt_model/problem_type: invalid problem - non-square system with no costs');
        end
    end
    if om.is_mixed_integer() && ~strcmp(prob, 'NLEQ')
        prob = ['MI' prob];
    end
    om.prob_type = prob;    %% cache it
else
    prob = om.prob_type;    %% return cached type
end

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function num = meanwords(allSNumBi)
count = 0;
for i = 1:num(allSNumBi)

end

author: sidaw

repo: nbsvm

save_path: github-repos/MATLAB/sidaw-nbsvm


text:
function [powerterms_MR,powerterms_RZF,powerterms_MMMSE] = functionComputeULPowerLevels_impairments(H,Hhat,C,nbrOfRealizations,M,K,L,p,f,kappatUE,kapparBS)
%Compute and categorize the UL signal power for different receive combining
%schemes, under hardware impairments.
%
%INPUT:
%H                 = M x nbrOfRealizations x K x L x L matrix with the
%                    exact channel realizations
%Hhat              = M x nbrOfRealizations x K x L x L matrix with the MMSE
%                    channel estimates
%C                 = M x M x K x L x L matrix with estimation error
%                    correlation matrices when using MMSE estimation.
%nbrOfRealizations = Number of channel realizations
%M                 = Number of antennas per BS
%K                 = Number of UEs per cell
%L                 = Number of BSs and cells
%p                 = Uplink transmit power per UE (same for everyone)
%f                 = Pilot reuse factor
%kappatUE          = Hardware quality of the UEs' transmitters
%kapparBS          = Hardware quality of the BSs' receivers
%
%OUTPUT:
%powerterms_MR    = 6 x K x L matrix with signal and interference powers
%                   for each of the UEs with MR combining.
%                   powerterms_MR(:,k,l) is the vector for UE k in cell l,
%                   where the first element is the average desired signal
%                   power, the second element is the interference from UEs
%                   having the same pilot, the third term is the
%                   interference from UEs having different pilots, the
%                   fourth term is the transmitter distortion from other
%                   UEs, the fifth term is the receiver distortion from all
%                   UEs, and the last term is the self-distortion/interference
%powerterms_RZF   = Same as powerterms_MR but with RZF combining
%powerterms_MMMSE = Same as powerterms_MR but with M-MMSE combining
%
%
%This Matlab function was developed to generate simulation results to:
%
%Emil Bjornson, Jakob Hoydis and Luca Sanguinetti (2017),
%"Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency",
%Foundations and Trends in Signal Processing: Vol. 11, No. 3-4,
%pp. 154-655. DOI: 10.1561/2000000093.
%
%For further information, visit: https://www.massivemimobook.com
%
%This is version 1.0 (Last edited: 2017-11-04)
%
%License: This code is licensed under the GPLv2 license. If you in any way
%use this code for research that results in publications, please cite our
%monograph as described above.

%Store identity matrices of different sizes
eyeK = eye(K);
eyeM = eye(M);

%Generate pilot pattern
if f == 1

    pilotPattern = ones(L,1);

elseif f == 2 %Only works in the running example with its 16 BSs

    pilotPattern = kron(ones(2,1),[1; 2; 1; 2; 2; 1; 2; 1]);

elseif f == 4 %Only works in the running example with its 16 BSs

    pilotPattern = kron(ones(2,1),[1; 2; 1; 2; 3; 4; 3; 4]);

elseif f == 16 %Only works in the running example with its 16 BSs

    pilotPattern = (1:L)';

end

%Compute sum of all estimation error correlation matrices at every BS
C_totM = reshape(p*sum(sum(C,3),4),[M M L]);

%Prepare to store simulation results for the signal gains
signal_MR = zeros(K,L);
signal_RZF = zeros(K,L);
signal_MMMSE = zeros(K,L);

%Prepare to store simulation results for the norms of combining vectors
vectornorm_MR = zeros(K,L);
vectornorm_RZF = zeros(K,L);
vectornorm_MMMSE = zeros(K,L);

%Prepare to store simulation results for the sum of interference powers
interference_MR = zeros(K,L);
interference_RZF = zeros(K,L);
interference_MMMSE = zeros(K,L);

%Prepare to store simulation results for the sum of interference powers
%from UEs having the same pilot signals, including the user itself
reuseinterf_MR = zeros(K,L);
reuseinterf_RZF = zeros(K,L);
reuseinterf_MMMSE = zeros(K,L);

%Prepare to store simulation results for the self-interference power
selfinterf_MR = zeros(K,L);
selfinterf_RZF = zeros(K,L);
selfinterf_MMMSE = zeros(K,L);

%Prepare to store simulation results for the sum of receiver distortion
receiverdist_MR = zeros(K,L);
receiverdist_RZF = zeros(K,L);
receiverdist_MMMSE = zeros(K,L);

%% Go through all channel realizations
for n = 1:nbrOfRealizations

    %Go through all cells
    for j = 1:L

        %Extract channel realizations from all UEs to BS j
        Hallj = reshape(H(:,n,:,:,j),[M K*L]);

        %Extract channel estimate realizations from all UEs to BS j
        Hhatallj = reshape(Hhat(:,n,:,:,j),[M K*L]);

        %Extract cells that use same pilots as cell j
        groupMembers = find(pilotPattern==pilotPattern(j))';

        %Compute three different combining schemes
        V_MR = Hhatallj(:,K*(j-1)+1:K*j);
        V_RZF = p*V_MR/(p*(V_MR'*V_MR)+eyeK);
        V_MMMSE = p*(p*(Hhatallj*Hhatallj')+C_totM(:,:,j)+eyeM)\V_MR;

        %Go through all UEs in cell j
        for k = 1:K

            %Extract channel realizations from UEs that cause pilot
            %contamination to UE k in cell j
            HjkPC = reshape(H(:,n,k,groupMembers,j),[M length(groupMembers)]);

            %Extract channel from BS j to its k:th UE
            Hjk = H(:,n,k,j,j);

            %%MR combining
            w = V_MR(:,k)/norm(V_MR(:,k))^2; %Extract combining vector
            wrep = repmat(w,[1 K*L]);

            %Use Monte Carlo simulations to compute the expectations that
            %appear in the six different components in Section 6.3.3
            signal_MR(k,j) = signal_MR(k,j) + (w'*Hjk)/nbrOfRealizations;
            vectornorm_MR(k,j) = vectornorm_MR(k,j) + norm(w).^2/nbrOfRealizations;
            interference_MR(k,j) = interference_MR(k,j) + p*sum(abs(w'*Hallj).^2)/nbrOfRealizations;
            receiverdist_MR(k,j) = receiverdist_MR(k,j) + p*sum(sum(abs(wrep.*Hallj).^2,1))/nbrOfRealizations;
            reuseinterf_MR(k,j) = reuseinterf_MR(k,j) + p*sum(abs(w'*HjkPC).^2)/nbrOfRealizations;
            selfinterf_MR(k,j) = selfinterf_MR(k,j) + p*sum(abs(w'*Hjk).^2)/nbrOfRealizations;

            %%RZF combining
            w = V_RZF(:,k); %Extract combining vector
            wrep = repmat(w,[1 K*L]);

            %Use Monte Carlo simulations to compute the expectations that
            %appear in the six different components in Section 6.3.3
            signal_RZF(k,j) = signal_RZF(k,j) + (w'*Hjk)/nbrOfRealizations;
            vectornorm_RZF(k,j) = vectornorm_RZF(k,j) + norm(w).^2/nbrOfRealizations;
            interference_RZF(k,j) = interference_RZF(k,j) + p*sum(abs(w'*Hallj).^2)/nbrOfRealizations;
            receiverdist_RZF(k,j) = receiverdist_RZF(k,j) + p*sum(sum(abs(wrep.*Hallj).^2,1))/nbrOfRealizations;
            reuseinterf_RZF(k,j) = reuseinterf_RZF(k,j) + p*sum(abs(w'*HjkPC).^2)/nbrOfRealizations;
            selfinterf_RZF(k,j) = selfinterf_RZF(k,j) + p*sum(abs(w'*Hjk).^2)/nbrOfRealizations;

            %%MMMSE combining
            w = V_MMMSE(:,k); %Extract combining vector
            wrep = repmat(w,[1 K*L]);

            %Use Monte Carlo simulations to compute the expectations that
            %appear in the six different components in Section 6.3.3
            signal_MMMSE(k,j) = signal_MMMSE(k,j) + (w'*Hjk)/nbrOfRealizations;
            vectornorm_MMMSE(k,j) = vectornorm_MMMSE(k,j) + norm(w).^2/nbrOfRealizations;
            interference_MMMSE(k,j) = interference_MMMSE(k,j) + p*sum(abs(w'*Hallj).^2)/nbrOfRealizations;
            receiverdist_MMMSE(k,j) = receiverdist_MMMSE(k,j) + p*sum(sum(abs(wrep.*Hallj).^2,1))/nbrOfRealizations;
            reuseinterf_MMMSE(k,j) = reuseinterf_MMMSE(k,j) + p*sum(abs(w'*HjkPC).^2)/nbrOfRealizations;
            selfinterf_MMMSE(k,j) = selfinterf_MMMSE(k,j) + p*sum(abs(w'*Hjk).^2)/nbrOfRealizations;

        end

    end

end

%Set the normalized noise variance
sigma2 = 1;

%Compute the desired signal term, as described in Section 6.3.3
signalpower_MR = kappatUE*kapparBS*mean(mean(p*abs(signal_MR).^2./vectornorm_MR))/sigma2;
signalpower_RZF = kappatUE*kapparBS*mean(mean(p*abs(signal_RZF).^2./vectornorm_RZF))/sigma2;
signalpower_MMMSE = kappatUE*kapparBS*mean(mean(p*abs(signal_MMMSE).^2./vectornorm_MMMSE))/sigma2;

%Compute the interference from UEs having the same pilot, as described in
%Section 6.3.3
reuseinterference_MR = kappatUE*kapparBS*mean(mean((reuseinterf_MR-selfinterf_MR)./vectornorm_MR))/sigma2;
reuseinterference_RZF = kappatUE*kapparBS*mean(mean((reuseinterf_RZF-selfinterf_RZF)./vectornorm_RZF))/sigma2;
reuseinterference_MMMSE = kappatUE*kapparBS*mean(mean((reuseinterf_MMMSE-selfinterf_MMMSE)./vectornorm_MMMSE))/sigma2;

%Compute the interference from UEs having different pilots, as described in
%Section 6.3.3
noreuseinterference_MR = kappatUE*kapparBS*mean(mean((interference_MR-reuseinterf_MR)./vectornorm_MR))/sigma2;
noreuseinterference_RZF = kappatUE*kapparBS*mean(mean((interference_RZF-reuseinterf_RZF)./vectornorm_RZF))/sigma2;
noreuseinterference_MMMSE = kappatUE*kapparBS*mean(mean((interference_MMMSE-reuseinterf_MMMSE)./vectornorm_MMMSE))/sigma2;

%Compute the transmitter distortion from other UEs, as described in
%Section 6.3.3
transmitDistortion_MR = (1-kappatUE)*kapparBS*mean(mean((interference_MR-selfinterf_MR)./vectornorm_MR))/sigma2;
transmitDistortion_RZF = (1-kappatUE)*kapparBS*mean(mean((interference_RZF-selfinterf_RZF)./vectornorm_RZF))/sigma2;
transmitDistortion_MMMSE = (1-kappatUE)*kapparBS*mean(mean((interference_MMMSE-selfinterf_MMMSE)./vectornorm_MMMSE))/sigma2;

%Compute the receiver distortion from all UEs, as described in
%Section 6.3.3
receiverDistortion_MR = (1-kapparBS)*mean(mean(receiverdist_MR./vectornorm_MR))/sigma2;
receiverDistortion_RZF = (1-kapparBS)*mean(mean(receiverdist_RZF./vectornorm_RZF))/sigma2;
receiverDistortion_MMMSE = (1-kapparBS)*mean(mean(receiverdist_MMMSE./vectornorm_MMMSE))/sigma2;

%Compute the self-distortion and self-interference, as described in
%Section 6.3.3
selfinterference_MR = kapparBS*mean(mean((selfinterf_MR-kappatUE*p*abs(signal_MR).^2)./vectornorm_MR))/sigma2;
selfinterference_RZF = kapparBS*mean(mean((selfinterf_RZF-kappatUE*p*abs(signal_RZF).^2)./vectornorm_RZF))/sigma2;
selfinterference_MMMSE = kapparBS*mean(mean((selfinterf_MMMSE-kappatUE*p*abs(signal_MMMSE).^2)./vectornorm_MMMSE))/sigma2;

%Prepare to output the power values
powerterms_MR = [signalpower_MR reuseinterference_MR noreuseinterference_MR transmitDistortion_MR receiverDistortion_MR selfinterference_MR];
powerterms_RZF = [signalpower_RZF reuseinterference_RZF noreuseinterference_RZF transmitDistortion_RZF receiverDistortion_RZF selfinterference_RZF];
powerterms_MMMSE = [signalpower_MMMSE reuseinterference_MMMSE noreuseinterference_MMMSE transmitDistortion_MMMSE receiverDistortion_MMMSE selfinterference_MMMSE];

author: emilbjornson

repo: massivemimobook


text:
function LP = setCplexParam(LP, solverParams,verbFlag)
% Sets the parameters of the IBM ILOG CPLEX object according to the structure `solverParams`
% The `solverParams` structure has to contain the same structure as the
% Cplex.Param structue in a Cplex object. But values can be set by directly
% setting the respective parameter instead of the `Cur` value which would
% be necessary if directly modifying the Cplex.Param structure.
% USAGE:
%
%    LP = setCplexParam(LP, solverParams, verbFlag);
%
% INPUTS:
%    LP:              IBM-ILOG Cplex object
%    solverParams:    parameter structure for Cplex. Check `LP.Param`. For example,
%                     `[solverParams.simplex.display, solverParams.tune.display, solverParams.barrier.display, solverParams.sifting.display, solverParams.conflict.display] = deal(0);`
%                     `[solverParams.simplex.tolerances.optimality, solverParams.simplex.tolerances.feasibility] = deal(1e-9,1e-8);`
%
%                     The full set of parameters can be obtained by calling 'Cplex().Param'
%    verbFlag:        true to show which parameter input is problematic if any (optional, default true)

if nargin < 3
    verbFlag = true;
end
% if isempty(solverParams)
%     solverParams = struct();
% end
if isempty(fieldnames(solverParams))
    return
end

% get all the end parameters and their paths from LP.Param
%(e.g., LP.Param.simplex.display will have 'display' in paramList and 'LP.Param.simplex' in paramPath)
[paramList, paramPath] = getParamList(LP.Param, 0);

% similarly get all parameters and paths from the user-supplied parameter structure
% (but not going to the bottom level. User-supplied parameter structure no
% need to have *.Cur field after every parameter. But still supported if there is .Cur)
[paramUserList, paramUserPath, addCur] = getParamList(solverParams, 1);

% whether an user-supplied parameter can be matched to LP.Param
paramIden = false(numel(paramUserList), 1);
for p = 1:numel(paramUserList)
    f = strcmpi(paramList,paramUserList{p});
    if sum(f) == 1
        % uniquely matched
        paramIden(p) = true;
        str = ['LP.Param.' paramPath{f} '.Cur = solverParams.' paramUserPath{p}];
        if addCur(p)  % add back .Cur if in the user-supplied parameters
            str = [str '.Cur'];
        end
        str = [str ';'];
        eval(str);
    elseif sum(f) > 1
        % not uniquely matched (for parameters end with 'display' etc.)
        % compare the path then
        if ismember(lower(paramUserPath{p}), paramPath)
            paramIden(p) = true;
            str = ['LP.Param.' lower(paramUserPath{p}) '.Cur = solverParams.' paramUserPath{p}];
            if addCur(p)  % add back .Cur if in the user-supplied parameters
                str = [str '.Cur'];
            end
            str = [str ';'];
            eval(str);
        else
            if verbFlag
                warning('*.%s cannot be uniquely identified as a valid cplex parameter. Ignore.', paramUserPath{p});
            end
        end
    else
        % not matched at all
        if verbFlag
            warning('*.%s cannot be identified as a valid cplex parameter. Ignore.', paramUserPath{p});
        end
    end
end
end

author: opencobra

repo: cobratoolbox


text:
function opt= opt_overwriteVoids(opt, fld, default_value);

if ~isfield(opt, fld),
  return;
end

for k= 1:length(opt),
  if isempty(opt(k).(fld)),
    opt(k).(fld)= default_value;
  end
end

author: bbci


text:
function determ = r8mat_gedet ( a, n, pivot )

%*****************************************************************************80
%
%% R8MAT_GEDET computes the determinant of a matrix factored by R8MAT_GEFA.
%
%  Discussion:
%
%    This is a modified version of the LINPACK routine DGEDI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, (Society for Industrial and Applied Mathematics),
%    3600 University City Science Center,
%    Philadelphia, PA, 19104-2688.
%    ISBN 0-89871-172-X
%
%  Parameters:
%
%    Input, real A(N,N), the LU factors computed by R8MAT_GEFA.
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, integer PIVOT(N), as computed by R8MAT_GEFA.
%
%    Output, real DETERM, the determinant of the matrix.
%
  diag = r8mat_diag_get_vector ( n, a );

  determ = prod ( diag(1:n) );

  for i = 1 : n
    if ( pivot(i) ~= i )
      determ = - determ;
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function plot_branches(P,cover,segment,fig,ms,segind,BO)

n = nargin;
if n < 7
    BO = 1000;
    if n < 6
        segind = 1;
        if n < 5
            ms = 1;
            if n == 3
                fig = 1;
            end
        end
    end
end

Bal = cover.ball;
Segs = segment.segments;
SChi = segment.ChildSegment;
SPar = segment.ParentSegment;

if iscell(Segs{1})
    ns = max(size(Segs));
    Seg = cell(ns,1);
    for i = 1:ns
        m = size(Segs{i},1);
        S = zeros(0);
        for j = 1:m
            s = Segs{i}(j);
            s = s{:};
            S = [S; s];
        end
        Seg{i} = S;
    end
else
    Seg = Segs;
end

% Color the segments with unique colors
col = rand(ns,3);
for i = 2:ns
    C = col(SPar(i),:);
    c = col(i,:);
    while sum(abs(C-c)) < 0.2
        c = rand(1,3);
    end
    col(i,:) = c;
end

segments = segind;
C = SChi{segind};
b = 0;
while ~isempty(C) && b <= BO
    b = b+1;
    segments = [segments; C];
    C = SChi{segind};
end

ns = length(segment);
figure(fig)
for i = 1:ns
    if i == 2
        hold on
    end
    S = vertcat(Bal{Seg{segments(i)}});
    plot3(P(S,1),P(S,2),P(S,3),'.','Color',col(segments(i),:),'Markersize',ms)
end
hold off

author: InverseTampere

repo: TreeQSM


text:
function [ s, t ] = segment_point_coords_2d ( p1, p2, p )

%*****************************************************************************80
%
%% SEGMENT_POINT_COORDS_2D: coordinates of a point on a line segment in 2D.
%
%  Discussion:
%
%    A line segment is the finite portion of a line that lies between
%    two points P1 and P2.
%
%    By the coordinates of a point P with respect to a line segment [P1,P2]
%    we mean numbers S and T such that S gives us the distance from the
%    point P to the nearest point PN on the line (not the line segment!),
%    and T gives us the position of PN relative to P1 and P2.
%
%    If S is zero, then P lies on the line.
%
%    If 0 <= T <= 1, then PN lies on the line segment.
%
%    If both conditions hold, then P lies on the line segment.
%
%    If E is the length of the line segment, then the distance of the
%    point to the line segment is:
%
%      sqrt ( S**2 +  T**2    * E**2 )     if      T <= 0;
%             S                            if 0 <= T <= 1
%      sqrt ( S**2 + (T-1)**2 * E**2 )     if 1 <= T
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real P1(2,1), P2(2,1), the endpoints of the line segment.
%
%    Input, real P(2,1), the point to be considered.
%
%    Output, real S, the distance of P to the nearest point PN
%    on the line through P1 and P2.  (S will always be nonnegative.)
%
%    Output, real T, the relative position of the point PN
%    to the points P1 and P2.
%

%
%  If the line segment is actually a point, then the answer is easy.
%
  if ( p1(1:2,1) == p2(1:2,1) )

    t = 0.0;

  else

    bot = sum ( ( p2(1:2,1) - p1(1:2,1) ).^2 );

    t = ( p(1:2,1)  - p1(1:2,1) )' * ( p2(1:2,1) - p1(1:2,1) ) / bot;

  end

  pn(1:2,1) = p1(1:2,1) + t * ( p2(1:2,1) - p1(1:2,1) );

  s = sqrt ( sum ( ( p(1:2,1) - pn(1:2,1) ).^2 ) );

  return
end


text:
function detector = acfTrainBadacostTrees( varargin )
% Train aggregate channel features object detector.
%
% -------------------------------------------------------------------------
%  MODIFIED BY Jose M. Buenaposada
%  NOTE: Same as P.Dollar's implementation but changing to multiclass
%  classifiers for detection. For example, the detector can be trained
%  taking into account the car orientation. One of the classes is used as
%  negative (background).
% -------------------------------------------------------------------------
%
% Train aggregate channel features (ACF) object detector as described in:
%  P. Doll�r, R. Appel, S. Belongie and P. Perona
%   "Fast Feature Pyramids for Object Detection", PAMI 2014.
% The ACF detector is fast (30 fps on a single core) and achieves top
% accuracy on rigid object detection. Please see acfReadme.m for details.
%
% Takes a set of parameters opts (described in detail below) and trains a
% detector from start to finish including performing multiple rounds of
% bootstrapping if need be. The return is a struct 'detector' for use with
% acfDetect.m which fully defines a sliding window detector. Training is
% fast (on the INRIA pedestrian dataset training takes ~10 minutes on a
% single core or ~3m using four cores). Taking advantage of parallel
% training requires launching matlabpool (see help for matlabpool). The
% trained detector may be altered in certain ways via acfModify(). Calling
% opts=acfTrain() returns all default options.
%
% (1) Specifying features and model: The channel features are defined by
% 'pPyramid'. See chnsCompute.m and chnsPyramid.m for more details. The
% channels may be convolved by a set 'filters' to remove local correlations
% (see our NIPS14 paper on LDCF), improving accuracy but slowing detection.
% If 'filters'=[wFilter,nFilter] these are automatically computed. The
% model dimensions ('modelDs') define the window height and width. The
% padded dimensions ('modelDsPad') define the extended region around object
% candidates that are used for classification. For example, for 100 pixel
% tall pedestrians, typically a 128 pixel tall region is used to make a
% decision. 'pNms' controls non-maximal suppression (see bbNms.m), 'stride'
% controls the window stride, and 'cascThr' and 'cascCal' are the threshold
% and calibration used for the constant soft cascades. Typically, set
% 'cascThr' to -1 and adjust 'cascCal' until the desired recall is reached
% (setting 'cascCal' shifts the final scores output by the detector by the
% given amount). Training alternates between sampling (bootstrapping) and
% training an AdaBoost classifier (clf). 'nWeak' determines the number of
% training stages and number of trees after each stage, e.g. nWeak=[32 128
% 512 2048] defines four stages with the final clf having 2048 trees.
% 'pBoost' specifies parameters for AdaBoost, and 'pBoost.pTree' are the
% decision tree parameters, see adaBoostTrain.m for details. Finally,
% 'seed' is the random seed used and makes results reproducible and 'name'
% defines the location for storing the detector and log file.
%
% (2) Specifying training data location and amount: The training data can
% take on a number of different forms. The positives can be specified using
% either a dir of pre-cropped windows ('posWinDir') or dirs of full images
% ('posImgDir') and ground truth labels ('posGtDir'). The negatives can by
% specified using a dir of pre-cropped windows ('negWinDir'), a dir of full
% images without any positives and from which negatives can be sampled
% ('negImgDir'), and finally if neither 'negWinDir' or 'negImgDir' are
% given negatives are sampled from the images in 'posImgDir' (avoiding the
% positives). For the pre-cropped windows all images must have size at
% least modelDsPad and have the object (of size exactly modelDs) centered.
% 'imreadf' can be used to specify a custom function for loading an image,
% and 'imreadp' are custom additional parameters to imreadf. When sampling
% from full images, 'pLoad' determines how the ground truth is loaded and
% converted to a set of positive bbs (see bbGt>bbLoad). 'nPos' controls the
% total number of positives to sample for training (if nPos=inf the number
% of positives is limited by the training set). 'nNeg' controls the total
% number of negatives to sample and 'nPerNeg' limits the number of
% negatives to sample per image. 'nAccNeg' controls the maximum number of
% negatives that can accumulate over multiple stages of bootstrapping.
% Define 'pJitter' to jitter the positives (see jitterImage.m) and thus
% artificially increase the number of positive training windows. Finally if
% 'winsSave' is true cropped windows are saved to disk as a mat file.
%
% USAGE
%  detector = acfTrainBadacostTrees.m ( opts )
%  opts = acfTrainTrainBadacostTrees.m ()
%
% INPUTS
%  opts       - parameters (struct or name/value pairs)
%   (1) features and model:
%   .aRatioType - ['mean'] per class computation bb aspect ratio
%                 computation (see computePerClassAspectRatios)
%   .aRatioFixedWidth [0]
%   .pPyramid   - [{}] params for creating pyramid (see chnsPyramid)
%   .filters    - [] [wxwxnChnsxnFilter] filters or [wFilter,nFilter]
%   .modelDs    - [] model height+width without padding (eg [100 41])
%   .modelDsPad - [] model height+width with padding (eg [128 64])
%   .pNms       - [..] params for non-maximal suppression (see bbNms.m)
%   .stride     - [4] spatial stride between detection windows
%   .cascCal    - [.1] cascade calibration (affects speed/accuracy)
%   .nWeak      - [128] vector defining number weak clfs per stage
%   .pBoost     - [..] parameters for boosting (see badacostTrain.m)
%   .seed       - [0] seed for random stream (for reproducibility)
%   .name       - [''] name to prepend to clf and log filenames
%   (2) training data location and amount:
%   .posGtDir   - [''] dir containing ground truth
%   .posImgDir  - [''] dir containing full positive images
%   .negImgDir  - [''] dir containing full negative images
%   .posWinDir  - [''] dir containing cropped positive windows
%   .negWinDir  - [''] dir containing cropped negative windows
%   .imreadf    - [@imread] optional custom function for reading images
%   .imreadp    - [{}] optional custom parameters for imreadf
%   .pLoad      - [..] params for bbGt>bbLoad (see bbGt)
%   .nPos       - [inf] max number of pos windows to sample
%   .nNeg       - [5000] max number of neg windows to sample
%   .nPerNeg    - [25]  max number of neg windows to sample per image
%   .nAccNeg    - [10000] max number of neg windows to accumulate
%   .pJitter    - [{}] params for jittering pos windows (see jitterImage)
%   (3) save training results:
%   .winsSave   - [0] if true save cropped windows at each stage to disk
%   .savePath   - ['.'] path to save results in.
%
% OUTPUTS
%  detector   - trained object detector
%   .opts       - input parameters used for model training
%   .clf        - learned boosted tree classifier (see badacostTrain)
%   .info       - info about channels (see chnsCompute.m)
%
% EXAMPLE
%
% See also acfDetectBadacostTrees, acfTestBadacostTrees, chnsCompute,
% chnsPyramid, badacostTrain, bbGt, bbNms, jitterImage
%
% Piotr's Computer Vision Matlab Toolbox      Version NEW
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

% initialize opts struct
opts = initializeOpts( varargin{:} );
if(nargin==0), detector=opts; return; end

% load or initialize detector and begin logging
nm=fullfile(opts.savePath, [opts.name 'Detector.mat']); t=exist(nm,'file');
if(t), if(nargout), t=load(nm); detector=t.detector; end; return; end
t=fileparts(nm); if(~isempty(t) && ~exist(t,'dir')), mkdir(t); end
detector = struct( 'opts',opts, 'clf',[], 'info',[] );
startTrain=clock; nm=fullfile(opts.savePath, [opts.name 'Log.txt']);
if(exist(nm,'file')), diary(nm); diary('off'); delete(nm); end; diary(nm);
RandStream.setGlobalStream(RandStream('mrg32k3a','Seed',opts.seed));

% iterate bootstraping and training
for stage = 0:numel(opts.nWeak)-1
  diary('on'); fprintf([repmat('-',[1 75]) '\n']);
  fprintf('Training stage %i\n',stage); startStage=clock;

  % sample positives and compute info about channels
  if( stage==0 )
    [Is1,Ls1,IsOrig1] = sampleWins( detector, stage, 1 );
    t=ndims(Is1); if(t==3), t=Is1(:,:,1); else t=Is1(:,:,:,1); end
    t=chnsCompute(t,opts.pPyramid.pChns); detector.info=t.info;
  end

  % compute local decorrelation filters
  if( stage==0 && length(opts.filters)==2 )
    fs = opts.filters; opts.filters = [];
    X1 = chnsCompute1( IsOrig1, opts );
    fs = chnsCorrelation( X1, fs(1), fs(2) );
    opts.filters = fs; detector.opts.filters = fs;
  end

  % compute lambdas
  if( stage==0 && isempty(opts.pPyramid.lambdas) )
    fprintf('Computing lambdas... '); start=clock;
    ds=size(IsOrig1); ds(1:end-1)=1; IsOrig1=mat2cell2(IsOrig1,ds);
    ls=chnsScaling(opts.pPyramid.pChns,IsOrig1,0);
    ls=round(ls*10^5)/10^5; detector.opts.pPyramid.lambdas=ls;
    fprintf('done (time=%.0fs).\n',etime(clock,start));
  end

  % compute features for positives
  if( stage==0 )
    X1 = chnsCompute1( Is1, opts );
    X1 = reshape(X1,[],size(X1,4))';
    clear Is1 IsOrig1 ls fs ds t;
  end

  % sample negatives and compute features
  [Is0, ~] = sampleWins( detector, stage, 0 );
  X0 = chnsCompute1( Is0, opts ); clear Is0;
  X0 = reshape(X0,[],size(X0,4))';

  % accumulate negatives from previous stages
  if( stage>0 )
    n0=size(X0p,1); n1=max(opts.nNeg,opts.nAccNeg)-size(X0,1);
    if(n0>n1 && n1>0), X0p=X0p(randSample(n0,n1),:); end
    if(n0>0 && n1>0), X0=[X0p; X0]; end %#ok<AGROW>
  end; X0p=X0;

%   save_file = fullfile(opts.savePath, sprintf('%s_STAGE_%d_NEGATIVE_DATA.mat', opts.name, stage));
%   if ~exist(save_file, 'file')
%     save(save_file, 'X0p');
%   else
%     load(save_file);
%   end

  % train boosted BAdaCost clf and calibrate the cascade threshold.
  save_file_detector = fullfile(opts.savePath, sprintf('%s_STAGE_%d_Detector.mat', opts.name, stage));
  if ~exist(save_file_detector, 'file')
    detector.opts.pBoost = opts.pBoost;
    detector.opts.pBoost.nWeak = opts.nWeak(stage+1);
    [detector.clf, cost_curve] = badacostWithTreesTrain(X0, X1, Ls1, detector.opts.pBoost);

    % Plot and save cost curve figure
    h = figure;
    plot(cost_curve);
    xlabel('#weak learners');
    ylabel('cost');
    cost_curve_file = fullfile(opts.savePath, sprintf('%s_STAGE_%d_COST_CURVE', opts.name, stage));
    save([cost_curve_file '.txt'], '-ascii', 'cost_curve');
    figure(h);
    saveas(h, [cost_curve_file '.eps'], 'epsc');
%    saveas(h, [cost_curve_file '.fig'], 'fig');

%    [detector.opts.cascThr, h_calCascFig] = badacostCalibrateCascade(X0, X1, detector.clf, opts.cascCal);
    [detector.opts.cascThr, h_calCascFig] = badacostCalibrateCascade(X0, X1, detector.clf);
    % Watch out!!! This is faster but you can miss detections!!
%    detector.opts.cascThr=-2;

    % Compute the Training Confusion Matrix:
    xx = [X0; X1];
    yy = [ones(1,size(X0,1)), Ls1+1];
    yy_predicted = badacostWithTreesApply(xx, detector.clf);

    ConfMatrix = confusionmat(yy, yy_predicted);
    disp('Confusion Matrix:');
    disp(ConfMatrix);
    conf_mat_file = fullfile(opts.savePath, sprintf('%s_STAGE_%d_CONFUSION_MATRIX', opts.name, stage));
    save([conf_mat_file '.mat'], 'ConfMatrix');
    save([conf_mat_file '.txt'], '-ascii', 'ConfMatrix');

    % Plot cascade calibration figure
    casc_calib_file = fullfile(opts.savePath, sprintf('%s_STAGE_%d_CASCADE_CALIBRATION', opts.name, stage));
    figure(h_calCascFig);
    saveas(h_calCascFig,[casc_calib_file '.eps'], 'epsc');
%    saveas(h_calCascFig,[casc_calib_file '.fig'], 'fig');

    % Save trained detector
    save(save_file_detector, 'detector');
  else
    load(save_file_detector);

    %% Uncomment if you want to recompute the training confusion matrix.
    % Compute the Training Confusion Matrix:
    xx = [X0; X1];
    yy = [ones(1,size(X0,1)), Ls1+1];
    yy_predicted = badacostWithTreesApply(xx, detector.clf);

    ConfMatrix = confusionmat(yy, yy_predicted);
    disp('Confusion Matrix:');
    disp(ConfMatrix);
    conf_mat_file = fullfile(opts.savePath, sprintf('%s_STAGE_%d_CONFUSION_MATRIX', opts.name, stage));
    save([conf_mat_file '.mat'], 'ConfMatrix');
    save([conf_mat_file '.txt'], '-ascii', 'ConfMatrix');
  end

  % update log
  fprintf('Done training stage %i (time=%.0fs).\n',...
    stage,etime(clock,startStage)); diary('off');
end

% JMBUENA: Compute the window size for every subclass of the
% positive metaclass. We are going to keep this window size
% within the detector classifier. We will train the detector with
% fixed big window (selecting the features in this window) but we
% finally return the "per class" best fit window.
detector.clf.aRatio = computePerClassAspectRatios(detector.opts.posImgDir, ...
  detector.opts.posGtDir, detector.opts.pLoad, opts.aRatioType);
detector.clf.aRatioFixedWidth = opts.aRatioFixedWidth;

% save detector
save(fullfile(opts.savePath, [opts.name 'Detector.mat']),'detector');

% finalize logging
diary('on'); fprintf([repmat('-',[1 75]) '\n']);
fprintf('Done training (time=%.0fs).\n',...
  etime(clock,startTrain)); diary('off');
end

function opts = initializeOpts( varargin )
% Initialize opts struct.
dfs= {'aRatioType', 'mean', ...
  'pPyramid',{}, 'filters', [], ...
  'modelDs',[100 41], 'modelDsPad',[128 64], ...
  'pNms',struct(), 'stride',4, 'cascCal',.1, ...
  'nWeak',128, 'pBoost', {}, 'seed',0, 'name','', ...,
  'aRatioFixedWidth', 0, ...
  'posGtDir','', ...
  'posImgDir','', 'negImgDir','', 'posWinDir','', 'negWinDir','', ...
  'imreadf',@imread, 'imreadp',{}, 'pLoad',{}, 'nPos',inf, 'nNeg',5000, ...
  'nPerNeg',25, 'nAccNeg',10000, 'pJitter',{}, 'winsSave', 0, 'savePath', '.'};
opts = getPrmDflt(varargin,dfs,1);
% fill in remaining parameters
p=chnsPyramid([],opts.pPyramid); p=p.pPyramid;
p.minDs=opts.modelDs; shrink=p.pChns.shrink;
opts.modelDsPad=ceil(opts.modelDsPad/shrink)*shrink;
p.pad=ceil((opts.modelDsPad-opts.modelDs)/shrink/2)*shrink;
p=chnsPyramid([],p); p=p.pPyramid; p.complete=1;
p.pChns.complete=1; opts.pPyramid=p;
% initialize pNms, pBoost, pBoost.pTree, and pLoad
dfs={ 'type','maxg', 'overlap',.65, 'ovrDnm','min' };
opts.pNms=getPrmDflt(opts.pNms,dfs,-1);
dfs={ 'Cost',{}, 'stopAtNegWeight', 1, 'shrinkage', {}, 'resampling', {}, 'nWeak',0, ...
      'minDepth', 1, 'maxDepth', 4, 'verbose', 16, 'quantized', 0, ...
      'use_rus', 0, 'fracFtrs', 1, ...
      'variable_depth', 0};
opts.pBoost=getPrmDflt(opts.pBoost,dfs,1);
%dfs={'nBins',256,'maxDepth',2,'minWeight',.01,'fracFtrs',1,'nThreads',1e5};
%opts.pBoost.pTree=getPrmDflt(opts.pBoost.pTree,dfs,1);
opts.pLoad=getPrmDflt(opts.pLoad,{'squarify',{0,1}},-1);
opts.pLoad.squarify{2}=opts.modelDs(2)/opts.modelDs(1);
end

function [Is,Ls,IsOrig] = sampleWins( detector, stage, positive )
% Load or sample windows for training detector.
opts=detector.opts; start=clock;
if( positive ), n=opts.nPos; else n=opts.nNeg; end
% JMBUENA FIXME! How to deal with positive subclass labels from cropped images?
if( positive ), crDir=opts.posWinDir; else crDir=opts.negWinDir; end
if( exist(crDir,'dir') && stage==0 )
  % if window directory is specified simply load windows
  fs=bbGt('getFiles',{crDir}); nImg=length(fs); assert(nImg>0);
  if(nImg>n), fs=fs(:,randSample(nImg,n)); else n=nImg; end
  for i=1:n, fs{i}=[{opts.imreadf},fs(i),opts.imreadp]; end
  Is=cell(1,n); parfor i=1:n, Is{i}=feval(fs{i}{:}); end
else
  % sample windows from full images using sampleWins1()
  hasGt=positive||isempty(opts.negImgDir); fs={opts.negImgDir};
  if(hasGt), fs={opts.posImgDir,opts.posGtDir}; end
  fs=bbGt('getFiles',fs); nImg=size(fs,2); assert(nImg>0);
  if(~isinf(n)), fs=fs(:,randperm(nImg)); end;
  Is=cell(nImg*1000,1); Ls=cell(nImg*1000,1);
  diary('off'); tid=ticStatus('Sampling windows',1,30); k=0; i=0; batch=64;
  while( i<nImg && k<n )
    batch=min(batch,nImg-i); Is1=cell(1,batch); lbls=cell(1, batch);
    parfor j=1:batch, ij=i+j;
      I = feval(opts.imreadf,fs{1,ij},opts.imreadp{:}); %#ok<PFBNS>
      gt=[]; objs=[]; if(hasGt), [objs,gt]=bbGt('bbLoad',fs{2,ij},opts.pLoad); end
      [Is1{j}, lbls{j}]  = sampleWins1( I, objs, gt, detector, stage, positive );
    end
    Is1=[Is1{:}]; k1=length(Is1); Is(k+1:k+k1)=Is1;
    lbls=[lbls{:}]; Ls(k+1:k+k1)=lbls; k=k+k1;
    if(k>n), inds=randSample(k,n); Is=Is(inds); Ls=Ls(inds); k=n; end
    i=i+batch; tocStatus(tid,max(i/nImg,k/n));
  end
  Is=Is(1:k); Ls=Ls(1:k); diary('on');
  fprintf('Sampled %i windows from %i images.\n',k,i);
end
% optionally jitter positive windows
if(length(Is)<2), Is={}; return; end
nd=ndims(Is{1})+1; Is=cat(nd,Is{:}); Ls=cat(2,Ls{:}); IsOrig=Is;
if( positive && isstruct(opts.pJitter) )
  opts.pJitter.hasChn=(nd==4); Is=jitterImage(Is,opts.pJitter);
  ds=size(Is);
  Ls=repmat(Ls, [ds(nd),1]);
  Ls=reshape(Ls,[1,ds(nd)*length(Ls)]);
  ds(nd)=ds(nd)*ds(nd+1); Is=reshape(Is,ds(1:nd));
end
% make sure dims are divisible by shrink and not smaller than modelDsPad
ds=size(Is); cr=rem(ds(1:2),opts.pPyramid.pChns.shrink); s=floor(cr/2)+1;
e=ceil(cr/2); Is=Is(s(1):end-e(1),s(2):end-e(2),:,:); ds=size(Is);
if(any(ds(1:2)<opts.modelDsPad)), error('Windows too small.'); end
% optionally save windows to disk and update log
nm=fullfile(opts.savePath, [opts.name 'Is' int2str(positive) 'Stage' int2str(stage)]);
if( opts.winsSave ), save(nm,'Is','-v7.3'); end
fprintf('Done sampling windows (time=%.0fs).\n',etime(clock,start));
diary('off'); diary('on');
end

function [Is,lbls] = sampleWins1( I, objs, gt, detector, stage, positive )
% Sample windows from I given its ground truth gt.
opts=detector.opts; shrink=opts.pPyramid.pChns.shrink;
modelDs=opts.modelDs; modelDsPad=opts.modelDsPad;
if( positive ), bbs=gt; bbs=bbs(gt(:,5)==0,:); objs=objs(gt(:,5)==0,:); else
  if( stage==0 )
    % generate candidate bounding boxes in a grid
    [h,w,~]=size(I); h1=modelDs(1); w1=modelDs(2);
    n=opts.nPerNeg; ny=sqrt(n*h/w); nx=n/ny; ny=ceil(ny); nx=ceil(nx);
    [xs,ys]=meshgrid(linspace(1,w-w1,nx),linspace(1,h-h1,ny));
    bbs=[xs(:) ys(:)]; bbs(:,3)=w1; bbs(:,4)=h1; bbs=bbs(1:n,:);
  else
    % run detector to generate candidate bounding boxes
    bbs=acfDetectBadacost(I,detector); [~,ord]=sort(bbs(:,5),'descend');
    bbs=bbs(ord(1:min(end,opts.nPerNeg)),1:4);
  end
  if( ~isempty(gt) )
    % discard any candidate negative bb that matches the gt
    n=size(bbs,1); keep=false(1,n);
    for i=1:n, keep(i)=all(bbGt('compOas',bbs(i,:),gt,gt(:,5))<.1); end
    bbs=bbs(keep,:);
    objs=[];
  end
end
% grow bbs to a large padded size and finally crop windows
modelDsBig=max(8*shrink,modelDsPad)+max(2,ceil(64/shrink))*shrink;
r=modelDs(2)/modelDs(1); assert(all(abs(bbs(:,3)./bbs(:,4)-r)<1e-5));
r=modelDsBig./modelDs; bbs=bbApply('resize',bbs,r(1),r(2));
Is=bbApply('crop',I,bbs,'replicate',modelDsBig([2 1]));
if isempty(Is); lbls=Is; else
  lbls=cell(1,length(Is));
  if( positive )
    lbls(:)={objs(:).subclass};
  else
    lbls(:)={-1};
  end
end
end

function chns = chnsCompute1( Is, opts )
% Compute single scale channels of dimensions modelDsPad.
if(isempty(Is)), chns=[]; return; end
fprintf('Extracting features... '); start=clock; fs=opts.filters;
pChns=opts.pPyramid.pChns; smooth=opts.pPyramid.smooth;
dsTar=opts.modelDsPad/pChns.shrink; ds=size(Is); ds(1:end-1)=1;
Is=squeeze(mat2cell2(Is,ds)); n=length(Is); chns=cell(1,n);
parfor i=1:n
  C=chnsCompute(Is{i},pChns); C=convTri(cat(3,C.data{:}),smooth);
  if(~isempty(fs)), C=repmat(C,[1 1 size(fs,4)]);
    for j=1:size(C,3), C(:,:,j)=conv2(C(:,:,j),fs(:,:,j),'same'); end; end
  if(~isempty(fs)), C=imResample(C,.5); shr=2; else shr=1; end
  ds=size(C); cr=ds(1:2)-dsTar/shr; s=floor(cr/2)+1; e=ceil(cr/2);
  C=C(s(1):end-e(1),s(2):end-e(2),:); chns{i}=C;
end; chns=cat(4,chns{:});
fprintf('done (time=%.0fs).\n',etime(clock,start));
end

function filters = chnsCorrelation( chns, wFilter, nFilter )
% Compute filters capturing local correlations for each channel.
fprintf('Computing correlations... '); start=clock;
[~,~,m,n]=size(chns); w=wFilter; wp=w*2-1;
filters=zeros(w,w,m,nFilter,'single');
for i=1:m
  % compute local auto-scorrelation using Wiener-Khinchin theorem
  mus=squeeze(mean(mean(chns(:,:,i,:)))); sig=cell(1,n);
  parfor j=1:n
    T=fftshift(ifft2(abs(fft2(chns(:,:,i,j)-mean(mus))).^2));
    sig{j}=T(floor(end/2)+1-w+(1:wp),floor(end/2)+1-w+(1:wp));
  end
  sig=double(mean(cat(4,sig{mus>1/50}),4));
  sig=reshape(full(convmtx2(sig,w,w)),wp+w-1,wp+w-1,[]);
  sig=reshape(sig(w:wp,w:wp,:),w^2,w^2); sig=(sig+sig')/2;
  % compute filters for each channel from sig (sorted by eigenvalue)
  [fs,D]=eig(sig); fs=reshape(fs,w,w,[]);
  [~,ord]=sort(diag(D),'descend');
  fs=flipdim(flipdim(fs,1),2); %#ok<DFLIPDIM>
  filters(:,:,i,:)=fs(:,:,ord(1:nFilter));
end
fprintf('done (time=%.0fs).\n',etime(clock,start));
end

author: MengyangPu

repo: EDTER

save_path: github-repos/MATLAB/MengyangPu-EDTER


text:
function [h, class] = LMplot(varargin)
% Visualizes the polygons in an image.
%
% LMplot(annotation, img)
% or
% LMplot(database, ndx, HOMEIMAGES)
%
% Example:
%   [annotation, img] = LMread(filename, HOMEIMAGES)
%   LMplot(annotation, img)
%
%   thumbnail = size of the image
%   The plot uses only 7 colors, therefore, some times, different object.
%   Classes might have the same color assigned.
%
% If an object has the field 'confidence', the thickness of the bounding
% box will be equal to the confidence.
%
% If the object has the filed 'detection', then when the state is 'false'
% it will show the outline in green, and when the state is 'correct' it
% will show the outline in red.

switch length(varargin)
   case 1
       annotation = varargin{1};
       img = -1;
   case 2
       annotation = varargin{1};
       img = varargin{2};
   case 3
       D = varargin{1};
       ndx = varargin{2};

       if ndx>length(D)
           return
       end

       HOMEIMAGES = varargin{3};

       annotation = D(ndx).annotation;
       %img = imread(fullfile(HOMEIMAGES, annotation.folder, annotation.filename));
       img = LMimread(D, ndx, HOMEIMAGES); % Load image
   otherwise
       error('Too many input arguments.')
end

% Define colors
colors = 'rgbcmyw';
colors = hsv(15);

% Draw image
%figure;
if size(img,1)>1
    if isfloat(img)
        img = img-min(img(:));
        img = 255*img/max(img(:));
    end
    image(uint8(img)); axis('off'); axis('equal'); hold on
    if size(img,3) == 1
        colormap(gray(256))
    end
else
    axis('off'); axis('equal'); axis('ij'); hold on
end

% Draw each object (only non deleted ones)
h = []; class = [];
if isfield(annotation, 'object')
   Nobjects = length(annotation.object); n=0;
   for i = 1:Nobjects
       n = n+1;
       class{n} = annotation.object(i).name; % get object name
       col = colors(mod(sum(double(class{n})),15)+1, :);
       [X,Y] = getLMpolygon(annotation.object(i).polygon);
       X = double(X(:)); Y = double(Y(:));

       LineWidth = 4;
       if isfield(annotation.object(i), 'parts') && isfield(annotation.object(i).parts, 'ispartof')
           if ~isempty(annotation.object(i).parts.ispartof)
               class{n} = [class{n} ' (part)'];
               LineWidth = 1.5;
           end
       end

       if isfield(annotation.object(i), 'confidence')
           p = annotation.object(i).confidence;
           class{n} = [class{n} sprintf('[%2.2f]',p)];

           LineWidth = round(annotation.object(i).confidence);
           LineWidth = min(max([2 LineWidth]), 8);

           if isfield(annotation.object(i), 'detection')
               if strcmp(annotation.object(i).detection, 'false')
                   LineWidth = 1; %col = 'g';
               else
                   LineWidth = 4; %col = 'r';
               end
           end

           dim = get(gca,'position');
           if dim(3)>.5
               % If the plot occupies the entire figure, then we can plot the
               % confidence as a string.
               hold on
               xx = mean(X);
               %yy = mean(Y);
               yy = min(Y);
               ht=text(xx,yy, ...
                   sprintf('%5.2f', annotation.object(i).confidence), ...
                   'horizontalAlignment', 'center', 'verticalAlignment', 'bottom');
               %sprintf('%d %5.3f', i, annotation.object(i).confidence), ...
               %'horizontalAlignment', 'center', 'verticalAlignment', 'bottom');

               set(ht, 'color', col, 'fontsize', 10);
           else
               xx = mean(X);
               %yy = mean(Y);
               yy = min(Y);
               ht=text(xx,yy, class{n}, 'horizontalAlignment', 'center', 'verticalAlignment', 'bottom');
               set(ht, 'color', col, 'fontsize', 10);
           end
       end

       plot([X; X(1)],[Y; Y(1)], 'LineWidth', LineWidth, 'color', [0 0 0]); hold on
       h(n) = plot([X; X(1)],[Y; Y(1)], 'LineWidth', LineWidth/2, 'color', col);
       hold on
   end

   if nargout == 0
       legend(h, class);
   end
   %drawnow
end

author: CSAILVision

repo: LabelMeToolbox


text:
function buttonIcon = createButtonLabel(string,varargin)
% createButtonLabel: creates a tightly cropped RGB image (button-label) of a simple text string
%
%          This is very useful for labeling uicontrols (pushbuttons), with the syntax:
%
%                           set(hndl,'cdata',buttonIcon)
%
%          All valid Parameter-Value pairs, INCLUDING TEXT ROTATION, are
%          supported. Note that this function requires the Image Processing
%          Toolbox, and that it triggers the creation of a temporary
%          figure, which will be momentarily visible during
%          button-label creation.
%
%          Allows easy creation of multi-line text for labeling
%          pushbuttons.
%
% SYNTAX:  BUTTONICON = CREATEBUTTONLABEL(STRING,PVs,FIGOPT);
%
% USAGE: buttonIcon = createButtonLabel(string)
%              Create a tightly cropped RGB image of a text object of the
%              simple string indicated in input STRING. Uses default
%              properties (fontname, fontsize, fontweight, fontcolor,
%              bgcolor, rotation, etc.) for text object.
%
%         buttonIcon = createButtonLabel(string,PVs)
%              Provide additional Parameter-Value pairs, either in a cell
%              array of strings, or in a structure of parameters. All PV
%              pairs valid for text objects are supported.
%
%         buttonIcon = createButtonLabel(string,PVs,figOpt)
%              Provides optional argument for controlling the (requisite)
%              temporary figure. By default, the temporary figure will be
%              (visibly) created and destroyed within this function.
%              However, if you will be calling CREATEBUTTONLABEL multiply (labeling
%              tabs in a tabPanel, for instance), it might be less
%              distracting to reuse the same figure, and to destroy the
%              figure after the last call. This parameter will allow you to
%              do so.
%
% ARGUMENTS:
%
%          BUTTONICON: Output RGB image
%
%          STRING: The string to be captured as a tightly-cropped RGB
%              Image.
%
%          PVs (OPTIONAL): Cell array or structure of any parameter-value pairs
%                 valid for TEXT object. (See 'Text Properties' for details.)
%                 Particularly useful: {'rotation',90} or {'rotation',270}.
%
%          FIGOPT (OPTIONAL): Control parameter for temporary figure. Valid options:
%              1: REUSE/DESTROY (DEFAULT): reuse existing temporary figure, or create
%                 one if it doesn't exist. Destroy figure on completion.
%              2: REUSE existing temporary figure, or create one if it
%                 doesn't exist. Do not destroy figure.
%
% EXAMPLES:
%
%         1) Create cdata to label a vertically-oriented pushbutton, apply
%              structure of text parameters, destroy figure.
%
%         txt.rotation = 270;
%         txt.color = 'r';
%         txt.backgroundcolor = 'y';
%         buttonIcon = createButtonLabel('PRESS',txt);
%
%         2) Create cdata, apply PV pairs, leave figure open for re-use.
%
%         buttonIcon = createButtonLabel('PUSH',{'fontname','Helvetica','fontweight','bold'},2);
%
%         3) Multi-line text on a pushbutton:
%         buttonIcon =  createButtonLabel({'PUSH', 'ME!'},...
%                    {'fontname','Helvetica','fontweight',...
%                     'bold','horizontalAlignment','center',...
%                     'color','r'});
%
% SO: Why is this useful? Try this (and then try it without createButtonLabel!):
%
%         4)
%             x = createButtonLabel('PRESS',...
%                 {'rotation',270,'fontsize',18,'fontweight','b',...
%                 'color','c','backgroundcolor','y'});
%             figure;uicontrol('style','pushbutton','units','pixels',...
%                 'pos',[10 10 60 100],'callback','disp(''I was pushed!'')',...
%                 'backgroundcolor','y','cdata',x);

% Created by Brett Shoelson
% brett.shoelson@mathworks.com
% 05/14/07
%
% Modifications:
% 02/29/08 Force figure into 'normal' (not docked) state
% 06/22/09 Modified at the suggestion of Dave Forstot to allow
%    the creation of multi-line text buttons. (Thanks, Dave.)

% Validate and parse inputs
string = parseValidateInputs(string,varargin{:});

% Depending on figOpt parameter, branch to use/reuse/destroy temporary button-label figure
ExistFig = findobj('tag','tmpcreateButtonLabelfig');
if ~isempty(ExistFig)
    figure(ExistFig);
    clf;
else
    figure('units','pixels','position',[100 100 500 500],'tag','tmpcreateButtonLabelfig','windowstyle','normal');
end

ax = axes('units','pixels','position',[0 0 500 500],'visible','on','xtick',[],'ytick',[],'color','w');
tmp = text(250,250,string,'units','pixels','fontunits','pixels');

% Apply any input Parameter-Value pairs
applyPVs(tmp,PVs);
% Temporarily change fontcolor to black for detection, then return it to
% intended color
currColor = get(tmp,'color');
set(tmp,'color','k');
% Create button-label
F = getframe(ax);
im = im2double(F.cdata(:,:,1));
im = im < 1;
[y,x] = find(im);
set(tmp,'color',currColor);
set(ax,'color',get(tmp,'backgroundcolor'));
buttonIcon = getframe(ax,[min(x(:))-1 500-max(y(:))-1 max(x(:))-min(x(:))+2 max(y(:))-min(y(:))+2]);
buttonIcon = buttonIcon.cdata;

if figOpt == 1
    delete(findobj('tag','tmpcreateButtonLabelfig'));
end

% NESTED FUNCTIONS
    function string = parseValidateInputs(string,varargin)
        error(nargchk(1,3,nargin))
        error(nargoutchk(0,1,nargout))
        if ~(ischar(string) || iscellstr(string))
            try
                string = num2str(string);
            catch
                error('CREATEBUTTONLABEL: Requires at a minimum a single input, which must be a string.')
            end
        end
        % Defaults
        figOpt = 1;
        PVs = [];
        if ~isempty(varargin)
            if ~isempty(varargin{1})
                PVs = varargin{1};
            end
            if nargin == 3
                figOpt = varargin{2};
            end
            if ~ismember(figOpt,1:2)
                error('CREATEBUTTONLABEL: Invalid option for figOpt parameter.')
            end
        end
    end

    function applyPVs(obj,pvarray)
        if isempty(pvarray)
            return;
        end
        if isstruct(pvarray)
            set(obj,pvarray);
        else %Cell
            for ii = 1:2:numel(pvarray)
                set(obj,pvarray{ii},pvarray{ii+1});
            end
        end
    end

end

author: Sable

repo: mcbench-benchmarks


text:
function l = l1dn ( n, h )

%*****************************************************************************80
%
%% L1DN stores the 1D DN Laplacian as a full matrix.
%
%  Discussion:
%
%    The N grid points are assumed to be evenly spaced by H.
%
%    For N = 5, the discrete Laplacian with left Dirichlet and right
%    Neumann condition on [0,6] has the matrix form L:
%
%       2 -1  0  0  0
%      -1  2 -1  0  0
%       0 -1  2 -1  0
%       0  0 -1  2 -1
%       0  0  0 -1  1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 October 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of points.
%    N must be at least 3.
%
%    Input, real H, the spacing between points.
%
%    Output, real L(N,N), the Laplacian matrix.
%
  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'L1DN - Fatal error!\n' );
    fprintf ( 1, '  N < 3.\n' );
    error ( 'L1DN - Fatal error!' );
  end

  l = zeros ( n, n );

  i = 1;
  l(1,1) =  2.0 / h / h;
  l(1,2) = -1.0 / h / h;

  for i = 2 : n - 1
    l(i,i-1) = -1.0 / h / h;
    l(i,i) =    2.0 / h / h;
    l(i,i+1) = -1.0 / h / h;
  end

  i = n;
  l(n,n-1) = -1.0 / h / h;
  l(n,n) =    1.0 / h / h;

  return
end


text:
function drawRawPnts(SenFig,Raw)

% DRAWRAWPNTS  Draw raw points.
%   DRAWRAWPNTS(SENFIG,RAW)  redraws the 2D points in RAW on figure SenFig.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

if ~isempty(Raw.data.points.coord)
    set(SenFig.raw.points,...
        'xdata',Raw.data.points.coord(1,:),...
        'ydata',Raw.data.points.coord(2,:))
end

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
addpath ../../examples/imagenet;
load('./train_val_test_split.mat') ;
imdb.imageDir = '/home/zzd/Image-Text-Embedding/dataset/Flickr30k-prepare/flickr30k-images-256/'; %please change it to your full path

filename = dir([imdb.imageDir,'*jpg']);
cell_data = {filename(1:end).name}.';
full_path = cellfun(@(x) [imdb.imageDir,x],cell_data,'UniformOutput',false);
imdb.images.data = full_path;
imdb.images.set = set;

count = 0;
for i=1:numel(set)
    if(set(i)==1)
        count = count +1;
        imdb.images.label(i) = single(count);
    else
        imdb.images.label(i) = 0;
    end
end

% Compute image statistics (mean, RGB covariances, etc.)
train = find(imdb.images.set == 1) ;
images = fullfile(imdb.images.data(train(1:end))) ;
[averageImage, rgbMean, rgbCovariance] = getImageStats(images, ...
    'imageSize', [256 256], ...
    'numThreads', 16, ...
    'gpus', 3) ;
%imdb.averageImage = averageImage;
imdb.rgbMean = rgbMean;
imdb.rgbCovariance =rgbCovariance;

save('url_data.mat','imdb');

author: layumi

repo: Image-Text-Embedding


text:
function tMat = flatLevelTSeries(flat,scan,recomputeFlag);
% tMat = flatLevelTSeries(flat,[scan,recomputeFlag]);
%
% For Flat multi-level views:
%
% Create a 4D array of size x by y by slices by time,
% with the values from the tSeries interpolated across
% the (unmasked) flat mesh for each slice/time point.
% Also compute the mean image across levels for each
% hemisphere/time point.
%
% flat: flat view
%
% tSeries: tSeries for the flat view. If not
% entered, gets it from the view.
%
% scan: number of scan for tSeries. Defaults to
% current scan.
%
% recomputeFlag: if 1, will recompute even if a saved
% file already exists with the data. (Otherwise, will
% check and load existing files rather than recomputing.)
%
% Will also save a file 'interpolatedTSeries' in
% the scan's tSeries directory with the 4D matrix,
% as well as the appropriate coords/indices for
% mapping back and forth between the measured
% points on the flat map (in the normal tSeries
% files) and the interpolated image (same format
% as anat/amp/co/ph/map).
%
% ras, 10/2004
tic

if ieNotDefined('scan')
    scan = getCurScan(flat);
end

if ieNotDefined('recomputeFlag')
    recomputeFlag = 0;
end

% params
imageSize = flat.ui.imSize;
nSlices = 2 + sum(flat.numLevels);
nFrames = numFrames(flat,scan);
numLevels = flat.numLevels;
mask = flat.ui.mask;
savePath = fullfile(tSeriesDir(flat),sprintf('Scan%i',scan),'interpolatedTSeries.mat');

% check if the file already exists
% -- if so, load it and skip the computation:
if recomputeFlag==0 & exist(savePath,'file')
    fprintf('Loading %s ...\n',savePath);
    load(savePath,'tMat');
    return
end

% combine coords across hemispheres,
% get indices for going the reverse direction
coords = flat.coords;
indices = flat.indices;

% initialize the 4D array -- if there are memory
% constraints, this will hopefully give a hint:
% tMat = zeros(imageSize(1),imageSize(2),nSlices,nFrames);

% (There are --- making it a uint16 to reduce the load):
for i = 1:nFrames
    tMat(:,:,:,i) = uint16(zeros(imageSize(1),imageSize(2),nSlices));
end

%% get the min and max value for uint8 scaling ;(
%minVal = min(tSeries(find(tSeries)))
%maxVal = max(tSeries(find(tSeries)))

skipSlices = [];

waitHandle = mrvWaitbar(0,'Computing mean tSeries across gray levels...');
set(waitHandle,'Units','Normalized');
pos = get(waitHandle,'Position');
set(waitHandle,'Position',pos+[0 .1 0 0]);

for slice = 1:nSlices
    % load the tSeries slice
    tData = loadtSeries(flat,scan,slice);

        % need to figure out the hemisphere for this slice
        if slice==1 | (slice>2 & slice<=2+numLevels(1))
            h = 1;
        else
            h = 2;
        end

        % get coordinates from this slice
        subCoords = coords{slice}(1:2,:);

    if size(subCoords,2) < 10
        skipSlices = [skipSlices slice];
    else

        % interpolate each frame in tMat
        sliceWait = mrvWaitbar(0,'Interpolating slice across frames...');
        for t = 1:nFrames
            vals = tData(t,:);

            % The operator .' is the NON-CONJUGATE transpose.  Very important.
            tMat(:,:,slice,t) = myGriddata(subCoords,vals.',mask(:,:,h));

            mrvWaitbar(t/nFrames,sliceWait);
        end
    close(sliceWait);

    end

    mrvWaitbar(slice/nSlices,waitHandle);
end

close(waitHandle);

% % remove skipped slices
% tMat = tMat(:,:,setdiff(1:nSlices,skipSlices),:);

% save the 4D matrix
save(savePath,'tMat','coords','indices'); %,'minVal','maxVal');
fprintf('Saved %s. \n',savePath);

toc

return

% % % if arg out is specified, sample
% % % the mean tSeries at the data points
% % % measured (specified in coords/indices):
%  if nargout > 0
%  	waitHandle = mrvWaitbar(0,'Grabbing the mean tSeries...');
%
%  	for h = 1:2
%          % find coordinates of mean tSeries for this hemi
%          subCoords = find(coords(3,:)==h);
%
%         % get the indices for measured points for this hemi
%         measured = indices(find(indices(:,:,h)));
%
%         % error check: subCoords and measured should
%         % represent the same voxels:
%         [size(subCoords) size(measured)]
%         if size(subCoords)~=size(measured)
%             error('Logical error -- size mismatch b/w subCoords & measured.')
%         end
%
%         % I'll attempt to do a (fast) one-shot re-indexing,
%         % across time, but this requires that the indices be
%         % specified correctly. Replicate the 2D row/column
%         % specifications to run across frames, then get a
%         % single index into the big 4D matrix:
%         [x y] = ind2sub([imageSize(1) imageSize(2)],measured);
%         x = repmat(x,[1 nFrames]);
%         y = repmat(y,[1 nFrames]);
%         z = repmat(h,size(x));
%         t = repmat(1:nFrames,[length(measured) 1])';
%         t = t(:);
%         subIndices = sub2ind(size(tMat),x,y,z,t);
%
%         % assign the tSeries from the measured points
%         % (first convert back from uint8 to double)
%         subData = tMat(subIndices);
%         subData = minVal + maxVal .* double(subData);
%         tSeries(:,subCoords) = subData; % check for inaccuracies
%
%         mrvWaitbar(h/2,waitHandle);
%     end
%
%     close(waitHandle);
% end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ X ] = special_add_se3( X,S )

    s_theta=S(1:3);
    s_p=S(4:6);

    sizeS=size(S,1);
    NumberOfLandmarks=(sizeS-6)/3;

    Exps=so3_exp(s_theta);

    X.position= Exps*X.position+jaco_r(-s_theta)*s_p;

    if NumberOfLandmarks>=1
        s_landmarksMatrix=reshape(S(7:end),3,NumberOfLandmarks);
        X.landmarks(1:3,:)=X.landmarks(1:3,:)+s_landmarksMatrix;
    end

    X.orientation=Exps*X.orientation;

end

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
function template = PTKCreateTemplateForSegment(segment_mask, context, reporting)
    % PTKCreateTemplateForSegment. Function for creating a template image
    %     for a pulmonary segment regon
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    template = segment_mask.Copy;
end

author: tomdoel

repo: pulmonarytoolkit