text:
function inside = triangle_contains_point_2d ( t, p )

%*****************************************************************************80
%
%% TRIANGLE_CONTAINS_POINT_2D finds if a point is inside a triangle in 2D.
%
%  Discussion:
%
%    The routine assumes that the vertices are given in counter-clockwise
%    order.
%
%    The routine determines if a point P is "to the right of" each of the lines
%    that bound the triangle.  It does this by computing the cross product
%    of vectors from a vertex to its next vertex, and to P.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the triangle vertices.
%    The vertices should be given in counter clockwise order.
%
%    Input, real P(2), the point to be checked.
%
%    Output, logical INSIDE, is TRUE if the point is inside
%    the triangle or on its boundary.
%
  dim_num = 2;

  inside = 0;

  if ( 0.0 < ( p(1)   - t(1,1) ) * ( t(2,2) - t(2,1) ) ...
           - ( t(1,2) - t(1,1) ) * ( p(2)   - t(2,1) ) );
    return
  end

  if ( 0.0 < ( p(1)   - t(1,2) ) * ( t(2,3) - t(2,2) ) ...
           - ( t(1,3) - t(1,2) ) * ( p(2)   - t(2,2) ) );
    return
  end

  if ( 0.0 < ( p(1)   - t(1,3) ) * ( t(2,1) - t(2,3) ) ...
           - ( t(1,1) - t(1,3) ) * ( p(2)   - t(2,3) ) );
    return
  end

  inside = 1;

  return
end


text:
function [yd] = cm2yd(cm)
% Convert length from centimeters to yards.
% Chad A. Greene 2012
yd = cm*0.01093613298338;

author: Sable

repo: mcbench-benchmarks


text:
function [x,fval,exitflag,info] = opti_cplex(H,f,A,rl,ru,lb,ub,xint,sos,qc,x0,opts)
%OPTI_CPLEX Solve a LP/MILP/QP/MIQP/QCQP/MIQCQP using CPLEX
%
%   min 0.5*x'*H*x + f'*x      subject to:     rl <= A*x <= ru
%    x                                         qrl <= x'Q'x + l'x <= qru
%                                              lb <= x <= ub
%                                              for i = 1..n: xi in Z
%                                              for j = 1..m: xj in {0,1}
%
%   x = opti_cplex([],f,A,rl,ru,lb,ub,xint) solves a LP/MILP where f is the
%   objective vector, A,rl,ru are the linear constraints, lb,ub are the
%   bounds and xint is a string of integer variables ('C', 'I', 'B').
%
%   x = opti_cplex(H,f,A,rl,ru,lb,ub,xint) solves a QP/MIQP where H is the
%   objective matrix, and the remainder of the arguments are as above.
%
%   x = opti_cplex(H,...,xint,sos) sos is a structure with fields type,
%   index, weight for SOS.
%
%   x = opti_cplex(H,...,sos,qc) qc is structure with fields Q, l, qrl and
%   qru for quadratic constraints.
%
%   x = opti_cplex(H,...,qc,x0) supples an initial guess of the solution.
%
%   x = opti_cplex(H,f,...,x0,opts) uses opts to pass CPLEX parameters to
%   the solver. For example to set the node limit:
%       opts.mip.limits.nodes.Cur = 1000;
%
%   [x,fval,exitflag,info] = opti_cplex(...) returns the objective value at
%   the solution, together with the solver exitflag, and an information
%   structure.

%   Copyright (C) 2013 Jonathan Currie (I2C2)

t = tic;

% Handle missing arguments
if nargin < 12, opts = []; end
if nargin < 11, x0 = []; end
if nargin < 10, qc = []; end
if nargin < 9, sos = []; end
if nargin < 8, xint = []; end
if nargin < 7, ub = []; end
if nargin < 6, lb = []; end
if nargin < 5, error('You must supply at least 5 arguments to opti_cplex'); end

% Set Defaults (not returned from Cplex if infeasible)
x = []; fval = [];

%Check cplex exists
if(isempty(which('Cplex.p')))
    error('Cannot find Cplex the class! Please ensure the MATLAB Cplex interface is installed on your computer!');
end

%Build CPLEX Model
cp = Cplex('OPTI');
%Fill in common properties
cp.Model.obj = f;
if(isempty(lb))
    cp.Model.lb = -Inf(size(f));
else
    cp.Model.lb = lb;
end
if(isempty(ub))
    cp.Model.ub = Inf(size(f));
else
    cp.Model.ub = ub;
end
if(isempty(A))
    cp.Model.A = zeros(0,length(f));
    cp.Model.lhs = zeros(0,0);
    cp.Model.rhs = zeros(0,0);
else
    cp.Model.A = A;
    if(isempty(rl))
        cp.Model.lhs = -Inf(size(A,1),1);
    else
        cp.Model.lhs = rl;
    end
    if(isempty(ru))
        cp.Model.rhs = Inf(size(A,1),1);
    else
        cp.Model.rhs = ru;
    end
end
%Check and fill in remainder of optional args
if(~isempty(x0)),   cp.Start.x = x0; end
if(~isempty(H)),    cp.Model.Q = H; end
if(~isempty(xint)), cp.Model.ctype = upper(xint); end
if(~isempty(qc))
    %Check we have all fields
    if(~isstruct(qc) || ~isfield(qc,'Q') || ~isfield(qc,'l') || ~isfield(qc,'qrl') || ~isfield(qc,'qru'))
        error('Quadratic constraints must be a structure with fields ''Q'', ''l'', ''qrl'' and ''qru''');
    end
    if(iscell(qc.Q)) %multiple constraints
        if(length(qc.Q) ~= size(qc.l,2) || length(qc.qrl) ~= length(qc.Q) || length(qc.qrl) ~= length(qc.qru))
            error('Error with quadratic constraint dimensions');
        end
        for i = 1:length(qc.Q)
            addQC(cp,qc.Q{i},qc.l(:,i),qc.qrl(i),qc.qru(i));
        end
    else %single row constraint
        addQC(cp,qc.Q,qc.l,qc.qrl,qc.qru);
    end
end
if(~isempty(sos))
    %Check we have all fields
    if(~isstruct(sos) || ~isfield(sos,'type') || ~isfield(sos,'index') || ~isfield(sos,'weight'))
        error('SOS constraints must be a structure with fields ''type'', ''index'' and ''weight''');
    end
    %Ensure rows for Cplex
    if(iscell(sos.index) && size(sos.index,1) > size(sos.index,2))
        sos.index = sos.index';
    end
    if(iscell(sos.weight) && size(sos.weight,1) > size(sos.weight,2))
        sos.weight = sos.weight';
    end
    cp.addSOSs(sos.type,sos.index,sos.weight);
end

%Get objective constant term
if(~isempty(opts) && isfield(opts,'objbias') && ~isempty(opts.objbias))
    objbias = opts.objbias;
    opts = rmfield(opts,'objbias');
else
    objbias = 0;
end

%Assign User supplied options to Cplex
if(~isempty(opts))
    on = fieldnames(opts);
    if(length(on) == 1 && strcmpi('Param',on{1})) %remove Param if present
        opts = opts.Param;
    end
    %Bit of a hack to be able to set Cplex parameters, must be an easier way!
    fn = nestedfieldnames(opts);
    %Get cell arrays of each field
    cn = regexp(fn,'\.','split');
    %For each field, assign to param struct
    param = cp.Param;
    for i = 1:length(cn)
        try
            param = setfield(param,cn{i}{:},getfield(opts,cn{i}{:}));
        catch ME
            error('Error setting a Cplex parameter [%s]. Ensure you are setting the .Cur field of a valid parameter!\n\nError: %s',fn{i},ME.message);
        end
    end
    %Reassign parameters
    cp.Param = param;
end

%OPTI Default Parameters
cp.Param.output.clonelog.Cur = 0;

%Solve problem
sol = cp.solve();

%Assign Outputs
if(isfield(sol,'x')), x = sol.x; end
if(isfield(sol,'objval')), fval = sol.objval; end
if(isfield(sol,'itcnt')), info.Iterations = sol.itcnt; end
if(isfield(sol,'mipitcnt')), info.Nodes = sol.mipitcnt; end
if(isfield(sol,'miprelgap')), info.RelGap = sol.miprelgap; end
info.Time = toc(t); ex = '';
if(isfield(sol,'method'))
    switch(sol.method)
        case 0, ex = 'Automatic';
        case 1, ex = 'Primal Simplex';
        case 2, ex = 'Dual Simplex';
        case 3, ex = 'Network Simplex';
        case 4, ex = 'Barrier';
        case 5, ex = 'Sifting';
        case 6, ex = 'Concurrent Dual and Barrier';
        case 11, ex = 'Feasible Optimal';
        case 12, ex = 'MIP Solver';
        case 13, ex = 'Robust Algorithm';
    end
