text:
function options = dsdpset(varargin)
%DSDPSET  Create or alter the options for Optimization with DSDP
%
% options = dsdpset('param1',value1,'param2',value2,...) creates an DSDP
% options structure with the parameters 'param' set to their corresponding
% values in 'value'. Parameters not specified will be set to the DSDP
% default.
%
% options = dsdpset(oldopts,'param1',value1,...) creates a copy of the old
% options 'oldopts' and then fills in (or writes over) the parameters
% specified by 'param' and 'value'.
%
% options = dsdpset() creates an options structure with all fields set to
% dsdpset defaults.
%
% dsdpset() prints a list of all possible fields and their function.
%
% See supplied DSDP Documentation for further details of these options.

%   Copyright (C) 2013 Jonathan Currie (I2C2)

% Print out possible values of properties.
if (nargin == 0) && (nargout == 0)
    printfields();
    return
end
%Names and Defaults
Names = {'drho';'rpos';'r0';'penalty';'rho';'dbound';'gaptol';'rtol';'mu0';'maxtrust';...
         'steptol';'ptol';'pnormtol';'reuse';'zbar';'dlbound';'ybound';'fixed'};
Defaults = {1;0;-1;1e8;5;1e20;1e-7;1e-6;-1;1e10;0.05;1e-4;1e30;4;[];[];[];[]};

%Enter and check user args
try
    options = opticheckset(Names,Defaults,@checkfield,varargin{:});
catch ME
    throw(ME);
end

function checkfield(field,value)
%Check a field contains correct data type
switch lower(field)
    %Scalar double
    case {'r0','penalty','rho','mu0','dbound','pnormtol','zbar','dlbound','ybound'}
        err = opticheckval.checkScalarDbl(value,field);
    %Scalar non negative double
    case {'rtol','ptol','maxtrust','gaptol','steptol'}
        err = opticheckval.checkScalarNonNeg(value,field);
    %Scalar 0/1
    case {'rpos','drho'}
        err = opticheckval.checkScalar01(value,field);
    %Scalar non-negative integer
    case 'reuse'
        err = opticheckval.checkScalarIntNonNeg(value,field);
    %Double column vector
    case 'fixed'
        err = opticheckval.checkDualCol(value,field);
    otherwise
        err = MException('OPTI:SetFieldError','Unrecognized parameter name ''%s''.', field);
end
if(~isempty(err)), throw(err); end

function printfields()
%Print out fields with defaults
fprintf('          rpos: [ Use Penalty Parameter to enforce Feasibility {0}, 1 ] \n');
fprintf('          drho: [ Use Dynamic Strategy to choose Rho: 0, {1} ] \n');
fprintf('            r0: [ Set initial value for variable r in Dual {-1} ] \n');
fprintf('       penalty: [ Penalty Parameter Gamma {1e8} \n');
fprintf('           rho: [ Potential Parmater: {5} ] \n');
fprintf('           mu0: [ Barrier Parameter: {-1} ] \n');
fprintf('          rtol: [ Classify Dual as feasible if r is less than this tolerance: {1e-6} ] \n');
fprintf('          ptol: [ Classify Primal as feasible if infeasibility is less than this tolerance: {1e-4} ] \n');
fprintf('      maxtrust: [ Maximum Trust Radius on Step Direction: {1e10} ] \n');
fprintf('        gaptol: [ Convergence Gap Tolerance: {1e-7} ] \n');
fprintf('       steptol: [ Terminate solver if step length in Dual is below this tolerance: {0.05} ] \n');
fprintf('        dbound: [ Terminate solver if Dual Objective greater than this Value {1e20} ] \n');
fprintf('      pnormtol: [ Terminate the solver when duality gap is sufficiently small and pnorm is less than this quantity: {1e30} ] \n');
fprintf('         reuse: [ Number of times the Hessian of the Barrier Function will be reused {4} ] \n');
fprintf('          zbar: [ Upper bound on objective value at the solution {[]} ] \n');
fprintf('       dlbound: [ Lower bound on dual value {[]} ] \n');
fprintf('        ybound: [ Bound on the dual variables y {[]} ] \n');
fprintf('         fixed: [ Matrix of fixed variables, column 1 variable indices, column 2, fixed values {[]} ] \n');
fprintf('\n');

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function ref = physical_to_reference_tet4 ( t, n, phy )

%*****************************************************************************80
%
%% PHYSICAL_TO_REFERENCE_TET4 maps physical points to reference points.
%
%  Discussion:
%
%    Given the vertices of an order 4 physical tetrahedron and a point
%    (X,Y,Z) in the physical tetrahedron, the routine computes the value
%    of the corresponding point (R,S,T) in the reference tetrahedron.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(3,4), the coordinates of the vertices of the
%    physical tetrahedron.  The vertices are assumed to be the images of
%    (1,0,0), (0,1,0), (0,0,1) and (0,0,0) respectively.
%
%    Input, integer N, the number of points to transform.
%
%    Input, real PHY(3,N), the coordinates of physical points
%    to be transformed.
%
%    Output, real REF(3,N), the coordinates of the corresponding
%    points in the reference tetrahedron.
%
  a(1:3,1:3) = t(1:3,1:3);
  for i = 1 : 3
    a(i,1:3) = a(i,1:3) - t(i,4);
  end

  for i = 1 : 3
    ref(i,1:n) = phy(i,1:n) - t(i,4);
  end

  ref(1:3,1:n) = a(1:3,1:3) \ ref(1:3,1:n);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [pvalue Vlstar Vl] = WhiteRealityCheck( alternative , flag, benchmark , n , display)

%% This file contains the White reality check for data snooping
%  -   This can be used to test whether atleast one of the created models
%  outperforms a benchmark model
%  - This can also be used to test whether you have found atleast one profitable
%  trading strategy, or test whether you have found atleast one trading
%  strategy that outperforms the benchmark.

%  The H0 hypothesis is always  that you have not found an outperforming
%  strategy or model

% input:
% 'Alternative' is a matrix with  Returns ==> [-1 , +Inf ]
%  or residuals of predictions with a model ==> [ -Inf ; +Inf]
% Important, the rows stand for new observation the columns for the models
%  'Flag'  indicates which loss function you want to use
%  Flag = 1 test for model superiority vs benchmark by mean squared error
%  Flag = 2 test for trading return superiorty vs benchmark
%  Flag = 3 test for model superiority vs benchmark by absolute error
% 'Benchmark' contains a vector with - (1) The residuals from the predictions
% of the benchmark model , (2) The returns of your benchmark trading
% strategy, (3) Special case where Benchmark is the single number 0 where
% your benchmark is zero.
%  'n' is the number of bootstrapped series you want to create aka the
%  number of simulations. In academic papers this is usually  set
%  to atleast 500 for reasonable results
% 'display' is a number [ 0, 1] if 0 then display is set off if 1 then
% display is on

% Examples of possibles testing can be found in help.txt file

[r,~]=size(alternative);

%  condition whether you have a benchmark or you want to test vs 0
if benchmark ~=0|| numel(benchmark)>1
    mat = repmat(benchmark,1,c);
else
    mat = benchmark; % mat = 0 benchmark is zero
end

% %  loss functions
if flag ==1
     f = - alternative.^2 + mat.^2;
elseif flag==2
    f = log(1+alternative)-log(1+mat);
elseif flag ==3
    f = - abs (alternative) + abs(mat);
end

%  input for average block size for P&R bootstrap (geometric distribution)
input = inputdlg('What do you want as average block size for the Politis Romano stationary bootstrap?');
blockparam = 1/(str2num(input{1})+1);

% Actual Politis Romano bootstrap as described in Politis& Romano (1994)
fstarroof =  PolitisRomanoBootstrap( alternative , n , blockparam,display,flag,mat);

froof = mean(f,1);
Vl = max(sqrt(r)*froof);

delta = fstarroof - repmat(froof,n,1);
Vlstar = max(sqrt(r)*delta,[],2);

Vlstar = sort(Vlstar);

% 'Vl' and 'Vlstar' are the same as in the paper of
% Sullivan, Timmermann and White (1999) or White(2000)

better = Vlstar>Vl;

pvalue = sum(better)/n;

% compare Vl and Vlstar to get the p-value

author: Sable

repo: mcbench-benchmarks


text:
function cor_fn = samplingPanoBox(cor_id, corn, edg, edg2, im_h, im_w, options)

    corn = im2double(corn);
    edg = im2double(edg);

    score_fin = -Inf;
    score_fin_t = -Inf;
    search_grid = 5;
    line = cor_id(1:2:end,1)';
    line_can = repmat(line, size(line,2)*(2*search_grid+1),1);
    score_ln = interp2(corn,cor_id(:,1),cor_id(:,2));
    score_ln = log(score_ln);
    score_ln = [score_ln(1) + score_ln(2); score_ln(3) + score_ln(4);...
               score_ln(5) + score_ln(6); score_ln(7) + score_ln(8)];
    [~,score_id] = sort(score_ln);

    line_can(1:(2*search_grid+1),score_id(1)) = line(score_id(1))-search_grid:line(score_id(1))+search_grid;
    line_can((2*search_grid+1)+1:(2*search_grid+1)*2, score_id(2)) = line(score_id(2))-search_grid:line(score_id(2))+search_grid;
    line_can((2*search_grid+1)*2+1:(2*search_grid+1)*3, score_id(3)) = line(score_id(3))-search_grid:line(score_id(3))+search_grid;
    line_can((2*search_grid+1)*3+1:(2*search_grid+1)*4, score_id(4)) = line(score_id(4))-search_grid:line(score_id(4))+search_grid;
    edg2 = im2double(edg2);

    % no sampling
    %line_can = line;
    for line_n = 1:size(line_can,1)
        %disp(line_n)
        cor_ini = cor_id;
        cor_ini(1:2:end,1) = line_can(line_n,:);
        cor_ini(2:2:end,1) = line_can(line_n,:);
        [cor_ini, score_all] = sample_opt_pano_joint(cor_ini, im_w, im_h, corn, edg, edg2, options);
        if score_all > score_fin_t
            score_fin_t = score_all;
            cor_fn_t = cor_ini;
        end
        if mod(line_n,search_grid*2+1) == 0
            line_id = line_n/(search_grid*2+1);
            line_can(line_n+1:end,score_id(line_id)) = cor_fn_t(2*score_id(line_id),1);
            if score_fin_t > score_fin
                score_fin = score_fin_t;
                cor_fn = cor_fn_t;
            end
            score_fin_t = -Inf;
        end
    end
    % no sampling
    %if score_fin_t > score_fin
    %    score_fin = score_fin_t;
    %    cor_fn = cor_fn_t;
    %end
