text:
function [x,y,s,exitCode]=secondOrderConeInteriorPoint(A,b,c,blockLengths,params)
%%SECONDORDERCONEINTERIORPOINT Perform second order cone programming
%           (SOCP) using an interior point algorithm. This solves the
%           primal optimization problem:
%             minimize sum_{i=1}^n c_i'*x_i
%           subject to sum_{i=1}^n A_i*x_i=b
%             and x_i(1)-norm(x_i(2:end))>=0 (second order cone constraint)
%           which is equivalent to the dual optimization problem
%            maximize b'*y
%           subject to sum_{i=1}^n A_i'*y_i+s_i=c_i
%            and s_i(1)-norm(s_i(2:end))>=0
%           Many types of optimization problems can be formulated as SOCP
%           problems. This function only supports optimization in the real
%           domain.
%
%INPUTS: A A real mXsum(blockLengths) matrix containing coefficients for
%          all of the constraints on the primal. A must have full row rank.
%          With regard to the A_i in the problem formulation,
%          A=[A_1,a_2,...A_n]
%        b An mX1 real vector the coefficients of the dual objective
%          function to maximize.
%        c A sum(blockLengths)X1 vector containing the coefficients of the
%          primal objective function to minimize. With regard to the c_i in
%          the problem formulation, c=[c_1;c_2;...;c_n].
% blockLengths An nX1 vector containing the length of each of the cone
%          constraints. That is the first primal cone constraint is
%          x(1)-norm(x(2:blockLengths(1)))>=0. The second is
%          x(blockLengths(1)+1)-norm(x((blockLengths(1)+1):(blockLengths(1)+blockLengths(2)-1)))
%          and so on when considering the full stacked x vector (all of the
%          x_i stacked).
%   params An optional structure of parameters controlling the behaviour of
%          the algorithm. Possible members of the structure are:
%          gammaScal
%          'delta' A onstant parameter affecting the scale of a line search
%                  in the algorithm. 0<delta<1. The default if omitted is
%                  3/4.
%          'mu0'   The initial value of a positive real valued function
%                  affecting the smoothing function in the algorithm and
%                  thus the convergence. The default value if omitted is
%                  1/100.
%          'sigma' A constant parameter affecting the comvergence of the
%                  algorithm with delta. 0<sigma<1. The default if omitted
%                  is 1/4.
%          'maxIter' The maximum possible number of iterations to perform.
%                  The default if omitted is
%          'epsNormH' A tolerance value for declaring the convergence of
%                  the algorithm. This is the maximum value of norm(H) (H
%                  is defined in Equation 6 of [1]) below which convergence
%                  is declared. The default if omitted is 1e-6. Setting to
%                  zero will mean convergence is declared only when
%                  the minProgress contraint ocurs or the maximum number of
%                  iterations elapses.
%          'minProgress' A value such that the algorithm terminates if the
%                  change in normH is less than this value during multiple
%                  iterations. This allows termination when further
%                  iterations make little change. The default if omitted or
%                  an empty matrix is passed is epsNormH/(2*maxIter).
%                  Setting this to zero will make termination depend only
%                  on epsNormH and maxIter.
%
%OUTPUTS: x The sum(blockLengths)X1 set of stacked vectors solving the
%           primal optimization problem. x=[x_1;...;x_n] or an empty matrix
%           if an error occurred.
%         y The mX1 vector solving the dual optimization problem or an
%           empty matrix if an error occurred.
%         s The sum(blockLengths)X1 vector in the dual optimization
%           equality constraint or an empty matrix if an error occurred.
%  exitCode A value indicating how the algorithm terminated. Possible
%           values are
%           0 The accuracy constraint on norm(H) is satisfied (the
%             algorithm converged).
%           1 The maximum number of iterations elapsed without convergence.
%           2 Termination occurred because insufficient progress was made
%             (this depends on the params.minProgress value).
%           3 Termination occurred due to numerical problems.
%
%This function implements the algorithm described in [1]. The
%implementation slightly differs from the paper in that it can handle
%multiple cone constraints. The implementation does not take advantage of
%any sparsity in A.
%
%REFERENCES:
%[1] X. Chi and S. Liu, "A one-step smoothing Newton method for second-
%    order cone programming," Journal of Computational and Applied
%    Mathematics, vol. 223, no. 1, pp. 114-123, 1 Jan. 2009.
%
%February 2016 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

%ni(i) is the size of the ith block.
%Implementing algorithm 3.1 in the paper....

numBlocks=length(blockLengths);%The number of constraints.
m=size(A,1);
k=size(A,2);%The length of x and b.

%Get a cell array where blockIndices{i} holds the indices for the elements
%of the ith block.
blockIndices=makeIndexList(blockLengths);

%Default parameters and initialization before step 1.
gammaScal=0.95;
delta=3/4;%A constant.
sigma=1/4;%A constant
mu0=1/100;%A positive value.
maxIter=100;
epsNormH=1e-6;%Convergence criterion on the norm of H.
%Termination criterion for taking a minimum stepsize
minProgress=epsNormH/(2*maxIter);

%If the user provided any parameters
if(nargin>4&&~isempty(params))
    if(isfield(params,'gammaScal'))
        gammaScal=params.gammaScal;
    end

    if(isfield(params,'delta'))
        delta=params.delta;
    end

    if(isfield(params,'sigma'))
        sigma=params.sigma;
    end

    if(isfield(params,'mu0'))
        mu0=params.mu0;
    end

    if(isfield(params,'maxIter'))
        maxIter=params.maxIter;
    end

    if(isfield(params,'epsNormH'))
        epsNormH=params.epsNormH;
    end

    if(isfield(params,'minProgress'))
        minProgress=params.minProgress;
    end
end

%Initialization before step 1
mu=mu0*ones(numBlocks,1);
zBar=[zeros(k+m,1);mu];%A constant

%We want initial x values that satisfy the cone constraints. We could set
%then all zero, but that would put it right on the boundary, so we will
%make all of the first elements of the cone constraints 1 and the rest
%zero, which puts it a bit away from the boundary.
x=zeros(k,1);
curIdx=1;
for curBlock=1:numBlocks
   x(curIdx)=1;
   curIdx=curIdx+blockLengths(curBlock);
end

%Initial dual variables.
y=zeros(m,1);%Dual variables initially just zero.
s=c-A'*y;%s must satisfy A'*y+s=c.
H=HFun(x,s,mu,A,b,blockIndices);
normH=norm(H);
normHPrev=Inf;

lkPrev=0;

%gammaVal satisfies both of the constraints for gamma: It is less than 1
%(as gammaScal is less than 1) and it is less than 1/normH.
gammaVal=gammaScal*min(1,1/normH);

if(any(~isreal(H)|~isfinite(H)))
    x=[];
    y=[];
    s=[];
    exitCode=3;%Numerical errors.
    return;
end

for curIter=1:maxIter
%%Step 1
    if(normH<epsNormH)%If it converged.
        exitCode=0;
        return;
    elseif(normHPrev-normH<minProgress)
        %If normH did not sufficiently change between iterations, then
        %stop.
        exitCode=2;
        return;
    end

    %Rho is defined between equations 7 and 8.
    rho=gammaVal*normH*min(1,normH);