end
info.Algorithm = ['Cplex: ' ex];
if(isfield(sol,'statusstring'))
    info.Status = sol.statusstring;
    info.Status(1) = upper(info.Status(1));
else
    info.Status = 'Unknown';
end
if(isfield(sol,'status'))
    switch(sol.status)
        case {1,101,102,24}
            exitflag = 1; %ok
        case {25,10,11,12,105,107,111,112}
            exitflag = 0; %time/iterations/nodes
        case {3,5,103,115,106}
            exitflag = -1; %infeasible
        case {4,118,119}
            exitflag = -2; %unbounded
        case {13}
            exitflag = -5;
        otherwise
            exitflag = -6;
    end
else
    exitflag = -6;
end

%Add constant objective term
if(~isempty(fval))
    fval = fval + objbias;
end

%Assign Lambda
if(isfield(cp.Solution,'dual') && ~isempty(cp.Solution.dual))
    info.Lambda.lin = cp.Solution.dual; %not splitting up into ineq and eq as we are using row constraints
    rc = cp.Solution.reducedcost;
    zu = rc; zu(zu > 0) = 0;
    zl = rc; zl(zl < 0) = 0;
    info.Lambda.upper = zu;
    info.Lambda.lower = zl;
end

function addQC(cp,Q,l,rl,ru)
if(isinf(rl))
    cp.addQCs(l,Q,'L',ru);
elseif(isinf(ru))
    cp.addQCs(l,Q,'G',rl);
elseif(rl == ru)
    error('Cplex does not support quadratic equalities!');
else
    cp.addQCs(l,Q,'L',ru);
    cp.addQCs(l,Q,'G',rl);
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [soln,eqn,info] = DivCurl3PDWG(node,elem,bdFlag,pde,option,varargin)
%% DivCurl3PDWG div-curl system equation:
%        primal-dual weak-Galerkin in 3-D to solve the following system.
%       div(epsilon*u)=f  in \Omega,
%       curl(u) = g       in \Omega,
%       Dirichlet boundary condition epsilon*u \cdot n =g_n on \Gamma_0,
%
% The code is vectorized with no loops. Data structure follows iFEM tradition.
%
% Reference: A New Numerical Method for Div-Curl Systems with Low Regularity Assumptions
% S. Cao, C. Wang, J. Wang, https://arxiv.org/abs/2101.03466
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

if ~exist('option','var'), option = []; end
if ~isfield(option,'inverseh'), option.inverseh = true; end
if ~isfield(option,'dquadorder'), option.dquadorder = 3; end

[lambda,w] = quadpts3(option.dquadorder); % quadrature order is 1 or 2
nQuad = size(lambda,1);

if ~exist('bdFlag','var')
    bdFlag = setboundary3(node,elem,'Dirichlet');
end

%% Mesh related geometric quantities

% local edge to face indices, first two edges' local indices on a face
loc_e2f = [4 5; 2 3; 1 3; 1 2];
% local face to vertices indices, positive oriented
loc_f2v = [2 3 4; 1 4 3; 1 2 4; 1 3 2];

loc_e2v = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];

[elem2face,face] = dof3face(elem);
NT = size(elem,1);
NF = size(face,1);
elem2intdof = reshape(1:3*NT, 3, NT)';

[Dphi,volume] = gradbasis3(node,elem);

hK = (6*volume).^(1/3);

elem2ve = zeros(NT,3,6);
for e = 1:6 % six edges
   elem2ve(:,:,e) = node(elem(:,loc_e2v(e,2)),:)-node(elem(:,loc_e2v(e,1)),:);
end

elem2ve = elem2ve./repmat(sqrt(sum(elem2ve.^2,2)),1,3); % normalize for basis

faceNormal = cross(node(face(:,2),:) - node(face(:,1),:),...
    node(face(:,3),:) - node(face(:,2),:),2);
faceArea = 0.5*sqrt(sum(faceNormal.^2,2));
faceArea2elem = faceArea(elem2face);
faceArea2elemTotal = sum(faceArea2elem,2);

center = (node(elem(:,1),:) + node(elem(:,2),:) + ...
    node(elem(:,3),:) + node(elem(:,4),:))/4;

%% permeability tensor
t = cputime;  % record assembling time

if ~isfield(pde,'Eps'), pde.Eps = []; end

if ~isempty(pde.Eps) && isnumeric(pde.Eps)
    Eps2elem = repmat(pde.Eps, [1,1,NT]);   % Eps is a constant matrix
    Eps2elem = permute(Eps2elem,[3,1,2]); % switch the element idx to the 1st dim
end

if ~isempty(pde.Eps) && ~isnumeric(pde.Eps)
    Eps2elem = zeros(NT,3,3);  %#ok<PREALL>
    % Eps2elem is a NTx3x3 array so that Eps2elem(i,:,:) is the tensor at the
    % center of i-th element
    Epstmp = arrayfun(@(rowidx) pde.Eps(center(rowidx,:)), ...
        1:size(center,1), 'UniformOutput',0);
    Eps2elem = cat(3,Epstmp{:});
    Eps2elem = permute(Eps2elem,[3,1,2]); % switch the element idx to the 1st dim
end

%% compute eps \nabla_w phi
EpsDphi = zeros(NT,3,4);
for i = 1:4
    EpsDphi(:,:,i) = sum(bsxfun(@times, Eps2elem, Dphi(:,:,i)), 2);
end

%% assembling
% DoF order: # of Dofs
% \lambda_h: NT + NF, local 1+4
% u_h: 3*NT
% s_h: s_0 + s_b;  NT + NF (simply-connected case), local 1+4
% q_h: q_0 + q_b; 3*NT + 2*NF (for H(curl) vector field)
% assembling indices following locally within each group

%% equation 1 for u_h and lambda_h

blockUlambda = sparse(4*NT+NF,4*NT+NF);  %#ok<NASGU>
% (u_h, \epsilon\nabla_w phi)
% + \sum h_T^{-1} <\lambda_0 - \lambda_b, \varphi_0 - \varphi_b>_{boundary T}

%% (u_h, \epsilon\nabla_w phi_{face})
% locally the basis for P_0(T)^3 is edge [1, 2], [1, 3], and [1, 4]
% size (NT+NF, 3*NT)
C = sparse(NF,3*NT); %(u_h, \epsilon\nabla_w phi_{face})

for j = 1:3 % constant vector space spanned by 3 edge vectors
    for i = 1:4 % face 1 to 4
        % Cij = (q_j, eps \nabla_w \phi_{F_i})
        % \nabla_w \phi_{F_i} = -3\nabla \phi_i

        Cij = dot(elem2ve(:,:,j),-3*EpsDphi(:,:,i),2).*volume;

        ii = double(elem2face(:,i));
        jj = elem2intdof(:,j);
        C = C + sparse(ii,jj,Cij,NF,3*NT);
    end
end

Ugradphi = [sparse(NT, 3*NT); C];

%% \sum h_T^{-1} <\lambda_0 - \lambda_b, \varphi_0 - \varphi_b>_{boundary T}
% this is essetially a mass matrix stab
Stab_Mh = sparse(NT+NF,NT+NF); %#ok<NASGU>

L0ij = faceArea2elemTotal./hK;
Lambda00 = sparse(1:NT, 1:NT, L0ij, NT, NT); %% <\lambda_0, \lambda_0>

Lambda0b = sparse(NT, NF); %% <\lambda_b, \lambda_0>
Lambdabb = sparse(NF, NF); %% <\lambda_b, \lambda_b>
for j = 1:4 % face 1 to 4
    jj = double(elem2face(:,j));
    Lbj = faceArea2elem(:,j)./hK; % = <\lambda_{F_j}, \lambda_0>
    Lambda0b = Lambda0b + sparse(1:NT, jj, Lbj, NT, NF);
    Lambdabb = Lambdabb + sparse(jj, jj, Lbj, NF, NF);
end

