text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1. use vectorization to generate face and face2elem data structure
%% this vectorization code does not work for uniform tets

% tic
%
% matlabTri = 1;
% N = size(mesh.p,1);
% node = [mesh.p;mesh.eIntP];
% allCutElem = mesh.t(mesh.tLoc<0,:);
% isCutElem = (mesh.tLoc<0);
% vSign = [mesh.pLoc;zeros(size(mesh.eIntP,1),1)];
%
% isInterfaceNode = false(size(node,1), 1);
% isInterfaceNode(allCutElem(:)) = true; % include the vertices of allCutElem
% isInterfaceNode(N+1:end) = true;    % and the cut points and the aux points
% interfaceNode = node(isInterfaceNode,:);
% matlabversion = version();
% if matlabTri == 1
%     DT = delaunayTriangulation(interfaceNode);
%     tetElem = DT.ConnectivityList;
% elseif matlabTri == 0
%     constrant1 = [mesh.e(mesh.eLoc<0,1),N-mesh.eLoc(mesh.eLoc<0)];
%     constrant2 = [N-mesh.eLoc(mesh.eLoc<0),mesh.e(mesh.eLoc<0,2)];
%     constrant = [constrant1;constrant2];
%     DT = constrainedDelaunayTetGen(interfaceNode,constrant);
%     tetElem = DT.ConnectivityList;
% end
% [tetElem, ortidx, volume] = fixorder3(interfaceNode, tetElem);
%
% localidx2globalidx = find(isInterfaceNode); % tetElem points to interfaceNode
% tetElem = localidx2globalidx(tetElem);  % map to the global index of node
% bc = (node(tetElem(:, 1), :) + node(tetElem(:, 2), :) ...
%     + node(tetElem(:, 3), :) + node(tetElem(:, 4), :))/4.0;
% idx2cube = FindElemId(bc, mesh); % Get the cube index of each tetElem
% idx2cubeBad = (idx2cube<=0);
% idx2cube = idx2cube(~idx2cubeBad);
% tetElem = tetElem(isCutElem(idx2cube), :); % Just keep the tets in cut elem
% volume = volume(isCutElem(idx2cube), :);
% idx2cube = idx2cube(isCutElem(idx2cube));  % keep cube idx in cut elements
% % volumetet0 =  accumarray(idx2cube,volume);
% % volumetet0 = volumetet0(volumetet0>0);
% %[tetElemOld, ortidxOld, volumeOld] = fixorder3(interfaceNode, mesh.t(mesh.tLoc<0,:));
%
% % There might exist some bad tet elements whose vertices are in
% % different cubes. We just need to keep the tet in every cut elem. So
% % here we get rid of them. (but do not handle different tets)
% NT = size(tetElem,1);
% X = reshape(node(tetElem,1), NT, 4);
% Y = reshape(node(tetElem,2), NT, 4);
% Z = reshape(node(tetElem,3), NT, 4);
% isBadTet =  (max(X,[],2) - min(X,[],2)) > 2*h - eps ... % d > 2h means
%           | (max(Y,[],2) - min(Y,[],2)) > 2*h - eps ... % in differnt cube
%           | (max(Z,[],2) - min(Z,[],2)) > 2*h - eps;
% tetElem = tetElem(~isBadTet,:);
%
% toc

% plot to check
%showmesh3(node,tetElem);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% this approach is based on forlopp which is slow
tic

Np = size(mesh.p,1);
node = [mesh.p;mesh.eIntP];
allCutElem = mesh.t(mesh.tLoc<0,:);
isCutElem = (mesh.tLoc<0);
vSign = [mesh.pLoc;zeros(size(mesh.eIntP,1),1)];
isInterfaceNode = false(size(node,1), 1);
isInterfaceNode(allCutElem(:)) = true; % include the vertices of allCutElem
numCutElemV1 = sum(isInterfaceNode);
isInterfaceNode(Np+1:end) = true;    % and the cut points and the aux points
allCutElemReoderTmp = zeros(size(node,1),1);
numCutElemV2 = sum(isInterfaceNode);
allCutElemReoderTmp(isInterfaceNode) = 1:numCutElemV2;
allCutElemReoder = allCutElemReoderTmp(allCutElem);
interfaceNode = node(isInterfaceNode,:);
ntI = -min(mesh.tLoc);
intID = find(mesh.tLoc<0);

tetElem = zeros(12*ntI,4);
tetElemLoc = zeros(12*ntI,1);
idx2cube = zeros(12*ntI,1);
tetcount = 0;