%%Step 2
    gradH=HFunGrad(x,s,mu,A,blockIndices);
    [HGradSqrt,didFail]=chol(gradH,'lower');
    if(didFail)%Using pinv for stability if poorly conditioned.
        deltaZ=lsqminnorm(gradH,(rho*zBar-H));
    else%Otherwise, solve using forward and backward substitution.
        opts.LT=true;
        opts.UT=false;
        %Solve using forward substitution.
        temp=linsolve(HGradSqrt,rho*zBar-H,opts);
        opts.UT=true;
        opts.LT=false;
        %Solve using backward substitution.
        deltaZ=linsolve(HGradSqrt',temp,opts);
    end
    deltaX=deltaZ(1:k);
    deltaY=deltaZ((k+1):(k+m));
    deltaMu=deltaZ((k+m+1):end);

%Steps 3 and 4
    %Search for the smallest allowable integer lk >=0. In this instance,
    %rather than just searching from 0 and going upwards, we start at the
    %last found value. In large problems, this can make a difference as lk
    %does not necessarily change much and is often not just 0.

    lk=lkPrev;
    lambda=delta^lk;
    xNew=x+lambda*deltaX;
    yNew=y+lambda*deltaY;
    muNew=mu+lambda*deltaMu;
    sNew=c-A'*yNew;
    HNew=HFun(xNew,sNew,muNew,A,b,blockIndices);
    normHNew=norm(HNew);
    boolValPrev=normHNew<=(1-sigma*(1-gammaVal*mu0)*lambda)*normH;

    while(1)
        if(boolValPrev==true)
            if(lk==0)
                break;
            end
            lk=lk-1;
        else
            lk=lk+1;
        end
        lambda=delta^lk;
        xNew=x+lambda*deltaX;
        yNew=y+lambda*deltaY;
        muNew=mu+lambda*deltaMu;
        sNew=c-A'*yNew;
        HNew=HFun(xNew,sNew,muNew,A,b,blockIndices);
        normHNew=norm(HNew);
        boolVal=normHNew<=(1-sigma*(1-gammaVal*mu0)*lambda)*normH;

        if(boolVal==true && boolValPrev==false)
            break;
        end
        boolValPrev=boolVal;
    end

    %We subtract one, because the loop above wants to try to bound the
    %correct lk value.
    lkPrev=max(lk-1,0);

    x=xNew;
    y=yNew;
    mu=muNew;
    s=sNew;
    H=HNew;
    normHPrev=normH;
    normH=normHNew;
end

%Maximum number of iterations completed.
exitCode=1;

end

function blockIndices=makeIndexList(blockLengths)
%%MAKEINDEXLIST This takes the vector of block lentghs and returns a cell
%               array where the ith element contains the indices of the
%               components in the ith block.

numBlocks=length(blockLengths);

blockIndices=cell(numBlocks,1);
curIdx=1;
for curBlock=1:numBlocks
    blockIndices{curBlock}=curIdx:(curIdx+blockLengths(curBlock)-1);
    curIdx=curIdx+blockLengths(curBlock);
end
end

function H=HFun(x,s,mu,A,b,blockIndices)
%%HFUN This implements the function H(z) from Equation 6. Note that s is
%%used instead of y, because c-A'*y is supposed to be passed as the second
%%argument to phiFun and that is equal to s.

H=[b-A*x;
   phiFun(x,s,mu,blockIndices);
   mu];
end

function gradG=HFunGrad(x,s,mu,A,blockIndices)
%%GRADH This function implements the gradient of H(z) in Equation 10. It
%       has been modified to handle more than one constraint.

m=size(A,1);
k=size(A,2);
numBlocks=length(blockIndices);

M=zeros(k,k);
N=zeros(k,k);
P=zeros(k,numBlocks);

%Allocate space and zero.
gradG=zeros(m+k+numBlocks,m+k+numBlocks);

gradG(1:m,1:k)=-A;%The first m rows.

for curBlock=1:numBlocks
    sel=blockIndices{curBlock};
    xCur=x(sel);
    sCur=s(sel);
    muCur=mu(curBlock);

    %e is defined toward the beginning of section 2.
    blockLength=length(sel);
    e=zeros(blockLength,1);
    e(1)=1;

    w1Bar=xCur+muCur*sCur;%s=c-A'*y
    w2Bar=muCur*xCur+sCur;
    wBar=vectorSqrt(vectorSquared(w1Bar)+vectorSquared(w2Bar)+2*muCur^2*e);

    %Lw is defined at the beginning of Section 2. Lw should be positive
    %definite, but is often pooly conditioned, so we use a pseudoinverse.
    Lw=Lx(wBar);
    [LwChol,cholDidFail]=chol(Lw,'lower');
    if(cholDidFail)
        LwInv=pinv(Lx(wBar));
    else
        opts.LT=true;
        opts.UT=false;
        %Solve using forward substitution.
        temp=linsolve(LwChol,eye(size(Lw)),opts);
        opts.UT=true;
        opts.LT=false;
        %Solve using backward substitution.
        LwInv=linsolve(LwChol',temp,opts);
    end

    Lw1=Lx(w1Bar);
    Lw2=Lx(w2Bar);

    M(sel,sel)=(1+muCur)*eye(blockLength,blockLength)-LwInv*(Lw1+muCur*Lw2);
    N(sel,sel)=(1+muCur)*eye(blockLength,blockLength)-LwInv*(muCur*Lw1+Lw2);
    P(sel,curBlock)=xCur+sCur-LwInv*(Lw1*sCur+Lw2*xCur+2*muCur*e);
end

gradG((m+1):(m+k),1:k)=M;
gradG((m+1):(m+k),(k+1):(k+m))=-N*A';
gradG((m+1):(m+k),(k+m+1):(k+m+numBlocks))=P;
gradG((m+k+1):(m+k+numBlocks),(m+k+1):(m+k+numBlocks))=eye(numBlocks,numBlocks);

end

function val=Lx(x)
%%LX Create the arrow-shaped matrix Lx that is given toward the beginning
%    of section 2. x must be a single set of cone variables.

    x0=x(1);
    xBar=x(2:end);
    xBarLen=length(xBar);

    val=[x0,xBar';
         xBar,x0*eye(xBarLen,xBarLen)];
end

function phi=phiFun(x,s,mu,blockIndices)
%%PHIBAR This implements phi(x,s,mu) from Equation 4. It has been modified
%%to deal with multiple blocks. in this instance, the value for each block
%%is stacked.

numBlocks=length(blockIndices);

phi=zeros(length(x),1);

%jordan_sqrt(jordan_square
for curBlock=1:numBlocks
    sel=blockIndices{curBlock};

    %e is defined toward the beginning of section 2.
    numIdx=length(sel);
    e=zeros(numIdx,1);
    e(1)=1;

    muCur=mu(curBlock);
    xCur=x(sel);
    sCur=s(sel);

    phi(sel)=(1+muCur)*(xCur+sCur)-vectorSqrt(vectorSquared(xCur+muCur*sCur)+vectorSquared(muCur*xCur+sCur)+2*muCur^2*e);
end

end

function xSquared=vectorSquared(x)
%%VECTORSQUARED In Section 2, a definition of the square a second order
%cone vector in terms of a spectral factorization is given. It is
%%implemented in this function.

[lambda,u]=vectorSpectralFact(x);

xSquared=lambda(1)^2*u(:,1)+lambda(2)^2*u(:,2);

end

function sqrtX=vectorSqrt(x)
%%VECTORSQRT In Section 2, a definition of the square root of a second
%%order cone vector in terms of a spectral factorization is given. It is
%%implemented in this function.

[lambda,u]=vectorSpectralFact(x);

%sqrtX should be real.
sqrtX=sqrt(lambda(1))*u(:,1)+sqrt(lambda(2))*u(:,2);

end

function [lambda,u]=vectorSpectralFact(x)
%%SPECTRALFACT In section 2, a spectral factorization of a second order
%       cone vector is defined. The vector is factored into spectral values
%       lambda and spectral vectors u. This function computes the
%       factorization.

x0=x(1);
x1=x(2:end);
x1Mag=norm(x1);

%lambda should be positive. However, finite precision errors might push it
%negative.
lambda=abs([x0-x1Mag;x0+x1Mag]);

x1Norm=x1/x1Mag;
if(~all(isfinite(x1Norm)))
    x1Norm=zeros(size(x1));
    x1Norm(1)=1;
end

u=(1/2)*[[1;-x1Norm],[1;x1Norm]];
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function check = fisk_check ( a, b, c )

%*****************************************************************************80
%
%% FISK_CHECK checks the parameters of the Fisk PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, the parameters of the PDF.
%    0.0 < B,
%    0.0 < C.
%
%    Output, logical CHECK, is true if the parameters are legal.
%
  if ( b <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FISK_CHECK - Fatal error!\n' );
    fprintf ( 1, '  B <= 0.\n' );
    check = 0;
    return
  end

  if ( c <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FISK_CHECK - Fatal error!\n' );
    fprintf ( 1, '  C <= 0.\n' );
    check = 0;
    return
  end

  check = 1;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ Y ] = ATMAV( X, h )
% ATMED  - The fuzzy filter with Triangular function with MEDian center
% x     - input (noisy image)
% y     - output (de-noised image)
% F     - window function, weighting function

if nargin < 2
    h=3;    % local search window size: h=3; h=5; h=7; h=9; h=11;
    f=1;    % padding value for window: f=1; f=2; f=3; f=4; f=5;
else
    f=(h-1)/2;    % padding value for window: f=1; f=2; f=3; f=4; f=5;
end
F = zeros(h*h,1);
Y  = zeros(size(X));

X = double(padarray(X,[f f],'symmetric'));
[m n ~] = size(X);

for i=1+f:1:m-f
   for j=1+f:1:n-f
       x=reshape(X(i-f:i+f, j-f:j+f),[],1);
       xmin = min(x);
       xmav = mean(x);
       xmax = max(x);

       F(:,:) = 0;
       if (xmav-xmin==0)||(xmax-xmav==0)
           F(:,:) = 1;
       else
           ind1 = find((x>=xmin)&(x<=xmav));
           F(ind1) = 1-(xmav-x(ind1))/(xmav-xmin);

           ind2 = find((x>=xmav)&(x<=xmax));
           F(ind2) = 1-(x(ind2)-xmav)/(xmax-xmav);
       end

       Y(i-f,j-f) = sum(sum(F.*x))/sum(sum(F));
       clear xmax xmin xmav ind1 ind2;
   end
end
end

author: Sable

repo: mcbench-benchmarks


text:
function [h,miss,stds] = plotRoc( D, varargin )
% Function for display of rocs (receiver operator characteristic curves).
%
% Display roc curves. Consistent usage ensures uniform look for rocs. The
% input D should have n rows, each of which is of the form:
%  D = [falsePosRate truePosRate]
% D is generated, for example, by scanning a detection threshold over n
% values from 0 (so first entry is [1 1]) to 1 (so last entry is [0 0]).
% Alternatively D can be a cell vector of rocs, in which case an average
% ROC will be shown with error bars. Plots missRate (which is just 1 minus
% the truePosRate) on the y-axis versus the falsePosRate on the x-axis.
%
% USAGE
%  [h,miss,stds] = plotRoc( D, prm )
%
% INPUTS
%  D    - [nx2] n data points along roc [falsePosRate truePosRate]
%         typically ranges from [1 1] to [0 0] (or may be reversed)
%  prm  - [] param struct
%   .color    - ['g'] color for curve
%   .lineSt   - ['-'] linestyle (see LineSpec)
%   .lineWd   - [4] curve width
%   .logx     - [0] use logarithmic scale for x-axis
%   .logy     - [0] use logarithmic scale for y-axis
%   .marker   - [''] marker type (see LineSpec)
%   .mrkrSiz  - [12] marker size
%   .nMarker  - [5] number of markers (regularly spaced) to display
%   .lims     - [0 1 0 1] axes limits
%   .smooth   - [0] if T compute lower envelop of roc to smooth staircase
%   .fpTarget - [] return miss rates at given fp values (and draw lines)
%   .xLbl     - ['false positive rate'] label for x-axis
%   .yLbl     - ['miss rate'] label for y-axis
%
% OUTPUTS
%  h    - plot handle for use in legend only
%  miss - average miss rates at fpTarget reference values
%  stds - standard deviation of miss rates at fpTarget reference values
%
% EXAMPLE
%  k=2; x=0:.0001:1; data1 = [1-x; (1-x.^k).^(1/k)]';
%  k=3; x=0:.0001:1; data2 = [1-x; (1-x.^k).^(1/k)]';
%  hs(1)=plotRoc(data1,struct('color','g','marker','s'));
%  hs(2)=plotRoc(data2,struct('color','b','lineSt','--'));
%  legend( hs, {'roc1','roc2'} ); xlabel('fp'); ylabel('fn');
%
% See also
%
% Piotr's Computer Vision Matlab Toolbox      Version 3.02
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

% get params
[color,lineSt,lineWd,logx,logy,marker,mrkrSiz,nMarker,lims,smooth, ...
  fpTarget,xLbl,yLbl] = getPrmDflt( varargin, {'color' 'g' 'lineSt' '-' ...
  'lineWd' 4 'logx' 0 'logy' 0 'marker' '' 'mrkrSiz' 12 'nMarker' 5 ...
  'lims' [] 'smooth' 0 'fpTarget' [] 'xLbl' 'false positive rate' ...
  'yLbl' 'miss rate' } );
if( isempty(lims) ); lims=[logx*1e-5 1 logy*1e-5 1]; end

% ensure descending fp rate, change to miss rate, optionally 'nicefy' roc
if(~iscell(D)), D={D}; end; nD=length(D);
for j=1:nD, assert(size(D{j},2)==2); end
for j=1:nD, if(D{j}(1,2)<D{j}(end,2)), D{j}=flipud(D{j}); end; end
for j=1:nD, D{j}(:,2)=1-D{j}(:,2); assert(all(D{j}(:,2)>=0)); end
if(smooth), for j=1:nD, D{j}=smoothRoc(D{j}); end; end

% plot: (1) h for legend only, (2) markers, (3) error bars, (4) roc curves
hold on; axis(lims); xlabel(xLbl); ylabel(yLbl);
prmMrkr = {'MarkerSize',mrkrSiz,'MarkerFaceColor',color};
prmClr={'Color',color}; prmPlot = [prmClr,{'LineWidth',lineWd}];
h = plot( 2, 0, [lineSt marker], prmMrkr{:}, prmPlot{:} ); %(1)
DQ = quantizeRocs( D, nMarker, logx, lims ); DQm=mean(DQ,3);
if(~isempty(marker))
  plot(DQm(:,1),DQm(:,2),marker,prmClr{:},prmMrkr{:} ); end %(2)
if(nD>1), DQs=std(DQ,0,3);
  errorbar(DQm(:,1),DQm(:,2),DQs(:,2),'.',prmClr{:}); end %(3)
if(nD==1), DQ=D{1}; else DQ=quantizeRocs(D,100,logx,lims); end
DQm = mean(DQ,3); plot( DQm(:,1), DQm(:,2), lineSt, prmPlot{:} ); %(4)

% plot line at given fp rate
m=length(fpTarget); miss=zeros(1,m); stds=miss;
if( m>0 )
  assert( min(DQm(:,1))<=min(fpTarget) ); DQs=std(DQ,0,3);
  for i=1:m, j=find(DQm(:,1)<=fpTarget(i)); j=j(1);
    miss(i)=DQm(j,2); stds(i)=DQs(j,2); end
  fp=min(fpTarget); plot([fp fp],lims(3:4),'Color',.7*[1 1 1]);
  fp=max(fpTarget); plot([fp fp],lims(3:4),'Color',.7*[1 1 1]);
end

% set log axes
if( logx==1 )
  ticks=10.^(-8:8);
  set(gca,'XScale','log','XTick',ticks);
end
if( logy==1 )
  ticks=[.001 .002 .005 .01 .02 .05 .1 .2 .5 1];
  set(gca,'YScale','log','YTick',ticks);
end
if( logx==1 || logy==1 ), grid on;
  set(gca,'XMinorGrid','off','XMinorTic','off');
  set(gca,'YMinorGrid','off','YMinorTic','off');
end

end

function DQ = quantizeRocs( Ds, nPnts, logx, lims )
% estimate miss rate at each target fp rate
nD=length(Ds); DQ=zeros(nPnts,2,nD);
if(logx==1), fps=logspace(log10(lims(1)),log10(lims(2)),nPnts);
else fps=linspace(lims(1),lims(2),nPnts); end; fps=flipud(fps');
for j=1:nD, D=[Ds{j}; 0 1]; k=1; fp=D(k,1);
  for i=1:nPnts
    while( k<size(D,1) && fp>=fps(i) ), k=k+1; fp=D(k,1); end
    k0=max(k-1,1); fp0=D(k0,1); assert(fp0>=fp);
    if(fp0==fp), r=.5; else r=(fps(i)-fp)/(fp0-fp); end
    DQ(i,1,j)=fps(i); DQ(i,2,j)=r*D(k0,2)+(1-r)*D(k,2);
  end
end
end

function D1 = smoothRoc( D )
D1 = zeros(size(D));
n = size(D,1); cnt=0;
for i=1:n
  isAnkle = (i==1) || (i==n);
  if( ~isAnkle )
    dP=D1(cnt,:); dC=D(i,:); dN=D(i+1,:);
    isAnkle = (dC(1)~=dP(1)) && (dC(2)~=dN(2));
  end
  if(isAnkle); cnt=cnt+1; D1(cnt,:)=D(i,:); end
end
D1=D1(1:cnt,:);
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
% Euclidean projection on a hyperplane
% Section 8.1.1, Boyd & Vandenberghe "Convex Optimization"
% Joelle Skaf - 10/04/05
%
% The projection of x0 on a hyperplane C = {x | a'*x = b} is given by
%           minimize || x - x0 ||^2
%               s.t.    a'*x = b
% It is also given by P_C(x0) = x0 + (b - a'*x0)*a/||a||^2

% Input data
randn('seed',0);
n  = 10;
a  = randn(n,1);
b  = randn(1);
x0 = randn(n,1);

% Analytical solution
fprintf(1,'Computing the analytical solution ...');
pc_x0 = x0 + (b - a'*x0)*a/norm(a)^2;
fprintf(1,'Done! \n');

% Solution via QP
fprintf(1,'Computing the optimal solution by solving a QP ...');

cvx_begin quiet
    variable x(n)
    minimize ( square_pos(norm(x - x0)) )
    a'*x == b;
cvx_end

fprintf(1,'Done! \n');

% Verification
disp('--------------------------------------------------------------------------------');
disp('Verifying that p_C(x0) and x_star belong to the hyperplane C: ');
disp(['a^T*p_C(x0) - b = ' num2str(a'*pc_x0 - b)]);
disp(['a^T*x_star - b  = ' num2str(a'*x - b)]);
disp('Computing the distance between x0 and the hyperplane in each case');
disp(['||x0 - p_C(x0)|| = ' num2str(norm(x0 - pc_x0))]);
disp(['||x0 - x_star || = ' num2str(norm(x0 - x))]);
disp('Verifying that the analytical solution and the solution obtained via QP are equal: ');
[pc_x0 x]

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
%returns the first component of the frequency representation of the signal
classdef SpectralEntropy < Algorithm

    methods (Access = public)

        function obj = SpectralEntropy()
            obj.name = 'SpectralEntropy';
            obj.inputPort = DataType.kSignal;
            obj.outputPort = DataType.kFeature;
        end

        %receives a power spectrum
        function result = compute(~,powerSpectrum)

            %Normalization
            sumPower = sum(powerSpectrum + 1e-12);
            powerSpectrum = powerSpectrum / sumPower;

            %entropy calculation
            result = -sum(powerSpectrum .* log2(powerSpectrum+eps));
        end

        function metrics = computeMetrics(~,input)
            n = size(input,1);
            flops = 21 * n;
            memory = n * 4;
            outputSize = Constants.kFeatureBytes;
            metrics = Metric(flops,memory,outputSize);
        end
    end
end

author: avenix

repo: WDK

save_path: github-repos/MATLAB/avenix-WDK


text:
%% Copyright (C) 2014, 2016, 2019 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy is free software; you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published
%% by the Free Software Foundation; either version 3 of the License,
%% or (at your option) any later version.
%%
%% This software is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty
%% of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See
%% the GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public
%% License along with this software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defmethod @@sym fliplr (@var{A})
%% Flip a symbolic matrix horizontally.
%%
%% Example:
%% @example
%% @group
%% A = sym([1 2 pi; 4 5 6]);
%% fliplr (A)
%%   @result{} (sym 2×3 matrix)
%%       ⎡π  2  1⎤
%%       ⎢       ⎥
%%       ⎣6  5  4⎦
%% @end group
%% @end example
%%
%% @seealso{@@sym/flipud, @@sym/reshape}
%% @end defmethod

function B = fliplr (A)

  cmd = { 'A, = _ins'
          'if A is None or not A.is_Matrix:'
          '    A = sp.Matrix([A])'
          'return A[:, ::-1]' };

  B = pycall_sympy__ (cmd, sym(A));

end

%!test
%! % simple
%! syms x
%! A = [x 2; sym(pi) x];
%! B = [2 x; x sym(pi)];
%! assert (isequal (fliplr(A), B))

%!test
%! % simple, odd # cols
%! syms x
%! A = [x 2 sym(pi); x 1 2];
%! B = [sym(pi) 2 x; 2 1 x];
%! assert (isequal (fliplr(A), B))

%!test
%! % scalar
%! syms x
%! assert (isequal (fliplr(x), x))

repo: octsympy


text:
function [ X,i ] = gmm_sample(nbSamples,Priors,Mu,Sigma )
%GMM_SAMPLE Samples points from a Gaussian Mixture Model
%
%   input -----------------------------------------------------------------
%
%       o nbSamples : (1 x 1), number of samples to return
%
%       o Priors    : (1 x K), gmm weights
%
%       o Mu        : (D x K), gmm means
%
%       o Sigma     : (D x D x K), gmm covariance
%
%   output ----------------------------------------------------------------
%
%       o X         : (D x N), N samples of dimension D
%
%       o i         : (N x 1), index of Gaussian center from which the data
%                              points were sampled form.

D = size(Mu,1);
X = zeros(D,nbSamples);

i = discretesample(Priors', nbSamples);

if D == 1

    for v=1:nbSamples

        X(v) = normrnd(Mu(i(v)), sqrt(Sigma(i(v))) );

    end
else

    for v = 1:nbSamples
        j = i(v);
        X(:,v) = mvnrnd(Mu(:,j),Sigma(:,:,j));
    end

end

end

author: epfl-lasa


text:
% calculate non-negative leaset squares
function H = calc_nls_nmf(X, W, lambda)
    R = size(W,2);
    %H = inv(W'*W + lambda * eye(R)) * W'* X;
    H = (W'*W + lambda * eye(R)) \ W'* X;
    H = max(H, 1e-16);
end

author: hiroyuki-kasai

repo: NMFLibrary


function [price, stdErr] = Price_MC_Lookback_func(Spath, call, M, mult, disc)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% About: Calculates Looback option price, given the simulatd paths
% Author: Justin Lars Kirkby
%
% -----------------
% Params
% -----------------
% Contract Types
% Floating Strike (Lookback) Put:  max{S_m: 0<=m<=M} - S_T
% Floating Strike (Lookback) Call: S_T - min{S_m: 0<=m<=M}
%
% Spath = paths of underlying, dimension N_sim x M+1, where M = number of time steps (since includes S_0)
% call = 1 for call (else put)
% M = number of monitoring points, e.g. 252 for "daily" monitoring
% mult = time partitioning multiplier to reduce bias (e.g. mult = 2 or 5)
% S_0 = initial price
% disc = discount factor (e.g. exp(-r*T))
% T = Time (in years)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N_sim = size(Spath,1);  % number of paths

M_mult = M*mult;  %time partitioning to reduce bias

if call ~= 1
    curr_max = zeros(N_sim, 1);
    for n = 1:N_sim
        curr_max(n) = max(Spath(n, 1:mult:M_mult+1));
    end
else  % find_min
    curr_min = zeros(N_sim, 1);
    for n = 1:N_sim
        curr_min(n) = min(Spath(n, 1:mult:M_mult+1));
    end
end

if call ==1  % Floating Strike (Lookback) Call: S_T - min{S_m: 0<=m<=M}
   payoffs  = Spath(:, M_mult+1) - curr_min;

else  % Floating Strike (Lookback) Put:  max{S_m: 0<=m<=M} - S_T
    payoffs  = curr_max - Spath(:,M_mult+1);
end
price = disc*mean(payoffs);
stdErr = disc*std(payoffs) / sqrt(N_sim);

end


text:
function [c, ceq] = ma_surfaceVelConstraint(stateLog, eventID, lbHorzVel, ubHorzVel, bodyIDApply, celBodyData, maData)
%ma_semiMajorAxisConstraint Summary of this function goes here
%   Detailed explanation goes here

    normFact = 1;

    if(ischar(eventID) && strcmpi(eventID,'final'))
        eventNum = max(stateLog(:,13));
    else
        [~, eventNum] = getEventByID(eventID, maData.script);
    end

    eventLog = stateLog(stateLog(:,13)==eventNum,:);
    finalEntry = eventLog(end,:);

    bodyID = finalEntry(8);

    if(bodyID == bodyIDApply || bodyIDApply==-1)
        horzVel = ma_GALongLatAltTasks(finalEntry, 'horzVel', celBodyData);

        if(lbHorzVel == ubHorzVel)
            c = [0 0];
            ceq(1) = horzVel - ubHorzVel;
        else
            c(1) = lbHorzVel - horzVel;
            c(2) = horzVel - ubHorzVel;
            ceq = [0];
        end
        c = c/normFact;
        ceq = ceq/normFact;
    else
        c = [0 0];
        ceq = [0];
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [hdr] = read_wdq_header(filename)

% READ_WDQ_HEADER reads header information from wdq files
%
% Use as
%  [hdr] = read_wdq_header(filename)
%

% Copyright (C) 2011, Jan-Mathijs Schoffelen
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

% information about how to interpret the file are taken from the document
% 'CODAS data storage format'

fid = fopen_or_error(filename, 'r');
%tmp = fread(fid, 110); % the first 110 bytes are always interpretable equally
hdr = [];
hdr.nchan         = fread(fid, 1, 'uint16=>uint16');
% hdr.nchan         = hdr.nchan - bitand(hdr.nchan, 32);
hdr.nchan         = double(bitand(hdr.nchan, 255)); % only the lowest bits seem to be coding the number of channels
i2                = fread(fid, 1, 'uint16');
hdr.offset        = fread(fid, 1, 'uint8');
hdr.nbyteschanhdr = fread(fid, 1, 'uint8');
hdr.nbyteshdr     = fread(fid, 1, 'int16');
hdr.nbytesdat     = fread(fid, 1, 'uint32'); %FIXME only for unpacked files
hdr.nbytestrailer1 = fread(fid, 1, 'uint32');
hdr.nbytestrailer2 = fread(fid, 1, 'uint16');
i9  = fread(fid, 1, 'uint16'); % not used
i10 = fread(fid, 1, 'uint16'); % not used
i11 = fread(fid, 1, 'int16');
i12 = fread(fid, 4, 'int8');
i13 = fread(fid, 1, 'double');

% FIXME there seems to be a discrepancy here, according to the
% documentation it should be the next (commented) line. according to the
% file I have at hand it should be the uncommented line.
%hdr.fsample        = double(hdr.nchan)./i13;
hdr.fsample = 1./i13;
i14 = fread(fid, 1, 'uint32');
i15 = fread(fid, 1, 'uint32');
i16 = fread(fid, 1, 'uint32');
i17 = fread(fid, 1, 'uint32');
i18 = fread(fid, 1, 'uint32');
i19 = fread(fid, 2, 'uint16');
i20 = fread(fid, 1, 'int16');
i21 = fread(fid, 1, 'int16');
i22 = fread(fid, 1, 'int8');
i23 = fread(fid, 1, 'int8');
i24 = fread(fid, 1, 'int8');
i25 = fread(fid, 1, 'int8');
i26 = fread(fid, 32, 'int8');
i27 = fread(fid, 1, 'uint16');
i28 = fread(fid, 1, 'uint8');
empty = fread(fid, 1, 'uint8');
i29 = fread(fid, 1, 'int8');
i30 = fread(fid, 1, 'int8');
i31 = fread(fid, 1, 'int8');
empty = fread(fid, 1, 'uint8');
i32 = fread(fid, 1, 'int8');
i33 = fread(fid, 1, 'int8');

% interpret relevant stuff from i27
% other stuff may or may not be relevant
hdr.ispacked  = bitand(i27, 2^15)>0;
hdr.lowchan   = bitand(i27, 2^9) >0;

% get channel info
for k = 1:(hdr.nbyteshdr-hdr.offset-2)/hdr.nbyteschanhdr
  fseek(fid, hdr.offset + (k-1)*hdr.nbyteschanhdr, 'bof');
  chanhdr(k).scalingslope     = fread(fid, 1, 'float');
  chanhdr(k).scalingintercept = fread(fid, 1, 'float');
  chanhdr(k).scale            = fread(fid, 1, 'double');
  chanhdr(k).intercept        = fread(fid, 1, 'double');
  chanhdr(k).unit             = fread(fid, 1, 'uint8=>char');
  fread(fid, 1, 'uint8');
  fread(fid, 1, 'uint8'); % may be relevant for packed files
  chanhdr(k).channr = fread(fid, 1, 'uint8');
  chanhdr(k).label  = num2str(k);
end
hdr.chanhdr = chanhdr;

fseek(fid, hdr.nbyteshdr + hdr.nbytesdat, 'bof');
x = fread(fid, hdr.nbytestrailer1, 'uint8');
hdr.x = x; % what shall we do with this?

fseek(fid, hdr.nbyteshdr + hdr.nbytesdat + hdr.nbytestrailer1, 'bof');
annot = fread(fid, hdr.nbytestrailer2, 'uint8');
sel   = find(annot==0);
indx  = find(annot>0, 1, 'first');
sel   = [0;sel];
for k = 1:numel(sel)-1
  hdr.chanhdr(k).annot = char(annot(sel(k)+1:sel(k+1)-1))';
end

fclose(fid);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
classdef TestPlotLargeCylinderTethaedra < handle

    properties (Access = private)
        backgroundMesh
        boundaryMesh
        levelSet
        unfittedMesh
    end

    methods (Access = public)

        function obj = TestPlotLargeCylinderTethaedra()
            obj.createBackgroundMesh();
            obj.createBoundaryMesh();
            obj.createLevelSet();
            obj.createUnfittedMesh();
            obj.plotUnfittedMesh();
        end

    end

    methods (Access = private)

        function createBackgroundMesh(obj)
            x = linspace(0,1,10);
            y = linspace(0,1,10);
            z = linspace(0,2,20);
            [X,Y,Z] = meshgrid(x,y,z);
            coord  = [X(:) Y(:) Z(:)];
            d = delaunayTriangulation(coord);
            s.connec = d.ConnectivityList;
            s.coord  = coord;
            obj.backgroundMesh = Mesh(s);
        end

        function createBoundaryMesh(obj)
            s.backgroundMesh = obj.backgroundMesh;
            s.dimension = 1:obj.backgroundMesh.ndim;
            bC = BoundaryMeshCreatorFromRectangularBox(s);
            bM = bC.create();
            obj.boundaryMesh = bM;
        end

        function createLevelSet(obj)
            s.type = 'cylinder';
            s.fracRadius = 1.1;
            s.coord      = obj.backgroundMesh.coord;
            s.ndim       = obj.backgroundMesh.ndim;
            lsCreator = LevelSetCreator.create(s);
            obj.levelSet = lsCreator.getValue();
        end

        function createUnfittedMesh(obj)
            s.backgroundMesh = obj.backgroundMesh;
            s.boundaryMesh   = obj.boundaryMesh;
            uM = UnfittedMesh(s);
            uM.compute(obj.levelSet);
            obj.unfittedMesh = uM;
        end

        function plotUnfittedMesh(obj)
            figure();
            obj.unfittedMesh.plot();
            view([1 1 1])
        end

    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [ theta,f1,m,W_new,beta,grt ] = armijo_theta_fullcycle_v3( Ltheta,dir,f0,phi0, grad,grt,W,K,M,Pt,omega,Hd,Hr,G)
    m=0;
    rhom=0.95;
    rho0=1/Ltheta*100;
    len=-real(grad'*dir);
    sig=0.4;
    while(1)
        rho=rho0*rhom^m;
        phi=phi0-rho*grad;
        x=exp(1j.*phi);
        [f1,W_new,beta,grt] = par_fun_B(x,W,K,M,Pt,omega,Hd,Hr,G);
        if (f1-f0)>=sig*rho*len
            break
        end
        if (rho)<1/Ltheta/10
            break
        end
        m=m+1;
    end
    theta=x;
end

author: guohuayan


text:
function [D, j, key] = LMquery(D, fieldName, content, method)
%
% This function provides basic query functionalities. It allows you to
% locate images with certain objects.
%
% [Dout, j, key] = LMquery(D, fieldName, content, method)
%
%    D:         database struct
%    fieldName: can be any of the fields in the struct at the object level.
%    content:   a string describing the desired field you are looking for.
%               you can use the characters '+' and '-' within the content string in
%               order to perform 'and' and 'not' operations.
%    method:    when selected 'exact', the content must be match exactly
%             Methods:
%               [] = default = substring matching (air => chair)
%               'exact' = strings should match exactly
%               'word'  = the field should match words (air ~=> chair)
%
% Output:
%   Dout: new struct matching the query
%   j: indices to images that fit the query. Note Dout is not D(j). D(j)
%   contains all the objects in the original struct
%   key: pointers to images and objects in Dout with respect to D.
%       The first image in Dout is:
%         D(key(1).jimg).annotation.object(key(1).jobj)
%      [key(:).jimg] is the same as j
%
%
% The function LMquery searchs for images in the database.
% Queries can be done on any field of the struct array.
% The fields for an image are:
%    filename
%    folder
%    source.sourceImage
%    source.sourceAnnotation
%    object(:).name
%    object(:).deleted
%    object(:).verified
%    object(:).date
%    object(:).polygon.username
%    object(:).polygon.pt(:).x
%    object(:).polygon.pt(:).y
%    object(:).viewpoint.azimuth
%
% Look at the results of these two lines (HOMEIMAGES should be set with the path of the database):
% LMdbshowobjects(LMquery(D, 'object.name', 'plate+license'), HOMEIMAGES);
% LMdbshowobjects(LMquery(D, 'object.name', 'plate-license'), HOMEIMAGES);
%
% Performs an OR function. This will select objects that belong to one of this groups: 1) side views
% of cars, 2) buildings 3) roads or 4) trees.
% [D,j] = LMquery(D, 'object.name', 'car+side,building,road,tree');
%
% To make AND queries of different objects within an image, you can chain several queries.
% For instance, to get images that have cars, buildings and roads:
%    [Dfoo,j1] = LMquery(D, 'object.name', 'building');
%    [Dfoo,j2] = LMquery(D, 'object.name', 'road');
%    j = intersect(j1,j2);
%    LMshow(LMquery(D(j), 'object.name', 'car,building,road'));
%
% You can query other fields. The next lines allow you to select from the
% database, all the polygons that have been verified and that are not
% marked as deleted. This is something you might want to do often:
%    D = LMquery(D, 'object.verified', '1');
%    D = LMquery(D, 'object.deleted', '0');
%
% or D = LMquery(LMquery(D, 'object.verified', '1'), 'object.deleted', '0');
%
% You can also query based on the date of annotation:
%    LMshow(LMquery(database, 'date', '13-Jan'));
% [in the case of dates, D is an exception about the meaning of the
% character '-'. Here is just a separation between day and month].
%
% More examples:
% Queries for objects
% [Dj,j] = LMquery(D, 'object.name', 'building');
% LMdbshowscenes(D(j), HOMEIMAGES); % this shows all the objects in the images that contain buildings
% LMdbshowscenes(Dj, HOMEIMAGES); % this shows only the buildings
%
% Queries for images in specific folders.
% [Dj,j] = LMquery(D, 'folder', '05june05_static_indoor');
% LMdbshowscenes(D(j), HOMEIMAGES); % this shows all the objects
%
% The next example shows the annotated objects from images that come from the web:
% [Dj,j] = LMquery(D, 'folder', 'web');
% LMdbshowscenes(D(j), HOMEIMAGES); % this shows all the objects
%
% look for a specific image file:
% [Dj,j] = LMquery(D, 'filename', 'p1010843.jpg');
% LMdbshowscenes(D(j), HOMEIMAGES); % this shows all the objects
%
% look for objects annotated by one user:
% [Dj,j] = LMquery(D, 'object.polygon.username', 'atb');
% LMdbshowscenes(Dj, HOMEIMAGES); % this shows all the objects annotated by one user
%
% look for objects by viewpoint:
% [D,j] = LMquery(D, 'object.viewpoint.azimuth', '0', 'exact');
%

if (nargin < 4)
    method = '';
end

% Parse query
%if length(strfind(fieldName, 'date'))>0
%    query{1} = {['+' lower(content)]};
%else
%    query = parseContent(content);
%end

% Parse query
if ~isempty(strfind(fieldName, 'date')) || strcmp(method,'exact')
    query{1} = {['+' lower(content)]};
%elseif strcmp(method,'exact')
%    query{1} = {['+' lower(content)]};
else
    query = parseContent(content);
end

% Parse queried field
fieldName = ['annotation.' fieldName];
fieldName = parseContent(fieldName, '.');

% Locate the images that verify the query
Nimages = length(D);
queriedImages = logical(zeros(Nimages,1));

k=0;
for i = 1:Nimages
    v = D(i);

    [vn, jc] = queryfields(v, fieldName, query, method);

    if ~isempty(vn)
        queriedImages(i)=1;
        D(i) = vn;

        k = k + 1;
        key(k).jimg = i;
        key(k).jobj = jc;
    end
end

disp(sprintf('%d matches out of %d', sum(queriedImages), length(D)))

D = D(queriedImages);
j = find(queriedImages);

function [vn, jc] = queryfields(v, fieldName, query, method)

vn = [];
jc = [];
depthfield = length(fieldName);

if depthfield == 1
    % End recursion
    if isfield(v, fieldName{1}{1}(2:end))
        jc = findobject({v.(fieldName{1}{1}(2:end))}, query, method);
        vn = v(jc);
    else
        if strcmp(query{1}(1), '-')
            jc = 1:length(v);
            vn = v;
        end
    end
else
    if isfield(v, fieldName{1}{1}(2:end))
        % Recursive call
        N = length(v);
        q = logical(zeros(N,1));
        for i = 1:N
            vi = getfield(v, {i}, fieldName{1}{1}(2:end));
            [vi, jc] = queryfields(vi, fieldName(2:end), query, method);
            q(i) = ~isempty(vi);
            v(i).(fieldName{1}{1}(2:end)) = vi;
        end
        jq = find(q);
        vn = v(jq);
    else
        if strcmp(query{1}(1), '-')
            jc = 1:length(v);
            vn = v;
        end
    end
end

author: CSAILVision

repo: LabelMeToolbox


text:
% compile mex file
mex mex/perform_front_propagation_2d.cpp mex/perform_front_propagation_2d_mex.cpp mex/fheap/fib.cpp
mex mex/perform_front_propagation_3d.cpp mex/perform_front_propagation_3d_mex.cpp  mex/fheap/fib.cpp
mex mex/perform_circular_front_propagation_2d.cpp mex/perform_front_propagation_2d.cpp mex/fheap/fib.cpp

% Alec: these do not actually exist
% % anisotropic FM
% mex mex/anisotropic-fm//perform_front_propagation_anisotropic.cpp
% mex mex/anisotropic-fm-feth/fm2dAniso.cpp

% compiling skeleton
mex mex/skeleton.cpp
% compiling distance transform
mex mex/eucdist2.c

disp('Compiling perform_front_propagation_mesh, might time some time.');
rep = 'mex/';
files =  { ...
    'perform_front_propagation_mesh.cpp', ...
    'gw/gw_core/GW_Config.cpp',           ...
    'gw/gw_core/GW_FaceIterator.cpp',     ...
    'gw/gw_core/GW_SmartCounter.cpp',     ...
    'gw/gw_core/GW_VertexIterator.cpp',   ...
    'gw/gw_core/GW_Face.cpp',             ...
    'gw/gw_core/GW_Mesh.cpp',             ...
    'gw/gw_core/GW_Vertex.cpp',       ...
    'gw/gw_geodesic/GW_GeodesicFace.cpp', ...
    'gw/gw_geodesic/GW_GeodesicMesh.cpp',     ...
    'gw/gw_geodesic/GW_GeodesicPath.cpp',         ...
    'gw/gw_geodesic/GW_GeodesicPoint.cpp',            ...
    'gw/gw_geodesic/GW_TriangularInterpolation_Cubic.cpp', ...
    'gw/gw_geodesic/GW_GeodesicVertex.cpp',                    ...
    'gw/gw_geodesic/GW_TriangularInterpolation_Linear.cpp',      ...
    'gw/gw_geodesic/GW_TriangularInterpolation_Quadratic.cpp',  ...
};
str = 'mex ';
for i=1:length(files)
    str = [str rep files{i} ' '];
end
eval(str);

author: alecjacobson

repo: gptoolbox


text:
function [bboxLeft3M,bboxRight3M] = getBreastMask(scan3M)
% AI 1/16/18

%Get no. of slices,rows,  cols
nSlices = size(scan3M,3);
nCol = size(scan3M,2);
nRow = size(scan3M,1);
airTh = 400;
tol = 0;
filt = gausswin(13,1.5);
filt = filt./sum(filt(:));

%Loop over slices
bboxLeft3M = false(size(scan3M));
bboxRight3M = false(size(scan3M));
fVoxV = nan(1,nSlices);
bsizeV = nan(1,nSlices);

for i = 1: nSlices

    sliceM = scan3M(:,:,i);

    %Get index of last voxel in each column with inensity > threshold
    th_maskM = sliceM>airTh;
    th_maskM = bwareaopen(th_maskM,50);

    fVoxV(i) = nnz(th_maskM(:))/numel(th_maskM(:));
    rowIdxV = nan(nCol,1);
    for j = 1 : nCol
        idx = find(th_maskM(:,j),1,'last');
        if ~isempty(idx)
            rowIdxV(j) = idx;
        end
    end
    colIdxV = [1:nCol].';
    % startCol = find(~isnan(rowIdxV),1);
    % endCol = find(~isnan(rowIdxV),1,'last');

    try

        %Smooth outline
        outlineV = smooth(colIdxV,rowIdxV,'rlowess');

        %Find peaks
        outSmoothV = conv(outlineV,filt,'same');
        [~,boundsV] = findpeaks(outSmoothV,'MinPeakProminence',10,'SortStr','descend');
        boundsV = sort(boundsV(1:2));

        %AI 2/10/20
        %Added to handle cases where one side is not visible
        if abs(outSmoothV(boundsV(2))-outSmoothV(boundsV(1)))>100
            %skip slice
            emptyMaskM = false(size(th_maskM));
            bboxLeft3M(:,:,i) = emptyMaskM;
            bboxRight3M(:,:,i) = emptyMaskM;
            fVoxV(i) = 0;
            bsizeV(i)=0;
        else
            %Find valley between peaks
            flipSigV = nRow-outlineV(boundsV(1):boundsV(2));
            [~,locV] = findpeaks(flipSigV,'MinPeakProminence',30,'SortStr','descend');
            if isempty(locV)
                [~,locV] = findpeaks(flipSigV,'SortStr','descend');
            end
            loc = locV(1) + boundsV(1);
            rowLim = outlineV(loc);
            rowLim = rowLim-tol;

            % --- For testing ----
            %Display
            %   figure,scatter(colIdxV,rowIdxV,'bo'),axis([1 nCol 1 nRow]);
            %   hold on
            %   plot(colIdxV,outlineV,'--r')
            %---------------------

            outlineV(outlineV<rowLim) = rowLim;
            bsizeV(i) = max(outlineV)- rowLim;

            %Create mask
            xV = [flipud(colIdxV);colIdxV];
            yV = [repmat(rowLim,nCol,1);outlineV];
            maskM = poly2mask(xV,yV,nRow,nCol) & th_maskM;
            maskM = bwareaopen(maskM,100);
            leftMaskM = maskM;
            leftMaskM(:,loc+1:end) = 0;
            rightMaskM = maskM;
            rightMaskM(:,1:loc) = 0;

            %Morphological processing
            leftMaskM = imfill(leftMaskM,'holes');
            leftMaskM = imclose(leftMaskM,strel('disk',4));
            rightMaskM = imfill(rightMaskM,'holes');
            rightMaskM = imclose(rightMaskM,strel('disk',4));

            bboxLeft3M(:,:,i) = leftMaskM;
            bboxRight3M(:,:,i) = rightMaskM;
        end

    catch

        fVoxV(i) = 0;
        bsizeV(i) = 0;

    end
end

%Plot size across slices
sliceV = 1:length(bsizeV);
% Smooth to remove sudden spikes
filtV = gausswin(64,10.5);
filtV = filtV./sum(filtV);
smooth_bsizeV = medfilt1(bsizeV,10);
smooth_bsizeV = conv(smooth_bsizeV,filtV,'same');
% Get dir of change in size
diffV = diff(smooth_bsizeV);
signV = [0 sign(diffV)];
signV = movmean(signV,10);             %Smooth

%Identify slices to retain
startsV = find([1,diff(signV)]);       %Starts of runs
runsV = diff(find([1,diff(signV),1])); %Lengths of runs
[runsV,orderV] = sort(runsV,'descend');
startsV = startsV(orderV);

% Find start slice
count = 0;
incFlag = 0;
while ~incFlag
    count = count+1;
    start1Idx = startsV(count);
    incFlag = signV(start1Idx)>0;
end
end1Idx = start1Idx+runsV(count)-1;
% Find end slice
run2IdxV = startsV>end1Idx;
startsV = startsV(run2IdxV);
runsV = runsV(run2IdxV);
count = 0;
decFlag = 0;
while ~decFlag
    count = count+1;
    start2Idx = startsV(count);
    decFlag = signV(start2Idx)<0;
end
end2Idx = start2Idx+runsV(count)-1;
minSizIdx = find(smooth_bsizeV(start2Idx:end2Idx)>= smooth_bsizeV(start1Idx)*.8,1,'last');
end2Idx = start2Idx + minSizIdx - 1;
keepSlicesV = start1Idx:end2Idx;
skipIdxV = ~ismember(sliceV,keepSlicesV);

% --- For testing (plot of bsizeV)----
% h = figure('visible','On');
% plot(sliceV,bsizeV,'-kx','lineWidth',1);
% axis([1 numel(sliceV),0,max(smooth_bsizeV)])
% hold on
% plot(sliceV,smooth_bsizeV,'--bx','lineWidth',1);
% line([start1Idx,start1Idx],[0,max(smooth_bsizeV)],'Color','r','LineStyle','--');
% line([end2Idx,end2Idx],[0,max(smooth_bsizeV)],'Color','r','LineStyle','--');
% xlabel('Size');
% ylabel('Slice #');
%----------------------------------

%Skip selected slices
bboxLeft3M(:,:,skipIdxV) = false(size(scan3M,1),size(scan3M,2),nnz(skipIdxV));
bboxRight3M(:,:,skipIdxV) = false(size(scan3M,1),size(scan3M,2),nnz(skipIdxV));

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [E tempprot] = RemoveNegMeas(Epn, protocol)
% Removes negative or zero measurements from a set and returns a protocol
% with the corresponding elements removed to exclude the negative
% measurements from the fitting.
%
% Epn is the full set of measurements; E is that with the negative ones
% removed.
%
% protocol is the full protocol; tempprot is that with the entries
% corresponding to negative measurements removed.
%
% author: Daniel C Alexander (d.alexander@ucl.ac.uk)
%         Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

tempprot = protocol;
E = Epn;
if(min(Epn)<=0)
    posonly = find(Epn>0);
    E = Epn(posonly);
    nonposonly = setdiff([1:size(Epn,1)], posonly);
    tempprot.totalmeas = length(posonly);
    % find the b0 indices for the new protocol without the non-positive
    % measurements
    tempprot.b0_Indices = [];
    for i=1:length(protocol.b0_Indices)
        currentB = protocol.b0_Indices(i);
        index = find(nonposonly<=currentB, 1, 'last');
        if isempty(index)
            tempprot.b0_Indices = [tempprot.b0_Indices currentB];
        elseif (nonposonly(index) ~= currentB)
            currentB = currentB - index;
            tempprot.b0_Indices = [tempprot.b0_Indices currentB];
        end
    end
    tempprot.numZeros = length(tempprot.b0_Indices);
    if(strcmp(tempprot.pulseseq, 'PGSE') || strcmp(tempprot.pulseseq, 'STEAM'))
        tempprot.G = tempprot.G(posonly);
        tempprot.delta = tempprot.delta(posonly);
        tempprot.smalldel = tempprot.smalldel(posonly);
        tempprot.grad_dirs = tempprot.grad_dirs(posonly, :);
    else
        error('Need to adapt for other pulse sequences.');
    end
end

if length(tempprot.b0_Indices) == 0
    error('All b=0 measurements are negative');
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function varargout = drawEllipseAxes(varargin)
%DRAWELLIPSEAXES Draw the main axes of an ellipse as line segments.
%
%   drawEllipseAxes(ELLI)
%   drawEllipseAxes(..., STYLE)
%   drawEllipseAxes(..., NAME, VALUE)
%   Draw the axes of the ellipse given by ELLI onto the currrent axis.
%   STYLE specifies the drawing style using a short character array like
%   'b', 'k:', 'm-'...
%   More complex drawing style can be specified using plot-like parameter
%   name-value pairs.
%
%   drawEllipseAxes(AX, ELLI)
%   Specifies the axes to draw the ellipse on.
%
%   Example
%     elli = [50 50  40 20  30];
%     figure; hold on; axis equal; axis([0 100 0 100]);
%     drawEllipse(elli, 'LineWidth', 2, 'Color', 'b')
%     drawEllipseAxes(elli, 'k')
%
%   See also
%     ellipses2d, drawEllipse, drawEllipseArc
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2022-09-11, using Matlab 9.9.0.1570001 (R2020b) Update 4
% Copyright 2022 INRAE - BIA Research Unit - BIBS Platform (Nantes)

%% Extract input arguments

% extract handle of axis to draw on
if isAxisHandle(varargin{1})
    ax = varargin{1};
    varargin(1) = [];
else
    ax = gca;
end

% extract dawing style strings
styles = {};
for iElli = 1:length(varargin)
    if ischar(varargin{iElli})
        styles = varargin(iElli:end);
        varargin(iElli:end) = [];
        break;
    end
end

% retrieve ellipse parameters
ellipse = varargin{1};
x0 = ellipse(:, 1);
y0 = ellipse(:, 2);
a  = ellipse(:, 3);
b  = ellipse(:, 4);
theta = ellipse(:, 5);
nElli = length(x0);

%% Process drawing of a set of ellipses

% angular positions of edge extremities
ti = [0 pi  pi/2 3*pi/2];

% compute position of points to draw each ellipse
h = zeros(2 * nElli , 1);
for iElli = 1:nElli
    % pre-compute rotation angles (given in degrees)
    cot = cosd(theta(iElli));
    sit = sind(theta(iElli));

    % compute position of points used to draw current ellipse
    xt = x0(iElli) + a(iElli) * cos(ti) * cot - b(iElli) * sin(ti) * sit;
    yt = y0(iElli) + a(iElli) * cos(ti) * sit + b(iElli) * sin(ti) * cot;

    % stores handle to graphic object
    h(2 * iElli - 1) = plot(ax, xt(1:2), yt(1:2), styles{:});
    h(2 * iElli)     = plot(ax, xt(3:4), yt(3:4), styles{:});
end

% return handles if required
if nargout > 0
    varargout = {h};
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function model = modelSaturation(model,varargin)
% function model = modelSaturation(model,[threshold],[soft])
%
% replaces all values above some threshold with a specified value
% to model saturation effects
% default = 4 (4 x max HRF)
%
% a 3rd argument signals a soft thresholding, where
% the argument is the exponent (ex) in a thresholding model:
% y(y > thresh) = thresh - (thresh .^ ex) + (y(y > thresh) .^ ex);
%
% 3/30/01 Tor Wager
%

if nargin > 1, thresh = varargin{1};,else thresh = 4;,end

if nargin < 3
    model(model(:,:) > thresh) = thresh;
else
    if isempty(varargin{2}), ex = .5;, else, ex = varargin{2};,end
    model(model>thresh) = 2 - (2 .^ ex) + (model(model>thresh) .^ ex);
end

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function dat = getDateAndTime
% fucniton dat = getDateAndTime
% Simple fucntion that formats the date and time for use in file names.
% Example: dat = getDateAndTime;
%          dat > 20-Oct-2010_10h41m42s
%
% LMP - 20-Oct-2010
dateAndTime     = datestr(now);
dateAndTime(12) = '_';
dateAndTime(15) = 'h';
dateAndTime(18) = 'm';
dateAndTime(21) = 's';
dat = dateAndTime;
return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function tests = test_spm_create_vol
% Unit Tests for spm_create_vol
%__________________________________________________________________________
% Copyright (C) 2016 Wellcome Trust Centre for Neuroimaging

% $Id: test_spm_create_vol.m 6908 2016-10-21 11:29:25Z guillaume $

tests = functiontests(localfunctions);

function test_spm_create_vol_hdr(testCase)
V = struct(...
    'fname',   [tempname '.img'],...
    'dim',     [32 32 32],...
    'dt',      [spm_type('uint8') spm_platform('bigend')],...
    'mat',     eye(4),...
    'pinfo',   [1 0 0]',...
    'descrip', 'description');

V = spm_create_vol(V);
hdrfile = spm_file(V.fname,'ext','.hdr');
testCase.verifyEqual(exist(hdrfile,'file'),2);
spm_unlink(hdrfile);

function test_spm_create_vol_nii(testCase)
V = struct(...
    'fname',   [tempname '.nii'],...
    'dim',     [64 64 64],...
    'dt',      [spm_type('uint16') spm_platform('bigend')],...
    'mat',     eye(4),...
    'pinfo',   [1 0 0]',...
    'descrip', 'description');

V = spm_create_vol(V);
testCase.verifyEqual(exist(V.fname,'file'),2);
spm_unlink(V.fname);

function test_spm_create_vol_nii_min(testCase)
V = struct(...
    'fname', [tempname '.nii'],...
    'dim',   [64 64 64],...
    'mat',   eye(4));

V = spm_create_vol(V);
testCase.verifyEqual(exist(V.fname,'file'),2);
spm_unlink(V.fname);

function test_spm_create_vol_nii_multi(testCase)
V(1:4) = deal(struct(...
    'fname', '',...
    'dim',   [64 64 64],...
    'mat',   eye(4)));
for i=1:numel(V), V(i).fname = [tempname '.nii']; end

V = spm_create_vol(V);
for i=1:numel(V)
    testCase.verifyEqual(exist(V(i).fname,'file'),2);
    spm_unlink(V(i).fname);
end

author: spm


text:
function fnames = findfiles(extension,directory,recurse)
% Finds all files with the specified extension in the current directory and
% subdirectories (recursive). Returns a cell array with the fully specified
% file names.
%
%    files = findfiles( ext )
%       searches in the current directory and subdirectories.
%
%    files = findfiles( ext, directory)
%       starts search in directory specified.
%
%    files = findfiles( ext, directory, searchSubdirectories )
%       If searchSubdirectories == 0, the function will not descend
%       subdirectories. The default is that searchSubdirectories ~= 0.
%
% Note: if this function causes your recursion limit to be exceeded, then
% it is most probably trying to follow a symbolic link to a directory that
% has a symbolic link back to the directory it came from.
%
% Examples:
%
%    bmpfiles = findfiles('bmp')
%    dllfiles = findfiles('dll','C:\WinNT')
%    mfiles = findfiles('m',matlabroot)

% Copyright (C) 2001 Quester Tangent Corporation
% Author: Tony Christney, tchristney@questertangent.com
% $Id: findfiles.m 1.5 2001/05/25 01:54:17 tchristney Exp $

% $Log: findfiles.m $
% Revision 1.5  2001/05/25 01:54:17  tchristney
% Fixed bug where numMatches was not updated after recursive call, causing
% potential overwrites in the fnames cell array.
%
% Revision 1.4  2001/05/24 20:46:09  tchristney
% Recursion on/off support courtesy Brett Shoelson.
% <Brett.Shoelson@joslin.harvard.edu>
%
% Revision 1.2  2001/03/16 21:28:37  tchristney
% added some comments for external release, i.e. examples.
%
% Revision 1.1  2001/03/07 22:38:52  tchristney
% Initial revision
%

if nargin == 1
    directory = cd;
    recurse = (1==1);
elseif nargin == 2
    oldDir = cd;
    cd(directory);
    directory = cd;
    cd(oldDir);
    recurse = (1==1);
elseif nargin >= 3
    oldDir = cd;
    cd(directory);
    directory = cd;
    cd(oldDir);
end

d = dir(directory);

fnames = {};
numMatches = 0;
for i=1:length(d)
    % look for occurences of ['.' extension] in the file name
    extIndices = findstr(['.' extension],d(i).name);

    % if the file is not a directory, and the file has at least one occurence
    if ~d(i).isdir & ~isempty(extIndices)

        % then if the last occurrence is at the end of the file name,
        % add the file name to the list.
        if length(d(i).name) == (extIndices(length(extIndices)) + length(extension))
            numMatches = numMatches + 1;
            fnames{numMatches} = fullfile(directory,d(i).name);
        end
    % otherwise, descend directories appropriately.
    % note that this could result in a recursion limit error if it tries to
    % follow symbolic links that loop back on themselves...
    elseif recurse & d(i).isdir & ~strcmp(d(i).name,'.') & ~strcmp(d(i).name,'..')
        fnames = [fnames findfiles(extension,fullfile(directory,d(i).name),recurse)];
        numMatches = length(fnames);
    end
end

author: Sable


text:
function initializeTrackTypePanel( horizontalBar, leftBottom, panelSize, hdls, fontSize )
%INITIALIZETRACKTYPEPANEL Initialize track type panel with a dropdown list of different track type

track_type_readable_list = { 'All', 'None'};

is_none = @(tag) 0;
is_all = @(tag) 1;

filter_functions = {is_all, is_none};

trackTypePanel = uipanel(horizontalBar, ...
    'Units', 'normalized', ...
    'Position', [leftBottom panelSize], ...
    'BorderType', 'None');

uicontrol(trackTypePanel, ...
    'String', 'Track type', ...
    'Style', 'text',...
    'Units', 'normalized', ...
    'FontSize', fontSize, ...
    'Position', [0 .5 1 .5]);

uicontrol(trackTypePanel, ...
    'String', track_type_readable_list, ...
    'Style', 'popup',...
    'Units', 'normalized', ...
    'FontSize', fontSize, ...
    'Position', [0 0 1 .5], ...
    'Callback', {@trackTypeBoxCallback, hdls, filter_functions});

hdls.track_type_filter_function = is_all;
end

function trackTypeBoxCallback(hObject, ~, hdls, filter_functions)

track_type_index = get(hObject, 'Value');
switch track_type_index
    case 1 % ALL
        hdls.is_plotting_all = true;
        hdls.is_plotting_none = false;
    case 2 % NONE
        hdls.is_plotting_none = true;
        hdls.is_plotting_all = false;
    otherwise
        hdls.track_type_filter_function = filter_functions{track_type_index};
        hdls.is_plotting_all = false;
        hdls.is_plotting_none = false;
end

hdls.generate_data();

end

author: jfaghm

repo: OceanEddies

save_path: github-repos/MATLAB/jfaghm-OceanEddies


text:
function file_saveMatlab(file, dat, mrk, mnt, varargin)
%FILE_SAVEMATLAB - Save EEG data structures in Matlab format
%
% Synopsis:
%   file_saveMatlab(FILE, DAT, MRK, MNT, <OPT>)
%
% Arguments:
%   FILE:   CHAR    name of data file
%   DAT:    STRUCT  structure of continuous or epoched signals
%   MRK:    STRUCT  marker structure
%   MNT:    STRUCT  electrode montage structure
%
%   OPT: PROPLIST - Structure or property/value list of optional properties:
%   'Folder' - CHAR: Path to save the file. Default is the global variable BTB.MatDir
%                  unless FILE is an absolute path in which case it is ''.
%   'Channelwise' - BOOL: If true, signals are saved channelwise. This is an advantage
%                  for big files, because it allows to load selected
%                  channels. (default 1)
%   'Format' - CHAR: 'double', 'float', 'int16', or 'auto' (default).
%             In 'auto' mode, the function tries to find a lossless conversion
%             of the signals to INT16 (see property '.ResolutionList').
%             If this is possible '.format' is set to 'INT16', otherwise it is
%             set to 'DOUBLE'.
%   'Resolution' - CHAR|INT : Resolution of signals, when saving in format INT16.
%                 (Signals are divided by this factor before saving.) The resolution
%                 maybe selected for each channel individually, or globally for all
%                 channels. In the 'auto' mode, the function tries to find for each
%                 channel a lossless conversion to INT16 (see property
%                 'ResolutionList'). For all other channels the resolution producing
%                 least information loss is chosen (under the resolutions that avoid
%                 clipping). Possible values:
%                 'auto' (default),
%                 numerical scalar, or
%                 numerical vector of length 'number of channels'
%                 (i.e., length(DAT.clab)).
%   'ResolutionList' - DOUBLE: Vector of numerical values. These values are tested as
%                      resolutions to see whether lossless conversion to INT16 is
%                      possible. Default [1 0.5 0.1].
%   'Accuracy' - DOUBLE: used to define "losslessness" (default 10e-10)
%   'AddChannels' - BOOL: (true or false) Adds the channels in DAT to the
%                   existing MAT file (default 0)
%   'Vars': Additional variables that should be stored. 'opt.Vars' must be a
%           cell array of variable names, e.g., {'hdr', 'patient_info'}
%   'SaveParam': additonal parameters to be passes to the Matlab save function,
%           default {'-v7'}
%
%Description:
%   Saves data, marker, montage and (if requested) additional variables in
%   BTB format in a single matlab file.
%
% See also: file_loadBV file_loadMatlab

% Author(s): Benjamin Blankertz, Feb 2005

global BTB

props = {'Folder'            BTB.MatDir   'CHAR'
         'Channelwise'       1            '!BOOL'
         'Format'            'auto'       '!CHAR(double float int16 auto)'
         'Resolution'        'auto'       '!CHAR(auto)|!DOUBLE[1]|!DOUBLE[-]'
         'ResolutionList'   [1 0.5 0.1]   '!DOUBLE[-]'
         'Accuracy'          10e-10       '!DOUBLE[1]'
         'AddChannels'       0            '!BOOL'
         'Vars'              {}           'CHAR|CELL{CHAR}'
         'SaveParam'         {'-v7'}      'CELL'
        };

misc_checkType(dat,'STRUCT(x fs clab)');
misc_checkType(mrk,'STRUCT(time y)');
misc_checkType(mnt,'STRUCT(x y clab)');
opt= opt_proplistToStruct(varargin{:});
[opt, isdefault]= opt_setDefaults(opt, props);
opt_checkProplist(opt, props);

if ischar(opt.Vars),
  opt.Vars= {opt.Vars};
end
if fileutil_isAbsolutePath(file),
  opt.Folder = '';
end

nChans= length(dat.clab);
fullname= fullfile(opt.Folder, file);
if opt.AddChannels && ~exist([fullname '.mat'], 'file')
  warning('File does not exist: ignoring option ''AddChannels''.');
  opt.AddChannels= 0;
end

if opt.AddChannels,
  if ~opt.Channelwise,
    warning('AddChannels requested: forcing channelwise mode');
    opt.Channelwise= 1;
  end
  nfo_old= file_loadMatlab(fullname, 'Vars','nfo');
  if ~isdefault.format && opt.Format~=nfo_old.format,
    warning(sprintf('format mismatch: using %s', nfo_old.format));
  end
  opt.format= nfo_old.format;
end

if any(strcmpi(opt.Format, {'INT16','AUTO'}))
  %% Look for resolutions producing lossless conversion to INT16.
  if strcmpi(opt.Resolution, 'auto'),
    opt.Resolution= NaN*ones(1, nChans);
    for cc= 1:nChans,
      rr= 0;
      res_found= 0;
      while ~res_found && rr<length(opt.ResolutionList),
        rr= rr+1;
        X= dat.x(:,cc,:) / opt.ResolutionList(rr);
        X= X(:);
        if all( abs(X-round(X)) < opt.Accuracy ) && ...
              all(X>=-32768-opt.Accuracy) && all(X<=32767+opt.Accuracy),
          opt.Resolution(cc)= opt.ResolutionList(rr);
          res_found= 1;
        end
      end
      clear X;
    end
  end

  %% Expand global resolution.
  if length(opt.Resolution)==1,
    opt.Resolution= opt.Resolution*ones(1,nChans);
  end

  %% Check format of resolution.
  if ~all(isnumeric(opt.Resolution)) || length(opt.Resolution)~=nChans,
    error('property resolution has invalid format');
  end

  %% If for all channels lossless conversions were found, auto-select
  %% format INT16.
  if strcmpi(opt.Format, 'auto'),
    if all(~isnan(opt.Resolution)),
      opt.Format= 'INT16';
    else
      opt.Format= 'DOUBLE';
    end
  end
end

%% Check format of property format.
if ~ismember(upper(opt.Format), {'INT16','FLOAT','DOUBLE'},'legacy'),
  error('unknown format');
end

%% Select resolution for lossy conversion to INT16.
if strcmpi(opt.Format, 'INT16'),
  iChoose= find(isnan(opt.Resolution));
  for cc= 1:length(iChoose),
    ci= iChoose(cc);
    dat_ch= dat.x(:,ci,:);
    opt.Resolution(ci)= 1.000001*max(abs(dat_ch(:)))'/32767;
  end
end

%% Gather some summary information into structure 'nfo'.
nfo= struct_copyFields(dat,{'fs','clab'});
nfo.T= size(dat.x,1);
nfo.nEpochs= size(dat.x,3);
nfo.length= size(dat.x,1)*size(dat.x,3) / dat.fs;
nfo.format= opt.Format;
nfo.resolution= opt.Resolution;
nfo.file= fullname;
if isfield(mrk, 'time')
  nfo.nEvents= length(mrk.time);
else
  nfo.nEvents= 0;
end
if isfield(mrk, 'y'),
  nfo.nClasses= size(mrk.y,1);
else
  nfo.nClasses= 0;
end
if isfield(mrk, 'className')
  nfo.className= mrk.className;
else
  nfo.className= {};
end

%% if adding channels is requested, merge the nfo structures
if opt.AddChannels,
  nfo.clab= cat(2, nfo_old.clab, nfo.clab);
  nfo.resolution= cat(2, nfo_old.resolution, nfo.resolution);
end

%% Create directory if necessary
[filepath, filename]= fileparts(fullname);
if ~exist(filepath, 'dir'),
  [parentdir, newdir]=fileparts(filepath);
  [status,msg]= mkdir(parentdir, newdir);
  if status~=1,
    error(msg);
  end
  if isunix,
    unix(sprintf('chmod a-rwx,ug+rwx %s', filepath));
  end
end

saveParamNoVer= opt.SaveParam;
iVer= strmatch('-v', saveParamNoVer);
saveParamNoVer(iVer)= [];

if opt.AddChannels,
  %% update the nfo structure
  save(fullname, '-APPEND', 'nfo', saveParamNoVer{:});
  chan_offset= length(nfo_old.clab);
else
  save(fullname, 'mrk', 'mnt', 'nfo', opt.SaveParam{:});
  chan_offset= 0;
end

dat.file= fullname;
if opt.Channelwise,
  rhs= 'dat.x(:,cc,:)';
  switch(upper(opt.Format)),
   case 'INT16',
    evalstr= ['int16(round(' rhs '/opt.Resolution(cc)));'];
    dat.resolution= opt.Resolution;
   case 'FLOAT',
    evalstr= ['float(' rhs ');'];
   case 'DOUBLE',
    evalstr= [rhs ';'];
  end
  for cc= 1:nChans,
    varname= ['ch' int2str(chan_offset+cc)];
    eval([varname '= ' evalstr]);
    save(fullname, '-APPEND', varname, saveParamNoVer{:});
    clear(varname);
  end
  dat= rmfield(dat, 'x');
  %% the following field updates are needed for add_channels=1
  dat.clab= nfo.clab;
  if ~ismember(upper(opt.Format), {'FLOAT','DOUBLE'},'legacy'),
    dat.resolution= nfo.resolution;
  end
  save(fullname, '-APPEND', 'dat', saveParamNoVer{:});
else
  switch(upper(opt.Format)),
   case 'INT16',
    for cc= 1:nChans,
      dat.x(:,cc,:)= round( dat.x(:,cc,:) / opt.Resolution(cc) );
    end
    dat.x= int16(dat.x);
    dat.resolution= opt.Resolution;
   case 'FLOAT',
    dat.x= float(dat.x);
   case 'DOUBLE',
    %% nothing to do
  end
  save(fullname, '-APPEND', 'dat', saveParamNoVer{:});
end

%% Save additional variables, as requested.
if ~isempty(opt.Vars),
  vars= struct;
  for vv= 1:length(opt.Vars),
    vars.(opt.Vars{vv})= evalin('caller', opt.Vars{vv});
  end
  save(fullname, '-APPEND', '-STRUCT', 'vars', saveParamNoVer{:});
end

author: bbci


text:
function h = plotAngleDistribution(obj,varargin)
% plot axis distribution
%
% Syntax
%
%   plotAngleDistribution(mdf)
%   plotAngleDistribution(CS1,CS2)
%   plotAngleDistribution(grains.boundary.misorientation)
%
% Input
%  CS - @symmetry
%
% Options
%  resolution, xy*degree - resolution of the plots (given as angle)
%

[mtexFig,isNew] = newMtexFigure(varargin{:});
mtexFig.keepAspectRatio = false;

% compute angles
plotType = 'line';
if isa(obj,'symmetry')
  maxOmega = maxAngle(obj,varargin{:});
else
  maxOmega = maxAngle(obj.CS,obj.SS);
  if ~isa(obj,'SO3Fun'), plotType = 'bar'; end
end

% seach for existing bar plots and adjust bar center
h = findobj(mtexFig.gca,'type','bar','-or','type','hgGroup');
h = flipud(h(:));

unit = '%';
if ~isempty(h)

  midPoints = ensurecell(get(h,'XData'));
  midPoints= midPoints{1}*degree;
  bins = [2*midPoints(1)-midPoints(2),midPoints,2*midPoints(end)-midPoints(end-1)];
  bins = (bins(1:end-1) + bins(2:end))/2;
  density = ensurecell(get(h,'YData'));
  density = cellfun(@(x) x(:),density,'UniformOutput',false);
  density = horzcat(density{:});
  lg = ensurecell(get(h,'DisplayName'));

  if strcmp(plotType,'bar')
    delete(h); % remove old bars

    % add a new column
    density(:,end+1) = 0;

    % maybe we have to enlarge bins
    if maxOmega > max(bins)
      bins = 0:(bins(2)-bins(1)):maxOmega + 0.01;
      density(end+1:length(bins)-1,:) = 0;
    end
  else
    faktor = 100 / size(density,1);
  end

else

  if strcmp(plotType,'bar')

    % bin size given?
    if max(obj.angle) < maxOmega/2, maxOmega = max(obj.angle);end
    nbins = max(15,round(maxOmega/get_option(varargin,'resolution',5*degree)));

    % compute bins
    bins = linspace(-eps,maxOmega+0.01,nbins);
    density = zeros(nbins-1,1);
    lg = {};
  elseif check_option(varargin,'percent')
    faktor = 100;
  else
    faktor = 1;
    unit = 'mrd';
  end
end

% compute angle distribution
if isa(obj,'symmetry') || isa(obj,'SO3Fun')
  [density,omega] = calcAngleDistribution(obj,varargin{:});
else
  d = histcounts(obj.angle,bins).';
  midPoints = 0.5*(bins(1:end-1) + bins(2:end));
  density(:,end) = 100 * d ./ sum(d);
end

% plot angle distribution
if strcmp(plotType,'bar')

  h = optiondraw(bar(midPoints/degree,density,'parent',mtexFig.gca),varargin{:});
  xlim(mtexFig.gca,[0,max(bins)/degree])

  % update legend
  lg = [lg;{[obj.CS.mineral '-' obj.SS.mineral]}];
  for i=1:length(h)
    set(h(i),'DisplayName',lg{i});
  end

else

  h = optiondraw(plot(omega/degree,faktor * max(0,density),...
    'parent',mtexFig.gca),'LineWidth',2,varargin{:});

end

% finish
if isNew
  xlabel(mtexFig.gca,'Misorientation angle (degrees)');
  ylabel(mtexFig.gca,['Frequency (' unit ')']);
  drawNow(mtexFig,varargin{:});
end

if nargout == 0, clear h; end

author: mtex-toolbox

repo: mtex

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


text:
classdef hessenberg
% Functions to work with Hessenberg forms

methods(Static)

function [Q, R] = qr(A)
    % Givens method for computing the QR factorization A = QR of a Hessenberg matrix.
    % A is modified in the algorithm
    % GVL4: algorithm 5.2.5
    import spx.la.givens.rotation;
    [n,n] = size(A);
    % Space for saving the Q factor
    Q = eye(n, n);
    % iterate over each column from first to last but one.
    for j=1:n-1
        % We need to process only one pair of rows
        % The diagonal element and immediate sub-diagonal element.
        a = A(j,j);
        b = A(j+1,j);
        % Compute the needed rotation for making b 0.
        [c,s] = rotation(a,b);
        % Form the givens rotation matrix.
        G = [c s;-s c];
        % Rotate the two consecutive rows based submatrix
        A(j:j+1,j:n) = G'*A(j:j+1,j:n);
        % Rotate corresponding columns n Q matrix
        % by postmultiplication with G
        Q(:,j:j+1) = Q(:,j:j+1)*G;
    end
    % Return the triangular form
    R = A;
end % function

function [H, cs, ss] = qr_rq(H)
    % Takes a Hessenberg matrix H = QR and returns RQ
    % GVL4: algorithm 7.4.1
    import spx.la.givens.rotation;
    [n,n] = size(H);
    % Space for storing the c and s values
    cs = zeros(1, n);
    ss = zeros(1, n);
    % first the QR factorization
    % iterate over each column from first to last but one.
    for k=1:n-1
        % We need to process only one pair of rows
        % The diagonal element and immediate sub-diagonal element.
        a = H(k,k);
        b = H(k+1,k);
        % Compute the needed rotation for making b 0.
        [c,s] = rotation(a,b);
        % Store the rotation
        cs(k) = c;
        ss(k) = s;
        % Form the givens rotation matrix.
        G = [c s;-s c];
        % Rotate the two consecutive rows based submatrix
        H(k:k+1,k:n) = G'*H(k:k+1,k:n);
    end
    % Now the RQ formation
    for k=1:n-1
        G = [cs(k) ss(k); -ss(k) cs(k)];
        H(1:k+1,k:k+1) = H(1:k+1,k:k+1) * G;
    end
    % The product of these Givens rotations is also upper Hessenberg.
    % Product of two Hessenberg matrices is Hessenberg.
end % function

function [H, U0] = hess(A)
    % Reduction of a square matrix A into Hessenberg form
    % using Householder reflections
    % GVL4: algorithm 7.4.2
    [n,n] = size(A);
    import spx.la.house;
    % Iterate over columns of A
    for k=1:n-2
        % Compute the householder reflector for k-th column
        % covering only sub-diagonal elements
        [v,beta] = house.gen(A(k+1:n,k));
        % Update the submatrix of A by premultiplication with the Projector
        A(k+1:n,k:n) = A(k+1:n,k:n) - (beta*v)*(v'*A(k+1:n,k:n));
        % Postmultiply with the same projector without affecting
        % the 0's in the k-th column
        A(:,k+1:n) = A(:,k+1:n) - (A(:,k+1:n)*v)*(beta*v)';
        % Store the householder reflection vector
        A(k+2:n,k) = v(2:n-k);
    end
    % Extract the Hessenberg form
    H = triu(A,-1);
    if nargout > 1
        % Extract the orthogonal matrix U_0 such that
        % A = U_0 H U_0^T
        U0 = eye(n,n);
        U0(2:n,2:n) = house.q_back_accum_full(A(2:n,1:n-1));
    end
end % function

function [H, i1, i2] = backsearch(H, z)
    % Searches for a subdiagonal entry in H which is
    % negligible. From last row to first row.
    if nargin < 2
        z = size(H, 2);
    end
    i1 = z;
    i2 = z;
    % Compute the Frobenius norm of H
    h_norm = norm(H, 'fro');
    % The threshold for negligible entries
    thr = eps * h_norm
    % iterate over diagonal entries backwards
    while i1 >1
        % check for negligible subdiagonal entry.
        fprintf('%d, %d, %e\n', i1, i1-1, H(i1, i1-1));
        if (abs(H(i1, i1-1)) < thr )
            % Make this entry zero
            H(i1, i1-1) = 0;
            if (i1 == i2)
                i2 = i1 - 1;
                i1 = i1 - 1;
            else
                return;
            end
        else
            i1 = i1 - 1;
        end

    end
end

end % methods

end % classdef

author: indigits

repo: sparse-plex


text:
function [samples, energies, diagn] = metrop(f, x, options, gradf, varargin)
%METROP	Markov Chain Monte Carlo sampling with Metropolis algorithm.
%
%	Description
%	 SAMPLES = METROP(F, X, OPTIONS) uses the Metropolis algorithm to
%	sample from the distribution P ~ EXP(-F), where F is the first
%	argument to METROP.   The Markov chain starts at the point X and each
%	candidate state is picked from a Gaussian proposal distribution and
%	accepted or rejected according to the Metropolis criterion.
%
%	SAMPLES = METROP(F, X, OPTIONS, [], P1, P2, ...) allows additional
%	arguments to be passed to F().  The fourth argument is ignored, but
%	is included for compatibility with HMC and the optimisers.
%
%	[SAMPLES, ENERGIES, DIAGN] = METROP(F, X, OPTIONS) also returns a log
%	of the energy values (i.e. negative log probabilities) for the
%	samples in ENERGIES and DIAGN, a structure containing diagnostic
%	information (position and acceptance threshold) for each step of the
%	chain in DIAGN.POS and DIAGN.ACC respectively.  All candidate states
%	(including rejected ones) are stored in DIAGN.POS.
%
%	S = METROP('STATE') returns a state structure that contains the state
%	of the two random number generators RAND and RANDN. These are
%	contained in fields randstate,  randnstate.
%
%	METROP('STATE', S) resets the state to S.  If S is an integer, then
%	it is passed to RAND and RANDN. If S is a structure returned by
%	METROP('STATE') then it resets the generator to exactly the same
%	state.
%
%	The optional parameters in the OPTIONS vector have the following
%	interpretations.
%
%	OPTIONS(1) is set to 1 to display the energy values and rejection
%	threshold at each step of the Markov chain. If the value is 2, then
%	the position vectors at each step are also displayed.
%
%	OPTIONS(14) is the number of samples retained from the Markov chain;
%	default 100.
%
%	OPTIONS(15) is the number of samples omitted from the start of the
%	chain; default 0.
%
%	OPTIONS(18) is the variance of the proposal distribution; default 1.
%
%	See also
%	HMC
%

%	Copyright (c) Ian T Nabney (1996-2001)

if nargin <= 2
  if ~strcmp(f, 'state')
    error('Unknown argument to metrop');
  end
  switch nargin
    case 1
      % Return state of sampler
      samples = get_state(f);	% Function defined in this module
      return;
    case 2
      % Set the state of the sampler
      set_state(f, x);		% Function defined in this module
      return;
  end
end

if 0
seed = 42;
randn('state', seed);
rand('state', seed)
end

display = options(1);
if options(14) > 0
  nsamples = options(14);
else
  nsamples = 100;
end
if options(15) >= 0
  nomit = options(15);
else
  nomit = 0;
end
if options(18) > 0.0
  std_dev = sqrt(options(18));
else
  std_dev = 1.0;   % default
end
nparams = length(x);

% Set up string for evaluating potential function.
f = fcnchk(f, length(varargin));

samples = zeros(nsamples, nparams);		% Matrix of returned samples.
if nargout >= 2
  en_save = 1;
  energies = zeros(nsamples, 1);
else
  en_save = 0;
end
if nargout >= 3
  diagnostics = 1;
  diagn_pos = zeros(nsamples, nparams);
  diagn_acc = zeros(nsamples, 1);
else
  diagnostics = 0;
end

% Main loop.
n = - nomit + 1;
Eold = feval(f, x, varargin{:});	% Evaluate starting energy.
nreject = 0;				% Initialise count of rejected states.
while n <= nsamples

  xold = x;
  % Sample a new point from the proposal distribution
  x = xold + randn(1, nparams)*std_dev;
  %fprintf('netlab propose: xold = %5.3f,%5.3f, xnew = %5.3f,%5.3f\n',...
  %	xold(1), xold(2), x(1), x(2));

  % Now apply Metropolis algorithm.
  Enew = feval(f, x, varargin{:});	% Evaluate new energy.
  a = exp(Eold - Enew);			% Acceptance threshold.
  if (diagnostics & n > 0)
    diagn_pos(n,:) = x;
    diagn_acc(n,:) = a;
  end
  if (display > 1)
    fprintf(1, 'New position is\n');
    disp(x);
  end

  r = rand(1);
  %fprintf('netlab: n=%d, a=%f/%f=%5.3f (%5.3f), r=%5.3f\n',...
  %	  n, exp(-Enew), exp(-Eold), a, exp(-Enew)/exp(-Eold), r);
  if a > r	% Accept the new state.
    Eold = Enew;
    if (display > 0)
      fprintf(1, 'Finished step %4d  Threshold: %g\n', n, a);
    end
  else			% Reject the new state
    if n > 0
      nreject = nreject + 1;
    end
    x = xold;	% Reset position
    if (display > 0)
      fprintf(1, '  Sample rejected %4d.  Threshold: %g\n', n, a);
    end
  end
  if n > 0
    samples(n,:) = x;			% Store sample.
    if en_save
      energies(n) = Eold;		% Store energy.
    end
  end
  n = n + 1;
end

if (display > 0)
  fprintf(1, '\nFraction of samples rejected:  %g\n', ...
          nreject/(nsamples));
end

if diagnostics
  diagn.pos = diagn_pos;
  diagn.acc = diagn_acc;
end

% Return complete state of the sampler.
function state = get_state(f)

state.randstate = rand('state');
state.randnstate = randn('state');
return

% Set state of sampler, either from full state, or with an integer
function set_state(f, x)

if isnumeric(x)
  rand('state', x);
  randn('state', x);
else
  if ~isstruct(x)
    error('Second argument to metrop must be number or state structure');
  end
  if (~isfield(x, 'randstate') | ~isfield(x, 'randnstate'))
    error('Second argument to metrop must contain correct fields')
  end
  rand('state', x.randstate);
  randn('state', x.randnstate);
end
return

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
%% show our result
clc;
close all;
clear all;
addpath(genpath('.'))
load('SUNRGBDMetaUS.mat')
result_path = '../../metadata/sunrgbd/results_full';
%% change path
for imageId = 5 %480
disp(i);
bdb_3d = {};
bdb_2d = {};
disp(imageId);
bdb_result_path = fullfile(result_path, int2str(imageId), 'bdb_3d.mat');
layout_path = fullfile(result_path, int2str(imageId), 'layout.mat');
if ~exist(bdb_result_path, 'file')
    continue;
end
fig = figure;
r_path = fullfile(result_path, int2str(imageId), 'r_ex.mat');
load(r_path);
load(bdb_result_path);
data = SUNRGBDMeta(imageId);
[rgb,points3d,depthInpaint,imsize]=read3dPoints(data);
%%
for i = 1:size(bdb, 2)
    bdb_3d{i} = bdb{1, i};
end
roomLayout = load(layout_path);

% get 2d corners
for kk = 1:length(bdb_3d)
    bdb_2d_temp = get_corners_of_bb3d(bdb_3d{kk});
    bdb_2d_temp = bdb_2d_temp(:, [1, 3, 2]);
    bdb_2d_temp(:, 2) = - bdb_2d_temp(:, 2);
    bdb_2d_temp = (data.K*r_ex*bdb_2d_temp')';
    bdb_2d_corner(:, 1) = bdb_2d_temp(:, 1) ./ bdb_2d_temp(:, 3);
    bdb_2d_corner(:, 2) = bdb_2d_temp(:, 2) ./ bdb_2d_temp(:, 3);
    bdb_2d{kk} = bdb_2d_corner;
end

%% draw 2D
imshow(data.rgbpath);
hold on;
for kk =1:length(bdb_2d)
    draw_square_3d(bdb_2d{kk}, myObjectColor(kk), 5);
    x_max = min(bdb_2d{kk});
    text(double(x_max(1)), double(x_max(2)), num2str(kk), 'Color','red','FontSize',14);
end
hold off;
set(gca,'xtick',[])
set(gca,'ytick',[])
axis off;
% clf;
% print(fig, ['results/'  num2str(imageId) '_2d'], '-dpng')
% saveas(fig, ['results/'  num2str(imageId) '_2d.png']);
% close(fig);
pause(1);
%% draw 3D
fig = figure;
vis_point_cloud(points3d,rgb, 100000)
hold on;
for kk = 1:length(bdb_3d)
    vis_cube(bdb_3d{kk}, myObjectColor(kk), 2);
%     vis_cube(data.groundtruth3DBB(kk), myObjectColor(kk), 2);
end
Linewidth = 3;
maxhight = 1.2;
drawRoom(roomLayout.layout','b',Linewidth,maxhight);
hold off;
set(gca,'xtick',[])
set(gca,'ytick',[])
set(gca, 'ztick', [])
axis off;
% print(fig, ['results/'  num2str(imageId) '_3d'], '-dpng')
% saveas(fig, ['results/'  num2str(imageId) '_3d.png']);
% close(fig);
pause(1);
end

author: thusiyuan


text:
function []=uw_interp();
%UW_INTERP Interpolate grid using nearest neighbour
%
%   Andy Hooper May 2007
%
%   ============================================================================
%   01/2012 AH: Speed up read/write for triangle
%   01/2013 AH: Replace dsearch by dsearchn only for versions 2012 onwards
%   08/2014 DB: Suppress command line output
%   09/2015 AH: use matlab triangulation if triangle program not installed
%   ============================================================================

fprintf('Interpolating grid...\n')

uw=load('uw_grid','n_ps','n_ifg','nzix');

arch=computer('arch');
if strcmpi(arch(1:3),'win')
    use_triangle='n';
else
    tripath=system('which triangle >& /dev/null');
    if tripath==0
        use_triangle='y';
    else
        use_triangle='n';
    end
end

[y,x]=find(uw.nzix);
xy=[[1:uw.n_ps]',x,y];

if use_triangle=='y'
    nodename=['unwrap.1.node'];
    fid=fopen(nodename,'w');
    fprintf(fid,'%d 2 0 0\n',uw.n_ps);

    fprintf(fid,'%d %d %d\n',xy');
    fclose(fid);

    [a,b] = system('triangle -e unwrap.1.node > triangle.log');

    fid=fopen('unwrap.2.edge','r');
    header=str2num(fgetl(fid));
    N=header(1);
    edgs=fscanf(fid,'%d %d %d %d\n',[4,N])';
    fclose(fid);
    n_edge=size(edgs,1);
    if n_edge~=N
        error('missing lines in unwrap.2.edge')
    end

    fid=fopen('unwrap.2.ele','r');
    header=str2num(fgetl(fid));
    N=header(1);
    ele=fscanf(fid,'%d %d %d %d\n',[4,N])';
    fclose(fid);
    n_ele=size(ele,1);
    if n_ele~=N
        error('missing lines in unwrap.2.ele')
    end
else
    xy=double(xy);
    ele=delaunay(xy(:,2),xy(:,3));
    tr=triangulation(ele,xy(:,2),xy(:,3));
    edgs=edges(tr);
    n_edge=size(edgs,1);
    edgs=[[1:n_edge]',edgs];
    n_ele=size(ele,1);
    ele=[[1:n_ele]',ele];
end

z=[1:uw.n_ps];
[nrow,ncol]=size(uw.nzix);

[X,Y]=meshgrid(1:ncol,1:nrow);
matlab_version=version('-release');
if str2num(matlab_version(1:4))<2012
    Z=dsearch(x,y,ele(:,2:4),X,Y); % dsearch removed in MatlabR2012a
else
    Z=dsearchn([x,y],ele(:,2:4),[X(:),Y(:)]); %index from grid to pixel node
    Z = reshape(Z,nrow,ncol);
end
Zvec=Z(:);
grid_edges=[Zvec(1:end-nrow),Zvec(nrow+1:end)]; % col edges
Zvec=reshape(Z',nrow*ncol,1);
grid_edges=[grid_edges;[Zvec(1:end-ncol),Zvec(ncol+1:end)]]; % add row edges
[sort_edges,I_sort]=sort(grid_edges,2); % sort each edge to have lowest pixel node first
edge_sign=I_sort(:,2)-I_sort(:,1);
[alledges,I,J]=unique(sort_edges,'rows'); % grid_edges=alledges(J)
sameix=(alledges(:,1)==alledges(:,2));
alledges(sameix,:)=0; % set edges connecting identical nodes to (0,0)
[edgs,I2,J2]=unique(alledges,'rows');
n_edge=size(edgs,1)-1;
edgs=[[1:n_edge]',edgs(2:end,:)]; % drop (0,0)
gridedgeix=(J2(J)-1).*edge_sign; % index to edges
colix=reshape(gridedgeix(1:nrow*(ncol-1)),nrow,ncol-1);
rowix=reshape(gridedgeix(nrow*(ncol-1)+1:end),ncol,nrow-1)';

fprintf('   Number of unique edges in grid: %d\n',n_edge);

save('uw_interp','edgs','n_edge','rowix','colix','Z');

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function plot_Maps4(sFilename)
%
%
%
load(sFilename)

%% plot results mean z values
cmin=-4;cmax=4;
figure_w_normalized_uicontrolunits('Name','Z-Value','Position',[100 25 400 400]);
mData=nan(size(median(params.mZ,1, 'omitnan')));
mData=median(params.mZ,1, 'omitnan');
% mData=params.mZ;
pcolor(params.vX,params.vY,...
    reshape(mData,...
    size(params.vY,1),size(params.vX,1)));
h=colorbar;
% set(h,'XAxisLocation','Top','XTick',0.5,'XTickLabel', 'median(z)')
colormap(jet);
xlabel('longitude');
ylabel('latitude');
shading interp;
hold on;plot(params.mCatalog(:,1),params.mCatalog(:,2),'k.','MarkerSize',0.5);
plot_WiemerWyss1994
xlim([-117.1 -115.6]);
ylim([33 35.2]);

% figure_w_normalized_uicontrolunits('Name','Std Z-Value');
% mData=nan(size(mean(params.mZ,'omitnan')));
% mData=std(params.mZ, 'omitnan');
% pcolor(params.vX,params.vY,...
%     reshape(mData,...
%     size(params.vY,1),size(params.vX,1)));
% h=colorbar;
% % set(h,'XAxisLocation','Top','XTick',0.5,'XTickLabel', 'std(z)')
% colormap(jet);
% xlabel('longitude');
% ylabel('latitude');
% shading interp;
% hold on;plot(params.mCatalog(:,1),params.mCatalog(:,2),'k.','MarkerSize',0.5);

%
% figure_w_normalized_uicontrolunits('Name','Z-probability by cross-comparison');
% pcolor(params.vX,params.vY,...
%     reshape(calc_ProbColorbar2Value(1-params.vPcrZ),...
%     size(params.vY,1),size(params.vX,1)));
% plot_ProbColorbar2(cmin, cmax);
% % set(gca,'XAxisLocation','Top','XTick',0.5,'XTickLabel', 'P(z)')
%
% xlabel('longitude');
% ylabel('latitude');
% shading interp;
% hold on;plot(params.mCatalog(:,1),params.mCatalog(:,2),'k.','MarkerSize',0.5);
% % colorbar

figure_w_normalized_uicontrolunits('Name','Z-probability by overlap','Position',[100 25 400 400]);
pcolor(params.vX,params.vY,...
    reshape(calc_ProbColorbar2Value(1-params.vPolZ),...
    size(params.vY,1),size(params.vX,1)));
xlabel('longitude');
ylabel('latitude');
plot_ProbColorbar2(cmin, cmax);
% set(gca,'XAxisLocation','Top','XTick',0.5,'XTickLabel', 'P(z)')
shading interp;
hold on;plot(params.mCatalog(:,1),params.mCatalog(:,2),'k.','MarkerSize',0.5);
plot_WiemerWyss1994
xlim([-117.1 -115.6]);
ylim([33 35.2]);

%% plot results mean B-values
cmin=-4;cmax=4;
% figure_w_normalized_uicontrolunits('Name','Median Beta-Value');
% mData=nan(size(median(-params.mB, 'omitnan')));
% mData=median(params.mB, 'omitnan');
% pcolor(params.vX,params.vY,...
%     reshape(mData,...
%     size(params.vY,1),size(params.vX,1)));
% h=colorbar;
% % set(h,'XAxisLocation','Top','XTick',0.5,'XTickLabel', 'median(1-B)')
% colormap(jet);
% xlabel('longitude');
% ylabel('latitude');
% shading interp;
% hold on;plot(params.mCatalog(:,1),params.mCatalog(:,2),'k.','MarkerSize',0.5);
%
% figure_w_normalized_uicontrolunits('Name','Std Beta-Value');
% mData=nan(size(std(params.mB, 'omitnan')));
% mData=std(params.mB, 'omitnan');
% pcolor(params.vX,params.vY,...
%     reshape(mData,...
%     size(params.vY,1),size(params.vX,1)));
% h=colorbar;
% % set(h,'XAxisLocation','Top','XTick',0.5,'XTickLabel', 'std(1-B)')
% colormap(jet);
% xlabel('longitude');
% ylabel('latitude');
% shading interp;
% hold on;plot(params.mCatalog(:,1),params.mCatalog(:,2),'k.','MarkerSize',0.5);

%
% figure_w_normalized_uicontrolunits('Name','Beta-probability by cross-comparison');
% pcolor(params.vX,params.vY,...
%     reshape(calc_ProbColorbar2Value(params.vPcrB),...
%     size(params.vY,1),size(params.vX,1)));
% plot_ProbColorbar2(cmin, cmax);
% % set(gca,'XAxisLocation','Top','XTick',0.5,'XTickLabel', 'P(B)')
% xlabel('longitude');
% ylabel('latitude');
% shading interp;
% hold on;plot(params.mCatalog(:,1),params.mCatalog(:,2),'k.','MarkerSize',0.5);
% % colorbar

figure_w_normalized_uicontrolunits('Name','Beta-probability by overlap','Position',[100 25 400 400]);
pcolor(params.vX,params.vY,...
    reshape(calc_ProbColorbar2Value(params.vPolB),...
    size(params.vY,1),size(params.vX,1)));
plot_ProbColorbar2(cmin, cmax);
% set(gca,'XAxisLocation','Top','XTick',0.5,'XTickLabel', 'P(B)')
xlabel('longitude');
ylabel('latitude');
shading interp;
hold on;plot(params.mCatalog(:,1),params.mCatalog(:,2),'k.','MarkerSize',0.5);
plot_WiemerWyss1994
xlim([-117.1 -115.6]);
ylim([33 35.2]);

figure_w_normalized_uicontrolunits('Name','Resolution','Position',[100 25 400 400]);
pcolor(params.vX,params.vY,...
    reshape(params.mSamples_,...
    size(params.vY,1),size(params.vX,1)));
h=colorbar;
% set(h,'XAxisLocation','Top','XTick',0.5,'XTickLabel', 'N')
xlabel('longitude');
ylabel('latitude');
shading interp;
hold on;plot(params.mCatalog(:,1),params.mCatalog(:,2),'k.','MarkerSize',0.5);
plot_WiemerWyss1994
xlim([-117.1 -115.6]);
ylim([33 35.2]);

end

function plot_WiemerWyss1994
mLatLon=[[-117.1 33];
    [-115.6 33];
    [-115.6 35.2];
    [-117.1 35.2];
    [-117.1 33]];
hold on;plot(mLatLon(:,1),mLatLon(:,2),'k--');
mLanders=[-116.4 34.3];
hold on;plot(mLanders(1),mLanders(2),'ko','MarkerSize',20)
end

author: CelsoReyes


text:
function c=lpccc2cc(cc,np)
%LPCCC2PF Extrapolate complex cepstrum C=(CC)

%      Copyright (C) Mike Brookes 1998
%      Version: $Id: lpccc2cc.m,v 1.4 2007/05/04 07:01:38 dmb Exp $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

p=size(cc,2);
if nargin<2 np=p; end
if np<=p
   c=cc(:,1:np);
else
   c=lpcar2cc(lpccc2ar(cc),np);
end

author: decouples


text:
function wstruct = swing(wsize,color,ax)

% SWING  Create a graphics structure with a wing shape.
%   WSTRUCT = SWING(WSIZE,COLR,AX) creates a structure for a 'patch'
%   graphic object correspoinding to a wing of size WSIZE, color COLR, in
%   axes AX.
%
%   See also WING, PATCH.

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

if nargin < 3
    ax = gca;
end

wstruct = wing(wsize);
wstruct.handle = patch(...
    'parent'   ,ax,...
    'vertices' ,wstruct.vert,...
    'faces'    ,wstruct.faces,...
    'linewidth',2,...
    'facecolor','none',...
    'edgecolor',color,...
    'visible'  ,'on');

% ========== 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:
classdef LHSintegrator < handle

    properties (Access = protected)
        mesh
        fun
        quadrature
        quadratureOrder
    end

    methods (Access = public, Static)

        function obj = create(s)
            f = LHSintegratorFactory();
            obj = f.create(s);
        end

    end

    methods (Access = public)

        function obj = LHSintegrator(cParams)
            obj.init(cParams);
            obj.createQuadrature();
        end

    end

    methods (Access = protected)

        function init(obj, cParams)
            obj.fun      = cParams.fun;
            obj.mesh     = cParams.mesh;
            obj.setQuadratureOrder(cParams);
        end

        function setQuadratureOrder(obj, cParams)
            if isfield(cParams, 'quadratureOrder')
                obj.quadratureOrder = cParams.quadratureOrder;
            else
                obj.quadratureOrder = obj.fun.order;
            end
        end

        function createQuadrature(obj)
            quad = Quadrature.set(obj.mesh.type);
            quad.computeQuadrature(obj.quadratureOrder);
            obj.quadrature = quad;
        end

    end

end

author: SwanLab


text:
function r = name(a)

r = a.name;

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function [elem,idx,volume] = fixorder3(node,elem)
%% FIXORDER3 fix orientation of tetrahedron
%
%   elem = FIXORDER3(node,elem) computes signed volume of all tetrahedron
%   in the triangulation and switch the vertices such that all signed
%   volume is positive.
%
%   [elem,idx,volume] = FIXORDER3(node,elem) also outputs the index set of
%   elements whose area is negative and the absolute value of area.
%
% See also fixorientation, fixorder3
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

% compute signed volume of each tetrahedron
volume = simplexvolume(node,elem);
% find tetrahedron with negative volume and switch the vertices
idx = find(volume<0);
elem(idx,[2 3]) = elem(idx,[3 2]);
volume(idx) = -volume(idx);

repo: Aero-matlab


text:
function writeSTL(filename, V,F)
  % WRITESTL writes an STL file with vertex/face information. Uses single
  % precision.
  %
  % writeSTL(filename,V,F)
  %
  % Input:
  %  filename  path to .stl file
  %  V  #V by 3 list of vertices (Default units are in mm)
  %  F  #F by 3 list of triangle indices
  %
  stlwrite(filename,F,V);
end

author: alecjacobson