Stab_Mh = [Lambda00 -Lambda0b; -Lambda0b' Lambdabb];
% Stablization on M_h should have null space dimension 1 (global constant)

%% \sum h_T <s_0 - s_b, r_0 - r_b>_{boundary T}

Stab_Sh = sparse(NT+NF,NT+NF); %#ok<NASGU>

S0ij = faceArea2elemTotal.*hK;
S00 = sparse(1:NT, 1:NT, S0ij, NT, NT); %% <s_0, s_0>

Sb0 = sparse(NT, NF); %% <s_b, s_0>
Sbb = sparse(NF, NF); %% <s_b, s_b>
for j = 1:4 % face 1 to 4
    jj = double(elem2face(:,j));
    Sbj = faceArea2elem(:,j).*hK; % = <\lambda_{F_j}, \lambda_0>
    Sb0 = Sb0 + sparse(1:NT, jj, Sbj, NT, NF);
    Sbb = Sbb + sparse(jj, jj, Sbj, NF, NF);
end

Stab_Sh = [S00 -Sb0; -Sb0' Sbb];
% Stablization on s_h should have null space dimension 1 (global constant)

%% Equation 1 for u_h and q_h = {q_0, q_b}

%%  (u_h, \epsilon\nabla_w \times psi_{face})
% locally the basis for P_0(T)^3 is the edge vector representing
% the edge [1, 2], [1, 3], and [1, 4] (local indexing of iFEM tradition)
% locally the basis for the tangential component are
% the edge 1 and 2 vectors for the triangular face (local indexing on this face)
% size: (3*NT+2*NF, 3*NT)

D = sparse(2*NF,3*NT);
%(u_h, \nabla_w \times psi_{face}) with 2 basis vec at each face

for i = 1:4 % face 1 to 4
    for j = 1:3 % constant vector space spanned by 3 edge vectors
        for e = 1:2 % two basis vectors on each face
            % Dije = (q_j, eps \nabla_w \times \psi_{F_i,e})
            % \nabla_w \times [0, \psi_{F_i}] =
            % (rotation of \psi_{F_i} counter-clockwise by pi/2)*|F_i|/|K|
            % = 3 \psi_{F_i} \times \nabla phi_i

            curlpsiFe = 3*mycross(elem2ve(:,:,loc_e2f(i,e)), Dphi(:,:,i));
            Dije = dot(elem2ve(:,:,j),curlpsiFe,2).*volume;

            % k-th (global face indexing) face has 2 DoFs in W_h: 2k-1, and 2k
            ii = 2*double(elem2face(:,i)) - e + 1;
            jj = elem2intdof(:,j);
            D = D + sparse(ii,jj, Dije, 2*NF,3*NT);
        end
    end
end

Ucurlpsi = [sparse(3*NT, 3*NT); D];

%%  (\psi, \epsilon\nabla_w s_h)
Psigrads = sparse(3*NT+2*NF, NT+NF);

E = sparse(3*NT,NF);
for i = 1:3 % constant vector space spanned by 3 edge vectors W_{h,0}
    for j = 1:4 % face 1 to 4
        % Eij = (\psi_{K,i}, eps \nabla_w \phi_{F_j})
        % \nabla_w \phi_{F_j} = -3\nabla \phi_j

        Eij = dot(elem2ve(:,:,i),-3*EpsDphi(:,:,j),2).*volume;

        ii = elem2intdof(:,i);
        jj = double(elem2face(:,j));
        E = E + sparse(ii,jj,Eij,3*NT,NF);
    end
end

%%% E is C's transpose, the code above is for booking-keeping
E = C';
Psigrads(1:3*NT,NT+1:end) = E;

%% \sum h_T^{-1} <(q_0-q_b)\times n, (\psi_0-\psi_b)\times n>_{boundary T}

% \sum h_T^{-1} <q_0\times n, \psi_0 \times n>_{boundary T}
WS00 = sparse(3*NT, 3*NT);
for i = 1:3 % constant vector space spanned by 3 edge vectors W_{h,0}
    for j = 1:3 % constant vector space spanned by 3 edge vectors W_{h,0}
        S00ijface = 0;
        for f = 1:4
            psi0crossn = mycross(elem2ve(:,:,i), Dphi(:,:,f)); % test
            q0crossn = mycross(elem2ve(:,:,j), Dphi(:,:,f)); % trial

            S00ijface = S00ijface + 9*dot(q0crossn,psi0crossn,2)...
                .*volume.^2./faceArea2elem(:,f);
        end
        S00ijface = S00ijface./hK;
        ii = elem2intdof(:,i);
        jj = elem2intdof(:,j);
        WS00 = WS00 + sparse(ii,jj,S00ijface,3*NT, 3*NT);
    end
end

% \sum h_T^{-1} <q_b\times n, \psi_0 \times n>_{boundary T}
WSb0 = sparse(3*NT, 2*NF);
for i = 1:3 % constant vector space spanned by 3 edge vectors W_{h,0}
    for j = 1:4 % 4 face for W_{h,b} each basis has support only on F_j
        for e = 1:2 % each face has 2 DoFs
            psi0crossn = mycross(elem2ve(:,:,i), Dphi(:,:,j)); % test
            qbcrossn = mycross(elem2ve(:,:,loc_e2f(j,e)), Dphi(:,:,j)); % trial

            Sb0ijface = 9*dot(qbcrossn,psi0crossn,2)...
                .*volume.^2./faceArea2elem(:,j)./hK;

            ii = elem2intdof(:,i);
            % k-th (global face indexing) face has 2 DoFs
            jj = 2*double(elem2face(:,j)) - e + 1;
            WSb0 = WSb0 + sparse(ii,jj,Sb0ijface,3*NT,2*NF);
        end
    end
end

% \sum h_T^{-1} <q_b\times n, \psi_b \times n>_{boundary T}
WSbb = sparse(2*NF, 2*NF);
for f = 1:4
    % this term is non-zero only when both trial and test
    % are associated with the same face
    for i = 1:2 % each face has 2 DoFs
        for j = 1:2
            psibcrossn = mycross(elem2ve(:,:,loc_e2f(f,i)), Dphi(:,:,f));
            qbcrossn = mycross(elem2ve(:,:,loc_e2f(f,j)), Dphi(:,:,f));

            Sbbijface = 9*dot(qbcrossn,psibcrossn,2)...
                .*volume.^2./faceArea2elem(:,f)./hK;

            % k-th (global face indexing) face has 2 DoFs
            ii = 2*double(elem2face(:,f)) - i + 1;
            jj = 2*double(elem2face(:,f)) - j + 1;
            WSbb = WSbb + sparse(ii,jj,Sbbijface,2*NF,2*NF);
        end
    end
end

Stab_Wh = [WS00 -WSb0; -WSb0' WSbb];
% Stab_Wh should have null space dimension
% = span{gradient of continuous P_1 Langrange} = # node -1

%% overall stiffness matrix for the system
blockUlambda = [Stab_Mh Ugradphi; Ugradphi' sparse(3*NT,3*NT)];
blockOffDiag = sparse(4*NT+3*NF, 4*NT+NF);
blockOffDiag(NT+NF+1:end, NT+NF+1:end) = Ucurlpsi;
blockQs = [-Stab_Sh Psigrads'; Psigrads Stab_Wh];
bigA = [blockUlambda blockOffDiag'; blockOffDiag blockQs];

%% right hand sides

%% -(f,\varphi_0)
ft = zeros(NT,1);

for p = 1:nQuad
		% quadrature points in the x-y-z coordinate
		pxyz = lambda(p,1)*node(elem(:,1),:) ...
			 + lambda(p,2)*node(elem(:,2),:) ...
			 + lambda(p,3)*node(elem(:,3),:) ...
             + lambda(p,4)*node(elem(:,4),:);
		fp = pde.f(pxyz);
        ft = ft + w(p)*fp;
end

Fphi0 = -accumarray((1:NT)',ft.*volume,[NT 1]);

%% \sum_{T} <eps u\cdot n, \varphi_b>_{\partial T} near \partial \Omega
GNphib = zeros(NF,1); % <eps u\cdot n, \varphi_F>_{\partial T}
% update: face quadorder increased

if ~isfield(option,'gNquadorder')
    option.gNquadorder = 5;   % default order exact for linear gN
end
[lambdagN, weightgN] = quadpts(option.gNquadorder);               % linear bases
nQuadgN = size(lambdagN,1);

for i = 1:4 % 4 faces
    idxBdElem = find(bdFlag(:,i) > 0); % if an element has a boundary face

    if ~isempty(idxBdElem)
        bdNormal = -3*repmat(volume(idxBdElem,:),1,3).*Dphi(idxBdElem,:,:);
        % scaled normal with face area built-in
        loc_face = loc_f2v(i,:);

        EpsUdotN = zeros(size(idxBdElem,1),1);

        for pp = 1:nQuadgN
            % quadrature points in the x-y coordinate
            ppxyz = lambdagN(pp,1)*node(elem(idxBdElem,loc_face(1)),:) ...
                + lambdagN(pp,2)*node(elem(idxBdElem,loc_face(2)),:) ...
                + lambdagN(pp,3)*node(elem(idxBdElem,loc_face(3)),:);

            u_bd = pde.exactu(ppxyz);
            EpsU = sum(bsxfun(@times, Eps2elem(idxBdElem,:,:), u_bd), 2);
            EpsUdotN = EpsUdotN + ...
                weightgN(pp)*dot(squeeze(EpsU), bdNormal(:,:,i), 2);

        end

        GNphib = GNphib + accumarray(elem2face(idxBdElem,i), EpsUdotN, [NF 1]);
    end
end

%% (g, \psi_0)
Gpsi0 = zeros(3*NT,1);

for i = 1:3 % 3 vector basis on each element for \psi_0
    gtpsi0 = zeros(NT,1);
    for p = 1:nQuad
        % quadrature points in the x-y-z coordinate
        pxyz = lambda(p,1)*node(elem(:,1),:) ...
            + lambda(p,2)*node(elem(:,2),:) ...
            + lambda(p,3)*node(elem(:,3),:) ...
            + lambda(p,4)*node(elem(:,4),:);
        gp = pde.g(pxyz);

        gtpsi0 = gtpsi0 + w(p)*dot(gp, elem2ve(:,:,i), 2);
    end

    Gpsi0 = Gpsi0 + accumarray(elem2intdof(:,i), gtpsi0.*volume, [3*NT 1]);
end
%% RHS
Rhs = [Fphi0; GNphib]; % test: phi_0, phi_b
Rhs = [Rhs; zeros(3*NT, 1)]; % test: v
Rhs = [Rhs; zeros(NT, 1)]; % test: r_0
Rhs = [Rhs; zeros(NF, 1)]; % test: r_b
Rhs = [Rhs; Gpsi0]; % test: psi_0
Rhs = [Rhs; zeros(2*NF, 1)]; % test: psi_b

%% set up free dof
% DoF order: # of Dofs
% prefix with "idx": global indices
% the DoF variable itself is boolean (fastest on MATLAB)
% \lambda_h: NT + NF, local 1+4 (no BC)
% u_h: 3*NT (no BC)
% s_h: s_0 + s_b;  NT + NF (simply-connected case), s_b = 0 on Gamma_0
% q_h: q_0 + q_b; 3*NT + 2*NF (for H(curl) vector field) q_b\cross n =0
% assembling indices following locally within each group

isBdFace = false(NF,1);
isBdFace(elem2face(bdFlag(:) > 0)) = true;
idxBdFace = find(isBdFace);
idxBdFaceVec = [2*idxBdFace-1 2*idxBdFace]';
idxBdFaceVec = idxBdFaceVec(:);

freeDof = true(8*NT+4*NF,1);
idxBdSb = idxBdFace + 5*NT+NF; % boundary DoF indices for s_b
idxBdQb = idxBdFaceVec + 8*NT+2*NF; % boundary DoF indices for q_b

%%%%% fixing 1 in lambda_0 is enough, the others are optional %%%%%
%%% current config according to the paper: fixing 1 in lambda_0
%%% bd dofs of s and q
freeDof([1; idxBdSb; idxBdQb]) = false;

%%%%% fixing bd dofs for q
% freeDof([1; idxBdQb]) = false;

%%%%% fixing bd dofs for s
% freeDof([1; idxBdSb]) = false;

idxIntBdFace = [];
% old way of handling cavity (second betti number not zero)
if any(elem2face(bdFlag(:) == 2)) && false
    disp('Modifying cavity dofs')
    isIntBdFace = false(NF,1);
    isIntBdFace(elem2face(bdFlag(:) == 2)) = true;
    idxIntBdFace = find(isIntBdFace);
    idxIntBdSb = idxIntBdFace + 5*NT+NF;  % boundary DoF indices for s_b on cavity
    freeDof(idxIntBdSb) = true;
end

%% Record assembling time
assembleTime = cputime - t;
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 2
    fprintf('Time to assemble matrix equation %4.2g s\n',assembleTime);
end

%% direct solve
NDof = 8*NT+4*NF;

bigSoln = zeros(NDof,1);

[~,sys] = memory;
if sum(freeDof) > 2e6 && sys.PhysicalMemory.Available < 3e10
    % 2 million DoFs need about 27GB free memory to be solved
    % 3 million DoFs need about 39GB free memory
    fprintf('Number of DoF %d exceeds estimated max dofs allowed in memory, stop. \n',sum(freeDof));
    return
end
bigSoln(freeDof) = bigA(freeDof, freeDof)\Rhs(freeDof);

%% post-processing for cavity
% TO-DO: current code allows only 1 cavity
if any(elem2face(bdFlag(:) == 2))
    disp('Post-processing cavity dofs')
    intBdFaceVec = false(size(freeDof,1), 1);
    isIntBdFace = false(NF,1);
    isIntBdFace(elem2face(bdFlag(:) == 2)) = true;
    idxIntBdFace = find(isIntBdFace);
    idxIntBdSb = idxIntBdFace + 5*NT+NF;  % boundary DoF indices for s_b on cavity
    intBdFaceVec(idxIntBdSb) = true;

    AintBdFaceVec = bigA*intBdFaceVec;
    sbConst = dot((Rhs - bigA*bigSoln),AintBdFaceVec)/dot(AintBdFaceVec,AintBdFaceVec);
    fprintf('s_b on the interior boundary is %6.4g \n',sbConst);
    bigSoln = bigSoln + sbConst*intBdFaceVec;
end

%% computing error
errL2U = zeros(NT,1);

idxUDoF = (NT+NF+1:NT+NF+3*NT)';
uh = bigSoln(idxUDoF);
uhElem = zeros(NT,3);

idxq0DoF = (5*NT+2*NF+1:5*NT+2*NF+3*NT)';
qh0 = bigSoln(idxq0DoF);
qh0Elem = zeros(NT,3);

idxqbDoF = (8*NT+2*NF+1:8*NT+4*NF)';
qhb = bigSoln(idxqbDoF);

idxqDoF = (5*NT+2*NF+1:8*NT+4*NF)';
qh = bigSoln(idxqDoF);

idxlambdaDoF = (1:NT+NF)';
lambdah = bigSoln(idxlambdaDoF);

lambdah0 = bigSoln((1:NT)');
lambdahb = bigSoln((NT+1:NT+NF)');

idxsDoF = (4*NT+NF+1:5*NT+2*NF)';
sh = bigSoln(idxsDoF);
sh0 = bigSoln((4*NT+NF+1:5*NT+NF)');
shb = bigSoln((5*NT+NF+1:5*NT+2*NF)');

for k = 1:3 % three bases in each K
    uhElem = uhElem + repmat(uh(elem2intdof(:,k)),1,3).*elem2ve(:,:,k);
    qh0Elem = qh0Elem + repmat(qh0(elem2intdof(:,k)),1,3).*elem2ve(:,:,k);
end

[lambda,w] = quadpts3(option.dquadorder+2); % quadrature order is 1 or 2
nQuad = size(lambda,1);
for p = 1:nQuad
    pxyz = lambda(p,1)*node(elem(:,1),:) ...
        + lambda(p,2)*node(elem(:,2),:) ...
        + lambda(p,3)*node(elem(:,3),:) ...
        + lambda(p,4)*node(elem(:,4),:);

    up = pde.exactu(pxyz);
    EpsUminusUh = sum(bsxfun(@times, Eps2elem, up - uhElem), 2);
    errL2U = errL2U + w(p)*dot(squeeze(EpsUminusUh), up - uhElem, 2);
end

errL2U = sqrt(errL2U.*volume);
errL2UTotal = norm(errL2U);

% idxfreeFace = find(~isBdFace);
% idxFreeFaceVec = [2*idxfreeFace-1 2*idxfreeFace]';
% idxFreeFaceVec = idxFreeFaceVec(:);
% errorStabTotal = sqrt(qh(idxFreeFaceVec)'*...
%     Stab_Wh(idxFreeFaceVec,idxFreeFaceVec)*qh(idxFreeFaceVec)...
%     +lambdah'*Stab_Mh*lambdah);
errorStabTotal = sqrt(qh'*Stab_Wh*qh+lambdah'*Stab_Mh*lambdah);
errorsTotal = sqrt(sh'*Stab_Sh*sh);

%% Output information
if nargout == 1
    soln = bigSoln;
else
    soln = struct('u',uh,...
                  'uh2elem', uhElem,...
                  'q0',qh0,...
                  'qb',qhb,...
                  'lambda0',lambdah0,...
                  'lambdab',lambdahb,...
                  's0',sh0,...
                  'sb',shb);
    eqn = struct('A',bigA,'b',Rhs,...
                'face',face,...
                'freeDof',freeDof,...
                'bdDofScalar',idxBdFace,...
                'bdDofVec',idxBdFaceVec,...
                'bdDofCavity',idxIntBdFace,...
                'errorU',errL2UTotal, ...
                'errorqlambda', errorStabTotal,...
                'errorS', errorsTotal,...
                'errorUelem',errL2U,...
                'NFace',NF);
    info.assembleTime = assembleTime;
end

%% end of function
end

repo: ifem


text:
function [pred_boxes, scores] = fast_rcnn_im_detect(conf, caffe_net, im, boxes, max_rois_num_in_gpu)
% [pred_boxes, scores] = fast_rcnn_im_detect(conf, caffe_net, im, boxes, max_rois_num_in_gpu)
% --------------------------------------------------------
% Fast R-CNN
% Reimplementation based on Python Fast R-CNN (https://github.com/rbgirshick/fast-rcnn)
% Copyright (c) 2015, Shaoqing Ren
% Licensed under The MIT License [see LICENSE for details]
% --------------------------------------------------------

    [im_blob, rois_blob, ~] = get_blobs(conf, im, boxes);

    % When mapping from image ROIs to feature map ROIs, there's some aliasing
    % (some distinct image ROIs get mapped to the same feature ROI).
    % Here, we identify duplicate feature ROIs, so we only compute features
    % on the unique subset.
    [~, index, inv_index] = unique(rois_blob, 'rows');
    rois_blob = rois_blob(index, :);
    boxes = boxes(index, :);

    % permute data into caffe c++ memory, thus [num, channels, height, width]
    im_blob = im_blob(:, :, [3, 2, 1], :); % from rgb to brg
    im_blob = permute(im_blob, [2, 1, 3, 4]);
    im_blob = single(im_blob);
    rois_blob = rois_blob - 1; % to c's index (start from 0)
    rois_blob = permute(rois_blob, [3, 4, 2, 1]);
    rois_blob = single(rois_blob);

    total_rois = size(rois_blob, 4);
    total_scores = cell(ceil(total_rois / max_rois_num_in_gpu), 1);
    total_box_deltas = cell(ceil(total_rois / max_rois_num_in_gpu), 1);
    for i = 1:ceil(total_rois / max_rois_num_in_gpu)

        sub_ind_start = 1 + (i-1) * max_rois_num_in_gpu;
        sub_ind_end = min(total_rois, i * max_rois_num_in_gpu);
        sub_rois_blob = rois_blob(:, :, :, sub_ind_start:sub_ind_end);

        net_inputs = {im_blob, sub_rois_blob};

        % Reshape net's input blobs
        caffe_net.reshape_as_input(net_inputs);
        output_blobs = caffe_net.forward(net_inputs);

        if conf.test_binary
            % simulate binary logistic regression
            scores = caffe_net.blobs('cls_score').get_data();
            scores = squeeze(scores)';
            % Return scores as fg - bg
            scores = bsxfun(@minus, scores, scores(:, 1));
        else
            % use softmax estimated probabilities
            scores = output_blobs{2};
            scores = squeeze(scores)';
        end

        % Apply bounding-box regression deltas
        box_deltas = output_blobs{1};
        box_deltas = squeeze(box_deltas)';

        total_scores{i} = scores;
        total_box_deltas{i} = box_deltas;
    end

    scores = cell2mat(total_scores);
    box_deltas = cell2mat(total_box_deltas);

    pred_boxes = fast_rcnn_bbox_transform_inv(boxes, box_deltas);
    pred_boxes = clip_boxes(pred_boxes, size(im, 2), size(im, 1));

    % Map scores and predictions back to the original set of boxes
    scores = scores(inv_index, :);
    pred_boxes = pred_boxes(inv_index, :);

    % remove scores and boxes for back-ground
    pred_boxes = pred_boxes(:, 5:end);
    scores = scores(:, 2:end);
end

function [data_blob, rois_blob, im_scale_factors] = get_blobs(conf, im, rois)
    [data_blob, im_scale_factors] = get_image_blob(conf, im);
    rois_blob = get_rois_blob(conf, rois, im_scale_factors);
end

function [blob, im_scales] = get_image_blob(conf, im)
    [ims, im_scales] = arrayfun(@(x) prep_im_for_blob(im, conf.image_means, x, conf.test_max_size), conf.test_scales, 'UniformOutput', false);
    im_scales = cell2mat(im_scales);
    blob = im_list_to_blob(ims);
end

function [rois_blob] = get_rois_blob(conf, im_rois, im_scale_factors)
    [feat_rois, levels] = map_im_rois_to_feat_rois(conf, im_rois, im_scale_factors);
    rois_blob = single([levels, feat_rois]);
end

function [feat_rois, levels] = map_im_rois_to_feat_rois(conf, im_rois, scales)
    im_rois = single(im_rois);

    if length(scales) > 1
        widths = im_rois(:, 3) - im_rois(:, 1) + 1;
        heights = im_rois(:, 4) - im_rois(:, 2) + 1;

        areas = widths .* heights;
        scaled_areas = bsxfun(@times, areas(:), scales(:)'.^2);
        [~, levels] = min(abs(scaled_areas - 224.^2), [], 2);
    else
        levels = ones(size(im_rois, 1), 1);
    end

    feat_rois = round(bsxfun(@times, im_rois-1, scales(levels))) + 1;
end

function boxes = clip_boxes(boxes, im_width, im_height)
    % x1 >= 1 & <= im_width
    boxes(:, 1:4:end) = max(min(boxes(:, 1:4:end), im_width), 1);
    % y1 >= 1 & <= im_height
    boxes(:, 2:4:end) = max(min(boxes(:, 2:4:end), im_height), 1);
    % x2 >= 1 & <= im_width
    boxes(:, 3:4:end) = max(min(boxes(:, 3:4:end), im_width), 1);
    % y2 >= 1 & <= im_height
    boxes(:, 4:4:end) = max(min(boxes(:, 4:4:end), im_height), 1);
end

author: ShaoqingRen


text:
function r83p_print ( n, a, title )

%*****************************************************************************80
%
%% R83P_PRINT prints a R83P matrix.
%
%  Discussion:
%
%    The R83P storage format stores a periodic tridiagonal matrix as
%    a 3 by N array, in which each row corresponds to a diagonal, and
%    column locations are preserved.  The matrix value
%    A(1,N) is stored as the array entry A(3,N), and the matrix value
%    A(N,1) is stored as the array entry A(1,1).
%
%  Example:
%
%    Here is how a R83P matrix of order 5 would be stored:
%
%      A51 A12 A23 A34 A45
%      A11 A22 A33 A44 A55
%      A21 A32 A43 A54 A15
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 April 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, real A(3,N), the R83P matrix.
%
%    Input, string TITLE, a title.
%
  r83p_print_some ( n, a, 1, 1, n, n, title );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef GMSHComputer < handle

    properties (Access = public)
        computation
        variables
    end

    properties (Access = private)
        testName
    end

    methods (Access = public)
        function obj = GMSHComputer(cParams)
            obj.testName = cParams.testName;
        end

        function compute(obj)
            filePath = fullfile('tests','Source','ReadingFilesTests','ReadingFiles','testReadingGmsh.msh');
            reader = GmsReader(filePath);
            reader.read();
            readData = reader.getDataBase();
            readData.isElemInThisSet = readData.isElemInThisSet{1};
            obj.variables = readData;
            obj.computation = obj;
        end
    end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [ R, subgrad, data ] = risk2( data, W )
%RISK evaluates risk term.
%
% Synopsis:
%  [R,subgrad,data] = risk(data)
%  [R,subgrad,data] = risk(data,W)
%
% Description:
%  Let the risk term be defined as
%
%  R(W) = 1/m sum_{i=1}^m [ max_{y \in Y}( (L(y_i, y) + <w, \Psi(x_i, y)> ) - <w, \Psi(x_i, y+i)> ]
%
%  This function returns value R and subgradient SUBGRAD of the
%  risk R(W) at W.
%
% 10-08-10 Michal Uricar
% 11-07-11 Michal Uricar, corners dataset (checked-only) + loss function modification

    options = data.options;

    if (nargin < 2)
        W = buildPsi2(options, data.tmpData, data.Y{1});
        W = sparse(length(W), 1);
    end

    psi_xiyhat = sparse(length(W), data.nImages);
    psi_xiyi = sparse(length(W), data.nImages);
    suma_all = zeros(1, data.nImages);

    for i = 1 : data.nImages
%     parfor i = 1 : data.nImages
        [lbpdat lbp_sparse] = getPsiMat(data, i);
        GT = data.Y{i};

        % \hat{y_i} = argmax_{y \in Y} [ L(y_i, y) + <w, \Psi(x_i, y)> ]
        L = computeL(options, GT, data.kappa(i));       % with normalization
%         L = computeL(options, GT);                      % without normalization
        y_hat = argmax_mex(options, W, data.mapTable, lbp_sparse, L);

        % \Psi(x_i, \hat{y_i})
        psi_xiyhat(:, i) = buildPsi2(options, lbpdat, y_hat);
        % \Psi(x_i, y_i)
        psi_xiyi(:, i) = buildPsi2(options, lbpdat, GT);

        % sum_{i = 1}^{m} ( max_{y \in Y}( L(y_i, y) + <w, \Psi(x_i, y_i)> )> ) - <w, \Psi(x_i, y_i)> )
        % with normalization coefficient kappa
        suma_all(i) = 100 * data.kappa(i) * ...
                        sum(1/data.options.M * sqrt(sum( (y_hat - GT).^2 ))) ...
                        + W'*psi_xiyhat(:, i) - W'*psi_xiyi(:, i);
%         % without normalization coefficient kappa
%         suma_all(i) = sum(1/data.options.M * sqrt(sum( (y_hat - GT).^2 ))) ...
%                         + W'*psi_xiyhat(:, i) - W'*psi_xiyi(:, i);
    end;

    suma = sum(suma_all);

    % \hat{R}(w) = 1/m \sum_{i = 1}^m [ max_{y \in Y}( L(y_i, y) + <w, \Psi(x_i, y_i)> )> ) - <w, \Psi(x_i, y_i)> ]
    R = 1/data.nImages * suma;
    % \partial_w \hat{R}(w) = 1/m \sum_{i = 1}^m ( \Psi(x_i, \hat{y_i}) - \Psi(x_i, y_i) )
    subgrad = 1/data.nImages * sum(psi_xiyhat - psi_xiyi, 2);

end

author: uricamic

repo: flandmark

save_path: github-repos/MATLAB/uricamic-flandmark


text:
function [elem,idx,volume,bdFlag] = fixorder3(node,elem,bdFlag)
%% 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 volume is negative.
%
%   [elem,idx,volume,bdFlag] = FIXORDER3(node,elem,bdFlag) changes the bdFlag
%   for boundary conditions.
%
% See also fixorder
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if (nargin==2)
    bdFlag = [];
end
% compute volume and elemSign of each tetrahedron
[volume,elemSign] = simplexvolume(node,elem);
% find tetrahedron with negative volume and switch the vertices
idx = find(elemSign==-1);
elem(idx,[2 3]) = elem(idx,[3 2]);
if exist('bdFlag','var') && ~isempty(bdFlag)
    bdFlag(idx,[2 3]) = bdFlag(idx,[3 2]);
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: Doyen Sahoo
% Contributors: Bin LI, Steven C.H. Hoi
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ ] = configurationMenu( )

    % This will be used to add and delete algorithms to the configuration

    menuId = 8; % Configuration Menu
    displayMenu(menuId);
    prompt = 'Please enter your choice (1-6):';
    choice = input(prompt);
    switch(choice),
        case 1,
            % Adding a Strategy
            cd ../GUI/config/;
            disp('Instructions:');
            disp('1. Create new strategy, and save it in folder /Strategy');
            disp('2. If file name is ''newStrategy.m'', create and appropritely populate ''newStrategy_config.m'' in the same folder');
            disp('3. An example of this type of files is shown in ''template.m'' and ''template_config.m''');
            newStrategyFileName = input('Enter name of file that you want to add (e.g. ''newStrategy''):');
            configFile          = input('Enter name of configuration you want to change (e.g. ''config_myconfig''):');
            addStrategy( newStrategyFileName, configFile );
            cd ../../PGUI/;
            configurationMenu;
        case 2,
            % Removing a Strategy
            prompt = ('Name of config file do you want to update (e.g. ''config_myconfig''):');
            configFile = input(prompt);
            cd ../GUI/config/;
            load(configFile);

            % Display all Strategies
            disp('Strategies in current config:');
            nAlgo = length(algorithmList);
            for i = 1:1:nAlgo
                msg = strcat(num2str(i), '.');
                msg = strcat(msg, algorithmList{i});
                fprintf(msg);
                fprintf('\n');
            end

            prompt = ('Enter ID of Strategy to be removed:');
            strategyId = input(prompt);
            removeStrategy( strategyId, configFile );

            cd ../../PGUI/;
            configurationMenu;
        case 3,
            % Adding a dataset
            cd ../GUI/config/;
            disp('First, save new dataset in folder /Data, in the correct format.');
            prompt              = 'Enter name of file (e.g. ''newData''):';
            dataFileName        = input(prompt);
            prompt              = 'Enter Data Description (e.g. '' New Data - 1st Jan, 2014 to 1st Jan, 2015):';
            dataDescriptionName = input(prompt);
            prompt              = 'Enter frequency (e.g. 1 for daily data):';
            newDataFrequency    = input(prompt);
            prompt              = 'Enter name of config file do you want to update (e.g. ''config_myconfig''):';
            configFile          = input(prompt);
            addData( dataFileName, dataDescriptionName, newDataFrequency, configFile );
            cd ../../PGUI/;
            configurationMenu;
        case 4,
            % Removing a dataset
            cd ../GUI/config/;
            prompt = ('Name of config file do you want to update (e.g. ''config_myconfig''):');
            configFile = input(prompt);
            load(configFile);

            % Display all data sets
            disp('Datasets in current config:');
            nData = length(dataList);
            for i = 1:1:nData
                msg = strcat(num2str(i), '.');
                msg = strcat(msg, dataList{i});
                fprintf(msg);
                fprintf('\n');
            end

            prompt = ('Enter ID of Data to be removed:');
            dataId = input(prompt);
            removeData( dataId, configFile );
            cd ../../PGUI/;
            configurationMenu;
        case 5,
            prompt = ('Name of config file do you want to make active (e.g. ''config_myconfig''):');
            configFile = input(prompt);
            cd ../GUI/config/;
            load(configFile);
            save('config', 'algorithmList', 'algorithmName', 'algorithmParameters', 'dataFrequency', 'dataList', 'dataName', 'defaultParameters', 'windowRisk');
            disp('Update Configuration. Restart toolbox to read the updated configuration');
            cd ../../PGUI/;
            configurationMenu;
        case 6,
            disp('Exiting Configuration Manager --> to Home');
            homeMenu;
        otherwise,
            disp('ERROR: Please enter a valid input');
            configurationMenu;
    end
end

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function hP = plugin_board()
% Zynq Platform PCore
% Use Plugin API to create board plugin object

%   Copyright 2015 The MathWorks, Inc.

% Call the common board definition function
hP = AnalogDevicesAGCDemo.adrv9361z7035.common.plugin_board('fmc lvds agc', 'Rx');

author: analogdevicesinc


text:
function [ varargout ] = readLeCroyWaveJetWfm( varargin )
%readLeCroyWaveJetWfm - Read binary waveform (.wfm) files as saved by
% LeCroy WaveJet oscilloscopes.
%
% readLeCroyWaveJetWfm( filename )
%
%   Loads the specified file, shows some information about the file,
%   and displays the waveform(s).
%
% [ data ] = readLeCroyWaveJetWfm( filename )
% [ data nfo ] = readLeCroyWaveJetWfm( filename )
%
%   Load the specified file and returns its data.
%   data - Array with loaded waveform(s). If only one channel was recorded
%          this is a single column with the waveform. If more than one
%          channel was recorded this is a 2D array where each column
%          corresponds with a channel. If a lower number channel is not
%          recorded it will contain a column with NaN's. If a higher number
%          channel is not recorded, it's column is removed from the data
%          array. E.g. ch1: x, ch2: 0, ch3: x, ch4: 0 will result in a 2D
%          array with 3 column, where the second column contains NaN's.
%   nfo - structure with time axis information and text header
%         nfo.x0    - time of first sample (trigger will be at t = 0)
%         nfo.dx    - sampling time, i.e. time between two consecutive samples
%         nfo.notes - text field with complete oscilloscope setup.
%

% Copyright (c) 2008-2010, Paul Wagenaars
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
%
% Changelog:
% v0.1 (2008-06-29)
%   - initial release
%
% v0.2 (2008-08-??)
%   - in averaging mode the binary WFM waveform seems to have a half a
%     sample offset compared to the ASCII CSV file. The binary waveform is
%     now shifted half sample downwards.
%   - if averaging mode is turned on the data of inactive channels was
%     skipped incorrectly. This has been fixed.
%
% v0.3 (2010-03-02)
%   - made function compatible with GNU Octave.
%   - fixed waveforms with multiple waveform modes, e.g. when averaging and
%     equivalent sampling are specified at the same time.
%

% Define some constants
global RLWJW_XWAV_DEF_INF_NOACQ RLWJW_XWAV_DEF_INF_NORM RLWJW_XWAV_DEF_INF_PEAK
global RLWJW_XWAV_DEF_INF_EQU   RLWJW_XWAV_DEF_INF_AVG  RLWJW_XWAV_DEF_INF_ROLL
global RLWJW_XWAV_DEF_INF_ILV
RLWJW_XWAV_DEF_INF_NOACQ = uint16(0);  % 0x0000 /* Not acquired */
RLWJW_XWAV_DEF_INF_NORM = uint16(1);   % 0x0001 /* Normal wave */
RLWJW_XWAV_DEF_INF_PEAK = uint16(2);   % 0x0002 /* Peak wave */
RLWJW_XWAV_DEF_INF_EQU = uint16(4);    % 0x0004 /* Equ wave */
RLWJW_XWAV_DEF_INF_AVG = uint16(8);    % 0x0008 /* Average wave */
RLWJW_XWAV_DEF_INF_ROLL = uint16(16);  % 0x0010 /* Roll wave */
RLWJW_XWAV_DEF_INF_ILV = uint16(32);   % 0x0020 /* Interleave wave */

% Maximum memory enum
global MAX_MEM_ENUM
MAX_MEM_ENUM = 1024 * [0 0.5 ; 1 1 ; 2 10 ; 3 100 ; 4 500];

% V/div enum
global RLWJW_VDIV_ENUM
RLWJW_VDIV_ENUM = [10 5 2 1 500e-3 200e-3 100e-3 50e-3 20e-3 10e-3 5e-3 2e-3];

% time/div enum
global RLWJW_TDIV_ENUM
RLWJW_TDIV_ENUM = [50     20     10     5     2     1     500e-3  200e-3  100e-3 ...
             50e-3  20e-3  10e-3  5e-3  2e-3  1e-3  500e-6  200e-6  100e-6 ...
             50e-6  20e-6  10e-6  5e-6  2e-6  1e-6  500e-9  200e-9  100e-9 ...
             50e-9  20e-9  10e-9  5e-9  2e-9  1e-9  500e-12 200e-12 100e-12 ...
             50e-12 20e-12 10e-12 5e-12 2e-12 1e-12 500e-15 200e-15 100e-15];

% Display informatie about waveform(s)
infomode = (nargout == 0);

if nargin < 1
    error('No filename specified.');
else
    filename = varargin{1};
end

if infomode
    fprintf('Loading %s\n', filename);
end

% Open file
[fid, msg] = fopen(filename, 'r', 'b');
if fid == -1
    error('Unable to open %s: %s', filename, msg);
end

% Skip text header
fullheader = '';
while true
    line = fgetl(fid);
    fullheader = sprintf('%s%s\n', fullheader, line);
    if ~ischar(line)
        fclose(fid);
        error('Unable to find end of setup information. Is this a correct LeCroy WaveJet .wfm file?');
    elseif strcmp(line, '[END OF SETUP]');
        % The 'fgetl' function seems to consume one byte too much.
        % So, here we check the byte and rewind in necessary.
        fseek(fid, -1, 'cof');
        data = fread(fid, 1, 'uint8');
        if (data ~= 10) && (data ~= 13)
            fseek(fid, -1, 'cof');
        end
        break;
    end
end

% Read the XWAV_TDE_WAVEINF structure
wavinf = load_XWAV_TDE_WAVEINF(fid);

% Read waveform data
data = nan(wavinf.timebase.memoryLength + 1, 4);
dataStart = wavinf.timebase.validWfmStart;
memLen = wavinf.timebase.memoryLength;
bankSize = wavinf.timebase.bankSize;
validChannels = nan(1,4);
if bitand(wavinf.timebase.wfmMode, RLWJW_XWAV_DEF_INF_AVG)
    bytesPerSample = 2;
else
    bytesPerSample = 1;
end
sampleType = sprintf('int%d', 8 * bytesPerSample);
maxValue = 2^(8 * bytesPerSample);
for j = 1:4
    if wavinf.channels(j).isValid
        validChannels(j) = 1;
        singleChannel = fread(fid, bankSize, sampleType);
        if bitand(wavinf.timebase.wfmMode, RLWJW_XWAV_DEF_INF_AVG)
            % In average mode there is a half a sample offset compared to
            % the csv file. Here we correct that.
            singleChannel = singleChannel - 127;
        end
        singleChannel = singleChannel * 8 * wavinf.channels(j).vdiv / maxValue;
        singleChannel = singleChannel - wavinf.channels(j).offset;
        data(:, j) = singleChannel((5 + dataStart):(memLen + dataStart + 5));
    else
        validChannels(j) = 0;
        fseek(fid, bankSize * bytesPerSample, 'cof');
    end
end
fclose(fid);
clear singleChannel;

% Select only the "good" channels
if sum(validChannels) == 1
    data = data(:, validChannels == 1);
elseif sum(validChannels) < 4
    maxIndex = find(validChannels == 1, 1, 'last');
    data = data(:, 1:maxIndex);
end

% Create nfo structure
totalTime = 10.0 * wavinf.timebase.tdiv;
nfo.dx = totalTime / (wavinf.timebase.memoryLength);
nfo.x0 = -0.5 * totalTime - wavinf.timebase.delay + nfo.dx;
nfo.notes = fullheader;

if infomode
    showInfo(wavinf, nfo, data);
end

% Set the output arguments
if nargout >= 1
    varargout{1} = data;
end
if nargout >= 2
    varargout{2} = nfo;
end

end  % End of main function

%%
%% SUBFUNCTIONS
%%

%% load_XWAV_TDE_WAVEINF
% Loads a wave information structure from file.
% The structure (C code):
% typedef struct
% {
%   XWAV_TDE_WAVEINF_TBASE  TbaseInfo;  /* Time Base Info. */
%   XWAV_TDE_WAVEINF_CH     ChInfo[4];  /* Ch Info. */
% } XWAV_TDE_WAVEINF;
function [waveinfo] = load_XWAV_TDE_WAVEINF(fid)
    waveinfo.timebase = load_XWAV_TDE_WAVEINF_TBASE(fid);
    waveinfo.channels = load_XWAV_TDE_WAVEINF_CH(fid);
end

%% load_XWAV_TDE_WAVEINF_TBASE
% Loads a time base information structure from file.
% The structure (C code):
% typedef struct
% {
%   SVMNG_TDE_TRGDLY delay;        /* Trigger delay time?i-500s..750s?j*/
%   XWAV_TDE_TSTAMP  Tstamp;       /* Time stamp */
%   IWDP_TDE_BANKS   bankNr;       /* Current bank size setting in enum */
%   Uint             mleng;        /* Memory length */
%   Uint             BankSize;     /* Actual bank size (512..512*1024) */
%   Uint             TrgSkew;      /* Trigger skew(TDC): WDP internal value */
%   Uint             WavValidAdr;  /* Start address for valid waveform data */
%   Ushort           wvinf;        /* Bitfielf of waveform mode */
%   Ushort           AvgTimes;     /* Number of average (0,1?`256) */
%   char             tim;          /* time/div in enum */
%   char             spclk;        /* sampling clock in enum */
%   Uchar            RdmInitVal;   /* Random seed: WDP internal value */
%   Uchar            WavSkew;      /* Fraction of delay and trigger skew in a
%                                     sampling clock. 0..99 for normalwaveform,
%                                     0..49 for interleaved waveform,
%                                     (0..100)/(equivalent sampling) for EQU. */
%   Uchar            maxmem;       /* MAX MEMORY LENGTH	0:0.5k, 1:1k, 2:10k,
%                                     3:100k, 4:500k */
%   Ushort           MaskBank;     /* Mask for bank number */
% } XWAV_TDE_WAVEINF_TBASE;
%
% Sizes of the differnet bank numbers (bankNr):
% IWDP_ENM_BANKS = [8 16 32 64 128 256 512 1024 2*1024 4*1024 8*1024 ...
%      16*1024 32*1024 64*1024 128*1024 256*1024 512*1024 1048576 ...
%      2*1048576 4*1048576 Inf];
function [tbase] = load_XWAV_TDE_WAVEINF_TBASE(fid)
    global MAX_MEM_ENUM RLWJW_TDIV_ENUM
    tbase.delay = load_SVMNG_TDE_TRGDLY(fid);
    tbase.timeStamp = load_XWAV_TDE_TSTAMP(fid);
    tbase.bankNr = fread(fid, 1, 'uint32');
    tbase.memoryLength = fread(fid, 1, 'uint32');
    tbase.bankSize = fread(fid, 1, 'uint32');
    tbase.triggerSkew = fread(fid, 1, 'uint32');
    tbase.validWfmStart = fread(fid, 1, 'uint32');
    tbase.wfmMode = fread(fid, 1, 'uint16');
    tbase.nrAvgs = fread(fid, 1, 'uint16');
    tbase.tdivEnum = fread(fid, 1, 'uint8');
    tbase.tdiv = RLWJW_TDIV_ENUM(1 + tbase.tdivEnum);
    tbase.samplingClockEnum = fread(fid, 1, 'uint8');
    tbase.rndInit = fread(fid, 1, 'uint8');
    tbase.wavSkew = fread(fid, 1, 'uint8');
    tbase.maxMemory = MAX_MEM_ENUM(1 + fread(fid, 1, 'uint8'), :);
    tbase.bankMask = fread(fid, 1, 'uint16');
    fread(fid, 1, 'uint8');
end

%% load_SVMNG_TDE_TRGDLY
% // Trigger delay structure
% typedef struct {
%  Uint sgn;  // -1:-/0:+
%  Uint sec;  // 0?`750s     (1s?`750s)
%  Uint us;   // 0?`999999us (1us?`999ms)
%  Uint ps;   // 0?`999999ps (1ps?`999ns)
% } SVMNG_TDE_TRGDLY;
function [delay] = load_SVMNG_TDE_TRGDLY(fid)
    sgn = fread(fid, 1, 'uint32');
    sec = fread(fid, 1, 'uint32');
    us = fread(fid, 1, 'uint32');
    ps = fread(fid, 1, 'uint32');
    if sgn == 0
        delay = 1.0;
    else
        delay = -1.0;
    end
    delay = delay * (sec + 1e-6 * us + 1e-12 * ps);
end

%% load_XWAV_TDE_TSTAMP
% Load time stamp information.
%
% // Date structure
% typedef struct {
%   Ushort  year;
%   Uchar   month;
%   Uchar   day;
% } SVMNG_TDE_DATE;
%
% // Time structure
% typedef struct {
%   Uchar   hour;
%   Uchar   min;
%   Uchar   sec;
% } SVMNG_TDE_TIME;
%
% typedef struct {
%   SVMNG_TDE_DATE  date;   // DATE year:2000..2099, month:1..12, day:1..31
%   SVMNG_TDE_TIME  time;   // TIME hour:0..23, min:0..59, sec:0..59
%   Uchar           msec;   // 100ms
% } XWAV_TDE_TSTAMP;
function [tstamp] = load_XWAV_TDE_TSTAMP(fid)
    year = fread(fid, 1, 'uint16');
    month = fread(fid, 1, 'uint8');
    day = fread(fid, 1, 'uint8');
    hour = fread(fid, 1, 'uint8');
    min = fread(fid, 1, 'uint8');
    sec = fread(fid, 1, 'uint8');
    msec = 100 * fread(fid, 1, 'uint8');
    tstamp = [year month day hour min sec msec];
end

%% load_XWAV_TDE_WAVEINF_CH
% Load channel information for a single channel.
function [chs] = load_XWAV_TDE_WAVEINF_CH(fid)
    global RLWJW_VDIV_ENUM
    chs = struct('offset', {NaN NaN NaN NaN}, ...
        'vdivEnum', {NaN NaN NaN NaN}, ...
        'isValid', {false false false false});
    for i = 1:4
        chs(i).offset = fread(fid, 1, 'int32');
        chs(i).vdivEnum = fread(fid, 1, 'uint8');
        chs(i).vdiv = RLWJW_VDIV_ENUM(1 + chs(i).vdivEnum);
        chs(i).offset = chs(i).offset / 3200 * chs(i).vdiv;
        chs(i).isValid = (fread(fid, 1, 'uint8') > 0);
        chs(i).unknown = fread(fid, 1, 'uint16');
    end
end

%% Show information about the waveform(s)
function showInfo(wavinf, nfo, data)
    global RLWJW_XWAV_DEF_INF_NOACQ RLWJW_XWAV_DEF_INF_NORM
    global RLWJW_XWAV_DEF_INF_PEAK  RLWJW_XWAV_DEF_INF_EQU
    global RLWJW_XWAV_DEF_INF_AVG   RLWJW_XWAV_DEF_INF_ROLL
    global RLWJW_XWAV_DEF_INF_ILV
    fprintf('Time stamp: %04d-%02d-%02d %02d:%02d:%02d.%03d\n', ...
        wavinf.timebase.timeStamp);
    fprintf('Memory length: %d (max.mem %dk)\n', wavinf.timebase.memoryLength, ...
        wavinf.timebase.maxMemory(2)/1024);
    fprintf('Bank size: %d\n', wavinf.timebase.bankSize);
    fprintf('Waveform mode:');
    if wavinf.timebase.wfmMode == RLWJW_XWAV_DEF_INF_NOACQ
        fprintf(' not_acquired');
    else
        if bitand(wavinf.timebase.wfmMode, RLWJW_XWAV_DEF_INF_NORM)
            fprintf(' normal');
        end
        if bitand(wavinf.timebase.wfmMode, RLWJW_XWAV_DEF_INF_PEAK)
            fprintf(' peak');
        end
        if bitand(wavinf.timebase.wfmMode, RLWJW_XWAV_DEF_INF_EQU)
            fprintf(' equivalent_sampling');
        end
        if bitand(wavinf.timebase.wfmMode, RLWJW_XWAV_DEF_INF_AVG)
            fprintf(' averaged(%dx)', wavinf.timebase.nrAvgs);
        end
        if bitand(wavinf.timebase.wfmMode, RLWJW_XWAV_DEF_INF_ROLL)
            fprintf(' roll');
        end
        if bitand(wavinf.timebase.wfmMode, RLWJW_XWAV_DEF_INF_ILV)
            fprintf(' interleaved');
        end
    end
    fprintf('\n');
    fprintf('Hor. resolution: %g s/div (%d)\n', wavinf.timebase.tdiv, ...
        wavinf.timebase.tdivEnum);
    fprintf('Trigger delay: %g s\n', wavinf.timebase.delay);
    fprintf('Valid wfm start: %d\n', wavinf.timebase.validWfmStart);

    fprintf('                    Channel 1  Channel 2  Channel 3  Channel 4 \n');
    fprintf('offset         (V)  %9g  %9g  %9g  %9g\n', wavinf.channels.offset);
    fprintf('resolution (V/div)  %9g  %9g  %9g  %9g\n', wavinf.channels.vdiv);
    fprintf('channel on or off?  %9g  %9g  %9g  %9g\n', wavinf.channels.isValid);
    fprintf('unknown             %9d  %9d  %9d  %9d\n', wavinf.channels.unknown);

    xf = (wavinf.timebase.memoryLength) * nfo.dx + nfo.x0;
    t = linspace(nfo.x0, xf, wavinf.timebase.memoryLength + 1);
    figure;
    plot(t, data);
    grid on;
    xlabel('Time (s)');
end  % end of showInfo function


text:
function a = r8vec_normalize_l1 ( n, a )

%*****************************************************************************80
%
%% R8VEC_NORMALIZE_L1 normalizes an R8VEC to have unit sum.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, real A(N), the vector to be normalized.
%
%    Output, real A(N), the entries of A should have unit sum.  However,
%    if the input vector has zero sum, the routine halts.
%
  a_sum = sum ( a(1:n) );

  if ( a_sum == 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8VEC_NORMALIZE_L1 - Fatal error!\n' );
    fprintf ( 1, '  The vector entries sum to 0.\n' );
    error ( 'R8VEC_NORMALIZE_L1 - Fatal error!' );
  end

  a(1:n) = a(1:n) / a_sum;

  return
end