for i = 1:ntI
    tID = intID(i);
    t_e = mesh.t_e(tID,:);
    pLocK = mesh.pLoc(mesh.t(tID,:));
    vert0 = mesh.p(mesh.t(tID,:),:);
    nodeid = [allCutElemReoder(i,:),numCutElemV1-mesh.eLoc(t_e(mesh.eLoc(t_e)<0))'];
    vert2 = vert0(pLocK>0,:); % plus domain
    vert1 = vert0(pLocK<0,:); % minus domain
    intpt0 = mesh.eIntP(-mesh.eLoc(t_e(mesh.eLoc(t_e)<0)),:);
    id1 = [find(pLocK<0)', 5:4+size(intpt0,1)];
    id2 = [find(pLocK>0)', 5:4+size(intpt0,1)];
    p1 = [vert1;intpt0]; DT = delaunayTriangulation(p1); t1 = DT.ConnectivityList;
    p2 = [vert2;intpt0]; DT = delaunayTriangulation(p2); t2 = DT.ConnectivityList;
    tetElem(tetcount+1:tetcount+size(t1,1),:) = nodeid(id1(t1));
    idx2cube(tetcount+1:tetcount+size(t1,1)) = tID;
    tetElemLoc(tetcount+1:tetcount+size(t1,1)) = 1; % inside subdomain
    tetcount = tetcount+size(t1,1);
    tetElem(tetcount+1:tetcount+size(t2,1),:) = nodeid(id2(t2));
    idx2cube(tetcount+1:tetcount+size(t2,1)) = tID;
    tetElemLoc(tetcount+1:tetcount+size(t2,1)) = 2; % outside subdomain
    tetcount = tetcount+size(t2,1);
end
tetElem(tetcount+1:end,:) = [];
tetElemLoc(tetcount+1:end,:) = [];
idx2cube(tetcount+1:end,:) = [];
[tetElem, ortidx, volume] = fixorder3(interfaceNode, tetElem);
% there are some tets which are coplane, need to get rid of them
idPlane = (volume<=10^(-12));
tetElem(idPlane,:) = [];
volume(idPlane,:) = [];
tetElemLoc(idPlane,:) = [];
idx2cube(idPlane,:) = [];
localidx2globalidx = find(isInterfaceNode); % tetElem points to interfaceNode
tetElem = localidx2globalidx(tetElem);  % map to the global index of node
% there are some tets which are on the interface, but contained in the surrounding tets
% need to get rid of them (but the following algorithm may not be robust)
tetSign = vSign(tetElem);
idSliver = (sum(abs(tetSign),2)==0);
tetElem(idSliver,:) = [];
volume(idSliver,:) = [];
tetElemLoc(idSliver,:) = [];
idx2cube(idSliver,:) = [];
PolyVolume = accumarray([idx2cube,tetElemLoc],volume);
PolyVolume = PolyVolume(mesh.tLoc<0,:);
h = (PolyVolume(:,1) + PolyVolume(:,2)).^(1/3);
% the first component contains the volume of the polyhedron on the
% subdomain 1, the first component contains the volume of the polyhedron
% on the subdomain 2

toc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% try another approach for vectorization

% tic
%
% Np = size(mesh.p,1);
% Ncube = size(mesh.t,1)/6;
% node = [mesh.p;mesh.eIntP];
% allCutElem = mesh.t(mesh.tLoc<0,:);
% isCutElem = (mesh.tLoc<0);
% vSign = [mesh.pLoc;zeros(size(mesh.eIntP,1),1)];
% intID = find(mesh.tLoc<0);
% isInterfaceNode = false(size(node,1), 1);
% isInterfaceNode(allCutElem(:)) = true; % include the vertices of allCutElem
% isInterfaceNode(Np+1:end) = true;    % and the cut points and the aux points
%
% for j = 1:6
%
%     Indj = (j-1)*Ncube + 1:Ncube;
%     Indj = intersect(intID,Indj);
%     allCutElemj = mesh.t(Indj,:);
%     isInterfaceNodej = false(size(node,1), 1);
%     isInterfaceNodej(allCutElemj(:)) = true; % include the vertices of allCutElem
%     tmp = mesh.eLoc(mesh.t_e(Indj,:));
%     tmp = tmp(tmp<0);
%     %isInterfaceNodej(Np-tmp) = true;    % and the cut points and the aux points
%     interfaceNodej = node(isInterfaceNodej,:);
%
%     DT = delaunayTriangulation(interfaceNodej);
%     tetElemj = DT.ConnectivityList;
%     [tetElemj, ortidx, volume] = fixorder3(interfaceNodej, tetElemj);
%
%     localidx2globalidxj = find(isInterfaceNodej); % tetElem points to interfaceNode
%     tetElemj = localidx2globalidxj(tetElemj);  % map to the global index of node
%     bcj = (node(tetElemj(:, 1), :) + node(tetElemj(:, 2), :) ...
%         + node(tetElemj(:, 3), :) + node(tetElemj(:, 4), :))/4.0;
%     idx2cubej = FindElemId(bcj, mesh); % Get the cube index of each tetElem
%     idx2cubeBadj = (idx2cubej<=0);
%     idx2cubej = idx2cubej(~idx2cubeBadj);
%     tetElemj = tetElemj(~idx2cubeBadj,:);
%     [godElemj,godElemjIDa,godElemjIDb] = intersect(idx2cubej,Indj);
%     idx2cubej = idx2cubej(godElemjIDa);
%     tetElemj = tetElemj(godElemjIDa,:);
%
% %     tetElem = tetElem(isCutElem(idx2cube), :); % Just keep the tets in cut elem
% %     volume = volume(isCutElem(idx2cube), :);
% %     idx2cube = idx2cube(isCutElem(idx2cube));  % keep cube idx in cut elements
%
% end
%
% toc

%% Get triangular faces for interrior elements and interface
tic
localFace = [2 3 4; 1 4 3; 1 2 4; 1 3 2];
NT = size(tetElem,1);
tface = zeros(4*NT, 3);
tface2elem = zeros(4*NT, 1);
iface = zeros(4*NT, 3);
iface2elem = zeros(4*NT, 1);
% find the interior tet elements
isIntTet = min(vSign(tetElem),[], 2) == -1; % can not be == -1 as there is sliver
intTet = tetElem(isIntTet,:);
% find the corresponding cube indices
intIdx2cube = idx2cube(isIntTet);
% find triangular interface
T = auxstructure3(intTet);
neighbor = T.neighbor; % if a face is on the boundary, then the neighbor element index is itself
clear T;
tmp = (1:size(intTet, 1))';
ct = 0;
ci = 0;
for i = 1:4
    face = intTet(:, localFace(i,:));
    % find the triangle faces of polyhedron
    % 1. face and its neighbor associated to different cubes, and
    % 2. face is not on the boundary of all cut elements (which are squares)
    isPolyTriFace = ((neighbor(:, i) == tmp) | (intIdx2cube ~= intIdx2cube(neighbor(:,i)))) &...
        (sum(abs(vSign(face)), 2) > 0);% & (sum(abs(vSign(face)), 2) > 0);
    c2 = ct + sum(isPolyTriFace);
    tface((ct+1):c2,:) = face(isPolyTriFace,:);
    tface2elem((ct+1):c2,:) = intIdx2cube(isPolyTriFace); % the indices of the polyhedron
    ct = c2;
    % note that only interior elements are treated

    % find the triangle faces on interface with normal points to exterior
    % 1. face is on the boundary and
    % 2. all vertices are on the interface
    isInterface = (sum(abs(vSign(face)), 2) == 0);

    % add to interface
    c4 = ci + sum(isInterface);
    iface((ci+1):c4,:) = face(isInterface,:); % the interface tri faces.
    iface2elem((ci+1):c4,:) = intIdx2cube(isInterface); % the indices of the polyhedron
    ci = c4;
end
iface((ci+1):end,:) = [];
iface2elem((ci+1):end,:) = [];
c2old = c2;
% plot to check
%trisurf(tface(1:c2,:),node(:,1),node(:,2),node(:,3))
%trisurf(iface(1:ci,:),node(:,1),node(:,2),node(:,3))

% Find the triangular faces for exterior elements
extTet = tetElem(~isIntTet, :);
extIdx2cube = idx2cube(~isIntTet);

T = auxstructure3(extTet);
neighbor = T.neighbor;
clear T;
tmp = (1:size(extTet,1))';
for i = 1:4
    face = extTet(:, localFace(i,:));
    % find the triangle faces of polyhedron
    % 1. face and its neighbor associated to different cubes, and
    % 2. face is not on the boundary of all cut elements (which are squares)
    isPolyTriFace = ((neighbor(:, i) == tmp) | (extIdx2cube ~= extIdx2cube(neighbor(:,i)))) &...
        (sum(abs(vSign(face)), 2) > 0);
    c2 = ct + sum(isPolyTriFace);
    tface((ct+1):c2,:) = face(isPolyTriFace,:);
    tface2elem((ct+1):c2,:) = extIdx2cube(isPolyTriFace);
    % index of exterior polyhedron is append to the end of elem
    ct = c2;
end
toc

% plot to check
%trisurf(tface(c2old+1:c2,:),node(:,1),node(:,2),node(:,3))
%trisurf(iface,node(:,1),node(:,2),node(:,3))

tface((ct+1):end,:) = [];    % remove empty meomory
tface2elem((ct+1):end) = [];

face2elem = tface2elem;
face2elemLoc = [ones(c2old,1);2*ones(c2-c2old,1)]; % face2elemLoc is the same as face location
face = tface;

interfaceData.tface = tface;
interfaceData.tface2elem = tface2elem;
interfaceData.iface = iface;
interfaceData.iface2elem = iface2elem;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2.use face and face2elem data structure to compute ife functions

vSign = [mesh.pLoc;zeros(size(mesh.eIntP,1),1)];
node = [mesh.p;mesh.eIntP];
face = interfaceData.tface;
face2elem = interfaceData.tface2elem;
N = size(node, 1); Ndof = N;
NF = size(face,1); NC = size(mesh.t, 1);

NP = max(face2elem);
isCutPoly = false(NP, 1);
isCutPoly(face2elem) = true;
% this is the same as:
% isCutPoly = (mesh.tLoc<0);
% isCutPoly = isCutPoly(1:max(face2elem));
cutPolyIdx = zeros(NP, 1);
isExtPoly = false(NP,1);
isExtPoly(NC+1:NP) = true;

NP = sum(isCutPoly); % this new NP is max(-mesh.tLoc(mesh.tLoc<0));
cutPolyIdx(isCutPoly) = 1:NP;
% this is the same as:
% cutPolyIdx = zeros(size(mesh.tLoc,1), 1);
% cutPolyIdx(mesh.tLoc<0) = -mesh.tLoc(mesh.tLoc<0);
% cutPolyIdx = cutPolyIdxNew(1:max(face2elem));
face2elem = cutPolyIdx(face2elem);
% the old face2elem is the global index of interface elements
% the new one is only the index of interface elements
isExtPoly = isExtPoly(isCutPoly);
%%%clear isCutPoly cutPolyIdx

% prepare new data structure to compute projections
poly2node = sparse(face2elem(:)*ones(1,3), face(:), 1, NP,N);
poly2node = (poly2node > 0);
NV = poly2node*ones(N, 1);
%centroid = poly2node*node./[NV, NV, NV];
%KP = pde.A(centroid);

% generate local IFE basis functions which are not associated with any DoFs
% IFEbasis(:,i,:,:) for minus subdomain (i=1) or plus domain (i=2)
iface = interfaceData.iface;
tface = interfaceData.tface;
iface2elem = interfaceData.iface2elem;
[iface2elemReduce, uniID] = unique(iface2elem);
% length(uniID) should be NP
ifaceReduce = iface(uniID,:);
[IntFnormal,IntFarea,unitNormal,unittgt1,unittgt2] = facenormaltgt(node,ifaceReduce);
[Fnormal,Farea,FunitNormal,Funittgt1,Funittgt2] = facenormaltgt(node,tface);
% IFEbasis contains the associated IFE function on each face
IFEbasis = zeros(size(face2elem,1),3,3);
% for two pieces of an ife function, only the normal vector is different
% on the minus subdomain (1) ife = n; on the plus subdomain (2) ife = n*bm/bp
IFEbasis(:,1,:) = unitNormal(face2elem,:);
IFEbasis(:,2,:) = unittgt1(face2elem,:);
IFEbasis(:,3,:) = unittgt2(face2elem,:);
piece1tmp = (face2elemLoc==1);
IFEbasis(~piece1tmp,1,:) = IFEbasis(~piece1tmp,1,:)*bm/bp;
% generate boundary integral for computing projections
BIFE =  zeros(length(face2elem),3);
% BIFE contains  (beta*grad vi.n)*|f|/3 for i=1,2,3 (vi is the test function)
BIFE(:,1) = sum(squeeze(IFEbasis(:,1,:)).*Fnormal,2)/2;
BIFE(:,2) = sum(squeeze(IFEbasis(:,2,:)).*Fnormal,2)/2;
BIFE(:,3) = sum(squeeze(IFEbasis(:,3,:)).*Fnormal,2)/2;
% devided by 2 is because Fnormal is unitnormal*(2*area)
BIFE(piece1tmp,:) = bm*BIFE(piece1tmp,:)/3.;
BIFE(~piece1tmp,:) = bp*BIFE(~piece1tmp,:)/3;
% generate the matrix for solving IFE projections which is a diagonal matrix
VdiagTmp1 = (bm*PolyVolume(:,1) + bm^2/bp*PolyVolume(:,2)).^(-1);
VdiagTmp2 = (bm*PolyVolume(:,1) + bp*PolyVolume(:,2)).^(-1);
VdiagTmp3 = (bm*PolyVolume(:,1) + bp*PolyVolume(:,2)).^(-1);
Vdiag1 = VdiagTmp1(face2elem);
Vdiag2 = VdiagTmp2(face2elem);
Vdiag3 = VdiagTmp3(face2elem);
Vdiag = zeros(3*length(face2elem),1);
Vdiag(1:3:end-2) = Vdiag1;
Vdiag(2:3:end-1) = Vdiag2;
Vdiag(3:3:end) = Vdiag3;
Mdiag = spdiags(Vdiag,0,3*length(face2elem),3*length(face2elem));
BIFE = reshape(Mdiag*reshape(BIFE',[],1),3,[])';
% this updated new BIFE(i,:) contains (old)BIFE(i,:)*M^(-1) where M is the
% coefficient matrix associated with the element for this face.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute projections to IFE spaces and form the matrices
nnz = sum(NV.^2);
iiP = zeros(nnz, 1);
jjP = zeros(nnz, 1);
ssP = zeros(nnz, 1);
indexP = 0;
iiF = zeros(nnz, 1);
jjF = zeros(nnz, 1);
ssF = zeros(nnz, 1);
indexF = 0;
unv = unique(NV);
hF = h(face2elem);
b = zeros(Ndof, 1);

for kk = 1:length(unv) % group polys according to their # of nodes

    tic
    %% generate IFE functions including the IFE basis functions on each interface
    %% elements which are not associated with any DoFs, and projections of gradients
    %% and the IFE functions with matching face average
    nv = unv(kk);
    isCurrentFace = (NV(face2elem) == nv);
    CurrentFace = tface(isCurrentFace,:);
    CNF = size(CurrentFace,1);
    Currentfaceloc = face2elemLoc(isCurrentFace);
    Currentface2elem = face2elem(isCurrentFace);
    CurrentFnormal = FunitNormal(isCurrentFace,:);
    isCurrentElem = false(NP, 1);
    isCurrentElem(face2elem(isCurrentFace)) = true;
    CurrentPolyVolume = PolyVolume(isCurrentElem,:);
    % Current elem to node matrix
    currentPoly2node = poly2node(isCurrentElem, :);
    CNP = size(currentPoly2node, 1);
    currentPolyLocalIdx = zeros(NP, 1);
    currentPolyLocalIdx(isCurrentElem) = 1:CNP;

    % currentElem(i,:) contains the node index of the (current) i-th element
    [I, ~] = find(currentPoly2node');
    currentElem = reshape(I, nv, [])';
    % localIdx(i,j) contains the (current) i-th element having the
    % node(dof=1,2,3,4) on the j-th location (j is the global node index)
    localIdx = sparse(repmat((1:CNP)', 1, nv), currentElem, ones(CNP, 1)*(1:nv), CNP, N);
    clear currentPoly2node

    % Deal with current triangle face cases
    %isCTFace = isCurrentFace;
    tFace = face(isCurrentFace, :);
    NN = sum(isCurrentFace);
    subs1 = currentPolyLocalIdx(face2elem(isCurrentFace));
    subs2 = [ones(NN, 1); 2*ones(NN, 1); 3*ones(NN, 1)];
    val = BIFE(isCurrentFace,:);

    % compute the projection of gradients
    GP1 = zeros(CNP,3,nv);
    % coefficient of n t1 and t2 for the subelement on the subdomain 1
    for m = 1:3
        subs3 = full(localIdx((CurrentFace(:, m) - 1)*CNP + subs1));
        % subs1: new order (range in 1:CNP) of the element index
        % subs2: three components in the vector
        % subs3: the local index (DoF) (w.r.t. element) of the m-th node of this face
        GP1 = GP1 + accumarray([repmat(subs1,3,1), subs2, repmat(subs3, 3, 1)], val(:), [CNP, 3, nv]);
    end
    GP2 = GP1; GP2(:,1,:) = bm/bp*GP2(:,1,:);
    % generate vector basis associated with interface plane
    unitNormalCurrent = unitNormal(isCurrentElem,:);
    unittgt1Current = unittgt1(isCurrentElem,:);
    unittgt2Current = unittgt2(isCurrentElem,:);
    NTCurrent = zeros(CNP,3,3);
    NTCurrent(:,:,1) = unitNormalCurrent;
    NTCurrent(:,:,2) = unittgt1Current;
    NTCurrent(:,:,3) = unittgt2Current;
    GP1 = mytimes(NTCurrent,GP1);
    GP2 = mytimes(NTCurrent,GP2);
    % the updated GP contains three component values of the projection vector
    % generate IFE vectors on each face

    % generate some data structure on faces to compute IFE functions (not gradients)
    GP1tmp = zeros(NP,3,nv); GP2tmp = zeros(NP,3,nv);
    GP1tmp(isCurrentElem,:,:) = GP1; GP2tmp(isCurrentElem,:,:) = GP2;
    Gface = zeros(CNF,3,nv);
    Gface(Currentfaceloc==1,:,:) = GP1tmp(Currentface2elem(Currentfaceloc==1),:,:);
    Gface(Currentfaceloc==2,:,:) = GP2tmp(Currentface2elem(Currentfaceloc==2),:,:);
    clear GP1tmp GP2tmp
    % compute the IFE functions with the constant such that integral on
    % faces equal the virtual element function
    VPtmp1 = zeros(CNP,nv); % sum_i DoFi(u)|f|/3 on each face f
    VPtmp2 = zeros(CNP,nv); % grad*(xf-xm)|f|
    VPtmp3 = zeros(CNP,nv); % the surface area of each element repeated nv times
    FareaCurrent = Farea(isCurrentFace); %|f|/3
    for m = 1:3
        subs3 = full(localIdx((CurrentFace(:, m) - 1)*CNP + subs1));
        VPtmp1 = VPtmp1 + accumarray([subs1, subs3], FareaCurrent(:)/3, [CNP, nv]);
    end
    Xmf = (node(ifaceReduce(Currentface2elem,1),:) + node(ifaceReduce(Currentface2elem,2),:) +...
        node(ifaceReduce(Currentface2elem,3),:))/3; % the Xm point on the approximate interface plane
    Xf = (node(CurrentFace(:,1),:) + node(CurrentFace(:,2),:) +...
        node(CurrentFace(:,3),:))/3; % the centroid on each triangular face
    VPtmp2tmp = ((Xf - Xmf)).*FareaCurrent;
    for i = 1:nv
        tmp = sum(Gface(:,:,i).*VPtmp2tmp,2);
        VPtmp2(:,i) = accumarray(subs1, tmp,[CNP, 1]);
    end
    VPtmp3 = repmat(accumarray(subs1, FareaCurrent(:)),[1, nv]);
    IFEConst = (VPtmp1 - VPtmp2)./VPtmp3;
    clear currentPolyLocalIdx localIdx
    clear subs1 subs2 subs3
    toc

    %% form the stiffness matrices and stabilization matrices
    CurrentFace2DoF = zeros(NP,nv);
    CurrentFace2DoF(isCurrentElem,:) = currentElem;
    CurrentFace2DoF = CurrentFace2DoF(Currentface2elem,:);
    FunitNormalCurrent = FunitNormal(isCurrentFace,:);
    Funittgt1Current = Funittgt1(isCurrentFace,:);
    Funittgt2Current = Funittgt2(isCurrentFace,:);
    FNTCurrent = zeros(CNF,3,3);
    FNTCurrent(:,1,:) = Funittgt1Current;
    FNTCurrent(:,2,:) = Funittgt2Current;
    FNTCurrent(:,3,:) = FunitNormalCurrent;
    GfaceP = zeros(size(Gface));
    for n = 1:nv
        GfaceP(:,:,n) =  Gface(:,:,n) - sum(Gface(:,:,n).*CurrentFnormal,2).*CurrentFnormal;
    end
    GfaceP = mytimes(FNTCurrent,GfaceP);
    GfaceP = GfaceP(:,1:2,:);
    % compute the face gradients based on DoFs
    % generate vector basis associated with each face
    Fnode  = zeros(CNF,3,3);
    for l = 1:3
        Fnode(:,:,l) = node(CurrentFace(:,l),:);
    end
    FacePCoord = mytimes(FNTCurrent,Fnode);
    gradLambdatmp = zeros(CNF,2,3);  % barycentric coordinates
    gradLambda = zeros(CNF,2,nv);
    signtmp = ones(CNF,2); signtmp(:,1) = -1;
    gradLambdatmp(:,:,1) = (FacePCoord(:,[2,1],3) - FacePCoord(:,[2,1],2)).*signtmp./(2*Farea(isCurrentFace));
    gradLambdatmp(:,:,2) = (FacePCoord(:,[2,1],1) - FacePCoord(:,[2,1],3)).*signtmp./(2*Farea(isCurrentFace));
    gradLambdatmp(:,:,3) = -gradLambdatmp(:,:,1) - gradLambdatmp(:,:,2);
    for i = 1:nv
        for j = 1:3
            IDij = (CurrentFace(:,j) == CurrentFace2DoF(:,i));
            gradLambda(IDij,:,i) = gradLambdatmp(IDij,:,j);
        end
    end

%     for n = 1:nv
%         for m = 1:nv
%             Xn = node(currentElem(:, n), :) - centroid(isCurrentElem,:);
%             IminusP(:, n, m) = - 1/nv - dot(Xn, BP(:, :, m), 2)./volume(isCurrentElem); % stabilization
%             if(n == m)
%                 IminusP(:, n, m) = 1 + IminusP(:, n, m);
%             end
%         end
%     end

    for n = 1:nv
        for m = 1:nv
            iiP(indexP+1:indexP + CNP) = currentElem(:, n);
            jjP(indexP+1:indexP + CNP) = currentElem(:, m);
            ssP(indexP+1:indexP + CNP) = bm*dot(GP1(:,:,n), GP1(:,:,m),2).*CurrentPolyVolume(:,1)+...
                bp*dot(GP2(:,:,n), GP2(:,:,m),2).*CurrentPolyVolume(:,2);
            indexP = indexP + CNP;

            iiF(indexF+1:indexF + CNF) = CurrentFace2DoF(:, n);
            jjF(indexF+1:indexF + CNF) = CurrentFace2DoF(:, m);
            ssF(indexF+1:indexF + CNF) = hF(isCurrentFace).*...
                dot((gradLambda(:,:,m) - GfaceP(:,:,m)),(gradLambda(:,:,n) - GfaceP(:,:,n)),2).*...
                Farea(isCurrentFace);
            indexF = indexF + CNF;
        end
    end

    %% compute the RHS with only evaluating the integral at the centroid
    %% approach 1 (will give the optimal order but slightly increase the error)
%     currentElemSign = vSign(currentElem);
%     currentElemSign1 = currentElemSign<=0; nv1 = sum(currentElemSign1,2);
%     currentElemSign2 = currentElemSign>=0; nv2 = sum(currentElemSign2,2);
% %     nodeX = node(:,1); currentElemNodeX = nodeX(currentElem); clear nodeX
% %     nodeY = node(:,2); currentElemNodeY = nodeY(currentElem); clear nodeY
% %     nodeZ = node(:,1); currentElemNodeZ = nodeZ(currentElem); clear nodeZ
%     currentElemNodeX = reshape(node(currentElem',1),size(currentElem,2),[])';
%     currentElemNodeY = reshape(node(currentElem',2),size(currentElem,2),[])';
%     currentElemNodeZ = reshape(node(currentElem',3),size(currentElem,2),[])';
%     centroid1 = [sum(currentElemNodeX.*currentElemSign1,2)./nv1,...
%         sum(currentElemNodeY.*currentElemSign1,2)./nv1, ...
%         sum(currentElemNodeZ.*currentElemSign1,2)./nv1];
%     centroid2 = [sum(currentElemNodeX.*currentElemSign2,2)./nv2,...
%         sum(currentElemNodeY.*currentElemSign2,2)./nv2, ...
%         sum(currentElemNodeZ.*currentElemSign2,2)./nv2];
%     ft1 = CurrentPolyVolume(:,1).*pde.f(centroid1(:,1),centroid1(:,2),centroid1(:,3));
%     ft2 = CurrentPolyVolume(:,2).*pde.f(centroid2(:,1),centroid2(:,2),centroid2(:,3));
%
%     Xm = (node(ifaceReduce(isCurrentElem,1),:) + node(ifaceReduce(isCurrentElem,2),:) +...
%         node(ifaceReduce(isCurrentElem,3),:))/3;
%     ifeEvalcen1 = zeros(CNP,nv); ifeEvalcen2 = zeros(CNP,nv);
%     for n = 1:nv
%         ifeEvalcen1(:,n) = sum(GP1(:,:,n).*(centroid1-Xm),2) + IFEConst(:,n);
%         ifeEvalcen2(:,n) = sum(GP2(:,:,n).*(centroid2-Xm),2) + IFEConst(:,n);
%     end
%
%     Currentb = ft1.*ifeEvalcen1 + ft2.*ifeEvalcen2;
%     b = b + accumarray(currentElem(:), Currentb(:), [Ndof, 1]);

    %% approach 2
    idx2cubeNew = -mesh.tLoc(idx2cube);
    TetID = isCurrentElem(idx2cubeNew);
    currentvolume = volume(TetID);
    ElemDoF = zeros(sum(TetID),nv);
    ElemDoF(isCurrentElem,:) = currentElem;
    ElemDoF = ElemDoF(idx2cubeNew(TetID),:);
    X1 = node(tetElem(TetID,1),:);
    X2 = node(tetElem(TetID,2),:);
    X3 = node(tetElem(TetID,3),:);
    X4 = node(tetElem(TetID,4),:);
    ng = 1;
    [gx, gy, gz] = gaussPtetra(X1, X2, X3, X4, ng);
    gw = gaussWtetra(ng);
    Xm = (node(ifaceReduce(isCurrentElem,1),:) + node(ifaceReduce(isCurrentElem,2),:) +...
        node(ifaceReduce(isCurrentElem,3),:))/3;

    Bas1 = zeros(size(isCurrentElem,1),3,nv); Bas2 = Bas1;
    Bas1(isCurrentElem,:,:) = GP1;
    Bas2(isCurrentElem,:,:) = GP2;
    Bas1 = Bas1(idx2cubeNew(TetID),:,:); Bas2 = Bas2(idx2cubeNew(TetID),:,:);
    Bas = zeros(sum(TetID),3,nv);
    %%%
    IFEc0 = zeros(size(isCurrentElem,1),nv);
    IFEc0(isCurrentElem,:) = IFEConst;
    IFEc0 = IFEc0(idx2cubeNew(TetID),:);
    %%%
    Xmtet = zeros(size(isCurrentElem,1),3);
    Xmtet(isCurrentElem,:) = Xm;
    Xmtet = Xmtet(idx2cubeNew(TetID),:);
    %%%
    currentTetDoF = zeros(size(isCurrentElem,1),nv);
    currentTetDoF(isCurrentElem,:) = currentElem;
    currentTetDoF = currentTetDoF(idx2cubeNew(TetID),:);
    %%%
    piecetmp = tetElemLoc(TetID);
    Bas(piecetmp==1,:,:) = Bas1(piecetmp==1,:,:);
    Bas(piecetmp==2,:,:) = Bas2(piecetmp==2,:,:);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    % mass matrix

%     for i = 1:nv
%         for j = 1:nv
%             Basi = squeeze(Bas(:,:,i)); IFEc0i = IFEc0(:,i);
%             Basj = squeeze(Bas(:,:,j)); IFEc0j = IFEc0(:,j);
%             uhi = (gx-Xmtet(:,1)).*Basi(:,1) + (gy-Xmtet(:,2)).*Basi(:,2) +...
%                 (gz-Xmtet(:,3)).*Basi(:,3) + IFEc0i;
%             uhj = (gx-Xmtet(:,1)).*Basj(:,1) + (gy-Xmtet(:,2)).*Basj(:,2) +...
%                 (gz-Xmtet(:,3)).*Basj(:,3) + IFEc0j;
%             Currentb(:,i) = sum(gw*ft.*uhi.*uhj,2).*currentvolume;
%             bii(count+1:count+cnt) = currentTetDoF(:,i);
%             count = count+cnt;
%         end
%     end

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

    ft = pde.f(gx,gy,gz);
    Currentb = zeros(size(ft));
    cnt = size(ft,1);
    bii = nv*cnt;
    count = 0 ;
    for i = 1:nv
        Basi = squeeze(Bas(:,:,i)); IFEc0i = IFEc0(:,i);
        uhi = (gx-Xmtet(:,1)).*Basi(:,1) + (gy-Xmtet(:,2)).*Basi(:,2) +...
            (gz-Xmtet(:,3)).*Basi(:,3) + IFEc0i;
        Currentb(:,i) = sum(gw*ft.*uhi,2).*currentvolume;
        bii(count+1:count+cnt) = currentTetDoF(:,i);
        count = count+cnt;
    end
    b = b + sparse(bii,1,reshape(Currentb,[],1),Ndof,1);
    %b = b + accumarray([idx2cubeNew(TetID),currentTetDoF(:)], Currentb(:), [Ndof, 1]);

end
K = sparse(iiP, jjP, ssP, Ndof, Ndof);
S = sparse(iiF, jjF, ssF, Ndof, Ndof);

repo: ifem


text:
function [r, v] = orb2eci(mu, oev)

% convert classical orbital elements to eci state vector

% input

%  mu     = gravitational constant (km**3/sec**2)
%  oev(1) = semimajor axis (kilometers)
%  oev(2) = orbital eccentricity (non-dimensional)
%           (0 <= eccentricity < 1)
%  oev(3) = orbital inclination (radians)
%           (0 <= inclination <= pi)
%  oev(4) = argument of perigee (radians)
%           (0 <= argument of perigee <= 2 pi)
%  oev(5) = right ascension of ascending node (radians)
%           (0 <= raan <= 2 pi)
%  oev(6) = true anomaly (radians)
%           (0 <= true anomaly <= 2 pi)

% output

%  r = eci position vector (kilometers)
%  v = eci velocity vector (kilometers/second)

% Orbital Mechanics with Matlab

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

r = zeros(3, 1);
v = zeros(3, 1);

% unload orbital elements array

sma = oev(1);
ecc = oev(2);
inc = oev(3);
argper = oev(4);
raan = oev(5);
tanom = oev(6);

slr = sma * (1 - ecc * ecc);

rm = slr / (1 + ecc * cos(tanom));

arglat = argper + tanom;

sarglat = sin(arglat);
carglat = cos(arglat);

c4 = sqrt(mu / slr);
c5 = ecc * cos(argper) + carglat;
c6 = ecc * sin(argper) + sarglat;

sinc = sin(inc);
cinc = cos(inc);

sraan = sin(raan);
craan = cos(raan);

% position vector

r(1) = rm * (craan * carglat - sraan * cinc * sarglat);
r(2) = rm * (sraan * carglat + cinc * sarglat * craan);
r(3) = rm * sinc * sarglat;

% velocity vector

v(1) = -c4 * (craan * c6 + sraan * cinc * c5);
v(2) = -c4 * (sraan * c6 - craan * cinc * c5);
v(3) = c4 * c5 * sinc;

author: Sable

repo: mcbench-benchmarks


text:
%% L1QP_FeatureSign solves nonnegative quadradic programming
%% using Feature Sign.
%%
%%    min  0.5*x'*A*x+b'*x+\lambda*|x|
%%
%% [net,control]=NNQP_FeatureSign(net,A,b,control)
%%
%%
%%

function [x]=L1QP_FeatureSign_yang(lambda,A,b)

A = double(A);
b = double(b);

EPS = 1e-9;
x=zeros(size(A, 1), 1);           %coeff

grad=A*sparse(x)+b;
[ma mi]=max(abs(grad).*(x==0));

while true,

  if grad(mi)>lambda+EPS,
    x(mi)=(lambda-grad(mi))/A(mi,mi);
  elseif grad(mi)<-lambda-EPS,
    x(mi)=(-lambda-grad(mi))/A(mi,mi);
  else
    if all(x==0)
      break;
    end
  end

  while true,
    a=x~=0;   %active set
    Aa=A(a,a);
    ba=b(a);
    xa=x(a);

    %new b based on unchanged sign
    vect = -lambda*sign(xa)-ba;
    x_new= Aa\vect;
    idx = find(x_new);
    o_new=(vect(idx)/2 + ba(idx))'*x_new(idx) + lambda*sum(abs(x_new(idx)));

    %cost based on changing sign
    s=find(xa.*x_new<=0);
    if isempty(s)
      x(a)=x_new;
      loss=o_new;
      break;
    end
    x_min=x_new;
    o_min=o_new;
    d=x_new-xa;
    t=d./xa;
    for zd=s',
      x_s=xa-d/t(zd);
      x_s(zd)=0;  %make sure it's zero
%       o_s=L1QP_loss(net,Aa,ba,x_s);
      idx = find(x_s);
      o_s = (Aa(idx, idx)*x_s(idx)/2 + ba(idx))'*x_s(idx)+lambda*sum(abs(x_s(idx)));
      if o_s<o_min,
        x_min=x_s;
        o_min=o_s;
      end
    end

    x(a)=x_min;
    loss=o_min;
  end

  grad=A*sparse(x)+b;

  [ma mi]=max(abs(grad).*(x==0));
  if ma <= lambda+EPS,
    break;
  end
end

author: thomas-koehler

repo: SupER


text:
function  test_general()

    clc;
    clear;
    close all;

    %% Set algorithms
    if 0
        algorithms = gd_solver_list('ALL');
    else
        algorithms = gd_solver_list('BFGS');
        algorithms = {'L-BFGS-WOLFE','SD-STD','SD-BKT','Newton-STD','Newton-DAMP','Newton-CHOLESKY','NCG-BTK'};
        %algorithms = {'TR-NEWTON','TR-QUASI-NEWTON','L-BFGS-WOLFE'};
        algorithms = {'TR-DOGLEG-QUASI-NEWTON','TR-CAUCHY-QUASI-NEWTON','TR-TRUSTONE-QUASI-NEWTON','TR-CG-QUASI-NEWTON','L-BFGS-WOLFE'};
        algorithms = {'SD-BKT', 'NCG-BTK','TR-DOGLEG-QUASI-NEWTON','L-BFGS-BKT', 'BFGS-H-BKT','BB'};
    end

    %% set cost, gradient and hessian
    if 0
        % f = x(1)^2 + 2 * x(2)^2
        % This is equivalent to the quadratic case wheree A = [1,0; 0,2] and b = [0;0].
        d = 2;
        f = @(x) x(1)^2 + 2 * x(2)^2;
        g = @(x) [2 * x(1); 4 * x(2)];
        h = @(x) [2 0; 0 4];
        w_init = [2; 1];

    elseif 0
        % f = 100 * x(1)^4 + 0.01 * x(2)^4
        d = 2;
        f = @(x) 100 * x(1)^4 + 0.01 * x(2)^4;
        g = @(x) [400 * x(1)^3; 0.04 * x(2)^3];
        h = @(x) [1200 * x(1)^2, 0; 0 0.12 * x(2)^2];
        w_init = [1; 1];

    elseif 0
        % f = sqrt(1+(1).^2) + sqrt(1+x(2).^2)
        d = 2;
        f = @(x) sqrt(1+x(1).^2) + sqrt(1+x(2).^2);
        g = @(x) [x(1)/sqrt(x(1).^2+1); x(2)/sqrt(x(2).^2+1)];
        h = @(x) [1/(x(1).^2+1).^1.5,0; 0, 1/(x(2).^2+1).^1.5];
        w_init = [10; 10];

    elseif 1
        d = 2;
        f = @(x)7/5+(x(1)+2*x(2)+2*x(1)*x(2)-5*x(1)^2-5*x(2)^2)/(5* exp(x(1)^2+x(2)^2));
        g = @(x)[(1+2*x(2)-10*x(1)-2*x(1)*(x(1)+ 2*x(2)+2*x(1)*x(2)-5*x(1)^2 -5*x(2)^2))/(5* exp(x(1)^2+x(2)^2));
                (2+2*x(1)-10*x(2)-2*x(2)*(x(1)+2*x(2)+ 2*x(1)*x(2)-5*x(1)^2 -5*x(2)^2))/(5* exp(x(1)^2+x(2)^2))];
        h = @(x) [];

        w_init = [0;0.5];

    end

    %% define problem definitions
    problem = general(f, g, h, [], d, [], [], [], [], [], [], []);

    % Calculate the solution
    fminunc_options.TolFun = 1e-36;
    [w_opt,f_opt] = fminunc(f, w_init, fminunc_options);
    fprintf('%f, %f, %f\n', w_opt(1), w_opt(2), f_opt);

    % initialize
    w_list = cell(1);
    info_list = cell(1);

    %% perform algorithms
    for alg_idx=1:length(algorithms)
        fprintf('\n\n### [%02d] %s ###\n\n', alg_idx, algorithms{alg_idx});

        clear options;
        % general options for optimization algorithms
        options.w_init = w_init;
        options.tol_gnorm = 1e-10;
        options.max_iter = 100;
        options.verbose = true;
        options.f_opt = f_opt;
        options.store_w = true;

        switch algorithms{alg_idx}
            case {'TR-DOGLEG-NEWTON'}

                options.H_mode = 'NEWTON';
                [w_list{alg_idx}, info_list{alg_idx}] = tr(problem, options);

            case {'TR-DOGLEG-QUASI-NEWTON'}

                options.H_mode = 'QUASI-NEWTON';
                options.subprob_solver = 'DOGLEG';
                [w_list{alg_idx}, info_list{alg_idx}] = tr(problem, options);

            case {'TR-CAUCHY-QUASI-NEWTON'}

                options.H_mode = 'QUASI-NEWTON';
                options.subprob_solver = 'CAUCHY';
                [w_list{alg_idx}, info_list{alg_idx}] = tr(problem, options);

            case {'TR-TRUSTONE-QUASI-NEWTON'}

                options.H_mode = 'QUASI-NEWTON';
                options.subprob_solver = 'TRUSTONE';
                [w_list{alg_idx}, info_list{alg_idx}] = tr(problem, options);

            case {'TR-CG-QUASI-NEWTON'}

                options.H_mode = 'QUASI-NEWTON';
                options.subprob_solver = 'CG';
                [w_list{alg_idx}, info_list{alg_idx}] = tr(problem, options);

%             case {'TR-Lanczos-QUASI-NEWTON'}
%
%                 options.H_mode = 'QUASI-NEWTON';
%                 options.subprob_solver = 'Lanczos';
%                 [w_list{alg_idx}, info_list{alg_idx}] = tr(problem, options);

            case {'SD-STD'}

                options.step_alg = 'fix';
                options.step_init = 1;
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'SD-BKT'}

                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'SD-EXACT'}

                options.step_alg = 'exact';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'SD-WOLFE'}

                options.step_alg = 'strong_wolfe';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'SD-SCALE-EXACT'}

                options.sub_mode = 'SCALING';
                options.step_alg = 'exact';
                [w_list{alg_idx}, info_list{alg_idx}] = sd(problem, options);

            case {'Newton-STD'}

                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            case {'Newton-DAMP'}

                options.sub_mode = 'DAMPED';
                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            case {'Newton-CHOLESKY'}

                options.sub_mode = 'CHOLESKY';
                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = newton(problem, options);

            case {'CG-PRELIM'}

                options.sub_mode = 'PRELIM';
                options.step_alg = 'exact';
                %options.beta_alg = 'PR';
                [w_list{alg_idx}, info_list{alg_idx}] = cg(problem, options);

            case {'CG-BKT'}

                options.sub_mode = 'STANDARD';
                options.step_alg = 'backtracking';
                %options.beta_alg = 'PR';
                [w_list{alg_idx}, info_list{alg_idx}] = cg(problem, options);

            case {'CG-EXACT'}

                options.sub_mode = 'STANDARD';
                options.step_alg = 'exact';
                %options.beta_alg = 'PR';
                [w_list{alg_idx}, info_list{alg_idx}] = cg(problem, options);

            case {'CG-PRECON-EXACT'}

                options.sub_mode = 'PRECON';
                % diagonal scaling
                options.M = diag(diag(A));
                options.step_alg = 'exact';
                options.beta_alg = 'PR';

                [w_list{alg_idx}, info_list{alg_idx}] = cg(problem, options);

            case {'NCG-BTK'}

                options.sub_mode = 'STANDARD';
                options.step_alg = 'backtracking';
                options.beta_alg = 'PR';
                [w_list{alg_idx}, info_list{alg_idx}] = ncg(problem, options);

            case {'NCG-WOLFE'}

                options.sub_mode = 'STANDARD';
                options.step_alg = 'strong_wolfe';
                options.beta_alg = 'PR';
                [w_list{alg_idx}, info_list{alg_idx}] = ncg(problem, options);

            case {'BFGS-H-BKT'}

                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'BFGS-H-EXACT'}

                options.step_alg = 'exact';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'BFGS-B-BKT'}

                options.step_alg = 'backtracking';
                options.update_mode = 'B';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'BFGS-B-EXACT'}

                options.step_alg = 'exact';
                options.update_mode = 'B';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'DAMPED-BFGS-BKT'}

                options.step_alg = 'backtracking';
                options.update_mode = 'Damping';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'DAMPED-BFGS-EXACT'}

                options.step_alg = 'exact';
                options.update_mode = 'Damping';
                [w_list{alg_idx}, info_list{alg_idx}] = bfgs(problem, options);

            case {'L-BFGS-BKT'}

                options.step_alg = 'backtracking';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'L-BFGS-EXACT'}

                options.step_alg = 'exact';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'L-BFGS-WOLFE'}

                options.step_alg = 'strong_wolfe';
                [w_list{alg_idx}, info_list{alg_idx}] = lbfgs(problem, options);

            case {'BB'}

                options.step_alg = 'exact';
                [w_list{alg_idx}, info_list{alg_idx}] = bb(problem, options);

            case {'SGD'}

                options.batch_size = 1;
                options.step = 0.1 * options.batch_size;
                %options.step_alg = 'decay';
                options.step_alg = 'fix';

                [w_list{alg_idx}, info_list{alg_idx}] = sgd(problem, options);

            otherwise
                warn_str = [algorithms{alg_idx}, ' is not supported.'];
                warning(warn_str);
                w_list{alg_idx} = '';
                info_list{alg_idx} = '';
        end

    end

    %% plot all
    close all;

    % display iter vs cost/gnorm
    display_graph('iter','optimality_gap', algorithms, w_list, info_list);
    display_graph('time','optimality_gap', algorithms, w_list, info_list);
    %display_graph('iter','gnorm', algorithms, w_list, info_list);

    % draw convergence sequence
    w_history = cell(1);
    cost_history = cell(1);
    for alg_idx=1:length(algorithms)
        w_history{alg_idx} = info_list{alg_idx}.w;
        cost_history{alg_idx} = info_list{alg_idx}.cost;
    end
    draw_convergence_sequence(problem, w_opt, algorithms, w_history, cost_history);

end

author: hiroyuki-kasai

repo: SGDLibrary


text:
function datr = ksFilter(buff, ops)

if ~isfield(ops, 'b1')
    if isfield(ops,'fslow')&&ops.fslow<ops.fs/2
        [b1, a1] = butter(3, [ops.fshigh/ops.fs,ops.fslow/ops.fs]*2, 'bandpass');
    else
        [b1, a1] = butter(3, ops.fshigh/ops.fs*2, 'high');
    end
else
    b1 = ops.b1; a1 = ops.a1;
end

if ops.GPU
    dataRAW = gpuArray(buff);
else
    dataRAW = buff;
end
dataRAW = dataRAW';
dataRAW = single(dataRAW);
dataRAW = dataRAW(:, ops.chanMap.chanMap);

% subtract the mean from each channel
dataRAW = dataRAW - mean(dataRAW, 1);

datr = filter(b1, a1, dataRAW);
datr = flipud(datr);
datr = filter(b1, a1, datr);
datr = flipud(datr);

% CAR, common average referencing by median
if getOr(ops, 'CAR', 1)
    datr = datr - median(datr, 2);
end

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function [varargout]=repop(varargin)
% Replicating arithemetical and logical operators.
%
% [Z]=repop(X,operator,Y [,options])
%
% Does element by element operations on X and Y where non-same sized
% dimensions are implicity wrapped round to match the size of the larger
% to give a result matrix Z with size max(size(X),size(Y));
%
% What this means is that if you give a [Nx1] and a [1xM] input you get
% a [NxM] output, or if you give it a [Nx1xQ] and a [NxMx1xP] you get
% [NxMxQxP] output etc..
% Note that the size of the input *completely and uniquely* determines the
% size of the output.
%
% In general this is at least 2x faster than the equivalent matlab code
% using repmats and has the advantage of requiring no additional memory.
%
% Example Usage:
%     X = randn(10000,10);                  % example signal with data in rows
%     stdX = repop(X,'-',mean(X,1));        % subtract the mean vector
%     stdX = repop(stdX,'/',std(stdX,0,1)); % divide by std-deviation
%
% Operator can be one of:
%
% Arthemetical -- returns a double matrix
%   '+','.+',plus   - Implicitly repmatted elementwise addition
%   '-','.-',minus  - Implicitly repmatted elementwise addition
%   '*','.*',times  - Implicitly repmatted elementwise multiplication
%   '^','.^',power  - Implicitly repmatted elementwise raise X to power Y
%   '\','.\',ldivide- Implicitly repmatted elementwise divide Y by X
%   '/','./',rdivide- Implicitly repmatted elementwise divide X by Y
%   'min'           - Implicitly repmatted elementwise min of X by Y
%   'max'           - Implicitly repmatted elementwise max of X by Y
%
% Relational -- returns a logical matrix
% N.B. for complex inputs the <,>,<=,>= operators are based upon abs(x)
% (not real(x) as in matlab)
%   '==',eq         - Implicitly repmatted elementwise equality
%   '~=',ne         - Implicitly repmatted elementwise dis-equality
%   '<' ,lt         - Implicitly repmatted elementwise less than
%   '>' ,gt         - Implicitly repmatted elementwise greater than
%   '<=',le         - Implicitly repmatted elementwise less than equal
%   '>=',ge         - Implicitly repmatted elementwise greater than equal
%
% N.B. the operator can go in any of the 3 argument positions, i.e.
% these are all valid: repop(X,'-',Y), repop(X,Y,'-'), repop('-',X,Y)
%
% The optional final argument is a string of single letter switches
% consisting off
%  'm'  -- allows replication of non-unit dimensions if the larger
%          dimensions size is an integer multiple of the smaller ones
%  'n'  -- allow replication of non-unit dimensions in *all* cases
%  'i'  -- perform "inplace" operation.  This means that we use a
%          *dangerous* matlab *hack* to perform the operation *without*
%          allocating new memory for the output, but by simply overwriting
%          the memory used for X in the input.  Thus the following code is
%          a memory (and time) efficient way to increment X.
%              X = repop(X,'+',1,'i');
%
%
% Class support of input P:
%     float: double, single
%
% SEE ALSO:   repop_testcases rplus rminus rtimes rpower rldivide rrdivide
%             req rne rlt rgt rle rge
%
% Copyright 2006-     by Jason D.R. Farquhar (jdrf@zepler.org)
% Permission is granted for anyone to copy, use, or modify this
% software and accompanying documents for any uncommercial
% purposes, provided this copyright notice is retained, and note is
% made of any changes that have been made. This software and
% documents are distributed without any warranty, express or
% implied
%
% Inspired by code from Douglas M. Schwarz and & Aki Vehtari.

% The rest of this code is a mex-hiding mechanism which compilies the mex if
% this runs and recursivly calls itself.
% Based upon code from: http://theoval.sys.uea.ac.uk/matlab
cwd  = pwd; % store the current working directory
name = mfilename('fullpath'); % get the directory where we are
% find out what directory it is defined in
name(name=='\')='/'; % deal with dos'isms
dir=name(1:max(find(name == '/')-1)); % dir is everything before final '/'
try % try changing to that directory
   cd(dir);
catch   % this should never happen, but just in case!
   cd(cwd);
   error(['unable to locate directory containing ''' name '.m''']);
end

try % try recompiling the MEX file
   fprintf(['Compiling ' mfilename ' for first use\n']);
   mex('ddrepop.c','dsrepop.c','sdrepop.c','ssrepop.c','repop_util.c','repop_mex.c','mxInfo.c','mxInfo_mex.c','-O','-output',mfilename);
   fprintf('done\n');
catch
   % this may well happen happen, get back to current working directory!
   cd(cwd);
   error('unable to compile MEX version of ''%s''%s\n%s%s', name, ...
         ', please make sure your', 'MEX compiler is set up correctly', ...
         ' (try ''mex -setup'').');
end

cd(cwd); % change back to the current working directory
rehash;  % refresh the function and file system caches

% recursively invoke MEX version using the same input and output arguments
[varargout{1:nargout}] = feval(mfilename, varargin{:});

% bye bye...

return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text: shallPrint = false;

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function problem = prepare_speed_problem(speedpath,speedimgpath,imgidx,nrOutliers)
load([speedpath,'keypoints-2.mat']);
imgs         = dir([speedimgpath,'*.jpg']);
nrImgs       = length(imgs);
model_3d     = double( model_3d' * 0.00000586 );
nrLandmarks  = 11;
camera_K     = [3003.41, 0, 960;
                0, 3003.41, 600;
                0, 0, 1];
camera_Kinv  = inv(camera_K);
img          = imread([speedimgpath,imgs(imgidx).name]);
keypoints    = double( squeeze(prediction_2d(imgidx,:,:))' );
keypoints(:,1:nrOutliers) = [1920 * rand(1,nrOutliers); 1200 * rand(1,nrOutliers)];
keypoints_2d = camera_Kinv * [keypoints;ones(1,nrLandmarks)];
keypoints_2d = keypoints_2d(1:2,:);
R_gt         = quat2rotm( gt_pose(imgidx,1:4) );
t_gt         = gt_pose(imgidx,5:end)';

problem.N    = nrLandmarks;
problem.X    = model_3d;
problem.x    = keypoints_2d;
problem.FOV  = 90;
problem.translationBound = 25.0;
problem.noiseBoundSq     = (0.2)^2;
problem.noiseBound       = sqrt(problem.noiseBoundSq);
problem.depthBound       = 0.0;
problem.pixNoiseBound    = 1e-2;
problem.R_gt             = R_gt;
problem.t_gt             = t_gt;
problem.img              = img;

problem.keypoints_pix    = keypoints;

problem.camera_K         = camera_K;

fprintf('SPEED: N: %d, num outliers: %d, noiseBound: %g.\n',...
    problem.N,nrOutliers,problem.noiseBound);

if ~check_translation(problem.t_gt,problem.translationBound,problem.FOV)
    error('Problem assumption wrong.')
end

end

author: MIT-SPARK

repo: CertifiablyRobustPerception


text:
function F = subsref(F,Y)
%subsref           Overloaded indexing

F = subsref(lmi(F),Y);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
%% Discretisation of Orientation Space
% represents a set of orientations
%%
%

author: mtex-toolbox

repo: mtex

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


text:
function lambda = givens_eigenvalues ( n )

%*****************************************************************************80
%
%% GIVENS_EIGENVALUES returns the eigenvalues of the GIVENS matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );

  for i = 1 : n
    angle = ( 2 * i - 1 ) * pi / ( 4 * n );
    lambda(i,1) = 0.5 / ( cos ( angle ) )^2;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [v_0_vec, v_vec, a, e, p] = LambertBattin_v2(r_0_vec, r_vec, delta_t, t_m, N)
%
% USAGE = [v_0_vec, v_vec, a, e, p] = LambertBattin_v2(r_0_vec, r_vec, delta_t, t_m, N)
%

% Function constants
% mu = 3.986e5;
mu = 398600.4418; % Sun gravitational constant (km^3/sec^2)
tol = 1e-6;
iter_max = 500;

r_0 = norm(r_0_vec);
r = norm(r_vec);

cos_delta_nu = dot(r_0_vec, r_vec) / (r_0 * r);
sin_delta_nu = t_m * sqrt(1 - cos_delta_nu^2);
delta_nu = atan2(sin_delta_nu, cos_delta_nu);
if delta_nu < 0

    delta_nu = 2 * pi + delta_nu;

end

c = sqrt(r_0^2 + r^2 - 2 * r_0 * r * cos_delta_nu);
s = (r_0 + r + c) / 2;

epsilon = (r - r_0) / r_0;
tan_sq_2w = (epsilon^2 / 4) / (sqrt(r / r_0) + (r / r_0) * (2 + sqrt(r / r_0)));
r_0p = sqrt(r_0 * r) * ((cos(delta_nu / 4))^2 + tan_sq_2w);

sin_sq_delta_nu_4 = (sin(delta_nu / 4))^2;
cos_sq_delta_nu_4 = (cos(delta_nu / 4))^2;
cos_delta_nu_2 = cos(delta_nu / 2);
if (delta_nu > 0) && (delta_nu < pi)

    l = (sin_sq_delta_nu_4 + tan_sq_2w) / (sin_sq_delta_nu_4 + tan_sq_2w + cos_delta_nu_2);

elseif (delta_nu > pi) && (delta_nu < 2*pi)

    l = (cos_sq_delta_nu_4 + tan_sq_2w - cos_delta_nu_2) / (cos_sq_delta_nu_4 + tan_sq_2w);

else

    error('Delta_nu out of range');

end

m = (mu * delta_t^2) / (8 * r_0p^3);

if N == 0

    % Initial estimate for x
    x_new = l;

    loop_flag = 1;
    i_iter = 1;
    while loop_flag == 1

        x_old = x_new;
        [x_new, y] = original_sub(x_old, l, m, N);

        if abs(x_new - x_old) < tol

            loop_flag = 0;

        elseif i_iter < iter_max

            i_iter = i_iter + 1;

        else

            loop_flag = 0;
            x_new = NaN;
            y = NaN;

        end

    end

    x_f = x_new;
    y_f = y;

else % N > 0

    loop_flag_R = 1;
    loop_flag_L = 1;
    i_iter_R = 1;
    i_iter_L = 1;

    % Initial estimate for x_R
    x_new = l;

    while loop_flag_R == 1

        x_old = x_new;
        [x_new, y] = reverse_sub(x_old, l, m, N);

        if y < 1

            % In divergent region for the reversed successive substitution
            % Begin original successive substitution to find x_L
            x_new = l;
            x_old = x_new;
            [x_new, y] = original_sub(x_old, l, m, N);

            while loop_flag_L == 1

                x_old = x_new;
                [x_new, y] = original_sub(x_old, l, m, N);

                if abs(x_new - x_old) < tol

                    loop_flag_L = 0;

                elseif i_iter_R < iter_max

                    i_iter_R = i_iter_R + 1;

                else

                    loop_flag_L = 0;
                    x_new = NaN;
                    y = NaN;

                end

            end

            y_L = y;
            x_L = x_new;

            % Reset initial guess for x_R
            x_new = x_new / 2;
            x_old = x_new - 2 * tol;

        end

        if abs(x_new - x_old) < tol

            loop_flag_R = 0;

        elseif i_iter_R < iter_max

            i_iter_R = i_iter_R + 1;

        else

            loop_flag_R = 0;
            x_new = NaN;
            y = NaN;

        end

    end

    y_R = y;
    x_R = x_new;

    % Only continue if x_L and y_L were not found in above iteration
    if loop_flag_L == 1

        % Initial estimate for x_L & y
        x_new = 2 * x_R;

        while loop_flag_L == 1

            x_old = x_new;
            [x_new, y] = original_sub(x_old, l, m, N);

            if abs(x_new - x_old) < tol

                loop_flag_L = 0;

            elseif i_iter_L < iter_max

                i_iter_L = i_iter_L + 1;

            else

                loop_flag_L = 0;
                x_new = NaN;
                y = NaN;

            end

        end

        y_L = y;
        x_L = x_new;

    end

    y_f = [y_R; y_L];
    x_f = [x_R; x_L];

end

if N == 0

    n = 1;

else

    n = 2;

end

a = zeros(n,1);
e = zeros(n,1);
p = zeros(n,1);
v_vec = zeros(n,3);
v_0_vec = zeros(n,3);
for i = 1:n

    a(i) = mu * delta_t^2 / (16 * r_0p^2 * x_f(i) * y_f(i)^2);
    e(i) = sqrt((epsilon^2 + 4 * (r / r_0) * (sin(delta_nu / 2))^2 * ((l - x_f(i)) / (l + x_f(i)))^2) / (epsilon^2 + 4 * (r / r_0) * (sin(delta_nu / 2))^2));
    p(i) = (4 * r_0p^2 * r_0 * r * y_f(i)^2 * (1 + x_f(i))^2 * (sin(delta_nu / 2))^2) / (mu * delta_t^2);
    if a(i) > 0

        f = 1 - (r / p(i)) * (1 - cos_delta_nu);
        g = (r * r_0 * sin_delta_nu) / sqrt(mu * p(i));
        g_dot = 1 - (r_0 / p(i)) * (1 - cos_delta_nu);

    else

        alpha_h = 2 * asinh(sqrt(s / (-2 * a(i))));
        beta_h = 2 * asinh(sqrt((s - c) / (-2 * a(i))));

        delta_H = alpha_h - beta_h;

        f = 1 - a(i) * (1 - cosh(delta_H)) / r;
        g = delta_t - sqrt(-a(i)^3 / mu) * (sinh(delta_H) - delta_H);
        g_dot = 1 - a(i) * (1 - cosh(delta_H)) / r;

    end

    v_0_vec(i,:) = (1 / g) .* (r_vec - f .* r_0_vec);
    v_vec(i,:) = (1 / g) .* (g_dot .* r_vec - r_0_vec);

end

%--------------------------------------------------------------------------

function [y] = cubic_solv(c1, c2, c3)

P = c2 - c1^2 / 3;
Q = c3 + (2 * c1^3 - 9 * c1 * c2) / 27;
D = Q^2 / 4 + P^3 / 27;

T = (-Q / 2 + sqrt(D))^(1/3) + (-Q/2 - sqrt(D))^(1/3);
y = T - c1 / 3;

%--------------------------------------------------------------------------

function [x_new, y] = original_sub(x_old, l, m, N)

E = 2 * atan(sqrt(x_old));
c = -m * (N * pi + E - sin(E)) / (4 * (tan(E / 2))^3);
y = cubic_solv(-1, 0, c);
x_new = sqrt(((1 - l) / 2)^2 + (m / y^2)) - ((1 + l) / 2);

%--------------------------------------------------------------------------

function [x_new, y] = reverse_sub(x_old, l, m, N)

y = sqrt(m / ((l + x_old) * (1 + x_old)));
E_0 = 2 * atan(sqrt(x_old));
E = y2E(E_0, y, m, N);
x_new = (tan(E / 2))^2;

%--------------------------------------------------------------------------

function [E] = y2E(E_0, y, m, N)

tol = 1e-6;
iter_max = 200;

q = 4 * (y^3 - y^2) / m;
h = (N * pi + E_0 - sin(E_0)) / (tan(E_0 / 2))^3 - q;

if h < 0

    loop_flag_h = 1;
    i_iter_h = 1;
    while loop_flag_h == 1

        E_0 = E_0 / 2;
        h = (N * pi + E_0 - sin(E_0)) / (tan(E_0 / 2))^3 - q;

        if h >= 0

            loop_flag_h = 0;

        elseif i_iter_h < iter_max

            i_iter_h = i_iter_h + 1;

        else

            loop_flag_h = 0;
            E_0 = NaN;

        end

    end

end

i_iter_E = 1;
loop_flag_E = 1;
while loop_flag_E == 1

    h = (N * pi + E_0 - sin(E_0)) / (tan(E_0 / 2))^3 - q;
    h_prime = (-3 * (cos(E_0 / 2))^2) * (N * pi + E_0 - sin(E_0)) / (2 * (sin(E_0 / 2))^4) + (1 - cos(E_0)) / (tan(E_0 / 2))^3;
    E = E_0 - h / h_prime;

    if abs(E - E_0) < tol

        loop_flag_E = 0;

    elseif i_iter_E < iter_max

        i_iter_E = i_iter_E + 1;
        E_0 = E;

    else

        loop_flag_E = 0;
        E = NaN;

    end

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function pass = test_partitionCombine( )
% Test the partition and combine methods.

tol = 1e3*chebfunpref().cheb2Prefs.chebfun2eps;

% Check that partitioning an empty diskfun gives two empty diskfuns
f = diskfun;
[feven,fodd] = partition(f);
pass(1) = isempty(feven) & isempty(fodd);

feven = diskfun;
fodd = diskfun;
f = combine(feven,fodd);
pass(2) = isempty(f);

% Check that an even diskfun is partitioned correctly.
f = diskfun(@(x,y) sin(pi*x.*y));  % Strictly even/pi-periodic
[feven,fodd] = partition(f);
pass(3) = isequal(feven,f);
pass(4) = isempty(fodd);

% Check that odd diskfun is partitioned correctly.
f = diskfun(@(x,y) sin(pi*x));  % Strictly odd/anti-periodic
[feven,fodd] = partition(f);
pass(5) = isequal(fodd,f);
pass(6) = isempty(feven);

% Check that the even and odd terms are partitioned correctly
fe = @(x,y) sin(pi*x.*y);  % Strictly even/pi-periodic
fo = @(x,y) sin(pi*x);  % Strictly odd/anti-periodic
f = diskfun(@(x,y) fe(x,y) + fo(x,y));
[feven,fodd] = partition(f);
pass(7) = norm(diskfun(fe)-feven) < tol;
pass(8) = norm(diskfun(fo)-fodd) < tol;

% Check that combine puts an even diskfun and empty spherfun back
% together.
feven = diskfun(@(x,y) sin(pi*x.*y));  % Strictly even/pi-periodic
fodd = diskfun;
f = combine(feven,fodd);
pass(9) = norm(feven-f) < tol;

% Check that combine puts an even diskfun and empty spherfun back
% together.
feven = diskfun;
fodd = diskfun(@(x,y) sin(pi*x));  % Strictly odd/anti-periodic
f = combine(feven,fodd);
pass(10) = norm(fodd-f) < tol;

% Check that combine puts non-empty odd and even diskfuns back together.
fe = @(x,y) sin(pi*x.*y);  % Strictly even/pi-periodic
fo = @(x,y) sin(pi*x);  % Strictly odd/anti-periodic
fcombine = diskfun(@(x,y) fe(x,y) + fo(x,y));
f = combine(diskfun(fe),diskfun(fo));
pass(11) = norm(fcombine-f) < tol;

% Check that combine cannot put together diskfuns that are not strictly
% odd or even.
try
    f = diskfun(@(x,y) sin(pi*x.*y) + sin(pi*x));
    g = combine(f,f);
    pass(12) = false;
catch ME
    pass(12) = strcmp(ME.identifier, 'CHEBFUN:DISKFUN:combine:parity');
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [p,t,pfun,tfun] = fixmesh(p,t,pfun,tfun)

%*****************************************************************************80
%
%  FIXMESH: Ensure that triangular mesh data is consistent.
%
%  [p,t,pfun,tfun] = fixmesh(p,t,pfun,tfun);
%
%  p     : Nx2 array of nodal XY coordinates, [x1,y1; x2,y2; etc]
%  t     : Mx3 array of triangles as indices, [n11,n12,n13; n21,n22,n23;
%          etc]
%  pfun  : (Optional) NxK array of nodal function values. Each column in
%          PFUN corresponds to a dependent function, PFUN(:,1) = F1(P),
%          PFUN(:,2) = F2(P) etc, defined at the nodes.
%  tfun  : (Optional) MxK array of triangle function values. Each column in
%          TFUN corresponds to a dependent function, TFUN(:,1) = F1(T),
%          TFUN(:,2) = F2(T) etc, defined on the triangles.
%
% The following checks are performed:
%
%  1. Nodes not refereneced in T are removed.
%  2. Duplicate nodes are removed.
%  3. Triangles are ordered counter-clockwise.
%  4. Triangles with an area less than 1.0e-10*eps*norm(A,'inf')
%     are removed
%
%  Author:
%
%    Darren Engwirda
%

TOL = 1.0e-10;

if (nargin<4)
   tfun = [];
   if (nargin<3)
      pfun = [];
      if nargin<2
         error('Wrong number of inputs');
      end
   end
elseif (nargin>4)
   error('Wrong number of inputs');
end
if (nargout>4)
   error('Wrong number of outputs');
end
if (numel(p)~=2*size(p,1))
   error('P must be an Nx2 array');
end
if (numel(t)~=3*size(t,1))
   error('T must be an Mx3 array');
end
if (any(t(:))<1) || (max(t(:))>size(p,1))
   error('Invalid T');
end
if ~isempty(pfun)
   if (size(pfun,1)~=size(p,1)) || (ndims(pfun)~=2)
      error('PFUN must be an NxK array');
   end
end
if ~isempty(tfun)
   if (size(tfun,1)~=size(t,1)) || (ndims(tfun)~=2)
      error('TFUN must be an Mxk array');
   end
end

% Remove duplicate nodes
[i,i,j] = unique(p,'rows');
if ~isempty(pfun)
   pfun = pfun(i,:);
end
p = p(i,:);
t = reshape(j(t),size(t));

% Triangle area
A = triarea(p,t);
Ai = A<0.0;
Aj = abs(A)>TOL*norm(A,'inf');

% Flip node numbering to give a counter-clockwise order
t(Ai,[1,2]) = t(Ai,[2,1]);

% Remove zero area triangles
t = t(Aj,:);
if ~isempty(tfun)
   tfun = tfun(Aj,:);
end

% Remove un-used nodes
[i,j,j] = unique(t(:));
if ~isempty(pfun)
   pfun = pfun(i,:);
end
p = p(i,:);
t = reshape(j,size(t));

end      % fixmesh()


text:
function a = rectangle_adj ( row_num, col_num, n )

%*****************************************************************************80
%
%% RECTANGLE_ADJ returns the RECTANGLE_ADJ matrix.
%
%  Discussion:
%
%    This is the adjacency matrix for a set of points arranged in
%    a ROW_NUM by COL_NUM grid.
%
%  Diagram:
%
%     1---5---9
%     |   |   |
%     2---6--10
%     |   |   |
%     3---7--11
%     |   |   |
%     4---8--12
%
%  Example:
%
%    ROW_NUM = 4
%    COL_NUM = 3
%
%    0  1  0  0    1  0  0  0    0  0  0  0
%    1  0  1  0    0  1  0  0    0  0  0  0
%    0  1  0  1    0  0  1  0    0  0  0  0
%    0  0  1  0    1  0  0  1    0  0  0  0
%
%    1  0  0  0    0  1  0  0    1  0  0  0
%    0  1  0  0    1  0  1  0    0  1  0  0
%    0  0  1  0    0  1  0  1    0  0  1  0
%    0  0  0  1    0  0  1  0    0  0  0  1
%
%    0  0  0  0    1  0  0  0    0  1  0  0
%    0  0  0  0    0  1  0  0    1  0  1  0
%    0  0  0  0    0  0  1  0    0  1  0  1
%    0  0  0  0    0  0  0  1    0  0  1  0
%
%  Properties:
%
%    A is integral, therefore det ( A ) is integral, and
%    det ( A ) * inverse ( A ) is integral.
%
%    A is a zero/one matrix.
%
%    A is block tridiagonal.
%
%    A is an adjacency matrix.
%
%    A is related to the "LIGHTS_OUT" matrix.
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ROW_NUM, COL_NUM, the number of rows and
%    columns in the rectangle.
%
%    Input, integer N, the order of A.
%
%    Output, real A(N,N), the matrix.
%
  a = zeros ( n, n );

  for j_block = 1 : col_num

    j = ( j_block - 1 ) * row_num + 1;

    for i_block = 1 : col_num

      i = ( i_block - 1 ) * row_num + 1;

      if ( j_block == i_block )
        a(i:i+row_num-1,j:j+row_num-1) = line_adj ( row_num );
      elseif ( abs ( j_block - i_block ) == 1 )
        a(i:i+row_num-1,j:j+row_num-1) = identity ( row_num, row_num );
      else
        a(i:i+row_num-1,j:j+row_num-1) = zero ( row_num, row_num );
      end

    end
  end

  return
end


text:
%% regression
d = 100;
beta = 1e-1;
X = rand(1,d);
w = randn;
b = randn;
t = w'*X+b+beta*randn(1,d);
x = linspace(min(X),max(X),d);   % test data

%% RVM regression by Mackay fix point update
[model,llh] = rvmRegFp(X,t);
plot(llh);
[y, sigma] = linRegPred(model,x,t);
figure
plotCurveBar(x,y,sigma);
hold on;
plot(X,t,'o');
hold off

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
function loser = updategrid(x, minHoleTime, minHoleSize, maxHoleTime, maxHoleSize)
global F
global alive
global snakeSize
global holeCountDown
global lastPos

res = size(F);
loser = 0;
for i = alive
    pixels = getpixels(x(i,:), snakeSize);
    index1 = size(F,1) * lastPos(:,i*2) + lastPos(:,i*2-1);
    index2 = size(F,1) * pixels(:,2) + pixels(:,1);
    F(index1) = 0;

    if any(pixels(:) < 1) || any(pixels(:,1) >= res(1)) || any(pixels(:,2) >= res(2))
        loser = i;
        alive(alive == i) = [];
        F(index1) = i;
        return
    end
    if sum(full(F(index2))) > 1
        alive(alive == i) = [];
        F(index1) = i;
        return
    end

    F(index2) = i;
    if holeCountDown(i,1) ~= 0
        F(index1) = i;
    else
        F(index1) = 0;
        holeCountDown(i,1) = 1;
        holeCountDown(i,2) = holeCountDown(i,2) - 1;
    end
    if holeCountDown(i,2) == 0
        F(index2) = i;
        holeCountDown(i,1) = round(rand() * (maxHoleTime - minHoleTime)) + minHoleTime;
        holeCountDown(i,2) = round(rand() * (maxHoleSize - minHoleSize)) + minHoleSize;
    end
    lastPos(:,i*2-1:i*2) = pixels;
end


text:
function [Rnew,Ki]=finvupdate(a,R,k_newvsold,knew)

k=k_newvsold;
k22=knew;
e=R'*R*k;
g= 1/   sqrt(k22-k'*R'*R*k);
Rnew=[ R,zeros(length(k),1);-g*e',g];
if(nargout==2)
    Ki=Rnew'*Rnew;
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function test_issue1568

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_databrowser

global ft_default
ft_default.representation = 'table';

%%

ntrials = 10;
nchans = 16;
fsample = 1000;
trialgap = 0;

data = [];
for i=1:nchans
  data.label{i} = sprintf('chan%d', i);
end
for i=1:ntrials
  data.time{i} = (1:fsample)/fsample;
  data.trial{i} = randn(nchans, fsample);
  data.sampleinfo(i,1) = (i-1)*fsample + 1 + (i*trialgap); % create some gaps between the trials
  data.sampleinfo(i,2) = data.sampleinfo(i,1) + fsample - 1;
end

if strcmp(ft_default.representation, 'numeric')
  data.trialinfo = [1:ntrials]';
else
  trialnum = arrayfun(@num2str, 1:ntrials, 'UniformOutput', false)';
  data.trialinfo = table(trialnum);
end

% this is for visualising the provenance
data.cfg.version.name = 'a';
data.cfg.previous.version.name = 'b';
data.cfg.previous.previous.version.name = 'c';
data.cfg.previous.previous.previous.version.name = 'd';
data.cfg.previous.previous.previous.previous.version.name = 'e';

%%

close all

cfg = [];
cfg.ylim = [-6 6];
% cfg.viewmode = 'vertical';

tmp1.begsample = 101;
tmp1.endsample = 110;
tmp1.channel = 'chan1';

tmp2.begsample = 121;
tmp2.endsample = 130;
tmp2.channel = 'chan2';

tmp3.begsample = 141;
tmp3.endsample = 150;
tmp3.channel = 'chan3';

% let's use three different types of artifacts, just like visual, eog, muscle, clip, etc.
cfg.artfctdef.clip.artifact = struct2table(tmp1);
cfg.artfctdef.threshold.artifact = struct2table(tmp2);
cfg.artfctdef.other.artifact = struct2table(tmp3);
cfg.artfctdef.someoption = 'no';

event(1).type     = 'ttl1';
event(1).value    = 1;
event(1).sample   = 701;
event(1).duration = 1;
event(1).offset   = 0;

event(2).type     = 'ttl2';
event(2).value    = 2;
event(2).sample   = 721;
event(2).duration = 10;
event(2).offset   = 0;

event(3).type     = 'ttl3';
event(3).value    = 3;
event(3).sample   = 741;
event(3).duration = 10;
event(3).offset   = [];

cfg.event = event;

cfg.plotartifactlabels = 'type';
cfg.ploteventlabels = 'type=value';
ft_databrowser(cfg, data);


text:
function varargout = mydategps (epoch)
    error(nargoutchk(0, 3, nargout, 'struct'));

    sec_per_hour = 60^2;
    sec_per_day  = 24 * sec_per_hour;
    sec_per_week = 7 * sec_per_day;
    %epoch0 = mydatenum([1980 01 13]);  % WRONG!
    epoch0 = mydatenum([1980 01 06]);
    % as per <http://facility.unavco.org/data/glossary.html#GPS week>

    sow = mydatesow(epoch);
    dow = floor(sow ./ sec_per_day);
    sod = sow - dow .* sec_per_day;
    week = (epoch - epoch0) ./ sec_per_week;
    week = floor(week);

    switch nargout
    case {0, 1}
        varargout = {sow};
    case 2
        varargout = {week, sow};
    case 3
        varargout = {week, dow, sod};
    end
end

%!test
%! % using SOPAC's date converter:
%! % <http://sopac.ucsd.edu/scripts/convertDate.cgi>
%! year = 2009;
%! mon = 10;
%! day = 28;
%! week = 1555;
%! dow = 3;
%! sod = 100*rand;
%! vec = [year, mon, day, 0, 0, sod];
%! num = mydatenum(vec);
%! [week2, dow2, sod2] = mydategps (num);
%! %[ [week2, dow2, sod2]; [week, dow, sod]; ...
%! %  [week2, dow2, sod2] - [week, dow, sod] ]  % DEBUG
%! myassert(week2, week);
%! myassert(dow2, dow);
%! myassert(sod2, sod, sqrt(eps));

%!test
%! % using SOPAC's date converter:
%! % <http://sopac.ucsd.edu/scripts/convertDate.cgi>
%! year = 2000;
%! mon = 01;
%! day = 01;
%! week = 1042;
%! dow = 6;
%! sod = 100*rand;
%! vec = [year, mon, day, 0, 0, sod];
%! num = mydatenum(vec);
%! [week2, dow2, sod2] = mydategps (num);
%! %[ [week2, dow2, sod2]; [week, dow, sod]; ...
%! %  [week2, dow2, sod2] - [week, dow, sod] ]  % DEBUG
%! myassert(week2, week);
%! myassert(dow2, dow);
%! myassert(sod2, sod, sqrt(eps));

%!test
%! % using SOPAC's date converter:
%! % <http://sopac.ucsd.edu/scripts/convertDate.cgi>
%! year = 2010;
%! mon = 06;
%! day = 05;
%! week = 1586;
%! dow = 6;
%! sod = 0;
%! vec = [year, mon, day, 0, 0, sod];
%! num = mydatenum(vec);
%! [week2, dow2, sod2] = mydategps (num);
%! %[ [week2, dow2, sod2]; [week, dow, sod]; ...
%! %  [week2, dow2, sod2] - [week, dow, sod] ]  % DEBUG
%! myassert(week2, week);
%! myassert(dow2, dow);
%! myassert(sod2, sod, sqrt(eps));


text:
% CSparse: a Concise Sparse matrix Package.
%
%   Matrices used in CSparse must in general be either sparse and real,
%   or dense vectors.  Ordering methods can accept any sparse matrix.
%
%   cs_add       - sparse matrix addition.
%   cs_amd       - approximate minimum degree ordering.
%   cs_chol      - sparse Cholesky factorization.
%   cs_cholsol   - solve A*x=b using a sparse Cholesky factorization.
%   cs_counts    - column counts for sparse Cholesky factor L.
%   cs_dmperm    - maximum matching or Dulmage-Mendelsohn permutation.
%   cs_dmsol     - x=A\b using the coarse Dulmage-Mendelsohn decomposition.
%   cs_dmspy     - plot the Dulmage-Mendelsohn decomposition of a matrix.
%   cs_droptol   - remove small entries from a sparse matrix.
%   cs_esep      - find an edge separator of a symmetric matrix A
%   cs_etree     - elimination tree of A or A'*A.
%   cs_gaxpy     - sparse matrix times vector.
%   cs_lsolve    - solve a sparse lower triangular system L*x=b.
%   cs_ltsolve   - solve a sparse upper triangular system L'*x=b.
%   cs_lu        - sparse LU factorization, with fill-reducing ordering.
%   cs_lusol     - solve Ax=b using LU factorization.
%   cs_make      - compiles CSparse for use in MATLAB.
%   cs_multiply  - sparse matrix multiply.
%   cs_nd        - generalized nested dissection ordering.
%   cs_nsep      - find a node separator of a symmetric matrix A.
%   cs_permute   - permute a sparse matrix.
%   cs_print     - print the contents of a sparse matrix.
%   cs_qr        - sparse QR factorization.
%   cs_qleft     - apply Householder vectors on the left.
%   cs_qright    - apply Householder vectors on the right.
%   cs_qrsol     - solve a sparse least-squares problem.
%   cs_randperm  - random permutation.
%   cs_sep       - convert an edge separator into a node separator.
%   cs_scc       - strongly-connected components of a square sparse matrix.
%   cs_scc2      - cs_scc, or connected components of a bipartite graph.
%   cs_sparse    - convert a triplet form into a sparse matrix.
%   cs_sqr       - symbolic sparse QR factorization.
%   cs_symperm   - symmetric permutation of a symmetric matrix.
%   cs_transpose - transpose a real sparse matrix.
%   cs_updown    - rank-1 update/downdate of a sparse Cholesky factorization.
%   cs_usolve    - solve a sparse upper triangular system U*x=b.
%   cs_utsolve   - solve a sparse lower triangular system U'*x=b.
%   cspy         - plot a matrix in color.
%   ccspy        - plot the connected components of a matrix.

% Example:
%   help cs_add

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

% helper function:
%   cs_must_compile - return 1 if source code f must be compiled, 0 otherwise


text:
function p_data = p07_data ( dim_num, data_num )

%*****************************************************************************80
%
%% P07_DATA returns the data for problem p07.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the spatial dimension of the dependent
%    variables.
%
%    Input, integer DATA_NUM, the number of data points.
%
%    Output, real P_DATA(DIM_NUM,DATA_NUM), the data.
%
  p_data = [ ...
   0.0, 1.0; ...
   1.0, 2.0; ...
   4.0, 2.0; ...
   5.0, 1.0 ]';

  return
end


text:
function sphere_llq_grid_line_count_test ( )

%*****************************************************************************80
%
%% SPHERE_LLQ_GRID_LINE_COUNT_TEST tests SPHERE_LLQ_GRID_LINE_COUNT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 April 2015
%
%  Author:
%
%    John Burkardt
%
  lat_num = 3;
  long_num = 4;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SPHERE_LLQ_GRID_LINE_COUNT_TEST\n' );
  fprintf ( 1, '  SPHERE_LLQ_GRID_LINE_COUNT counts the lines used for a\n' );
  fprintf ( 1, '  grid based on quadrilaterals defined by latitude and longitude\n' );
  fprintf ( 1, '  lines on a sphere in 3D.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '     LAT_NUM    LONG_NUM   LINE_NUM\n' );
  for lat_num = 1 : 2 : 17
    fprintf ( 1, '\n' );
    long_num = 1;
    for long_log = 1 : 4
      long_num = long_num * 2;
      line_num = sphere_llq_grid_line_count ( lat_num, long_num );
      fprintf ( 1, '  %8d  %8d  %8d\n', lat_num, long_num, line_num );
    end
  end

  return
end


text:
function []=PS_calc_ifg_std
% PS_CALC_IFG_STD() calculate std for each ifg
%
%   Andy Hooper, June 2006
%
%   ======================================================
%   09/2006 AH: small baselines added
%   02/2010 AH: More informative info displayed
%   ======================================================

fprintf('\nEstimating noise standard deviation (degrees)...\n')

small_baseline_flag=getparm('small_baseline_flag');

load psver
psname=['ps',num2str(psver)];
phname=['ph',num2str(psver)];
pmname=['pm',num2str(psver)];
bpname=['bp',num2str(psver)];
ifgstdname=['ifgstd',num2str(psver)];

ps=load(psname);
pm=load(pmname);
bp=load(bpname);

if exist([phname,'.mat'],'file')
    phin=load(phname);
    ph=phin.ph;
    clear phin
else
    ph=ps.ph;
end

n_ps=length(ps.xy);
master_ix=sum(ps.master_day>ps.day)+1;

if strcmpi(small_baseline_flag,'y')
    ph_diff=angle(ph.*conj(pm.ph_patch).*exp(-j*(repmat(pm.K_ps,1,ps.n_ifg).*bp.bperp_mat)));
else
    bperp_mat=[bp.bperp_mat(:,1:ps.master_ix-1),zeros(ps.n_ps,1,'single'),bp.bperp_mat(:,ps.master_ix:end)];
    ph_patch=[pm.ph_patch(:,1:master_ix-1),ones(n_ps,1),pm.ph_patch(:,master_ix:end)];
    ph_diff=angle(ph.*conj(ph_patch).*exp(-j*(repmat(pm.K_ps,1,ps.n_ifg).*bperp_mat+repmat(pm.C_ps,1,ps.n_ifg))));
end

ifg_std=[sqrt(sum(ph_diff.^2)/n_ps)*180/pi]';
if strcmpi(small_baseline_flag,'y')
  for i=1:ps.n_ifg
    fprintf('%3d %s_%s %3.2f\n',i,datestr(ps.ifgday(i,1)),datestr(ps.ifgday(i,2)),ifg_std(i))
  end
else
  for i=1:ps.n_ifg
    fprintf('%3d %s %3.2f\n',i,datestr(ps.day(i)),ifg_std(i))
  end
end
fprintf('\n')

save(ifgstdname,'ifg_std');

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function M = femUnk2Qud(fe,domain)
%+========================================================================+
%|                                                                        |
%|              OPENFEM - LIBRARY FOR FINITE ELEMENT METHOD               |
%|           openFem is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal & Francois Alouges (c) 2017-2018.          |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%|             francois.alouges@polytechnique.edu                         |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : femUnknown.m                                  |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal & François Alouges            |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Unknowns and reduction matrix for constrained |
%|  `---'  |                finite elements                               |
%+========================================================================+

% Unknown to degrees of freedom matrix
[~,P] = fe.unk;

% Surfacic restriction (trace)
if (size(fe.msh.elt,2) == size(domain.msh.elt,2) + 1)
    bound  = fe.msh.bnd;
    P      = restriction(fe,bound) * P;
    fe.msh = bound;
end

% Degrees of freedom to quadrature matrix
if strcmp(fe.typ(1),'P')
    M = femLagrangePn(fe,domain);
elseif strcmp(fe.typ,'RWG')
    M = femRaoWiltonGlisson(fe,domain);
elseif strcmp(fe.typ,'NED')
    M = femNedelec(fe,domain);
else
    error('fem.m : unavailable case')
end

% Unknown to quadrature
if iscell(M)
    M{1} = M{1} * P;
    M{2} = M{2} * P;
    M{3} = M{3} * P;
else
    M = M * P;
end
end


text:
filename='Gripping_tetrahedra_coarse';
 ptype = 'MACRO';
method = 'SIMPALL';
materialType = 'ISOTROPIC';
initial_case = 'full';
cost = {'nonadjoint_compliance'};
weights = [1, 0.1];
constraint = {'volumeConstraint'};
constraint_case = 'INEQUALITY';
optimizer = 'MMA'; incrementFactor = 1;
filterType = 'P1';

nsteps = 1;
Vfrac_final = 0.4;
optimality_final = 1e-3;
constr_final =1e-3;

Vfrac_initial = 1;
optimality_initial = 1e-3;
constr_initial = 1e-3;
TOL.rho_plus = 1;
TOL.rho_minus = 0;
TOL.E_plus = 1;
TOL.E_minus = 1e-3;
TOL.nu_plus = 1/3;
TOL.nu_minus = 1/3;


text:
function asa245_test03 ( )

%*****************************************************************************80
%
%% ASA245_TEST03 demonstrates the use of GAMMA_LN.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 September 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA245_TEST03:\n' );
  fprintf ( 1, '  GAMMALN computes the logarithm of the\n' );
  fprintf ( 1, '  Gamma function.\n' );
  fprintf ( 1, '  GAMMALN is a builtin function in MATLAB.\n' );
  fprintf ( 1, '  We compare the result to tabulated values.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '          X                     ' );
  fprintf ( 1, 'FX                        FX2\n' );
  fprintf ( 1, '                                ' );
  fprintf ( 1, '(Tabulated)               (GAMMALN)                DIFF\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = gamma_log_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fx2 = gammaln ( x );

    fprintf ( 1, '  %24.16f  %24.16e  %24.16e  %10.4e\n', ...
    x, fx, fx2, abs ( fx - fx2 ) );

  end

  return
end


text:
function [varargout] = wavefilter(wname, type)
%WAVEFILTER Create wavelet decomposition and reconstruction filters.
%   [VARARGOUT] = WAVEFILTER(WNAME, TYPE) returns the decomposition
%   and/or reconstruction filters used in the computation of the
%   forward and inverse FWT (fast wavelet transform).
%
%   EXAMPLES:
%     [ld, hd, lr, hr] = wavefilter('haar') Get the low and highpass
%                                           decomposition (ld, hd)
%                                           and reconstruction
%                                           (lr, hr) filters for
%                                           wavelet 'haar'.
%     [ld, hd] = wavefilter('haar','d')     Get decomposition filters
%                                           ld and hd.
%     [lr, hr] = wavefilter('haar','r')     Get reconstruction
%                                           filters lr and hr.
%
%   INPUTS:
%     WNAME             Wavelet Name
%     ---------------------------------------------------------
%     'haar' or 'db1'   Haar
%     'db4'             4th order Daubechies
%     'sym4'            4th order Symlets
%     'bior6.8'         Cohen-Daubechies-Feauveau biorthogonal
%     'jpeg9.7'         Antonini-Barlaud-Mathieu-Daubechies
%
%     TYPE              Filter Type
%     ---------------------------------------------------------
%     'd'               Decomposition filters
%     'r'               Reconstruction filters
%
%   See also WAVEFAST and WAVEBACK.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% Check the input and output arguments.
if (nargin == 1 && nargout ~= 4) || (nargin == 2 && nargout ~= 2)
   error('Invalid number of output arguments.');
end

if nargin == 1 && ~ischar(wname)
   error('WNAME must be a string.');
end

if nargin == 2 && ~ischar(type)
   error('TYPE must be a string.');
end

% Create filters for the requested wavelet.
switch lower(wname)
case {'haar', 'db1'}
   ld = [1 1]/sqrt(2);     hd = [-1 1]/sqrt(2);
   lr = ld;                hr = -hd;

case 'db4'
   ld = [-1.059740178499728e-002 3.288301166698295e-002 ...
         3.084138183598697e-002 -1.870348117188811e-001 ...
         -2.798376941698385e-002 6.308807679295904e-001 ...
         7.148465705525415e-001 2.303778133088552e-001];
   t = (0:7);
   hd = ld;    hd(end:-1:1) = cos(pi * t) .* ld;
   lr = ld;    lr(end:-1:1) = ld;
   hr = cos(pi * t) .* ld;

case 'sym4'
   ld = [-7.576571478927333e-002 -2.963552764599851e-002 ...
         4.976186676320155e-001 8.037387518059161e-001 ...
         2.978577956052774e-001 -9.921954357684722e-002 ...
         -1.260396726203783e-002 3.222310060404270e-002];
   t = (0:7);
   hd = ld;    hd(end:-1:1) = cos(pi * t) .* ld;
   lr = ld;    lr(end:-1:1) = ld;
   hr = cos(pi * t) .* ld;

case 'bior6.8'
   ld = [0 1.908831736481291e-003 -1.914286129088767e-003 ...
         -1.699063986760234e-002 1.193456527972926e-002 ...
         4.973290349094079e-002 -7.726317316720414e-002 ...
         -9.405920349573646e-002 4.207962846098268e-001 ...
         8.259229974584023e-001 4.207962846098268e-001 ...
         -9.405920349573646e-002 -7.726317316720414e-002 ...
         4.973290349094079e-002 1.193456527972926e-002 ...
         -1.699063986760234e-002 -1.914286129088767e-003 ...
         1.908831736481291e-003];
   hd = [0 0 0 1.442628250562444e-002 -1.446750489679015e-002 ...
         -7.872200106262882e-002 4.036797903033992e-002 ...
         4.178491091502746e-001 -7.589077294536542e-001 ...
         4.178491091502746e-001 4.036797903033992e-002 ...
         -7.872200106262882e-002 -1.446750489679015e-002 ...
         1.442628250562444e-002 0 0 0 0];
   t = (0:17);
   lr = cos(pi * (t + 1)) .* hd;
   hr = cos(pi * t) .* ld;

case 'jpeg9.7'
   ld = [0 0.02674875741080976 -0.01686411844287495 ...
         -0.07822326652898785 0.2668641184428723 ...
         0.6029490182363579 0.2668641184428723 ...
         -0.07822326652898785 -0.01686411844287495 ...
         0.02674875741080976];
   hd = [0 0.09127176311424948 -0.05754352622849957 ...
         -0.5912717631142470 1.115087052456994 ...
         -0.5912717631142470 -0.05754352622849957 ...
         0.09127176311424948 0 0];
   t = (0:9);
   lr = cos(pi * (t + 1)) .* hd;
   hr = cos(pi * t) .* ld;

otherwise
   error('Unrecognizable wavelet name (WNAME).');
end

% Output the requested filters.
if (nargin == 1)
   varargout(1:4) = {ld, hd, lr, hr};
else
   switch lower(type(1))
   case 'd'
      varargout = {ld, hd};
   case 'r'
      varargout = {lr, hr};
   otherwise
      error('Unrecognizable filter TYPE.');
   end
end

author: dipum

repo: dipum-toolbox

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


text:
function Image_out=Intensity_adjust(Image_in);
%This code extends the intensity limits of  an image from [low_in; high_in]
%to [0;255]. It is most convenient for applications where operations such
%as edge detection or black to white transformations are to be performed
%on a sequence of images. Threshold values are fixed for different images.
%
%   CLASS SUPPORT
%   -------------
%   The input is an RGB image.
%
%
%   Syntaxes:
%   Image_out=Intensity_adjust('sky.jpg');
%   The output is an RGB image of the same size as the initial image.
%   This code is written by:
%                         Nassim Khaled
%                         American University of Beirut
%                         Research and Graduate Assistant
%
%  Developed under Matlab 7                       Date:  April,2005

a=imread(Image_in);
s=rgb2gray(a);
b=double(s);
low_in=min(b(:))/255;
high_in=max(b(:))/255;
Image_out = imadjust(a,[low_in; high_in],[0;1]);


text:
function [tx, categoryFeatures] = getBoundaryClassifierFeatures(bndinfo, X, ind)
% X is the raw data
% ind is the set of indices for which the features should be computed
%
% tx:
%   Edge features (1-6)
%        1:  Pb
%        2:  Length / Perimeter
%        3:  Smoothness
%        4:  Angle
%      5-6:  Continuity
%      7-8:  Convexity (area and ratio) - not used
%        9:  Chain length
%   Region features (7-17)+3
%    10-11:  Area
%       12:  Color Mean Difference
%       13:  Color Entropy Difference
%       14:  Gradient Entropy Difference - not used
%    15-16:  Position (x,y)
%    17-18:  Extent Overlap (x, y)
%    an additional 10 features of position/overlap
%   Geometry features (16-39)+3
%    19-28:  Geometric Context Mean
%    29-33:  Geometric Context Difference
%    34   :  Geometric Context Sum Abs Difference
%    35-36:  Geometric Context Most Likely Label (G V or S)
%    37-40:  Depth under- and over-estimates for each side
%    31-43:  Depth, min1-min2, max1-max2, min(max12) - max(min12)
%    44-47:  Depthcol, each sp, diff, abs diff

ng = 5; % five geometric classes

ndata = numel(ind);

[imh, imw] = size(bndinfo.wseg);

tx = zeros([ndata 60], 'single');

if isempty(ind)
    return;
end

spLR = bndinfo.edges.spLR;
s1 = spLR(ind, 1);
s2 = spLR(ind, 2);
wseg = bndinfo.wseg;

categoryFeatures = [];
f = 0;

%% Edge features
tx(:, f+1) = X.edge.pb(ind);

perim = zeros(bndinfo.nseg, 1);
for k = 1:numel(X.edge.length)
    perim(spLR(k, 1)) = perim(spLR(k, 1)) + X.edge.length(k);
    perim(spLR(k, 2)) = perim(spLR(k, 2)) + X.edge.length(k);
end
minperim = min([perim(s1) perim(s2)], [], 2);
tx(:, f+2) = X.edge.length(ind) ./ minperim; % edge length / perim of smaller region

% juncts = bndinfo.edges.junctions(ind, :);
% jpos1 = bndinfo.junctions.position(juncts(:, 1), :);
% jpos2 = bndinfo.junctions.position(juncts(:, 2), :);
% directLength = abs(jpos2(:, 1)-jpos1(:,1)) + abs(jpos2(:, 2)-jpos1(:,2));
tx(:, f+3) = X.edge.smoothness(ind); % measure of smoothess

theta = X.edge.theta;
% discrete angle
tx(:, f+4) = max(ceil((mod(theta(ind),2*pi) / (pi*2) * 16 - 1E-10)),1);
categoryFeatures(end+1) = f+4;

% relative angle (continuity)
%theta = mod([theta ; theta+pi]/pi*180, 360);
theta1 = mod(X.edge.thetaStart*180/pi, 360);
theta2 = mod(X.edge.thetaEnd*180/pi, 360);
maxc = zeros(ndata, 2);
eadj = bndinfo.edges.adjacency;
ne = bndinfo.ne;
for k = 1:ndata
    ki = ind(k);
    ra = abs(theta2(ki)-theta1(eadj{ki}));
    ra = ra - 180*(ra>180);
    if isempty(ra), maxc(k,1) = 0;
    else maxc(k,1) = min(ra);
    end
    ra = mod(abs(theta2(ne+ki)-theta1(eadj{ne+ki})), 180+1E-5);
    if isempty(ra), maxc(k,2) = 0;
    else maxc(k,2) = min(ra);
    end
end
tx(:, f+(5:6)) = [min(maxc, [], 2) max(maxc, [], 2)];

area1 = X.region.area(s1);
area2 = X.region.area(s2);
%tx(:, f+7) = X.edge.convArea(ind) ./ min([area1 area2], [], 2);
%tx(:, f+8) = 0; %X.edge.convRatio;

%ind2 = (X.edge.edge2chain(ind)>0);
%tx(ind2, f+9) = X.edge.chainsize(X.edge.edge2chain(ind(ind2)));
tx(:, f+9) = X.edge.edge2chain(ind);

f = f + 9;

%% Region features

% area

tx(:, f+(1:2)) = [min([area1 area2], [], 2) max([area1 area2], [], 2)];

% color
tx(:, f+3) = sqrt(sum((X.region.colorMean(s1, :)-X.region.colorMean(s2, :)).^2, 2));

ch = X.region.colorHist+1E-10;
for k = 1:ndata
    h1 = ch(s1(k), :);  e1 = sum(-log(h1).*h1);
    h2 = ch(s2(k), :);  e2 = sum(-log(h2).*h2);
    e12 = (area1(k)*e1 + area2(k)*e2)/(area1(k)+area2(k));
    h3 = (area1(k)*h1 + area2(k)*h2)/(area1(k)+area2(k));
    e3 = sum(-log(h3).*h3);
    tx(k, f+4) = e3-e12;
end

% gradient
% ch = X.region.gradHist+1E-10;
% for k = 1:ndata
%     h1 = ch(s1(k), :);  e1 = sum(-log(h1).*h1);
%     h2 = ch(s2(k), :);  e2 = sum(-log(h2).*h2);
%     e12 = (area1(k)*e1 + area2(k)*e2)/(area1(k)+area2(k));
%     h3 = (area1(k)*h1 + area2(k)*h2)/(area1(k)+area2(k));
%     e3 = sum(-log(h3).*h3);
%     tx(k, f+5) = e3-e12;
% end

% position
tx(:, f+6) = (X.region.y(s1, 3))-(X.region.y(s2, 3)); % difference of tops
tx(:, f+7) = (X.region.y(s1, 1))-(X.region.y(s2, 1)); % difference of bottoms
tx(:, f+8) = (X.region.y(s1, 3))-(X.region.y(s2, 1)); % top1 - bottom2
tx(:, f+9) = (X.region.y(s1, 1))-(X.region.y(s2, 3)); % bottom1 - top2
tx(:, f+10) = X.region.y(s1, 3) - X.region.y(s1, 1); % top1 - bottom1
tx(:, f+11) = X.region.y(s2, 3) - X.region.y(s2, 1); % top2 - bottom2
tx(:, f+12) = (X.region.x(s1, 1))-(X.region.x(s2, 1)); % left1 - left2
tx(:, f+13) = (X.region.x(s1, 3))-(X.region.x(s2, 3)); % right1 - right2
tx(:, f+14) = X.region.x(s1, 3) - X.region.x(s1, 1); % right1 - left1
tx(:, f+15) = X.region.x(s2, 3) - X.region.x(s2, 1); % right2 - left2

%tx(:, f+6) = (X.region.x(s1, 2)-X.region.x(s2, 2)) / imw;
%tx(:, f+7) = (X.region.y(s1, 2)-X.region.y(s2, 2)) / imh;

% x alignment
x1 = X.region.x(s1, [1 3]);
x2 = X.region.x(s2, [1 3]);
tx(:, f+16) = (min([x1(:, 2) x2(:, 2)], [], 2)-max([x1(:, 1) x2(:, 1)], [], 2)) ./ ...
    (max([x1(:, 2) x2(:, 2)], [], 2)-min([x1(:, 1) x2(:, 1)], [], 2));

% determine whether regions are x-aligned at boundary
jpos1 = ceil(bndinfo.junctions.position(bndinfo.edges.junctions(ind, 1), :));
jpos2 = ceil(bndinfo.junctions.position(bndinfo.edges.junctions(ind, 2), :));
if jpos1(:, 1)>jpos2(:, 1) % make jpos1 the left-most junction
    tmp = jpos2;
    jpos2 = jpos1;
    jpos1 = tmp;
end
tx(:, f+17) = (jpos2(:, 1)-jpos1(:, 1))/imw; % boundary width
jpos1(:, 1) = max(jpos1(:, 1)-3, 1);  % go a little to left of left junction
jpos2(:, 1) = min(jpos2(:, 1)+3, imw); % go a little to right or right junction
jpos1(:, 2) = min(jpos1(:, 2), imh);
jpos2(:, 2) = min(jpos2(:, 2), imh);
js1 = wseg((jpos1(:, 1)-1)*imh + jpos1(:, 2)); % region slightly to left
js2 = wseg((jpos2(:, 1)-1)*imh + jpos2(:, 2)); % region slightly to right
tx(:, f+18) = (js1~=s1) & (js1~=s2) & (js2~=s1) & (js2~=s2); % whether x-aligned

% y overlap
y1 = X.region.y(s1, [1 3]);
y2 = X.region.y(s2, [1 3]);
tx(:, f+19) = (min([y1(:, 2) y2(:, 2)], [], 2)-max([y1(:, 1) y2(:, 1)], [], 2)) ./ ...
    (max([y1(:, 2) y2(:, 2)], [], 2)-min([y1(:, 1) y2(:, 1)], [], 2));

f = f + 19;

%% 3D Geometry features

% geometric context features
gc = X.region.geomContext;

tx(:, f+(1:ng)) = gc(s1, :);
tx(:, f+ng+(1:ng)) = gc(s2, :);
tx(:, f+2*ng+(1:ng)) = tx(:, f+(1:ng))-tx(:, f+ng+(1:ng));
tx(:, f+3*ng+1) = sum(abs(tx(:, f+2*ng+(1:ng))), 2)/2;

[maxval, maxlab] = max([gc(:, 1) sum(gc(:, 2:4), 2) gc(:, 5)], [], 2);
tx(:, f+3*ng+2) = (maxlab(s1)-1)*3+ maxlab(s2);
categoryFeatures(end+1) = f+3*ng+2;

f = f + 17;

% relative depth
tx(:, f+(1:4)) = [X.edge.depthmin(ind, 1:2)  X.edge.depthmax(ind, 1:2)];
tx(:, f+5) = X.edge.depthmin(ind, 1)-X.edge.depthmin(ind,2);
tx(:, f+6) = X.edge.depthmax(ind, 1)-X.edge.depthmax(ind,2);
tx(:, f+7) = min(X.edge.depthmax(ind, :), [], 2)- max(X.edge.depthmin(ind, :), [], 2);

tx(:, f+(8:9)) = X.region.depthcol([s1 s2]);
tx(:, f+10) = tx(:, f+8)-tx(:, f+9);
tx(:, f+11) = abs(tx(:, f+10));

%% Geometric T-junctions
% Ground-Vertical-Ground junctions
% In future make this faster by checking chains to see when g/v transitions
% to v/v
gvs1 = zeros(bndinfo.ne, 1);  gvs2 = zeros(bndinfo.ne, 1);
gvs1(ind) = maxlab(s1);
gvs2(ind) = maxlab(s2);
is_gv = (gvs1==1 & gvs2==2) | (gvs1==2 & gvs2==1);
is_vv = (gvs1==2 & gvs2==2);
ejuncts = bndinfo.edges.junctions;
juncts = cell(bndinfo.nj, 1);
jsize = zeros(size(juncts));
for k = 1:size(ejuncts, 1)
    j1 = ejuncts(k, 1);  j2 = ejuncts(k, 2);
    juncts{j1}(end+1) = k;
    juncts{j2}(end+1) = k;
    jsize(j1) = jsize(j1) + 1;
    jsize(j2) = jsize(j2) + 1;
end
juncts = cell2mat(juncts(jsize==3));
isGeomJunction = sum(is_gv(juncts), 2)==2 & sum(is_vv(juncts), 2)==1;
juncts = juncts(isGeomJunction, :);
nGeomJuncts = zeros(bndinfo.nseg, 1);
e2chain = X.edge.edge2chain;
spLR = [bndinfo.edges.spLR ; bndinfo.edges.spLR(:, [2 1])];
for k = 1:size(juncts, 1) % check that g/v --> v/v angle is less than g/v -->
    gve = juncts(k, is_gv(juncts(k, :)));
    % if ground on left, then reverse
    vve = juncts(k, is_vv(juncts(k, :)));
    if maxlab(spLR(gve(1), 1))==1, gve(1) = gve(1)+ne;  end
    if maxlab(spLR(gve(2), 1))==1, gve(2) = gve(2)+ne;  end
    if any(e2chain(vve)==e2chain(gve)) && (e2chain(vve)>0) % foreground on left
        chainind = sort(X.edge.chains{e2chain(vve)});
        chainind1 = chainind(chainind<=ne);
        chainind2 = chainind(chainind>ne)-ne;
        tmpind = ismember_sorted(ind, chainind1);
        tx(tmpind, f+12) = 1;
        tmpind = ismember_sorted(ind, chainind2);
        tx(tmpind, f+12) = 2;
        nGeomJuncts(spLR(vve, 1)) = nGeomJuncts(spLR(vve, 1))+1;
    elseif any(e2chain(vve+ne)==e2chain(gve)) && (e2chain(vve+ne)>0) % foreground on right
        chainind = sort(X.edge.chains{e2chain(vve+ne)});
        chainind1 = chainind(chainind<=ne);
        chainind2 = chainind(chainind>ne)-ne;
        tmpind = ismember_sorted(ind, chainind1);
        tx(tmpind, f+12) = 1;
        tmpind = ismember_sorted(ind, chainind2);
        tx(tmpind, f+12) = 2;
        nGeomJuncts(spLR(vve, 2)) = nGeomJuncts(spLR(vve, 2))+1;
    end
end
categoryFeatures(end+1) = f+12;

% add feature if segment has two or more geometric T-junctions
tmpind = find(tx(:, f+12)>0);
for k = tmpind'
    if tx(k, f+12)==1 && nGeomJuncts(spLR(ind(k), 1))>1
        tx(k, f+13) = 1;
    elseif tx(k, f+12)==2 && nGeomJuncts(spLR(ind(k), 2))>1
        tx(k, f+13) = 2;
    end
end
categoryFeatures(end+1) = f+13;

% Sky-Vertical-Vertical junctions (same as above but substitute sky for
% ground)
% In future make this faster by checking chains to see when s/v transitions
% to v/v
is_gv = (gvs1==3 & gvs2==2) | (gvs1==2 & gvs2==3);
is_vv = (gvs1==2 & gvs2==2);
ejuncts = bndinfo.edges.junctions;
juncts = cell(bndinfo.nj, 1);
jsize = zeros(size(juncts));
for k = 1:size(ejuncts, 1)
    j1 = ejuncts(k, 1);  j2 = ejuncts(k, 2);
    juncts{j1}(end+1) = k;
    juncts{j2}(end+1) = k;
    jsize(j1) = jsize(j1) + 1;
    jsize(j2) = jsize(j2) + 1;
end
juncts = cell2mat(juncts(jsize==3));
isGeomJunction = sum(is_gv(juncts), 2)==2 & sum(is_vv(juncts), 2)==1;
juncts = juncts(isGeomJunction, :);
nGeomJuncts = zeros(bndinfo.nseg, 1);
e2chain = X.edge.edge2chain;
for k = 1:size(juncts, 1) % check that g/v --> v/v angle is less than g/v -->
    gve = juncts(k, is_gv(juncts(k, :)));
    % if ground on left, then reverse
    vve = juncts(k, is_vv(juncts(k, :)));
    if maxlab(spLR(gve(1), 1))==3, gve(1) = gve(1)+ne;  end
    if maxlab(spLR(gve(2), 1))==3, gve(2) = gve(2)+ne;  end
    if any(e2chain(vve)==e2chain(gve)) && (e2chain(vve)>0) % foreground on left
        chainind = sort(X.edge.chains{e2chain(vve)});
        chainind1 = chainind(chainind<=ne);
        chainind2 = chainind(chainind>ne)-ne;
        tmpind = ismember_sorted(ind, chainind1);
        tx(tmpind, f+14) = 1;
        tmpind = ismember_sorted(ind, chainind2);
        tx(tmpind, f+14) = 2;
        nGeomJuncts(spLR(vve, 1)) = nGeomJuncts(spLR(vve, 1))+1;
    elseif any(e2chain(vve+ne)==e2chain(gve)) && (e2chain(vve+ne)>0) % foreground on right
        chainind = sort(X.edge.chains{e2chain(vve+ne)});
        chainind1 = chainind(chainind<=ne);
        chainind2 = chainind(chainind>ne)-ne;
        tmpind = ismember_sorted(ind, chainind1);
        tx(tmpind, f+14) = 1;
        tmpind = ismember_sorted(ind, chainind2);
        tx(tmpind, f+14) = 2;
        nGeomJuncts(spLR(vve, 2)) = nGeomJuncts(spLR(vve, 2))+1;
    end
end
categoryFeatures(end+1) = f+14;

% add feature if segment has two or more geometric T-junctions
tmpind = find(tx(:, f+14)>0);
for k = tmpind'
    if tx(k, f+14)==1 && nGeomJuncts(spLR(ind(k), 1))>1
        tx(k, f+15) = 1;
    elseif tx(k, f+14)==2 && nGeomJuncts(spLR(ind(k), 2))>1
        tx(k, f+15) = 2;
    end
end
categoryFeatures(end+1) = f+15;

%% Object features
if isfield(X, 'objectFeatures')
    tx = [tx X.objectFeatures(ind, :)];
end

author: Cloud-CV

repo: object-proposals