end

author: zouchuhang

repo: LayoutNet


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fbestloc.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% checks whether a new point in the shopping basket (found by local
% search) is better than the current best point; in that case it
% updates the best point and its function value
% Uses the following m-file:
% chrelerr.m

if fmi(nbasket0) < fbest
  fbest = fmi(nbasket0)
  xbest = xmin(:,nbasket0)
  chrelerr;
end

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function [eMat,f2eMat,f2eIsDirMat] = mapface2edge(vMat,fMat)
% MAPFACE2EDGE creates a mapping from faces to edges
%
% Input:
%   regular:
%       vMat: double[nVerts,3] - coordinates of vertices
%       fMat: double[nFace,3] - indices of face vertices in vMat
%
% Output:
%   eMat: double[nEdges,2] - contains indices of vertices corresponding
%       to each edge
%
%   f2eMat: double[nFaces,3] - contains indices of edges for
%       each face in this order (1-2, 2-3, 1-3)
%
%   f2eIsDirMat: logical[nFaces,3] - contains true if face
%       references edge in a direct order (i.e. 1-2 for instance)
%       and false if reference is in an opposite order
%
% $Author: Peter Gagarinov, PhD  <pgagarinov@gmail.com> $
% $Copyright: Peter Gagarinov, PhD,
%            Moscow State University,
%            Faculty of Computational Mathematics and Computer Science,
%            System Analysis Department 2011-2016 $
%
trObj = triangulation(fMat,vMat);
eMat=trObj.edges;
nEdges=size(eMat,1);
nFaces=size(fMat,1);
%% Build Face to Edges map and edge orientation map for each face
fMidCVec=trObj.edgeAttachments(eMat);
indShiftVec=cellfun('length',fMidCVec);
indF2EVec=zeros(sum(indShiftVec),1);
indF2EVec(1)=1;
indF2EVec(1+cumsum(indShiftVec(1:end-1)))=ones(nEdges-1,1);
indF2EVec=cumsum(indF2EVec);
indFVec=[fMidCVec{:}].';
%edges for each face are expected to be oriented as
%1-2, 2-3 , 3-1
indEdgeNumVec=...
    all(fMat(indFVec,[1,2])==eMat(indF2EVec,:),2)...
    -all(fMat(indFVec,[2,1])==eMat(indF2EVec,:),2)...
    +2*all(fMat(indFVec,[2,3])==eMat(indF2EVec,:),2)...
    -2*all(fMat(indFVec,[3,2])==eMat(indF2EVec,:),2)...
    +3*all(fMat(indFVec,[1,3])==eMat(indF2EVec,:),2)...
    -3*all(fMat(indFVec,[3,1])==eMat(indF2EVec,:),2);
%
[~,indSortVec]=sortrows([indFVec,abs(indEdgeNumVec)]);
indF2EVec=indF2EVec(indSortVec);
f2eMat=reshape(indF2EVec,3,nFaces).';
f2eIsDirMat=reshape(indEdgeNumVec(indSortVec)>0,3,nFaces).';

repo: OpenMAS


text:
function [x,xL,t,intersects]=findClosestPointLineEllipsoid(xc,A,x0,a,aType,epsConst)
%%FINDCLOSESTPOINTLINEELLIPSOID Given a 2D or 3D ellipsoid of the form
%           (x-xc)'*A*(x-xc)=1, and a line of the form x=x0+s*t, where t is
%           a scalar parameteric parameter, find the point on the ellipsoid
%           that is closest to the line as well as the point on the line
%           that is closest to the ellipsoid. If the line intersects the
%           ellipsoid, then more than once point can be returned. All
%           solutions and inputs are real.
%
%INPUTS: xc The xDimX1 center of the ellipsoid.
%         A The xDimXxDim symmetric positive definite matrix defining the
%           above equation for an ellipsoid.
%     x0, a These xDimX1 values define a line and the type of line is given
%           by aType. If aType=0, then the line is parameteric as x=x0+a*t
%           where t is a scalar value. If aType=1, then x0 and a are both
%           points on the line.
%     aType As mentioned above, this specifies how the line is
%           parameterized. The default if omitted or an empty matrix is
%           passed is 0.
%  epsConst In the rooting formulation of the solution, is is possible that
%           some candidate solutions are invalid (they do not satisfy the
%           constraint (x-xc)'*A*(x-xc)=1. Thus, the absolute value of a
%           transformed version of abs((x-xc)'*A*(x-xc)-1) must be
%           <=epsConst to be valid. The default if omitted or an empty
%           matrix is passed is 1e-9.
%
%OUTPUTS: x The xDimXnumSol set of closest points on the ellipsoid. If the
%           line and ellipsoid do not intersect, numSol=1.
%        xL The xDimXnumSol set of closest points on the line. These
%           correspond to those in x.
%         t The 1XnumSol set of parametric values such that the kth one is
%           xL(:,k)=x0+s*t(k).
% intersects If the line intersects the ellipsoid, this is true. Otherwise,
%           it is false.
%
%This function utilizes a change of coordiantes to solve the problem.
%We have to solve the optimization problem.
%minimize norm(x0+s*t-x)^2
%such that
%(x-xc)'*A*(x-xc)=1
%We perform an eigenvalue decomposition of A as A=U*D*U', where D is a
%diagonal matrix and U is in a unitary matrix (because A is symmetric and
%positive definite). Next, perform the change of variables
%x=U*y+xc
%s=U*c1        or c1=U'*s
%c=x0-xc=U*c2  or c2=U'*(x0-xc)
%c=U*c2
%Because U'*U=I, the optimization problem becomes:
%norm(c1*t-y+c2)^2
%Subject to
%y'*D*y=1
%Thus, we solve for y and transform back into x. The method of solution is
%Lagrangian relaxation. Chapter 4.1 or [1] goes over Lagrangian
%multipliers. The Lagrangian for the problem is
%L=norm(c1*t-y+c2)^2+lambda*(y'*D*y-1)
%Setting the derivatives of L with respect to t and y equal to 0, one gets
%a solution in terms of lambda (the Lagrangian parameter). One can then
%choose lambda to satisfy the constraint y'*D*y=1. In this instance, that
%can be formulated as involves solving a polynomial in terms of lambda.
%That is what this function does. Due to possible singularities in a ratio
%that is reduced to the polynomial, it is possible that invalid solutiosn
%will be found. Thus, the epsConst parameter is introduced to throw out
%possible solutions that do not sufficiently satisfy the condition. Of the
%multiple solutions found for lambda, only the one minimizing the cost
%function is kept.
%
%EXAMPLE 1:
%This example is done in 2D, so it can be easily plotted. Multiple lines
%are drawn. Two intersect an ellipse and two do not intersect the ellipse.
%The closest points on the ellipse are drawn in red. Those on the lines at
%black circles. Dashed black lines connect them.
% xc=[2;2];
% A=[0.5,-0.1;
%   -0.1, 0.2];
% figure(1)
% clf
% hold on
% drawEllipse(xc,A,1,'linewidth',2)
% axis([-2, 6, -2, 6])
% axis square
% t=linspace(-10,10,500);
% %Consider 4 lines.
% a=[1,0,0, -2;
%    2,0,-2, -1];
% b=[4,1,-2,3;
%    3,1,-1,-1;];
% numLines=size(a,2);
% for k=1:numLines
%     xy=bsxfun(@plus,a(:,k),bsxfun(@times,b(:,k),t));
%     plot(xy(1,:),xy(2,:),'linewidth',2)
%     [x,xL]=findClosestPointLineEllipsoid(xc,A,a(:,k),b(:,k));
%
%     scatter(x(1,:),x(2,:),200,'.r')
%     scatter(xL(1,:),xL(2,:),100,'ok')
%     if(size(xL,2)==1)
%         %If the line either grazes the unit circle or does not intercept
%         %the unit circle, connect the point on the line to the nearest
%         %point on the circle.
%         plot([x(1),xL(1)],[x(2),xL(2)],'--k')
%     end
% end
%
%EXAMPLE 2:
%This example is similar to the first one, but is in 3D.
% xc=[2;2;0];
% A=[0.5,-0.1, 0.15;
%   -0.1, 0.2, 0.15;
%    0.15,0.15,0.3];
% figure(1)
% clf
% hold on
% drawEllipse(xc,A,1,'edgecolor','none')
% light()
% %axis([-2, 6, -2, 6])
% %axis square
% t=linspace(-10,10,500);
% %Consider 4 lines.
% a=[1,0,0, -2;
%    2,0,-2, -1;
%    1,0,1, -1];
% b=[4,1,-2,3;
%    3,1,-1,-1;
%    1,2,0,0];
% numLines=size(a,2);
% for k=1:numLines
%     xyz=bsxfun(@plus,a(:,k),bsxfun(@times,b(:,k),t));
%     plot3(xyz(1,:),xyz(2,:),xyz(3,:),'linewidth',2)
%     [x,xL]=findClosestPointLineEllipsoid(xc,A,a(:,k),b(:,k));
%     scatter3(x(1,:),x(2,:),x(3,:),200,'.r')
%     scatter3(xL(1,:),xL(2,:),xL(3,:),100,'ok')
%     if(size(xL,2)==1)
%         %If the line either grazes the unit circle or does not intercept
%         %the unit circle, connect the point on the line to the nearest
%         %point on the circle.
%         plot3([x(1),xL(1)],[x(2),xL(2)],[x(3),xL(3)],'--k')
%     end
% end
% axis([-6, 6, -6, 6,-6,6])
% axis square
% view(-8,30)
%
%REFERENCES:
%[1] D. P. Bertsekas, Nonlinear Programming, 2nd ed. Belmont, MA:
%    Athena Science, 1999.
%
%October 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<6||isempty(epsConst))
    epsConst=1e-9;
end

if(nargin<5||isempty(aType))
    aType=0;
end

if(aType==1)
    s=a-x0;
else
    s=a;
end

%The special case of s being 0. This means that we just find the closest
%point on the ellipsoid to the point x0.
if(all(s==0))
    x=nearestPointOnEllipsoid(xc,A,x0,1);
    xL=x0;
    t=0;
    intersects=all(x==xL);
    return;
end

%Check the special case of the line intercepting the ellipsoid. In this
%instance,
[x,t]=findEllipsLineIntersect(xc,A,x0,s,0);
if(~isempty(x))
    t=t';
    xL=x;
    intersects=true;
    return;
end

%There is a line and it does not intercept the ellipsoid. Find the closest
%point.
intersects=false;

[~,D,U]=eig(A);
c1=U'*s;
c2=U'*(x0-xc);

[y,t]=solveDiagA(D,c2,c1,epsConst);

%Transform back to x
x=bsxfun(@plus,U*y,xc);
xL=bsxfun(@plus,x0,bsxfun(@times,s,t));

end

function [x,t]=solveDiagA(A,a,b,epsConst)

numDim=size(a,1);
switch(numDim)
    case 2
        a11=A(1,1);
        a22=A(2,2);
        a1=a(1);
        a2=a(2);
        b1=b(1);
        b2=b(2);

        firstTerm=-a11*a22*(b1^2+b2^2)*(a11*b1^2+a22*b2^2);
        sqrtTerm=sqrt(a11^3*a22^3*(a2*b1-a1*b2)^2*(b1^2+b2^2)^2*(a11*b1^2+a22*b2^2));
        denom=a11^2*a22^2*(b1^2+b2^2)^2;

        if(imag(sqrtTerm)~=0)
            %If there is some finite precision error.
            x=[];
            t=[];
            return
        end

        lambda(1)=(firstTerm+sqrtTerm)/denom;
        lambda(2)=(firstTerm-sqrtTerm)/denom;

        x=zeros(2,2);
        t=zeros(1,2);
        cost=zeros(2,1);
        constErr=zeros(2,1);
        for k=1:2
            denom=a22*b2^2+a11*(b1^2+a22*(b1^2+b2^2)*lambda(k));

            x(:,k)=[a22*b2*(a1*b2-a2*b1);a11*b1*(a2*b1-a1*b2)]/denom;
            t(k)=-(a2*a22*b2*(1+a11*lambda(k))+a1*a11*b1*(1+a22*lambda(k)))/denom;

            cost(k)=norm(a+b*t(k)-x);
            constErr(k)=abs(x(:,k)'*A*x(:,k)-1);
        end
        sel=(constErr<=epsConst);
        if(sum(sel)==0)
            %If nothing met the constraint.
            x=[];
            t=[];
            return
        end
        x=x(:,sel);
        t=t(sel);
        cost=cost(sel);

        [~,idx]=min(cost);

        x=x(:,idx);
        t=t(idx);
    case 3
        a11=A(1,1);
        a22=A(2,2);
        a33=A(3,3);

        a1=a(1);
        a2=a(2);
        a3=a(3);
        b1=b(1);
        b2=b(2);
        b3=b(3);

        c0=(a11*b1^2+a22*b2^2+a33*b3^2)*(a22*(-1+a3^2*a33)*b2^2-2*a2*a22*a3*a33*b2*b3+(-1+a2^2*a22)*a33*b3^2+a11*((-1+a2^2*a22+a3^2*a33)*b1^2-2*a1*b1*(a2*a22*b2+a3*a33*b3)+a1^2*(a22*b2^2+a33*b3^2)));
        c1=2*(a11*b1^2+a22*b2^2+a33*b3^2)*(-a11*a33*(b1^2+b3^2)-a22*a33*(b2^2+b3^2)+a11*a22*((-1+(a2^2+a3^2)*a33)*b1^2+(-1+(a1^2+a3^2)*a33)*b2^2-2*a2*a3*a33*b2*b3+(a1^2+a2^2)*a33*b3^2-2*a1*a33*b1*(a2*b2+a3*b3)));
        c2=(-a22^2*a33^2*(b2^2+b3^2)^2+a11^2*(-a33^2*(b1^2+b3^2)^2+a22^2*((-1+a3^2*a33)*(b1^2+b2^2)^2-2*a3*a33*(a1*b1+a2*b2)*(b1^2+b2^2)*b3+a33*(a1*b1+a2*b2)^2*b3^2)+a22*a33*((-4+a2^2*a33)*b1^4-2*a1*a2*a33*b1^3*b2+2*a1*a33*b1*b2*b3*(a3*b2-a2*b3)+b3^2*((-2+a3^2*a33)*b2^2-2*a2*a3*a33*b2*b3+a2^2*a33*b3^2)+b1^2*((-4+a1^2*a33)*b2^2-2*a2*a3*a33*b2*b3+2*(-2+a2^2*a33)*b3^2)))+a11*a22*a33*(-2*a33*(b1^2*b2^2+2*(b1^2+b2^2)*b3^2+2*b3^4)+a22*(-2*a1*a33*b1*(a2*b2+a3*b3)*(b2^2+b3^2)+(b2^2+b3^2)*((-4+a1^2*a33)*b2^2+a1^2*a33*b3^2)+b1^2*((-4+a2^2*a33)*b2^2+2*a2*a3*a33*b2*b3+(-2+a3^2*a33)*b3^2))));
        c3=-2*a11*a22*a33*(b1^2+b2^2+b3^2)*(a11*a22*(b1^2+b2^2)+a11*a33*(b1^2+b3^2)+a22*a33*(b2^2+b3^2));
        c4=-a11^2*(a22^2)*(a33^2)*((b1^2+b2^2+b3^2)^2);

        lambda=roots([c4;c3;c2;c1;c0]);
        lambda=lambda(imag(lambda)==0);
        numSol=length(lambda);

        if(numSol==0)
            %If there is some finite precision error. There should be two
            %solutions.
            x=[];
            t=[];
            return
        end

        x=zeros(3,numSol);
        t=zeros(1,numSol);
        cost=zeros(numSol,1);
        %Due to singularities, some solutions might not be correct. Thus,
        %we also evaluate the relative error for enforcing the constraint
        %and if the error is too high, those solutions are discarded.
        constErr=zeros(numSol,1);
        for k=1:numSol
            lambdaCur=lambda(k);
            denom=a22*b2^2+a33*b3^2+a22*a33*(b2^2+b3^2)*lambdaCur+a11*b1^2*(1+a22*lambdaCur)*(1+a33*lambdaCur)+a11*lambdaCur*(a33*b3^2+a22*(b2^2+a33*(b2^2+b3^2)*lambdaCur));
            x(:,k)=[-a3*a33*b1*b3*(1+a22*lambdaCur)-a2*a22*b1*b2*(1+a33*lambdaCur)+a1*(a22*b2^2+a33*b3^2+a22*a33*(b2^2+b3^2)*lambdaCur);
                     a11*b1*(a2*b1-a1*b2)+a33*b3*(-a3*b2+a2*b3)+a11*a33*(-b2*(a1*b1+a3*b3)+a2*(b1^2+b3^2))*lambdaCur;
                     a11*b1*(a3*b1-a1*b3)+a22*b2*(a3*b2-a2*b3)+a11*a22*(a3*(b1^2+b2^2)-(a1*b1+a2*b2)*b3)*lambdaCur]/denom;
            t(k)=-(a1*a11*b1*(1+a22*lambdaCur)*(1+a33*lambdaCur)+(1+a11*lambdaCur)*(a3*a33*b3*(1+a22*lambdaCur)+a2*a22*b2*(1+a33*lambdaCur)))/denom;
            cost(k)=norm(a+b*t(k)-x);
            constErr(k)=abs(x(:,k)'*A*x(:,k)-1);
        end
        sel=(constErr<=epsConst);
        if(sum(sel)==0)
            %If nothing met the constraint.
            x=[];
            t=[];
            return
        end
        x=x(:,sel);
        t=t(sel);
        cost=cost(sel);

        [~,idx]=min(cost);
        x=x(:,idx);
        t=t(idx);
    otherwise
        error('An unsupported dimensionality is being used.')
end
end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function mov=videoread(fname,varargin)
%VIDEOREAD Read an video file in a mannter compatible with Matlab's
%          built-in AVIREAD function, but using a videoReader object to do
%          the work instead.
%
%   MOV=VIDEOREAD(FNAME)
%     Read the video file FNAME into the Matlab movie structure MOV.  MOV
%     has two fields, "cdata" and "colormap".  Since all videoReader
%     plugins currently output 24-bit images in all cases, colormap will
%     always be blank.
%
%   MOV=VIDEOREAD(FNAME,VRARGS)
%     The same as MOV=VIDEOREAD(FNAME), but VRARGS are passed to the
%     videoReader constructor.  VRARGS should be a cell array of arguments.
%     E.g.
%         mov = videoread('tests/numbers.wmv3.avi', {'DirectShow'});
%
%   MOV=VIDEOREAD(FNAME,INDEX)
%     Reads only the frame(s) specified by INDEX.  INDEX can be a single
%     index or an array of indices into the video stream, where the first
%     frame is  number one.  IMPORTANT NOTES: If you use videoReader
%     directly, it assumes  frames start at 0, not 1.  We use 1-indexing
%     here in VIDEOREAD to maximize  compatibility with MathWork's AVIREAD.
%     Also, if you find yourself making several calls to VIDEOREAD to
%     extract different portions of the file and/or if the number of
%     frames is large, consider using videoReader directly instead of this
%     wrapper function.
%
%   MOV=VIDEOREAD(FNAME,VRARGS,INDEX)
%     The same as MOV=VIDEOREAD(FNAME,INDEX), but VRARGS are passed to the
%     videoReader constructor.  VRARGS should be a cell array of arguments.
%
%Copyright (c) 2006 Gerald Dalley
%See "MIT.txt" in the installation directory for licensing details (especially
%when using this library on GNU/Linux).

if nargin==1,
  vrargs = {};
elseif nargin==2,
  if (iscell(varargin{1})),
    vrargs = varargin{1};
  else
    vrargs = {};
    index  = varargin{1};
  end
elseif nargin==3,
  vrargs = varargin{1};
  index  = varargin{2};
end

vr = videoReader(fname, vrargs{:});
if ~(exist('index','var')),
  info = getinfo(vr);
  index = (0:info.numFrames-1); % 0-indexed
else
  index = index - 1; % 0-indexed
end

mov = struct;
currFrame = -1;
if (size(index)>0)
  % index supplied or deduced from file
  for i=1:length(index)
    f = index(i);
    if (f == currFrame+1)
      if ~next(vr), error('Could not advance to frame %d\n', f+1); end
    else
      if ~seek(vr, f), error('Could not seek to frame %d\n', f+1); end
    end
    currFrame = f;
    mov(i).cdata = getframe(vr);
    mov(i).colormap = [];
  end
else
  % no index was supplied and we could not deduce the number of frames
  % in the video.
  i=1;
  while next(vr)
    mov(i).cdata = getframe(vr);
    mov(i).colormap = [];
    i = i+1;
  end
end

vr = close(vr);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [xAL,xBL,doesClip]=clipLine2ConvexPolygon(xAL,xBL,v)
%%CLIPLINE2CONVEXPOLYGON Given the vertices of a convex polygon in 2D in
%       either clockwise or counterclockwise order, clip one or more lines
%       to the polygon. The lines are defined by two points on the lines,
%       though they extend beyond the two points.
%
%INPUTS: xAL, xBL These are 2 2XnumLines matrices where xAL(:,k) and
%                 xBL(:,k) are two points on the kth line.
%               v A 2XnumVertices set of vertices defining the perimeter of
%                 the convex polygon. The last vertex should not be a
%                 repeat of the first vertex.
%
%OUTPUTS: xAL, xBL These are 2XnumLines set of the starting and ending
%                  points of the lines clipped to the polygon. If the kth
%                  line is complete outside the polygon, then xAL(:,k) and
%                  xBL(:,k) will be NaNs.
%         doesClip A numLinesX1 boolean vector indicating whether or not
%                  each line is in or on the polygon at all.
%
%EXAMPLE:
%Draw a polygon and a few line segments in black. The line segments
%represent infinte lines. Clip the lines to the polygon and draw the
%clipped lines in red. One of the lines does not clip, so a black segment
%is drawn, but nothing is drawn in red.
% numVertices=5;
% v=zeros(2,numVertices);
% v(:,1)=[1;1];
% v(:,2)=[15;0];
% v(:,3)=[10;10];
% v(:,4)=[7;10];
% v(:,5)=[0;8];
%
% %Two line segments.
% xA=[[-3;10],[0;0],[4;6],[-2;1]];
% xB=[[12;10],[16;8],[2;7],[-3;3]];
% [xAT,xBT]=clipLine2ConvexPolygon(xA,xB,v);
%
% figure(1)
% clf
% hold on
% plot([v(1,:),v(1,1)],[v(2,:),v(2,1)],'linewidth',2)
% for k=1:size(xA,2)
%     plot([xA(1,k),xB(1,k)],[xA(2,k),xB(2,k)],'-k','linewidth',4)
%
%     if(~isnan(xAT(1,k)))
%         plot([xAT(1,k),xBT(1,k)],[xAT(2,k),xBT(2,k)],'-r','linewidth',2)
%     end
% end
%
%REFERENCES:
%[1] V. Skala, "An efficient algorithm for line clipping by convex
%    polygon," Computer and Graphics, vol. 17, no. 4, pp. 417-421, 1993.
%
%September 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

N=size(v,2);%N is the number of edges and vertices.

%If many lines are clipped at once, precomputing this outside the loop can
%be fastest.
siHat=zeros(2,N);
for k=1:(N-1)
    siHat(:,k)=v(:,k+1)-v(:,k);
end
siHat(:,N)=v(:,1)-v(:,N);

numPts=size(xAL,2);
doesClip=false(numPts,1);

for curPt=1:numPts
    xA=xAL(:,curPt);
    xB=xBL(:,curPt);

    si=bsxfun(@minus,v,xA);
    index=zeros(2,1);
    prevSpecialIdx=-2;

    k=0;
    i=N-1;
    j=0;
    s=xB-xA;
    xi=(si(2,i+1)*s(1)-si(1,i+1)*s(2));
    while((j<N)&&(k<2))
        eta=(si(2,j+1)*s(1)-si(1,j+1)*s(2));
        if(xi*eta<0)%There is an intersection
            index(k+1)=i+1;%The edge having the intersection.
            k=k+1;
            xi=eta;
        elseif(xi*eta==0)%xi and/or eta=0
            if(~(prevSpecialIdx==N-1&&i==0)&&~(prevSpecialIdx==i-1))
                if(eta==0)
                    if(xi~=0)%If not special case m.
                        index(k+1)=i+1;%The edge having the intersection.
                        k=k+1;
                        xi=eta;
                        prevSpecialIdx=i;
                    end
                else%xi==0 and eta is not zero.
                    index(k+1)=i+1;%The edge having the intersection.
                    k=k+1;
                    xi=eta;
                    prevSpecialIdx=i;
                end
            end
        end
        i=j;
        j=j+1;
    end

    if(k==0)
        doesClip(curPt)=false;%There is no intersection.
        xA=NaN;
        xB=NaN;
    else
        doesClip(curPt)=true;%There is an intersection.
        tMin=-Inf;
        tMax=Inf;
        i=index(1);
        t1=det([si(:,i),-siHat(:,i)])/det([s,-siHat(:,i)]);
        i=index(2);
        t2=det([si(:,i),-siHat(:,i)])/det([s,-siHat(:,i)]);

        if(t1<t2)
           temp=t2;
           t2=t1;
           t1=temp;
        end

        if(t2<tMax)
            xB=xA+s*t2;
        end

        if(t1>tMin)
            xA=xA+s*t1;
        end
    end

    xAL(:,curPt)=xA;
    xBL(:,curPt)=xB;
end

end

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


text:
function b = r8blt_to_r8ge ( n, ml, a )

%*****************************************************************************80
%
%% R8BLT_TO_R8GE copies a R8BLT matrix to a R8GE matrix.
%
%  Discussion:
%
%    The R8BLT storage format is used for a banded lower triangular matrix.
%    The matrix is assumed to be zero below the ML-th subdiagonal.
%    The matrix is stored in an ML+1 by N array, in which the diagonal
%    appears in the first row, followed by successive subdiagonals.
%    Columns are preserved.
%
%  Example:
%
%    N = 5, ML = 2
%
%    A11   0   0   0   0
%    A21 A22   0   0   0
%    A31 A32 A33   0   0
%      0 A42 A43 A44   0
%      0   0 A53 A54 A55
%                --- ---
%                    ---
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer ML, the lower bandwidth.
%
%    Input, real A(ML+1,N), the R8BLT matrix.
%
%    Output, real B(N,N), the R8GE matrix.
%
  for i = 1 : n
    for j = 1 : n
       if ( j <= i & i <= j + ml )
          b(i,j) = a(i-j+1,j);
       else
          b(i,j) = 0.0;
       end
    end
  end

  return
end


text:
function tocStatus( id, fracDone )
% Used to display the progress of a long process.
%
% For more information see ticStatus.
%
% USAGE
%  tocStatus( id, fracDone )
%
% INPUTS
%  id          - unique id of progress indicator
%  fracDone    - value in (0,1] indicating percent operation completed
%
% OUTPUTS
%
% EXAMPLE
%
% See also TICSTATUS
%
% Piotr's Image&Video Toolbox      Version 2.50
% Copyright 2012 Piotr Dollar.  [pdollar-at-caltech.edu]
% Please email me if you find bugs, or have suggestions or questions!
% Licensed under the Simplified BSD License [see external/bsd.txt]

global TT_STATUS TT_FREE_IDS

%%% error check
if( length(TT_STATUS)<id || TT_FREE_IDS(id)==1 )
  error('MATLAB:tocStatus:callTicstatusFirst', ...
    'You must call TICSTATUS before calling TOCSTATUS.');
end
[fracDone,er] = checkNumArgs( fracDone, [1 1], -1, 1 );
if(~isempty(er)), error(er); end
if( fracDone>1 ); error(['fracDone: ' num2str(fracDone) ' > 1'] ); end;

%%% get parameters
updateFreq  = TT_STATUS(id).updateFreq;
updateMinT  = TT_STATUS(id).updateMinT;
erasePrev   = TT_STATUS(id).erasePrev;
msg         = TT_STATUS(id).msg;
t0          = TT_STATUS(id).t0;
tLast       = TT_STATUS(id).tLast;
lenPrev     = TT_STATUS(id).lenPrev;

%%% update if enough time has passed
if( etime( clock, tLast )> updateFreq || (fracDone==1 && lenPrev>0) )
  tLast = clock;
  elptime = etime(clock,t0);
  fracDone = max( fracDone, .00001 );
  esttime = elptime/fracDone - elptime;
  if( lenPrev || (elptime/fracDone)>updateMinT )
    if( ~lenPrev ); fprintf('\n'); end

    % create display message
    fracdoneS = num2str(fracDone*100,'%.1f');
    if( elptime/fracDone < 600 )
      elptimeS  = num2str(elptime,'%.1f');
      esttimeS  = num2str(esttime,'%.1f');
      timetypeS = 's';
    else
      elptimeS  = num2str(elptime/60,'%.1f');
      esttimeS  = num2str(esttime/60,'%.1f');
      timetypeS = 'm';
    end
    if( ~isempty(msg) ); msg = [msg '   ']; end
    msg = [msg 'completed=' fracdoneS '%% [elapsed=' elptimeS ];
    msg = [msg timetypeS ' / remaining~=' esttimeS timetypeS ']' ];

    % erase previous display and create new display
    if( erasePrev ) % undo previous disp
      fprintf( repmat('\b', [1 lenPrev] ) ); end
    fprintf( msg );  % fprintf( [msg '\n'] );
    lenPrev = length( msg ) - 1; %note %% (+1 if using \n)
    TT_STATUS(id).tLast = tLast;
    TT_STATUS(id).lenPrev = lenPrev;
  end
end

%%% free id if done
if( fracDone==1 )
  if(lenPrev); fprintf('\n'); end
  TT_FREE_IDS(id) = 1;
end


text:
function [data,stats] = ExpPPHeadless(expdirs,jabfiles,varargin)

[paramsfile,statframelimits] = myparse(varargin,'paramsfile','ExpPP.params.txt',...
  'statframelimits',[-inf,inf]);

data = ExpPP.loadexps(expdirs,jabfiles);
data = ExpPP.ensureCustomGroupInit(data);
params = ReadParams(paramsfile);

statframelimits(end) = min(statframelimits(end),max([data.trxt1]));
computeparams = {'statframelimits' statframelimits};

[data,stats] = ExpPP.computeStats(data,params{:},computeparams{:});

% possibly not the cleanest way to do this, copy-pasted
function params = ReadParams(paramsfile)

fh = fopen(paramsfile);
if isequal(fh,-1)
  error('ExpPPHeadless:errRead','Cannot read postprocessing parameter file.');
end
out = textscan(fh,'%s%s%s%s%s','delimiter',',','collectoutput',true);
fclose(fh);
out = out{1};
assert(isequal(out(1,:),{'name' 'type' 'nanok' 'default' 'desc'}));
out = out(2:end,:);
tbl = struct();
tbl.name = out(:,1);
tbl.type = out(:,2);
tbl.nanok = str2double(out(:,3));
tbl.default = str2double(out(:,4));
tbl.desc = out(:,5);
data = [tbl.name num2cell(tbl.default) tbl.desc];
params = struct();
for i = 1:size(data,1)
  prm = data{i,1};
  val = data{i,2};

  assert(strcmp(prm,tbl.name{i}));

  % check val
  if isnan(val)
    if tbl.nanok(i)
      % none
    else
      warningNoTrace('ExpPPControl:nanVal',...
        'Parameter ''%s'' cannot be NaN. Using default value of %d.\n',...
        prm,tbl.default(i));
      val = tbl.default(i);
    end
  else
    switch tbl.type{i}
      case 'nat'
        if floor(val)~=val
          warningNoTrace('ExpPPControl:fracVal',...
            'Parameter ''%s'' should be an integer. Rounding to %d.\n',prm,round(val));
          val = round(val);
        end
        if val<0
          warningNoTrace('ExpPPControl:negVal',...
            'Parameter ''%s'' cannot be negative. Using 0.\n',prm);
          val = 0;
        end
      otherwise
        assert(false,'All types currently nat');
    end
  end

  % assign
  subparams = regexp(prm,'#','split');
  switch numel(subparams)
    case 1
      params.(subparams{1}) = val;
    case 2
      params.(subparams{1}).(subparams{2}) = val;
    otherwise
      assert(false,'Unexpected parameter.');
  end
end

params = struct2paramscell(params);

author: kristinbranson

repo: JAABA


text:
% CLASSIF_ERR Calculates the classification error.
%
% Usage
%    err = CLASSIF_ERR(labels, test_set, src)
%
% Input
%    labels (int): The predicted labels corresponding to the testing in-
%       stances.
%    test_set (int): The object indices of the testing instances.
%    src (struct): The source from which the objects were taken.
%
% Output
%    err (numeric): The classification error.
%
% See also
%    AFFINE_TEST, SVM_TEST

function err = classif_err(labels,test_set,src)
	truth = [src.objects(test_set).class];

	if isfield(src,'cluster')
		cluster = src.cluster;
	else
		cluster = 1:max(truth);
	end

	err = 1-sum(bsxfun(@eq,cluster(labels),cluster(truth)),2)/length(truth);
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function p = addFromProbedBinary(p)

% For a nonlinear operator, the envelope model is tighter on a smaller
% domain. Search for (binary d implies x >= L) and then compute bounds on
% f(x) on that domain, and add implies(d, L <= f(x) <= U)
if ~isempty(p.binary_variables) && ~isempty(p.F_struc)
    newCuts = zeros(0,size(p.F_struc,2));
    b = p.F_struc(:,1);
    A = p.F_struc(:,2:end);
    D = A(:,p.binary_variables);
    A(:,p.binary_variables) = 0;
    probedBoundsL = [];
    probedBoundsU = [];
    for i = p.K.f + 1:p.K.f + p.K.l
        [~,ix,valx] = find(A(i,:));
        if length(ix) == 1
            [~,id,vald] = find(D(i,:));
            if length(id) == 1
                if valx > 0
                    % b + valx*x + vald*d >= 0
                    % valx*x >= (-b-vald*x)
                    % d implies x >= (-b-vald)/valx
                    activated_lower = (-b(i)-vald)/valx;
                    probedBoundsL(end+1,1) = id;
                    probedBoundsL(end,2) = ix;
                    probedBoundsL(end,3) = activated_lower;
                elseif valx < 0
                    % b - (-valx)*x + vald*d >= 0
                    % b + vald*d >= -valx*x
                    % d implies (b+vald/(-valx)
                    activated_upper = (b(i)+vald)/(-valx);
                    probedBoundsU(end+1,1) = id;
                    probedBoundsU(end,2) = ix;
                    probedBoundsU(end,3) = activated_upper;
                end
            end
        end
    end
    if ~isempty(probedBoundsL) || ~isempty(probedBoundsU)
        for i = 1:length(p.binary_variables)
            iL = find(probedBoundsL(:,1) == i);
            iU = find(probedBoundsU(:,1) == i);
            if ~isempty(iL) || ~isempty(iU)
                ptemp = p;
                if ~isempty(iL)
                    for k = iL
                        ptemp.lb(probedBoundsL(k,2)) = probedBoundsL(k,3);
                    end
                end
                if ~isempty(iU)
                    for k = iU
                        ptemp.ub(probedBoundsU(k,2)) = probedBoundsU(k,3);
                    end
                end
                ptemp.lb(p.binary_variables(i)) = 1;
                pprobed = update_monomial_bounds(ptemp);
                pprobed = propagate_bounds_from_evaluations(pprobed);
                sL = find(pprobed.lb > ptemp.lb+1e-6);
                sU = find(pprobed.ub < ptemp.ub-1e-6);
                for j = 1:length(sL)
                    % {d implies x(s) >= q}
                    % x(s) >= q*d + L*(1-d)
                    % -L + x(s) + d*(L-q) >= 0
                    if ~ismember( sL(j),p.binary_variables)   % {d implies x(s) <= q}
                        q = pprobed.lb(sL(j));
                        L = p.lb(sL(j));
                        if ~isinf(q) && ~isinf(L)
                            newCuts(end+1,1) = -L;
                            newCuts(end,1+sL(j)) = 1;
                            newCuts(end,1+p.binary_variables(i)) = L-q;
                        end
                    end
                end
                for j = 1:length(sU)
                    if ~ismember( sU(j),p.binary_variables)   % {d implies x(s) <= q}
                        % x(s) <= q*d + U*(1-d)
                        % U - x(s) + d*(q-U) >= 0
                        q = pprobed.ub(sU(j));
                        U = p.ub(sU(j));
                        if ~isinf(q) && ~isinf(U)
                            newCuts(end+1,1) = U;
                            newCuts(end,1+sU(j)) = -1;
                            newCuts(end,1+p.binary_variables(i)) = q-U;
                        end
                    end
                end
            end
        end
    end
    if nnz(newCuts)>0
        p.F_struc = [p.F_struc(1:p.K.f+p.K.l,:);
            newCuts;
            p.F_struc(p.K.f+p.K.l+1:end,:)];
        p.K.l = p.K.l + size(newCuts,1);
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [p,t]=fixmesh(p,t)
%FIXMESH  Remove duplicated/unused nodes and fix element orientation.
%   [P,T]=FIXMESH(P,T)

%   Copyright (C) 2004-2006 Per-Olof Persson. See COPYRIGHT.TXT for details.

snap=max(max(p,[],1)-min(p,[],1),[],2)*1024*eps;
[foo,ix,jx]=unique(round(p/snap)*snap,'rows');
p=p(ix,:);
t=jx(t);

[pix,ix,jx]=unique(t);
t=reshape(jx,size(t));
p=p(pix,:);

flip=simpvol(p,t)<0;
t(flip,[1,2])=t(flip,[2,1]);

repo: ifem


text:
function mdsfig_3d(X,names,clus,linemat,cordir)
% ::
%
%    mdsfig_3d(X,names,clus,linemat,cordir)
%
% 3-D MDS figure
%
% Subfunction of mdsfig used if 3+ dims are available

numel=size(X,1); %%% define scaling %%%%%%%
a=X(:,1);b=X(:,2);c=X(:,3);
axa=ones(numel,1)*mean(a);
axb=ones(numel,1)*mean(b);
axc=ones(numel,1)*mean(c);
axa1=min(a):(max(a)-min(a))/(numel-1):max(a);
axb1=min(b):(max(b)-min(b))/(numel-1):max(b);
axc1=min(c):(max(c)-min(c))/(numel-1):max(c);

%%% plotting defaults %%%%
f1=figure('Color','w');
hold on;
set(gca,'FontSize',12);
colors = {'ro' 'go' 'bo' 'yo' 'co' 'mo' 'ko' 'r^' 'g^' 'b^' 'y^' 'c^' 'm^' 'k^'};
while max(clus) > length(colors), colors = [colors colors];,end

%draw central axes & scale
plot3(axa,axb,axc1,'color','k','linewidth',2);
plot3(axa,axb1,axc,'color','k','linewidth',2);
plot3(axa1,axb,axc,'color','k','linewidth',2);
axis([min(a) max(a) min(b) max(b) min(c) max(c)]);

%%% plot lines %%%%%%

drawlines(X,linemat,0,{'Positive' 'Negative'});

%%%% plot points %%%%
for j = 1:size(X,1);
    plot3(a(j),b(j),c(j),colors{clus(j)},'MarkerFaceColor',colors{clus(j)}(1),'Markersize',8);
end

if ~isempty(names);
for i = 1:size(X,1)
    text(X(i,1)+.025*[max(X(:,1))-min(X(:,1))],X(i,2)+.0,X(i,3)+.0,names{i},'Color','k','FontSize',12,'FontWeight','bold');
end
end

rotate=1;
if rotate
% rotate until closed
for v=1:45;
    if ishandle(f1);
        view(v,35);
        drawnow;
        pause (0.2)
        if v==44
            v==1;
            continue
        end
    else
        break
    end
end
end

% Second Figure -- different dims

tor_fig(2,2);hold off;
subplot(2,2,1);
mdsfig(X(:,[1 2]),names,clus,linemat);
subplot(2,2,2);
mdsfig(X(:,[1 3]),names,clus,linemat);
subplot(2,2,3);
mdsfig(X(:,[2 3]),names,clus,linemat);

subplot(2,2,4);
if size(X,2) > 3,
    mdsfig(X(:,[1 4]),names,clus,linemat);
else
    axis off
end

return

function [hhp,hhn] = drawlines(pc,sigmat,sigcol,legmat,varargin)

% sigcol is a scalar between zero and one, lower is 'more salient'
hold on

% linew
lw = 3;

color(1,:) = [0 0 0];   % first line color, 'positive'
color(2,:) = [0 .7 1];   % second line color, 'negative'
style{1} = '-';         % first line style
style{2} = '-';         % second line style

if length(varargin) > 0, color = varargin{1}; end
if length(varargin) > 1, style = varargin{2}; sigcol = 0; end

hhp=[];hhn=[];
for i=1:size(pc,1);
     for j=1:size(pc,1);
         if sigmat(i,j) > 0
             hhp = line([pc(i,1) pc(j,1)],[pc(i,2) pc(j,2)]);
             set(hhp,'Color',color(1,:) + [1 1 1] * abs(sigcol),'LineStyle',style{1},'LineWidth',lw - (lw*sigcol)-1)
         elseif sigmat(i,j) < 0
             hhn = line([pc(i,1) pc(j,1)],[pc(i,2) pc(j,2)]);
             set(hhn,'Color',color(2,:) + abs([sigcol sigcol 0]),'LineStyle',style{2},'LineWidth',lw - (lw*sigcol)-1)
         end
     end
end
legend ([hhp hhn],legmat);
return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function x = r8blt_sl ( n, ml, a, b, job )

%*****************************************************************************80
%
%% R8BLT_SL solves a R8BLT system.
%
%  Discussion:
%
%    The R8BLT storage format is appropriate for a banded lower triangular matrix.
%    The matrix is assumed to be zero below the ML-th subdiagonal.
%    The matrix is stored in an ML+1 by N array, in which the diagonal
%    appears in the first row, followed by successive subdiagonals.
%    Columns are preserved.
%
%    No factorization of the lower triangular matrix is required.
%
%  Example:
%
%    N = 5, ML = 2
%
%    A11   0   0   0   0
%    A21 A22   0   0   0
%    A31 A32 A33   0   0
%      0 A42 A43 A44   0
%      0   0 A53 A54 A55
%                --- ---
%                    ---
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 February 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, integer ML, the lower bandwidth.
%
%    Input, real A(ML+1,N), the R8BLT matrix.
%
%    Input, real B(N), the right hand side.
%
%    Input, integer JOB, is 0 to solve the untransposed system,
%    nonzero to solve the transposed system.
%
%    Output, real X(N), the solution vector.
%
  x(1:n) = b(1:n);

  if ( job == 0 )

    for j = 1 : n
      x(j) = x(j) / a(1,j);
      ihi = min ( j + ml, n );
      for i = j+1 : ihi
        x(i) = x(i) - a(i-j+1,j) * x(j);
      end
    end

  else

    for j = n : -1 : 1
      x(j) = x(j) / a(1,j);
      ilo = max ( j - ml, 1 );
      for i = ilo : j-1
        x(i) = x(i) - a(j-i+1,i) * x(j);
      end
    end

  end

  return
end


text:
classdef CEC2017_F2 < PROBLEM
% <single> <real> <constrained>
% CEC'2017 constrained optimization benchmark problem

%------------------------------- Reference --------------------------------
% G. Wu, R. Mallipeddi, and P. N. Suganthan, Problem definitions and
% evaluation criteria for the CEC 2017 competition on constrained real-
% parameter optimization, National University of Defense Technology, China,
% 2016.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    properties
        O;      % Optimal decision vector
        Mat;	% Rotation matrix
    end
    methods
        %% Default settings of the problem
        function Setting(obj)
            CallStack = dbstack('-completenames');
            load(fullfile(fileparts(CallStack(1).file),'CEC2017.mat'),'Data');
            obj.O = Data{2}.o;
            obj.M = 1;
            if isempty(obj.D) || obj.D < 30
                obj.D   = 10;
                obj.Mat = Data{2}.M_10;
            elseif obj.D < 50
                obj.D   = 30;
                obj.Mat = Data{2}.M_30;
            elseif obj.D < 100
                obj.D   = 50;
                obj.Mat = Data{2}.M_50;
            else
                obj.D   = 100;
                obj.Mat = Data{2}.M_100;
            end
            obj.lower    = zeros(1,obj.D) - 100;
            obj.upper    = zeros(1,obj.D) + 100;
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            PopObj = sum(cumsum(Z,2).^2,2);
        end
        %% Calculate constraint violations
        function PopCon = CalCon(obj,PopDec)
            Z = PopDec - repmat(obj.O(1:size(PopDec,2)),size(PopDec,1),1);
            Y = Z*obj.Mat';
            PopCon = sum(Y.^2-5000*cos(0.1*pi*Y)-4000,2);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function Model = verifyModelBorg(Model, varargin)
% Ensures that a model is in the correct format to be analyzed
% by any of the scripts in the `Tmodel` suite. It will add fields that are
% missing and expected for comparison. It will remove fields not in this list.
% Called by  `driveModelBorgifier`, calls `TmodelFields`, `fixNames`, `removeDuplicateNames`,
% `makeNamesUnique`, `buildRxnEquations`, `fixChemFormulas`, `orderModelFieldsBorg`, `organizeModelCool`.
%
% USAGE:
%
%    Model = verifyModelBorg(Model)
%
% INPUTS:
%    Vmodel:        Model from `readCbModel` or any of the `readModel` functions.
%
% OPTIONAL INPUTS:
%    'keepName':    Don't ask for verification of model name
%    'Verbose':     Print steps.
%
% OUTPUTS:
%    Model:         Model with additional fields and correct format.
%
% Please cite:
% `Sauls, J. T., & Buescher, J. M. (2014). Assimilating genome-scale
% metabolic reconstructions with modelBorgifier. Bioinformatics
% (Oxford, England), 30(7), 1036?8`. http://doi.org/10.1093/bioinformatics/btt747
%
% ..
%    Edit the above text to modify the response to help addMetInfo
%    Last Modified by GUIDE v2.5 06-Dec-2013 14:19:28
%    This file is published under Creative Commons BY-NC-SA.
%
%    Correspondance:
%    johntsauls@gmail.com
%
%    Developed at:
%    BRAIN Aktiengesellschaft
%    Microbial Production Technologies Unit
%    Quantitative Biology and Sequencing Platform
%    Darmstaeter Str. 34-36
%    64673 Zwingenberg, Germany
%    www.brain-biotech.de

askForModelName = true ; % Declare variables.
verbose = false ;
if ~isempty(varargin)
    if sum(strcmp('keepName', varargin))
        askForModelName = false ;
    end
    if sum(strcmp('Verbose', varargin))
        verbose = true ;
    end
end

nRxns = length(Model.rxns);
nMets = length(Model.mets);

% Current Field names in the model.
fieldNames = fieldnames(Model);

% Get desired field names.
fields = TmodelFields ;
Field.rxn = fields{1};
Field.rNum = fields{2};
Field.met = fields{3};
Field.mNum = fields{4};
Field.all = fields{5};

%% Remove fields that will not be used
% This is a hacky solution, and unused fields should be carried forward
rmFieldNames = false(size(fieldNames)) ;
for iField = 1:length(fieldNames)
    if strmatch(fieldNames{iField}, Field.all)
        continue
    else
        rmFieldNames(iField) = 1 ;
    end
end
Model = rmfield(Model, fieldNames(rmFieldNames)) ;

%% Pad fields with empty strings or zeros.
% Reaction related cell arrays.
for iField = 1:length(Field.rxn)
   fieldIndex = strcmp(Field.rxn{iField}, fieldNames);
   fieldIndex = find(fieldIndex);
   if isempty(fieldIndex) % Field does not currently exist
       if verbose
           fprintf(['Array .' Field.rxn{iField} ' not in Model. Adding.\n'])
       end
       Model.(Field.rxn{iField}) = cell(nRxns, 1);
       Model.(Field.rxn{iField})(:) = {''};
   else % Field exists, check if the length is correct
       vFieldLength = length(Model.(fieldNames{fieldIndex}));
       for icell = 1:vFieldLength
           while iscell(Model.(fieldNames{fieldIndex}){icell})
               Model.(fieldNames{fieldIndex}){icell} = [Model.(fieldNames{fieldIndex}){icell}{:}] ;
           end
       end
       if vFieldLength ~= nRxns
           vFieldLength = vFieldLength + 1 ;
           Model.(fieldNames{fieldIndex})(vFieldLength:nRxns) = {''};
       end
   end
end

% Reaction related double arrays.
for iField = 1:length(Field.rNum)
   fieldIndex = strcmp(Field.rNum{iField}, fieldNames);
   fieldIndex = find(fieldIndex);
   if isempty(fieldIndex) % Field does not currently exist
       if verbose
           fprintf(['Array .' Field.rNum{iField} ' not in Model. Adding.\n'])
       end
       Model.(Field.rNum{iField}) = zeros(nRxns, 1);
   else % Field exists, check if the length is correct
       vFieldLength = length(Model.(fieldNames{fieldIndex}));
       if vFieldLength ~= nRxns
           vFieldLength = vFieldLength + 1 ;
           Model.(fieldNames{fieldIndex})(vFieldLength:nRxns) = 0 ;
       end
   end
end

% Metabolite related cell arrays.
for iField = 1:length(Field.met)
   fieldIndex = strcmp(Field.met{iField}, fieldNames);
   fieldIndex = find(fieldIndex);
   if isempty(fieldIndex) % Field does not currently exist
       if verbose
           fprintf(['Array .' Field.met{iField} ' not in Model. Adding.\n'])
       end
       Model.(Field.met{iField}) = cell(nMets, 1);
       Model.(Field.met{iField})(:) = {''};
   else % Field exists, check if the length is correct
       vFieldLength = length(Model.(fieldNames{fieldIndex}));
       % weed out cells in cells
       for icell = 1:vFieldLength
           while iscell(Model.(fieldNames{fieldIndex}){icell})
               Model.(fieldNames{fieldIndex}){icell} = [Model.(fieldNames{fieldIndex}){icell}{:}] ;
           end
       end
       if vFieldLength ~= nMets
           vFieldLength = vFieldLength + 1 ;
           Model.(fieldNames{fieldIndex})(vFieldLength:nMets) = {''};
       end
   end
end

% Metabolite related double arrays.
for iField = 1:length(Field.mNum)
   fieldIndex = strcmp(Field.mNum{iField}, fieldNames);
   fieldIndex = find(fieldIndex);
   if isempty(fieldIndex) % Field does not currently exist
       if verbose
           fprintf(['Array .' Field.mNum{iField} ' not in Model. Adding.\n'])
       end
       Model.(Field.mNum{iField}) = zeros(nMets, 1);
   else % Field exists, check if the length is correct
       vFieldLength = length(Model.(fieldNames{fieldIndex}));
       if vFieldLength ~= nMets
           vFieldLength = vFieldLength + 1 ;
           Model.(fieldNames{fieldIndex})(vFieldLength:nMets) = 0 ;
       end
   end
end

%% Ensure reactions are forwards.
if verbose
    fprintf('Making sure reactions are all forwards\n')
end
% Find reverse reactions based on bounds.
revRxns = find(abs(Model.lb) > Model.ub);

% Do it.
for iRxn = 1:length(revRxns)
    % Reverse bounds.
    [Model.lb(revRxns(iRxn)), Model.ub(revRxns(iRxn))] = ...
        deal(-Model.ub(revRxns(iRxn)), -Model.lb(revRxns(iRxn)));

    % Change sign of stochiometrix matrix.
    metStoics = find(Model.S(:, revRxns(iRxn)));
    for iMet = 1:length(metStoics)
        if Model.S(metStoics(iMet), revRxns(iRxn)) > 0
            Model.S(metStoics(iMet), revRxns(iRxn)) = ...
                0 - Model.S(metStoics(iMet), revRxns(iRxn));
        else
            Model.S(metStoics(iMet), revRxns(iRxn)) = ...
                abs(Model.S(metStoics(iMet), revRxns(iRxn)));
        end
    end
end

%% Remove hyphens and obtuse characters from names.
if verbose
    fprintf('Fixing names of metabolites and reaction\n')
end
Model.mets = fixNames(Model.mets);
Model.metNames = fixNames(Model.metNames);
Model.metNames = removeDuplicateNames(Model.metNames);
Model.rxns = fixNames(Model.rxns);
Model.rxnNames = removeDuplicateNames(Model.rxnNames);

%% Make sure all reaction and metabolite names are unique.
if verbose
    fprintf('Checking if reaction IDs (.rxns) are unique.\n');
end
if length(Model.rxns) ~= length(unique(Model.rxns))
    fprintf('ERROR: Not all reactions are unique.\n')
    % Launch name correcting script.
    Model.rxns = makeNamesUnique(Model.rxns, Model.rxnNames);
end
if verbose
    fprintf('Checking if metabolite IDs (.mets) are unique.\n');
end
if length(Model.mets) ~= length(unique(Model.mets))
    fprintf('ERROR: Not all metabolites are unique.\n')
    Model.mets = makeNamesUnique(Model.mets, Model.metNames);
end

%% Make sure all metabolites have a compartment.
noComp = find(cellfun(@isempty, regexp(Model.mets, '\[\w\]$')));

if isempty(noComp)
    if verbose
        fprintf('All metabolites have comparment designation.\n')
    end
else
    % try to get compartment information from metNames
    nameComp = find(~cellfun(@isempty, regexp(Model.metNames, '\[\w\]$')));
    if ~isempty(nameComp)
        for inc = rowvector(nameComp)
            Model.mets{inc} = [Model.mets{inc} ...
                Model.metNames{inc}(strfind(Model.metNames{inc}, '['): ...
                strfind(Model.metNames{inc}, ']')) ];
            Model.metNames{inc} = ...
                Model.metNames{inc}(1:strfind(Model.metNames{inc}, '[')-1);
        end
    end

    noComp = find(cellfun(@isempty, regexp(Model.mets, '\[\w\]$')));

    if isempty(noComp)
        if verbose
            fprintf('All metabolites have comparment designation.\n')
        end
    else
        if verbose
            fprintf(['Some metabolites have no compartment designation. ' ...
            'Assigning as cytosolic.\n'])
        end
        for iMet = 1:length(noComp)
            Model.mets{noComp(iMet)} = [Model.mets{noComp(iMet)} '[c]'];
        end
    end
end

%% Give .rxnID or .metID names from .rxns and .mets.
Model.rxnID = Model.rxns ;
Model.metID = Model.mets ;

%% If model has SEED style reaction or metabolite IDs, add them to arrays.
% Reactions.
areSEEDids = ~cellfun(@isempty, regexp(Model.rxns, '^rxn\d{5}$'));
Model.rxnSEEDID(areSEEDids) = Model.rxns(areSEEDids);

% Metabolites.
noComp = Model.mets;
for iMet = 1:length(noComp)
    noComp{iMet} = noComp{iMet}(1:end - 3);
end
areSEEDids = ~cellfun(@isempty, regexp(noComp, '^cpd\d{5}$'));
Model.metSEEDID(areSEEDids) = noComp(areSEEDids);

%% Format chemical formulas
Model.metFormulas = fixChemFormulas(Model.metFormulas);

%% Rebuild reaction equations to ensure they use fixed names/abbreviations.
Model = buildRxnEquations(Model);

%% Ensure all vectors are column vectors.
colFields = {'rxn' 'rNum' 'met' 'mNum'};
for iName = 1:length(colFields)
    for iField = 1:length(Field.(colFields{iName}))
        Model.(Field.(colFields{iName}){iField}) = ...
            Model.(Field.(colFields{iName}){iField})(:);
    end
end

%% Make sure there is a model name and if not, ask for one.
if ~isfield(Model, 'description')
    needModelName = 1 ;
else
    if askForModelName
        fprintf(['Current model name is:' char(10) Model.description char(10)])
        answer = input('Keep name? (y/n): ', 's');
        if strcmpi(answer,'y') || strcmpi(answer, 'yes')
            fprintf('Keeping name\n')
            needModelName = 0 ;
        else
            needModelName = 1 ;
        end
    else
        needModelName = 0 ;
    end
end
try
    test.(Model.description) = true  ;
catch
    disp('Sorry, this model name is invalid. (Hint: it must start with a letter and only contain letters and numbers)')
    needModelName = 1 ;
end
while needModelName
    prompt = 'Input model name: ';
    modelName = input(prompt, 's');
    Model.description = modelName ;
    try
        test.(Model.description) = true  ;
        needModelName = 0 ;
    catch
        disp('Sorry, this model name is invalid. (Hint: it must start with a letter and only contain letters and numbers)')
    end
end

%% Reorder fields and organize model based on most common mets.
Model = orderModelFieldsBorg(Model);
Model = organizeModelCool(Model);

author: opencobra

repo: cobratoolbox


text:
function pass = test_coeffs_vals( pref )

% Grab some preferences
if ( nargin == 0 )
    pref = chebfunpref();
end
tol = 1e3*pref.techPrefs.chebfuneps;
j = 1;
%grab some coeffs
u = randnfunsphere(.3);
v = randnfunsphere(.4);
w = randnfunsphere(.3);
f = [u;v;w];

% test coeffs2:
[x, y, z] = coeffs2(f);
pass(j) = norm(spherefun.coeffs2spherefun(x)-u, inf) < tol;
pass(j+1) = norm(spherefun.coeffs2spherefun(y)-v, inf) < tol;
pass(j+2) = norm(spherefun.coeffs2spherefun(z)-w, inf) < 10*tol;
j = j+3;

%parameters
[x, y, z] = coeffs2(f, 50, 60);
pass(j) = norm(x-coeffs2(u, 50, 60), inf) < tol;
pass(j+1) = norm(y-coeffs2(v, 50, 60), inf) < tol;
pass(j+2) = norm(z-coeffs2(w, 50, 60), inf) < tol;
j = j+3;

% test coeffs2diskfunv:
f2 = spherefunv.coeffs2spherefunv(coeffs2(u), coeffs2(v),...
    coeffs2(w));
pass(j) = norm(f - f2) < tol;
j = j+1;

%test coeffs2vals:
[x, y, z] = coeffs2(f2);
[u, v, w] = spherefunv.coeffs2vals(x,y,z);
pass(j) = norm(spherefun.coeffs2vals(x)-u, 'inf') < tol;
pass(j+1) = norm(spherefun.coeffs2vals(y)-v, 'inf')< tol;
pass(j+2) = norm(spherefun.coeffs2vals(z)-w, 'inf')< tol;
j= j+2;

%testvals2coeffs:
[a,b,c] = spherefunv.vals2coeffs(u,v, w);
pass(j) = norm(x - a, 'inf');
pass(j+1) = norm(y - b, 'inf');
pass(j+2) = norm(z - c, 'inf');

if (nargout > 0)
    pass = all(pass(:));
end
end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%IM2OBJ Convert Matlab images or datafile to dataset object
%
%  B = IM2OBJ(IM,A)
%  B = IM2OBJ(IM,FEATSIZE)
%
% INPUT
%  IM        X*Y image, X*Y*C image, X*Y*K array of K images,
%            X*Y*C*K array of color images, or cell-array of images
%            The images may be given as a datafile.
%  A         Input dataset
%  FEATSIZE  Vector with desired feature sizes to solve ambiguities
%
% OUTPUT
%  B   Dataset with IM added
%
% DESCRIPTION
% Add standard Matlab images, as objects, to an existing dataset A. If A is
% not given, a new dataset is created. Images of type 'uint8' are converted
% to 'double' and divided by 256. The resulting feature size is X*Y or
% X*Y*C. The set of images IM may be given as a datafile.
%
% SEE ALSO (PRTools Guide)
% DATASETS, DATAFILES, IM2FEAT, FEATIM

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: im2obj.m,v 1.5 2008/07/28 08:59:47 duin Exp $

function a = im2obj (im,a)

	if iscell(im)
		imsize = size(im{1});
	else
		imsize = size(im);
	end

	nodataset = 0;
	if (nargin > 1)
		if (~isdataset(a)) & size(a,1) ~= 1
			error('second argument should be a dataset or a size vector');
		end
		if isdataset(a)
			featsize = getfeatsize(a);
		else
			featsize = a;
			nodataset = 1;
		end
	else % no information on feature size, assume most simple solution
		nodataset = 1;
		if length(imsize) == 2
			featsize = imsize;
		else
			featsize = imsize(1:end-1);
		end
	end
	if length(featsize) == length(imsize)
		nobj = 1;
	elseif length(featsize) == (length(imsize)-1)
		nobj = imsize(end);
	elseif length(featsize) == (length(imsize)-2) & imsize(3) == 1
		nobj = imsize(end);
	else
		wrongfeatsize;
	end

	if any(featsize ~= imsize(1:length(featsize)))
		wrongfeatsize;
	end

	if (isa(im,'cell'))

		% If IM is a cell array of images, unpack it and recursively call IM2OBJ
		% to add each image.

		im = im(:);																	% Reshape to 1D cell array.
		for i = 1:length(im)
			b = feval(mfilename,im{i});
			if ~isempty(a) & any(a.featsize ~= b.featsize)
				error('Images should have equal sizes')
			end
			a = [a; b];
		end

  elseif isdatafile(im)

		if (nargin < 2)
			a = prdataset([]);
		end
    testdatasize(im);
    a = [a; feval(mfilename,data2im(im))];

	else

		% If IM is an image or array of images, reshape it and add it in one go.

		% Convert to double, if necessary
		if (isa(im,'uint8'))
			im = double(im)/256;
		else
			im = double(im);
		end

		% ready for the real work, at last!

		if nobj == 1
			im = im(:)';
		else
			im = shiftdim(im,ndims(im)-1);
			im = reshape(im,nobj,prod(featsize));
		end

		if nodataset
			a = prdataset(im);
			a = setfeatsize(a,featsize);
		else
			a = [a; im];
		end

	end

return

function wrongfeatsize
	error('Desired feature size and size of supplied image array are inconsistent')
return

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function scansV = getScanAt(scanNum, xV, yV, zV, planC)
%"getScanAt"
%   Returns the scan values at the points defined by vectors xV, yV, zV, in the
%   requested scanNum.  Uses linear interpolation with zeros returned for
%   x,y,z points that are out of the scan bounds.
%
%   Any transM values are NOT included.  If the transformed scan is wanted,
%   pass xV,yV,zV through applyTransM before calling getScansAt.
%
%   If planC isn't specified, the global planC is used.
%
%LM: 02 May 03, JOD, first version.  Based on getDVH.m.
%    24 Feb 05, JRA, Modified header, now using some subfunctions.
%                    Also updated for new finterp3.
%
%Usage:
%   scansV = getScanAt(scanNum, xV, yV, zV, planC)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

indexS = planC{end};

[xVS, yVS, zVS] = getScanXYZVals(planC{indexS.scan}(scanNum));

%Use scan access function in case of remote variables.
% sA = getScanArray(planC{indexS.scan}(scanNum));

sA = planC{indexS.scan}(scanNum).scanArray;

%Interpolate to values in xV/yV/zV from the scanarray, 0 if out of bounds.
tol = 1e6*eps;
if numel(size(sA)) > 2
    scansV = finterp3(xV, yV, zV, sA, [xVS(1)-tol xVS(2)-xVS(1) xVS(end)+tol], [yVS(1)+tol yVS(2)-yVS(1) yVS(end)-tol], zVS, 0);
else
    scansV = finterp2(xVS,yVS,single(sA),xV,yV,0);
end
% scansV = interp3(xVS, yVS, zVS, sA , xV, yV, zV);

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ records, spectra, rawSpectra ] = hyperReadSpecpr( filename )
%HYPERREADSPECPR Reads USGS Specpr files.
%   hyperReadSpecpr reads USGS Specpr files.
%
% Usage
%   [ records, spectra, rawSpectra ] = hyperReadSpecpr( filename )
% Inputs
%   filename - Input filename
% Outputs
%   records - Individual records
%   spectra - The spectra post-processed
%   rawSpectra - The raw spectra
%
% References
%   http://speclab.cr.usgs.gov/specpr-format.html.
dbstop if error;
f = fopen(filename, 'r');
if (f == -1)
    error(sprintf('Failed to open file: %s'), filename);
end

% Ignore the first record.
r = uint32(fread(f, 1536, 'uint8'));

firstTime = 1;
records = {};
spectra = [];
rawSpectra = {};
numRawSpectra = 0;
numRecords = 0;

done = 0;
while (not(done))
    r = uint32(fread(f, 1536/4, 'uint32', 'ieee-be'));
    if (length(r) == 0)
        break;
    end
    numRecords
    r = swapbytes(r);
    r = typecast(r, 'uint8');
    % First two bits of file.  I am making this verbose here so it is clear what I
    % am doing.
    firstTwoBits = dec2bin(bitand(r(4), 3));

    %  Parse first two bits.
    if (firstTwoBits == '10')
        % This is a text record. Skip.
        numRecords = numRecords + 1;
    elseif (firstTwoBits == '0')
        % This is an actual (initial) data record.
        if (not(firstTime))
            numRecords = numRecords+1;
            records{record.irecno} = record;
        end
        record = [];
        data = [];
        firstTime = 0;
        iband = int32(zeros(2, 1));
        record.ititl = char(r(5:44)).';
        record.usernm = char(r(45:52)).';
        iscta = typecast(r(53:56), 'int32');
        isctb = typecast(r(57:60), 'int32');
        jdatea = typecast(r(61:64), 'int32');
        jdateb = typecast(r(65:68), 'int32');
        istb = typecast(r(69:72), 'int32');
        isra = typecast(r(73:76), 'int32');
        isdec = typecast(r(77:80), 'int32');
        record.itchan = swapbytes(typecast(r(81:84), 'int32'));
        irmas = typecast(r(85:88), 'int32');
        revs = typecast(r(89:92), 'int32');
        iband(1) = typecast(r(93:96), 'int32');
        iband(2) = typecast(r(97:100), 'int32');
        record.irwav = swapbytes(typecast(r(101:104), 'int32'));
        record.irespt = swapbytes(typecast(r(105:108), 'int32'));
        record.irecno = swapbytes(typecast(r(109:112), 'int32'));
        itpntr = typecast(r(113:116), 'int32');
        ihist = char(r(117:176)).';
        mhist = char(r(177:472)).';
        nruns = typecast(r(473:476), 'int32');
        siangl = typecast(r(477:480), 'int32');
        seangl = typecast(r(481:484), 'int32');
        sphase = typecast(r(485:488), 'int32');
        iwtrns = typecast(r(489:492), 'int32');
        itimch = typecast(r(493:496), 'int32');
        xnrm = typecast(r(497:500), 'int32');
        scatim = typecast(r(501:504), 'int32');
        timint = typecast(r(505:508), 'int32');
        tempd = typecast(r(509:512), 'int32');
        data = swapbytes(typecast(r(513:1536), 'single'));
        % Remove null data samples. Set to zero instead of -1.23e34.
        data(find(data < -1e34)) = 0;
        record.data = data;
    elseif (firstTwoBits == '1')
        % Continuation of data values.
        cData = swapbytes(typecast(r(5:1536), 'single'));
        cData(find(cData < -1e34)) = 0;
        data = [data; cData];
        record.data = data;
    else
        numRecords = numRecords + 1;
    end
end

% Convert to an array of signatures.
% Resample to model AVIRIS sensor
high = 2.40;
low = 0.4;
numBands = 224;
%d.data = sortrows(d.data, 1);
%[q, w, r ]= unique(d.data(:,1));
%d.data = d.data(w, :);
%lambda = d.data(:, 1);
%reflectance = d.data(:, 2);
s = length(records);
numSpectra = 0;
for q=1:s
    if (isempty(records{q}))
        continue;
    end
    if (records{q}.irwav == 0)
        continue;
    end

    if (not(isempty(strfind(records{q}.ititl, 'error'))))
        continue;
    end
    if (not(isempty(strfind(records{q}.ititl, 'Error'))))
        continue;
    end
    if (not(isempty(strfind(records{q}.ititl, 'Bandpass'))))
        continue;
    end
    if (not(isempty(strfind(records{q}.ititl, 'Wavelengths'))))
        continue;
    end
    % Find wavelengths
    if (isempty(records{records{q}.irwav}))
        continue;
    end
    lambdas = records{records{q}.irwav}.data;
    spectrum = records{q}.data;
    if (length(lambdas) ~= length(spectrum))
        fprintf('Error %d !!!\n', q);
        continue;
    end
    numRawSpectra = numRawSpectra + 1;
    rawSpectra{numRawSpectra}.name = records{q}.ititl;
    rawSpectra{numRawSpectra}.wavelengths = lambdas;
    rawSpectra{numRawSpectra}.reflectance = spectrum;
    goodIdx = find(lambdas > 0);
    lambdas = lambdas(goodIdx);
    spectrum = spectrum(goodIdx);
    % Ensure we have proper lower and upper bounds.
    if (lambdas(1) > low)
        %fprintf('Bad file: Lower wavelength value missing.');
        lambdas = [low; lambdas];
        spectrum = [spectrum(1); spectrum];
    end
    if (lambdas(end) < high)
        %fprintf('Bad file: Upper wavelength value missing.\n');
        %d.data = [];
        %fclose(fid);
        %return;
        lambdas = [lambdas; high];
        spectrum = [spectrum; spectrum(end)];
    end
    % Resample
    records{q}.ititl
    ts = timeseries(spectrum, lambdas);
    inc = (high-low) / (numBands-1);
    c = resample(ts, [low:inc:high], 'zoh');
    numSpectra = numSpectra + 1;
    spectra(numSpectra).data = [c.time c.data [1:numBands].'];
    spectra(numSpectra).name = records{q}.ititl;
end

return;

author: davidkun

repo: HyperSpectralToolbox