text:
%% DEMO_febio_0059_face_mask_loading
% Below is a demonstration for:
%
% * Building triangulated surface geometry for a face
% * Meshing the face using pentahedral elements
% * Building model of a tube
% * Defining the boundary conditions
% * Coding the febio structure
% * Running the model
% * Importing and visualizing results

%% Keywords
%
% * febio_spec version 3.0
% * febio, FEBio
% * face
% * contact, sliding, friction
% * pentahedral elements, penta6
% * static, solid
% * hyperelastic, Ogden
% * displacement logfile

%%

clear; close all; clc;

%%
% Plot settings
fontSize=15;
faceAlpha1=1;
faceAlpha2=0.3;
markerSize1=15;
markerSize2=10;
lineWidth=2;
cMap=spectral(250);

%% Control parameters

% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','temp');

% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=fullfile(savePath,[febioFebFileNamePart,'.txt']); %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement
febioLogFileName_force=[febioFebFileNamePart,'_force_out.txt']; %Log file name for exporting force
febioLogFileName_strainEnergy=[febioFebFileNamePart,'_energy_out.txt']; %Log file name for exporting strain energy density

% Geometry parameters
pointSpacingTissue=6;
faceTissueThickness=6;
pointSpacingMask=pointSpacingTissue/2;
maskRimWidth=5;
maskRimFilletRadius=6;
maskDiscRadius1=25;
maskDiscRadius2=maskDiscRadius1+4;
maskDiscOffset=25;
bezierTangency=0.1;

distInclude=40; %Distance from mask to include face in FEA

%Ray tracing parameters
optionStructRayTrace.tolEps        = 1e-6;
optionStructRayTrace.triSide       = 0;
optionStructRayTrace.rayType       = 'ray';
optionStructRayTrace.exclusionType = 'inclusive';
optionStructRayTrace.paired        = 0;

%Material parameters
c1_tissue=1e-3; %Shear-modulus-like parameter
m1_tissue=2; %Material parameter setting degree of non-linearity
k_tissue=c1_tissue*100; %Bulk modulus

c1_rim=c1_tissue*5; %Shear-modulus-like parameter
m1_rim=2; %Material parameter setting degree of non-linearity
k_rim=c1_rim*10; %Bulk modulus

% FEA control settings
numTimeSteps=15; %Number of time steps desired
max_refs=35; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=10; %Optimum number of iterations
max_retries=5; %Maximum number of retires
dtmin=(1/numTimeSteps)/100; %Minimum time step size
dtmax=(1/numTimeSteps); %Maximum time step size
symmetric_stiffness=0;
min_residual=1e-20;
runMode='external';

%Boundary condition parameters
initialOffset=0;
displacementMagnitude_z=-2-initialOffset; %Displacement applied

%Contact parameters
contactPenalty=10;
laugon=0;
minaug=1;
maxaug=10;
fric_coeff=0.5;

%% Load face geometry

testCase=1;
switch testCase
    case 1
        %Load surface model
        [Ff,Vf]=graphicsModels(9);

        %Surface markers
        V_markers=[65.51,49.94,217.14;... %Tip of the nose
            66.44,54.79,259.81;... %Nose between eyes
            53.81,115,263.39;... %Right eye outer corner
            126.5,46.62,269.75;... %Left eye outer corner
            85.67,69.44,194.1;... %Middle of mounth
            98.39,80.38,158]; %Bottom of chin
    case 2
        %Load surface model
        [Ff,Vf]=graphicsModels(13);

        %Surface markers
        V_markers=[3.50162,-181.107,-8.09110;...  %Tip of the nose
                   1.51627,-159.171,30.4679;... %Nose between eyes
                   -43.2473,-139.003,24.6407;... %Right eye outer corner
                   48.4245,-136.490,25.3913;... %Left eye outer corner
                   4.11905,-167.594,-36.9600;... %Middle of mounth
                   2.85290,-161.802,-73.4644]; %Bottom of chin
end

distEyes=sqrt(sum((V_markers(3,:)-V_markers(4,:)).^2,2));

%% Remeshing surface
optionStruct_remesh.pointSpacing=pointSpacingTissue; %Set desired point spacing
optionStruct_remesh.disp_on=0; % Turn off command window text display
[Ff,Vf]=ggremesh(Ff,Vf,optionStruct_remesh);

%%

ny=vecnormalize(V_markers(2,:)-V_markers(6,:));
nx=vecnormalize(V_markers(4,:)-V_markers(3,:));
nz=vecnormalize(cross(nx,ny));
nx=vecnormalize(cross(ny,nz));

Q=[nx;ny;nz]';

%%

% cFigure; hold on;
% gpatch(Ff,Vf,'w','none',1);
% % plotV(V_markers,'r.','MarkerSize',35);
% % text(V_markers(:,1)+4,V_markers(:,2),V_markers(:,3),{'1','2','3','4','5','6'},'FontSize',25);
% % quiverTriad(V_markers(1,:),Q,100);
% axisGeom; camlight headlight;
% colormap(spectral(250))
% gdrawnow;

%%

cFigure; hold on;
gpatch(Ff,Vf,'w','none',0.5);
plotV(V_markers,'r.','MarkerSize',35);
text(V_markers(:,1)+4,V_markers(:,2),V_markers(:,3),{'1','2','3','4','5','6'},'FontSize',25);
quiverTriad(V_markers(1,:),Q,100);
axisGeom; camlight headlight;
colormap(spectral(250))
gdrawnow;

%% Centre on nose and rotate to face face looking down Z-axis

Vf=Vf-V_markers(1,:);
Vf=Vf*Q;
V_markers=V_markers-V_markers(1,:);
V_markers=V_markers*Q;
nz=[0 0 1];%nz*Q;

%%

cFigure; hold on;
gpatch(Ff,Vf,'w','none',0.5);
plotV(V_markers,'r.','MarkerSize',35);
text(V_markers(:,1)+6,V_markers(:,2),V_markers(:,3)+15,{'1','2','3','4','5','6'},'FontSize',25);
axisGeom; camlight headlight;
colormap(spectral(250))
gdrawnow;

%% Construct mask rim curve

V1=V_markers(1,[1 2]);
V2=V_markers(2,[1 2]);
V3=V_markers(3,[1 2]);
V4=V_markers(4,[1 2]);
V5=V_markers(5,[1 2]);
V6=V_markers(6,[1 2]);

pp1=0.4*V1+0.6*V2;
pp2=0.6*V1+0.4*V3;
pp3=V3-[0 V3(2)]+[0 0.5*V5(2)+0.5*V1(2)];
pp4=[0.3*V2(1)+0.8*V3(1) 0.5*V5(2)+0.5*V6(2)];
pp5=V6;
pp6=[0.3*V2(1)+0.8*V4(1) 0.5*V5(2)+0.5*V6(2)];
pp7=V4-[0 V4(2)]+[0 0.5*V5(2)+0.5*V1(2)];
pp8=0.6*V1+0.4*V4;

V_rim_points=[pp1;pp2;pp3;pp4;pp5;pp6;pp7;pp8];
V_rim_points(:,3)=0;

[V_rim_points,indFaceIntersect]=traceToSurf(V_rim_points,-nz,Ff,Vf,optionStructRayTrace);
numRimControlPoints=size(V_rim_points,1);

V_rim_curve=evenlySpaceCurve(V_rim_points,pointSpacingMask,'pchip',1);
V_rim_curve=traceToSurf(V_rim_curve,-nz,Ff,Vf,optionStructRayTrace);
numPointsRimCurve=size(V_rim_curve,1);

Ne1=vecnormalize([V_rim_points(2:end,:); V_rim_points(1,:)]-V_rim_points(1:end,:));
Ne2=vecnormalize(V_rim_points - [V_rim_points(end,:); V_rim_points(1:end-1,:)]);
Ne=vecnormalize(0.5*Ne1+0.5*Ne2);

Nf=patchNormal(Ff,Vf); %Normal directions
Nff=Nf(indFaceIntersect(:,2),:);
N_rim_points=vecnormalize(cross(Nff,Ne));

V_rim_points1=V_rim_points-N_rim_points.*maskRimWidth/2;
V_rim_points2=V_rim_points+N_rim_points.*maskRimWidth/2;

%%

cFigure; hold on;
gpatch(Ff,Vf,'w','none',0.5);
plotV(V_markers,'r.','MarkerSize',25);
text(V_markers(:,1)+4,V_markers(:,2),V_markers(:,3),{'1','2','3','4','5','6'},'FontSize',25);

plotV(V_rim_points,'k.','MarkerSize',25);
plotV(V_rim_points1,'b.','MarkerSize',25);
plotV(V_rim_points2,'g.','MarkerSize',25);

quiverVec(V_rim_points,N_rim_points,maskRimWidth/2,'y');
quiverVec(V_rim_points,-N_rim_points,maskRimWidth/2,'y');
axisGeom; camlight headlight;
view(2);
gdrawnow;

%%

[~,indClose]=minDist(V_markers,Vf);
d=meshDistMarch(Ff,Vf,indClose([1 5]));

[~,indClose]=minDist(V_rim_curve,Vf);
d_rim_curve=meshDistMarch(Ff,Vf,indClose);
d_markers_max=max(d(indClose));

logicCloseVertices= d<=d_markers_max | d_rim_curve<distInclude;

logicCloseFaces= any(logicCloseVertices(Ff),2);
logicCloseFaces=triSurfLogicSharpFix(Ff,logicCloseFaces);

[Fs,Vs]=patchCleanUnused(Ff(logicCloseFaces,:),Vf);

ns=vecnormalize(mean(patchNormal(Fs,Vs)));

[Q]=pointSetPrincipalDir(Vs);
nz=Q(:,3)';
if dot(nz,ns)<1
    nz=-nz;
end

ny=vecnormalize(V_markers(2,:)-V_markers(1,:));
nx=cross(ny,nz);
ny=cross(nz,nx);
Q=[nx;ny;nz]';

Ffc=Ff(~logicCloseFaces,:);

%%

cFigure;  hold on;
gpatch(Ff,Vf,d,'k',0.5);
gpatch(Ff(logicCloseFaces,:),Vf,'none','b',1,2);
plotV(V_markers,'r.','MarkerSize',25);
plotV(V_rim_points,'k.','MarkerSize',15);
plotV(V_rim_curve,'k-','LineWidth',3);

axisGeom; camlight headlight;
colormap(spectral(250))
gdrawnow;

%%

cFigure;  hold on;
gpatch(Ff,Vf,'w','none',0.5);
% gpatch(Ff(logicCloseFaces,:),Vf,'none','b',1,2);
plotV(V_markers,'r.','MarkerSize',25);
plotV(V_rim_points,'k.','MarkerSize',15);
plotV(V_rim_curve,'k-','LineWidth',3);

axisGeom; camlight headlight;
colormap(spectral(250))
gdrawnow;

%%

cFigure;  hold on;
gpatch(Ff,Vf,'w','none');
gpatch(Fs,Vs,'w','b');
% patchNormPlot(Fs,Vs);
plotV(V_markers,'r.','MarkerSize',25);
% quiverTriad(V_markers(1,:),Q,50);
axisGeom;
camlight headlight;
% colormap(viridis(2)); icolorbar;
gdrawnow;

%%

Ebs=patchBoundary(Fs);
indBoundaryCurve=edgeListToCurve(Ebs);
indBoundaryCurve=indBoundaryCurve(1:end-1)';

[~,~,Ns]=patchNormal(Fs,Vs);

Vs2=Vs-Ns*faceTissueThickness;
Fs2=Fs;

Vsc=Vs2(indBoundaryCurve,:);%-Ns(indBoundaryCurve,:)*layerThickness;
[Fs2t,Vs2t]=regionTriMesh3D({Vsc},pointSpacingMask,0,'natural');
Vs2t_ori=Vs2t;
indBoundary=unique(patchBoundary(Fs2t));

[~,indMap]=minDist(Vsc,Vs2t(indBoundary,:));
indBoundaryCurve_2t=indBoundary(indMap);
indBoundaryCurve_2t=indBoundaryCurve_2t(:);

Vs2t=traceToSurf(Vs2t,-nz,Fs2,Vs2,optionStructRayTrace);

cParSmooth.n=3;
cParSmooth.Method='HC';
cParSmooth.RigidConstraints=indBoundaryCurve_2t;
Vs2t=patchSmooth(Fs2t,Vs2t,[],cParSmooth);

numNodesThickness=ceil(faceTissueThickness./pointSpacingTissue);
if numNodesThickness<2
    numNodesThickness=2;
end

cParLoft.numSteps=numNodesThickness;
cParLoft.closeLoopOpt=1;
cParLoft.patchType='tri';
[Fss,Vss,ind1,ind2]=polyLoftLinear(Vs(indBoundaryCurve,:),Vs2t(indBoundaryCurve_2t,:),cParLoft);

[Fb,Vb,Cb]=joinElementSets({Fs,Fs2t,Fss},{Vs,Vs2t,Vss});
[Fb,Vb]=mergeVertices(Fb,Vb);

%%

cFigure; hold on;
gpatch(Ffc,Vf,'w','none',0.5);
gpatch(Fb,Vb,Cb,'k',0.5);

% gpatch(Fs,Vs,'bw','k',1);
% % gpatch(Fs2t,Vs2t_ori,'gw','g',0.5);
% gpatch(Fs2t,Vs2t,'rw','k',1);
% gpatch(Fss,Vss,'gw','k',1);

% patchNormPlot(Fb,Vb);
axisGeom; camlight headlight;
colormap(spectral); icolorbar;
gdrawnow;

%%
% Get inner mesh point between top and bottom at nose
Pn=triSurfRayTrace(V_markers(1,:),-nz,Fb,Vb,optionStructRayTrace);
Pn=mean(Pn,1);

%%

cFigure; hold on;
%gpatch(Ffc,Vf,'w','none',0.5);
gpatch(Fb,Vb,Cb,'k',1);
plotV(Pn,'r.','MarkerSize',25);

% patchNormPlot(Fb,Vb);
axisGeom; camlight headlight;
colormap(spectral); icolorbar;
gdrawnow;

%%

%Create tetgen input structure
inputStruct.stringOpt='-pq1.2AaY'; %Options for tetgen
inputStruct.Faces=Fb; %Boundary faces
inputStruct.Nodes=Vb; %Nodes of boundary
inputStruct.faceBoundaryMarker=Cb;
inputStruct.regionPoints=Pn; %Interior points for regions
inputStruct.holePoints=[]; %Interior points for holes
inputStruct.regionA=tetVolMeanEst(Fb,Vb); %Desired tetrahedral volume for each region

% Mesh model using tetrahedral elements using tetGen
[meshOutput]=runTetGen(inputStruct); %Run tetGen

%%
% Access mesh output structure

E_face=meshOutput.elements; %The elements
V=meshOutput.nodes; %The vertices or nodes
F=meshOutput.faces; %Element faces (all)
CE=meshOutput.elementMaterialID; %Element material or region id
Fb=meshOutput.facesBoundary; %The boundary faces
Cb=meshOutput.boundaryMarker; %The boundary markers

%%
% Visualization

hf=cFigure;
subplot(1,2,1); hold on;
title('Input boundaries','FontSize',fontSize);
hp(1)=gpatch(Fb,V,Cb,'k',faceAlpha1);
hp(2)=plotV(Pn,'r.','MarkerSize',markerSize1);
legend(hp,{'Input mesh','Interior point(s)'},'Location','NorthWestOutside');
axisGeom(gca,fontSize); camlight headlight;
colormap(cMap); icolorbar;

hs=subplot(1,2,2); hold on;
title('Tetrahedral mesh','FontSize',fontSize);

% Visualizing using |meshView|
optionStruct.hFig=[hf,hs];

meshView(meshOutput,optionStruct);
hold on; plotV(Pn,'r.','MarkerSize',25);
axisGeom(gca,fontSize);
gdrawnow;

%%

pp9=V_markers(5,:);% (0.1*V_markers(1,:)+0.9*V_markers(5,:));

t=linspace(0.5*pi,2.5*pi,numPointsRimCurve+1)';
t=t(1:end-1);
Vcd1=maskDiscRadius1*[cos(t) sin(t) zeros(size(t))];
Vcd1=Vcd1+pp9;
Vcd1(:,3)=maskDiscOffset;

Vcd2=maskDiscRadius2*[cos(t) sin(t) zeros(size(t))];
Vcd2=Vcd2+pp9;
Vcd2(:,3)=maskDiscOffset;

%
cFigure; hold on;
gpatch(Fb,V,'w','none',0.5);
plotV(V_markers,'r.','MarkerSize',50);
text(V_markers(:,1)+4,V_markers(:,2),V_markers(:,3)+10,{'1','2','3','4','5','6'},'FontSize',25);

plotV(V_rim_points,'k.','MarkerSize',35);
plotV(V_rim_curve,'k-','LineWidth',6);
plotV(Vcd1,'b-','LineWidth',6);
plotV(Vcd2,'b-','LineWidth',6);
axisGeom;
camlight headlight;
view(2);
drawnow;

%%

V_rim_curve1=evenlySampleCurve(V_rim_points1,numPointsRimCurve,'pchip',1);
V_rim_curve2=evenlySampleCurve(V_rim_points2,numPointsRimCurve,'pchip',1);

V_rim_curve1=traceToSurf(V_rim_curve1,[0 0 -1],Fb(Cb==1,:),V,optionStructRayTrace);
V_rim_curve2=traceToSurf(V_rim_curve2,[0 0 -1],Fb(Cb==1,:),V,optionStructRayTrace);

%%
cFigure; hold on;
gpatch(Fb(Cb==1,:),Vb,'kw','none',0.5);

% patchNormPlot(Fm,Vm);
plotV(V_markers,'r.','MarkerSize',50);
text(V_markers(:,1)+2,V_markers(:,2),V_markers(:,3),{'1','2','3','4','5','6'},'FontSize',25);
plotV(V_rim_points,'k.','MarkerSize',35,'LineWidth',3);
plotV(V_rim_curve,'k-','LineWidth',3);

plotV(V_rim_points1,'g.','MarkerSize',35,'LineWidth',3);
plotV(V_rim_curve1,'g-','LineWidth',3);

plotV(V_rim_points2,'b.','MarkerSize',35,'LineWidth',3);
plotV(V_rim_curve2,'b-','LineWidth',3);

axisGeom; camlight headlight;
view(2);
drawnow;

%%

pointSpacingNow=mean(diff(pathLength(V_rim_curve1)));
numNodStrip=ceil(maskRimWidth./pointSpacingNow);
if numNodStrip<2
    numNodStrip=2;
end
if iseven(numNodStrip)
    numNodStrip=numNodStrip+1;
end

cParLoft.numSteps=numNodStrip;
cParLoft.closeLoopOpt=1;
cParLoft.patchType='tri';
[Fm,Vm,indStart_Vm,indEnd_Vm]=polyLoftLinear(V_rim_curve1,V_rim_curve2,cParLoft);
Fm=fliplr(Fm);
indStart_Vm=fliplr(indStart_Vm);
indEnd_Vm=fliplr(indEnd_Vm);

[~,~,Nm]=patchNormal(Fm,Vm);
Vm=traceToSurf(Vm,Nm,Fb(Cb==1,:),V,optionStructRayTrace);

%%
cFigure; hold on;
gpatch(Fb(Cb==1,:),Vb,'kw','none',0.5);
gpatch(Fm,Vm,'rw','k',1,1);
% patchNormPlot(Fm,Vm);
% plotV(V_markers,'r.','MarkerSize',50);
% text(V_markers(:,1)+4,V_markers(:,2),V_markers(:,3),{'1','2','3','4','5','6'},'FontSize',25);
plotV(V_rim_points,'k.','MarkerSize',35,'LineWidth',3);
plotV(Vm(indStart_Vm,:),'g-','LineWidth',3);
plotV(Vm(indEnd_Vm,:),'b-','LineWidth',3);
axisGeom; camlight headlight;
view(2);
drawnow;

%%

[~,~,Nm]=patchNormal(Fm,Vm);

pointSpacingNow=mean(diff(pathLength(Vm(indStart_Vm,:))));
nRim=ceil((pi/2*maskRimFilletRadius)/pointSpacingNow)+1;
if nRim<4
    nRim=4;
end

[Fr1,Vr1]=roundMesh(indStart_Vm,Vm,Nm,nRim,maskRimFilletRadius);
[Fr2,Vr2]=roundMesh(indEnd_Vm,Vm,Nm,nRim,maskRimFilletRadius);
indEnd_Vr1=size(Vr1)-numPointsRimCurve+1:1:size(Vr1);
indEnd_Vr2=fliplr(indEnd_Vr1);
[Fr1,Vr1]=quad2tri(Fr1,Vr1,'a');
[Fr2,Vr2]=quad2tri(Fr2,Vr2,'a');

%%

cFigure; hold on;
gpatch(Fb,V,'kw','none',0.5);
gpatch(Fm,Vm,'rw','k',1,1);
gpatch(Fr1,Vr1,'gw','k',1,1);
gpatch(Fr2,Vr2,'bw','k',1,1);
plotV(Vr1(indEnd_Vr1,:),'r-','LineWidth',3);
plotV(Vr2(indEnd_Vr2,:),'b-','LineWidth',3);
axisGeom; camlight headlight;
view(2);
drawnow;

%%

pointSpacingNow=mean(diff(pathLength(Vr1(indEnd_Vr1,:))));
numNodStripTop=ceil((maskRimWidth+2*maskRimFilletRadius)./pointSpacingNow);
if numNodStripTop<2
    numNodStripTop=2;
end
if iseven(numNodStripTop)
    numNodStripTop=numNodStripTop+1;
end

cParLoft.numSteps=numNodStripTop;
cParLoft.closeLoopOpt=1;
cParLoft.patchType='tri';
[Ft,Vt,indCurve1_Vt,indCurve2_Vt]=polyLoftLinear(Vr1(indEnd_Vr1,:),Vr2(indEnd_Vr2,:),cParLoft);
Ft=fliplr(Ft);
indCurve1_Vt=flipud(indCurve1_Vt(:));
% indCurve2_Vt=flipud(indCurve2_Vt(:));
[~,~,Nt]=patchNormal(Ft,Vt);

%%

cFigure; hold on;
gpatch(Fb,V,'w','none',1);
gpatch(Fm,Vm,'y','k',0,1);
gpatch(Ft,Vt,'gw','k',1,1);
gpatch(Fr1,Vr1,'rw','k',1,1);
gpatch(Fr2,Vr2,'bw','k',1,1);
plotV(Vr1(indEnd_Vr1,:),'r-','LineWidth',3);
plotV(Vr2(indEnd_Vr1,:),'b-','LineWidth',3);
axisGeom; camlight headlight;
view(2);
drawnow;

%% Join and merge rim surfaces

[F_rim,V_rim,C_rim]=joinElementSets({Fm,Fr1,Fr2,Ft},{Vm,Vr1,Vr2,Vt});
[F_rim,V_rim]=mergeVertices(F_rim,V_rim);
[F_rim,V_rim]=patchCleanUnused(F_rim,V_rim);

%%

cFigure; hold on;
gpatch(Fb(Cb==1,:),V,'w','none',1);
gpatch(F_rim,V_rim,C_rim,'none',1);
% patchNormPlot(F_rim,V_rim);
axisGeom; camlight headlight;
colormap spectral; icolorbar;
drawnow;

%%

V_loft1=Vt(indCurve1_Vt,:);
V_loft2=Vcd2;

V_loft3=Vt(indCurve2_Vt,:);
V_loft4=Vcd1;

[~,indMax]=max(V_loft4(:,1));
if indMax<numPointsRimCurve/2
    V_loft4=flipud(V_loft4);
end

[~,indMax]=max(V_loft2(:,1));
if indMax<numPointsRimCurve/2
    V_loft2=flipud(V_loft2);
end

N1=Nt(indCurve1_Vt,:);
N3=Nt(indCurve2_Vt,:);
N4=ones(size(V_loft4,1),1)*[0 0 1];

N2e=vecnormalize([V_loft2(2:end,:); V_loft2(1,:)]-V_loft2(1:end,:));
N2=vecnormalize(cross(N4,N2e));

[Fd1,Vd1,X,Y,Z]=bezierLoft(V_loft1,V_loft2,N1,N2,pointSpacingMask,bezierTangency);
[Fd1,Vd1]=quad2tri(Fd1,Vd1);
[Fd2,Vd2,X,Y,Z]=bezierLoft(V_loft3,V_loft4,N3,N4,pointSpacingMask,bezierTangency);
[Fd2,Vd2]=quad2tri(Fd2,Vd2);

pointSpacingNow=mean(diff(pathLength(V_loft2)));
n=ceil((maskDiscRadius2-maskDiscRadius1)./pointSpacingNow);
if n<2
    n=2;
end
if iseven(n)
    n=n+1;
end

cParLoft.numSteps=n;
cParLoft.closeLoopOpt=1;
cParLoft.patchType='tri';
[Fdt,Vdt]=polyLoftLinear(V_loft2,V_loft4,cParLoft);

[Fc,Vc]=regionTriMesh2D({V_loft4(:,[1 2])},pointSpacingNow,0,0);
Vc(:,3)=maskDiscOffset;

%%

cFigure; hold on;
gpatch(Fb(Cb==1,:),V,'kw','none',0.5);
gpatch(F_rim,V_rim,'kw','none',1);

gpatch(Fd1,Vd1,'rw','none',0.5);
gpatch(Fd2,Vd2,'bw','none',0.5);
gpatch(Fdt,Vdt,'gw','none',0.5);
gpatch(Fc,Vc,'yw','none',0.5);

plotV(V_loft1,'r-','LineWidth',3);
quiverVec(V_loft1,N1,5,'k');

plotV(V_loft4,'b-','LineWidth',2);
quiverVec(V_loft4,N4,5,'k');

plotV(V_loft3,'b-','LineWidth',3);
quiverVec(V_loft3,N3,5,'k');

plotV(V_loft2,'r-','LineWidth',2);

quiverVec(V_loft2,N2,5,'k');

axisGeom; camlight headlight;
drawnow;

%% Join and merge mask body components

[F_mask,V_mask,C_mask]=joinElementSets({Ft,Fd1,Fd2,Fdt,Fc},{Vt,Vd1,Vd2,Vdt,Vc});
[F_mask,V_mask]=mergeVertices(F_mask,V_mask);
[F_mask,V_mask]=patchCleanUnused(F_mask,V_mask);

%%
cFigure; hold on;
gpatch(Ff,Vf,'w','none',1);
% gpatch(Fb,V,'w','none',1);
gpatch(F_mask,V_mask,'bw','none',1);
gpatch(F_rim,V_rim,'kw','none',1);
axisGeom; camlight headlight;
view(2);
% colormap spectral; icolorbar;
drawnow;

%%
cFigure; hold on;
gpatch(Ff,Vf,'w','none',1);
% gpatch(Fb,V,'w','none',1);
% gpatch(F_mask,V_mask,'gw','none',0.5);
gpatch(F_rim,V_rim,'kw','none',1);
axisGeom; camlight headlight;
drawnow;

%%

cFigure; hold on;
% gpatch(Fp1,V,'w','none',1);
gpatch(F_mask,V_mask,'bw','none',1);
gpatch(F_rim,V_rim,'gw','none',1);
axisGeom;
camlight headlight;
view(2);
drawnow;

%%

cFigure; hold on;
% gpatch(Fb,V,'w','none',0.9);
gpatch(F_rim,V_rim,C_rim,'k',1);
% patchNormPlot(F_rim,V_rim);
axisGeom;
camlight headlight;
view(2);
colormap gjet; icolorbar;
drawnow;

%%

[V_regions]=getInnerPoint(F_rim,V_rim); % Define region points

%%

cFigure; hold on;
gpatch(F_rim,V_rim,C_rim,'k',1);
% plotV(V_regions,'r.','MarkerSize',markerSize1)
axisGeom;
camlight headlight; view(2);
drawnow;

%%

[regionA]=tetVolMeanEst(F_rim,V_rim); %Volume for regular tets

inputStruct.stringOpt='-pq1.2AaY';
inputStruct.Faces=F_rim;
inputStruct.Nodes=V_rim;
inputStruct.holePoints=[];
inputStruct.faceBoundaryMarker=C_rim; %Face boundary markers
inputStruct.regionPoints=V_regions; %region points
inputStruct.regionA=regionA;

% Mesh model using tetrahedral elements using tetGen
[meshOutput]=runTetGen(inputStruct); %Run tetGen

%%
% Access model element and patch data
Fb_rim=meshOutput.facesBoundary;
Cb_rim=meshOutput.boundaryMarker;
V_rim=meshOutput.nodes;
E_rim=meshOutput.elements;

% Visualizing mesh using |meshView|, see also |anim8|
meshView(meshOutput);

%% Join node sets
V_rim(:,3)=V_rim(:,3)+initialOffset;
Fb_rim=Fb_rim+size(V,1);
E_rim=E_rim+size(V,1);
V=[V;V_rim];

%%
cFigure;
gpatch(F,V,'w','k',1);
gpatch(Fb_rim,V,Cb_rim,'k',1);
axisGeom;
colormap gjet; icolorbar;
camlight headlight;
drawnow;

%% Define contact surfaces

% The rigid primary surface of the sphere
F_contact_primary=fliplr(Fb_rim(Cb_rim~=4,:));

% The deformable secondary surface of the slab
Fb_contact=Fb(Cb==1,:);
V_Fb_centre=patchCentre(Fb_contact,V);
D=minDist(V_Fb_centre,V_rim);
logicSecondary=D<=(2*pointSpacingTissue);
logicSecondary=triSurfLogicSharpFix(Fb_contact,logicSecondary,3);
F_contact_secondary=fliplr(Fb_contact(logicSecondary,:));

%%
% Visualize contact surfaces

cFigure; hold on;
title('Contact sets and normal directions','FontSize',fontSize);

gpatch(Fb,V,'kw','none',0.5);

hl(1)=gpatch(F_contact_primary,V,'gw','k',1);
patchNormPlot(F_contact_primary,V);
hl(2)=gpatch(F_contact_secondary,V,'rw','k',1);
patchNormPlot(F_contact_secondary,V);

legend(hl,{'Primary','secondary'});

axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Define boundary conditions

%Supported nodes
bcSupportList=unique(Fb(Cb==2,:));

%Prescribed displacement nodes
bcPrescribeList=unique(Fb_rim(Cb_rim==4,:));

%%
% Visualize BC's

hf=cFigure; hold on;
% title('Boundary conditions model','FontSize',fontSize);
gpatch(Fb,V,'kw','none',faceAlpha2);
gpatch(Fb_rim(Cb_rim~=4,:),V,'kw','none',faceAlpha2);
gpatch(Fb_rim(Cb_rim==4,:),V,'rw','none',1);
gpatch(Fb(Cb==2,:),V,'kw','none',1);
hl2(1)=plotV(V(bcPrescribeList,:),'r.','MarkerSize',markerSize2);
hl2(2)=plotV(V(bcSupportList,:),'k.','MarkerSize',markerSize2);
legend(hl2,{'BC prescribe','BC support'});
axisGeom(gca,fontSize); camlight headlight;
drawnow;

%% Defining the FEBio input structure
% See also |febioStructTemplate| and |febioStruct2xml| and the FEBio user
% manual.

%Get a template with default settings
[febio_spec]=febioStructTemplate;

%febio_spec version
febio_spec.ATTR.version='3.0';

%Module section
febio_spec.Module.ATTR.type='solid';

%Control section
febio_spec.Control.analysis='STATIC';
febio_spec.Control.time_steps=numTimeSteps;
febio_spec.Control.step_size=1/numTimeSteps;
febio_spec.Control.solver.max_refs=max_refs;
febio_spec.Control.solver.max_ups=max_ups;
febio_spec.Control.solver.symmetric_stiffness=symmetric_stiffness;
febio_spec.Control.time_stepper.dtmin=dtmin;
febio_spec.Control.time_stepper.dtmax=dtmax;
febio_spec.Control.time_stepper.max_retries=max_retries;
febio_spec.Control.time_stepper.opt_iter=opt_iter;

%Material section
materialName1='Material1';
febio_spec.Material.material{1}.ATTR.name=materialName1;
febio_spec.Material.material{1}.ATTR.type='Ogden unconstrained';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.c1=c1_tissue;
febio_spec.Material.material{1}.m1=m1_tissue;
febio_spec.Material.material{1}.c2=c1_tissue;
febio_spec.Material.material{1}.m2=-m1_tissue;
febio_spec.Material.material{1}.cp=k_tissue;

materialName2='Material2';
febio_spec.Material.material{2}.ATTR.name=materialName2;
febio_spec.Material.material{2}.ATTR.type='Ogden unconstrained';
febio_spec.Material.material{2}.ATTR.id=2;
febio_spec.Material.material{2}.c1=c1_rim;
febio_spec.Material.material{2}.m1=m1_rim;
febio_spec.Material.material{2}.c2=c1_rim;
febio_spec.Material.material{2}.m2=-m1_rim;
febio_spec.Material.material{2}.cp=k_rim;

% Mesh section
% -> Nodes
febio_spec.Mesh.Nodes{1}.ATTR.name='All'; %The node set name
febio_spec.Mesh.Nodes{1}.node.ATTR.id=(1:size(V,1))'; %The node id's
febio_spec.Mesh.Nodes{1}.node.VAL=V; %The nodel coordinates

% -> Elements
partName1='Part1';
febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of this part
febio_spec.Mesh.Elements{1}.ATTR.type='tet4'; %Element type
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E_face,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E_face; %The element matrix

partName2='Part2';
febio_spec.Mesh.Elements{2}.ATTR.name=partName2; %Name of this part
febio_spec.Mesh.Elements{2}.ATTR.type='tet4'; %Element type
febio_spec.Mesh.Elements{2}.elem.ATTR.id=size(E_face,1)+(1:1:size(E_rim,1))'; %Element id's
febio_spec.Mesh.Elements{2}.elem.VAL=E_rim; %The element matrix

% -> NodeSets
nodeSetName1='bcSupportList';
febio_spec.Mesh.NodeSet{1}.ATTR.name=nodeSetName1;
febio_spec.Mesh.NodeSet{1}.node.ATTR.id=bcSupportList(:);

nodeSetName2='bcPrescribeList';
febio_spec.Mesh.NodeSet{2}.ATTR.name=nodeSetName2;
febio_spec.Mesh.NodeSet{2}.node.ATTR.id=bcPrescribeList(:);

%MeshDomains section
febio_spec.MeshDomains.SolidDomain{1}.ATTR.name=partName1;
febio_spec.MeshDomains.SolidDomain{1}.ATTR.mat=materialName1;

febio_spec.MeshDomains.SolidDomain{2}.ATTR.name=partName2;
febio_spec.MeshDomains.SolidDomain{2}.ATTR.mat=materialName2;

% -> Surfaces
surfaceName1='contactSurface1';
febio_spec.Mesh.Surface{1}.ATTR.name=surfaceName1;
febio_spec.Mesh.Surface{1}.tri3.ATTR.id=(1:1:size(F_contact_primary,1))';
febio_spec.Mesh.Surface{1}.tri3.VAL=F_contact_primary;

surfaceName2='contactSurface2';
febio_spec.Mesh.Surface{2}.ATTR.name=surfaceName2;
febio_spec.Mesh.Surface{2}.tri3.ATTR.id=(1:1:size(F_contact_secondary,1))';
febio_spec.Mesh.Surface{2}.tri3.VAL=F_contact_secondary;

% -> Surface pairs
contactPairName='Contact1';
febio_spec.Mesh.SurfacePair{1}.ATTR.name=contactPairName;
febio_spec.Mesh.SurfacePair{1}.primary=surfaceName1;
febio_spec.Mesh.SurfacePair{1}.secondary=surfaceName2;

%Boundary condition section
% -> Fix boundary conditions
febio_spec.Boundary.bc{1}.ATTR.type='fix';
febio_spec.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{1}.dofs='x,y,z';

febio_spec.Boundary.bc{2}.ATTR.type='fix';
febio_spec.Boundary.bc{2}.ATTR.node_set=nodeSetName2;
febio_spec.Boundary.bc{2}.dofs='x,y';

febio_spec.Boundary.bc{3}.ATTR.type='prescribe';
febio_spec.Boundary.bc{3}.ATTR.node_set=nodeSetName2;
febio_spec.Boundary.bc{3}.dof='z';
febio_spec.Boundary.bc{3}.scale.ATTR.lc=1;
febio_spec.Boundary.bc{3}.scale.VAL=displacementMagnitude_z;
febio_spec.Boundary.bc{3}.relative=0;

%Contact section
febio_spec.Contact.contact{1}.ATTR.type='sliding-elastic';
febio_spec.Contact.contact{1}.ATTR.surface_pair=contactPairName;
febio_spec.Contact.contact{1}.two_pass=0;
febio_spec.Contact.contact{1}.laugon=laugon;
febio_spec.Contact.contact{1}.tolerance=0.2;
febio_spec.Contact.contact{1}.gaptol=0;
febio_spec.Contact.contact{1}.minaug=minaug;
febio_spec.Contact.contact{1}.maxaug=maxaug;
febio_spec.Contact.contact{1}.search_tol=0.01;
febio_spec.Contact.contact{1}.search_radius=0.1*sqrt(sum((max(V,[],1)-min(V,[],1)).^2,2));
febio_spec.Contact.contact{1}.symmetric_stiffness=0;
febio_spec.Contact.contact{1}.auto_penalty=1;
febio_spec.Contact.contact{1}.penalty=contactPenalty;
febio_spec.Contact.contact{1}.fric_coeff=fric_coeff;

%LoadData section
% -> load_controller
febio_spec.LoadData.load_controller{1}.ATTR.id=1;
febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{1}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{1}.points.point.VAL=[0 0; 1 1];

%Output section
% -> log file
febio_spec.Output.logfile.ATTR.file=febioLogFileName;
febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.node_data{1}.VAL=1:size(V,1);

febio_spec.Output.logfile.node_data{2}.ATTR.file=febioLogFileName_force;
febio_spec.Output.logfile.node_data{2}.ATTR.data='Rx;Ry;Rz';
febio_spec.Output.logfile.node_data{2}.ATTR.delim=',';
febio_spec.Output.logfile.node_data{2}.VAL=1:size(V,1);

febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_strainEnergy;
febio_spec.Output.logfile.element_data{1}.ATTR.data='sed';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.element_data{1}.VAL=1:size(E_face,1);

%% Quick viewing of the FEBio input file structure
% The |febView| function can be used to view the xml structure in a MATLAB
% figure window.

%%
% |febView(febio_spec); %Viewing the febio file|

%% Exporting the FEBio input file
% Exporting the febio_spec structure to an FEBio input file is done using
% the |febioStruct2xml| function.

febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode

%% Running the FEBio analysis
% To run the analysis defined by the created FEBio input file the
% |runMonitorFEBio| function is used. The input for this function is a
% structure defining job settings e.g. the FEBio input file name. The
% optional output runFlag informs the user if the analysis was run
% succesfully.

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.disp_log_on=1; %Display convergence information in the command window
febioAnalysis.runMode=runMode;%'internal';
febioAnalysis.t_check=0.25; %Time for checking log file (dont set too small)
febioAnalysis.maxtpi=1e99; %Max analysis time
febioAnalysis.maxLogCheckTime=10; %Max log file checking time

[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!

%% Import FEBio results

if runFlag==1 %i.e. a succesful run

    % Importing nodal displacements from a log file
    [time_mat, N_disp_mat,~]=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp)); %Nodal displacements
    time_mat=[0; time_mat(:)]; %Time

    N_disp_mat=N_disp_mat(:,2:end,:);
    sizImport=size(N_disp_mat);
    sizImport(3)=sizImport(3)+1;
    N_disp_mat_n=zeros(sizImport);
    N_disp_mat_n(:,:,2:end)=N_disp_mat;
    N_disp_mat=N_disp_mat_n;
    DN=N_disp_mat(:,:,end);

    V_def=V+DN;
    V_DEF=N_disp_mat+repmat(V,[1 1 size(N_disp_mat,3)]);
    X_DEF=V_DEF(:,1,:);
    Y_DEF=V_DEF(:,2,:);
    Z_DEF=V_DEF(:,3,:);

    C=sqrt(sum(DN(:,3).^2,2));

    %%
    % Importing element strain energies from a log file
    [~,E_energy,~]=importFEBio_logfile(fullfile(savePath,febioLogFileName_strainEnergy)); %Element strain energy

    %Remove nodal index column
    E_energy=E_energy(:,2:end,:);

    %Add initial state i.e. zero energy
    sizImport=size(E_energy);
    sizImport(3)=sizImport(3)+1;
    E_energy_mat_n=zeros(sizImport);
    E_energy_mat_n(:,:,2:end)=E_energy;
    E_energy=E_energy_mat_n;

    [FE_face,C_energy_face]=element2patch(E_face,E_energy(1:size(E_face,1),:,end),'tet4');
    [CV]=faceToVertexMeasure(FE_face,V,C_energy_face);

    %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations
    cMap_c=gjet(250);
    cMap=[linspacen([1 1 1],cMap_c(1,:),50)'; cMap_c];

%     cMap=cMap(4:end,:);

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);

    gpatch(Ffc,Vf,cMap(1,:),'none',1)
    hp1=gpatch(Fb(Cb==1,:),V_def,CV,'none',1); %Add graphics object to animate
    hp1.FaceColor='Interp';
    hp2=gpatch(Fb_rim,V_def,'w','none',0.25); %Add graphics object to animate
    hp3=gpatch(F_mask,V_mask,'w','none',0.25);

    axisGeom(gca,fontSize); camlight headlight;
    colormap(cMap); colorbar;
    caxis([0 max(C_energy_face)/10]);
    axis([min(X_DEF(:)) max(X_DEF(:)) min(Y_DEF(:)) max(Y_DEF(:)) min(Z_DEF(:)) max(Z_DEF(:))]);
    axis tight;

    % Set up animation features
    animStruct.Time=time_mat; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments
        DN=N_disp_mat(:,:,qt); %Current displacement
        V_def=V+DN; %Current nodal coordinates

        %         C=sqrt(sum(DN(:,3).^2,2)); %New color
        [FE_face,C_energy_face]=element2patch(E_face,E_energy(1:size(E_face,1),:,qt),'tet4');
        [CV]=faceToVertexMeasure(FE_face,V,C_energy_face);

        u=mean(DN(bcPrescribeList,:),1);
        V_mask_def=V_mask+u(ones(size(V_mask,1),1),:);

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp1 hp1 hp2 hp3]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData','Vertices','Vertices'}; %Properties of objects to animate
        animStruct.Set{qt}={V_def,CV,V_def,V_mask_def}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    drawnow;

end

%%
function [varargout]=traceToSurf(V1,N1,F2,V2,optionStructRayTrace)

if size(N1,1)==1
    N1=N1(ones(size(V1,1),1),:);
end

numPoints=size(V1,1);
indFacesIntersect=nan(size(V1,1),2);
for q=1:1:numPoints
    [P,indFaceIntersect,~,~]=triSurfRayTrace(V1(q,:),N1(q,:),F2,V2,optionStructRayTrace);
    if size(P,1)>1
        [~,indMin]=minDist(V1(q,:),P);
        %         [~,indMin]=min(d);
        P=P(indMin,:);
        indFaceIntersect=indFaceIntersect(indMin,:);
    end
    if ~isempty(P)
        V1(q,:)=P;
        indFacesIntersect(q,:)=indFaceIntersect;
    end
end

varargout{1}=V1;
varargout{2}=indFacesIntersect;

end

function [Fr,Vr]=roundMesh(indCurve,Vm,Nm,nc,stripRadius)

E=[indCurve(1:end)' [indCurve(2:end) indCurve(1)]'];
ind1=indCurve(1:end)';
ind2=[indCurve(2:end) indCurve(1)]';
ind3=[indCurve(end) indCurve(1:end-1)]';

N1f=Vm(ind2,:)-Vm(ind1,:);
N1b=Vm(ind1,:)-Vm(ind3,:);
Ne=vecnormalize((N1f+N1b)/2);

% Ne=vecnormalize(edgeVec(E,Vm));
Nf=-Nm(E(:,1),:);% -vecnormalize((Nm(E(:,1),:)+Nm(E(:,2),:))/2);
Ne2=vecnormalize(cross(Nf,Ne));

X=repmat(Vm(E(:,1),1),1,nc);
Y=repmat(Vm(E(:,1),2),1,nc);
Z=repmat(Vm(E(:,1),3),1,nc);

t=repmat(linspace(0,pi/2,nc),size(Z,1),1);

X=X+stripRadius.*sin(t).*repmat(Ne2(:,1),1,nc)-stripRadius.*cos(t).*repmat(Nf(:,1),1,nc)+stripRadius.*repmat(Nf(:,1),1,nc);
Y=Y+stripRadius.*sin(t).*repmat(Ne2(:,2),1,nc)-stripRadius.*cos(t).*repmat(Nf(:,2),1,nc)+stripRadius.*repmat(Nf(:,2),1,nc);
Z=Z+stripRadius.*sin(t).*repmat(Ne2(:,3),1,nc)-stripRadius.*cos(t).*repmat(Nf(:,3),1,nc)+stripRadius.*repmat(Nf(:,3),1,nc);

for q=2:1:size(X,2)
   v=evenlySampleCurve([X(:,q) Y(:,q) Z(:,q)],size(X,1),'pchip',1);
   X(:,q)=v(:,1);
   Y(:,q)=v(:,2);
   Z(:,q)=v(:,3);
end

[Fr,Vr]=grid2patch(X,Y,Z,[],[1 0 0]);

end
%%

function [F,V,X,Y,Z]=bezierLoft(P1,P4,N1,N4,pointSpacing,f)

D12=sqrt(sum((P1-P4).^2,2));
numPoints=ceil(max(D12)./pointSpacing);
if numPoints<2
    numPoints=2;
end

P2=P1+D12.*f.*N1;
P3=P4-D12.*f.*N4;

X=zeros(numPoints,size(P1,1));
Y=zeros(numPoints,size(P1,1));
Z=zeros(numPoints,size(P1,1));
for q=1:1:size(P1,1)
    p=[P1(q,:); P2(q,:); P3(q,:); P4(q,:)]; %Control points
    V_bezier=bezierCurve(p,numPoints*2); %Compute bezier curve
    V_bezier=evenlySampleCurve(V_bezier,numPoints,'pchip'); %resample evenly
    X(:,q)=V_bezier(:,1);
    Y(:,q)=V_bezier(:,2);
    Z(:,q)=V_bezier(:,3);
end

%Create quad patch data
[F,V] = surf2patch(X,Y,Z);
I=[(1:size(Z,1)-1)' (1:size(Z,1)-1)' (2:size(Z,1))' (2:size(Z,1))' ];
J=[size(Z,2).*ones(size(Z,1)-1,1) ones(size(Z,1)-1,1) ones(size(Z,1)-1,1) size(Z,2).*ones(size(Z,1)-1,1)];
F_sub=sub2ind(size(Z),I,J);
F=[F;F_sub];
F=fliplr(F);

end

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
classdef regressionTest < matlab.unittest.TestCase

    properties
        testDir = ''
        plotSolvers = []  % 1 to plot new run comparison by sln method
        openCompare = []  % 1 opens all new run vs. stored run plots for comparison of each solver
        regular
        regularCIC
        regularSS
        irregularCIC
        irregularSS
        OriginalDefault
    end

    methods(Access = 'public')

        function obj = regressionTest(plotSolvers, openCompare)
            arguments
                plotSolvers (1,1) double = 1
                openCompare (1,1) double = 1
            end
            % Assign arguments to test Class
            obj.plotSolvers = plotSolvers;
            obj.openCompare = openCompare;
            % Set test directory
            obj.testDir = fileparts(mfilename('fullpath'));
            % Save the visibility state at construction
            obj.OriginalDefault = get(0,'DefaultFigureVisible');
        end

    end

    methods (TestMethodSetup)
        function killPlots (~)
            set(0,'DefaultFigureVisible','off');
        end
    end

    methods(TestClassSetup)

        function runRegTest(testCase)
            cd(fullfile(testCase.testDir,       ...
                        'RegressionTests',      ...
                        'RegularWaves',         ...
                        'regular'))
            runLoadRegular;
            testCase.regular = load('regular.mat').("regular");
            savefig(fullfile('..', 'figReg'));
            cd(testCase.testDir);
        end

        function runRegCICTest(testCase)
            cd(fullfile(testCase.testDir,       ...
                        'RegressionTests',      ...
                        'RegularWaves',         ...
                        'regularCIC'))
            runLoadRegularCIC;
            testCase.regularCIC = load('regularCIC.mat').("regularCIC");
            savefig(fullfile('..', 'figRegCIC'));
            cd(testCase.testDir);
        end

        function runRegSSTest(testCase)
            cd(fullfile(testCase.testDir,       ...
                        'RegressionTests',      ...
                        'RegularWaves',         ...
                        'regularSS'))
            runLoadRegularSS;
            testCase.regularSS = load('regularSS.mat').("regularSS");
            savefig(fullfile('..', 'figRegSS'));
            cd(testCase.testDir);
        end

        function runIrregCICTest(testCase)
            cd(fullfile(testCase.testDir,   ...
                        'RegressionTests',  ...
                        'IrregularWaves',   ...
                        'irregularCIC'))
            runLoadIrregularCIC;
            testCase.irregularCIC = ...
                                load('irregularCIC.mat').("irregularCIC");
            savefig(fullfile('..', 'figIrregCIC'));
            cd(testCase.testDir);
        end

        function runIrregSSTest(testCase)
            cd(fullfile(testCase.testDir,   ...
                        'RegressionTests',  ...
                        'IrregularWaves',   ...
                        'irregularSS'))
            runLoadIrregularSS;
            testCase.irregularSS = load('irregularSS.mat').("irregularSS");
            savefig(fullfile('..', 'figIrregSS'));
            cd(testCase.testDir);
        end

    end

    methods(TestClassTeardown)

        function plotRegTests(testCase)
            % Plot Solver Comparisons
            if testCase.plotSolvers == 1
                cd(fullfile(testCase.testDir,   ...
                            'RegressionTests',  ...
                            'RegularWaves'));
                printPlotRegular;
                cd(fullfile(testCase.testDir,   ...
                            'RegressionTests',  ...
                            'IrregularWaves'));
                printPlotIrregular;
            end
            % Open new vs. org Comparisons
            if testCase.openCompare == 1
                cd(fullfile(testCase.testDir,   ...
                            'RegressionTests',  ...
                            'RegularWaves'));
                openfig('figReg.fig');
                openfig('figRegCIC.fig');
                openfig('figRegSS.fig');
                cd(fullfile(testCase.testDir,   ...
                            'RegressionTests',  ...
                            'IrregularWaves'));
                openfig('figIrregCIC.fig');
                openfig('figIrregSS.fig');
            end
            set(0,'DefaultFigureVisible',testCase.OriginalDefault);
            testCase.assertEqual(get(0,'DefaultFigureVisible'),     ...
                                 testCase.OriginalDefault);
        end
    end

    methods(Test)

        function body1_reg_disp_heave(testCase)
            % Body1 Displacement in Heave
            tol = 1e-10;
            org = testCase.regular.B1.WEC_Sim_org.heave;
            new = testCase.regular.B1.WEC_Sim_new.heave;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['Body1 Displacement in Heave, Diff = '     ...
                     num2str(max(abs(org-new))) '\n']);
        end

        function body2_reg_disp_heave(testCase)
            % Body2 Displacement in Heave
            tol = 1e-10;
            org = testCase.regular.B2.WEC_Sim_org.heave;
            new = testCase.regular.B2.WEC_Sim_new.heave;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['Body2 Displacement in Heave, Diff = '     ...
                     num2str(max(abs(org-new))) '\n']);
        end

        function bodyRel_reg_disp_heave(testCase)
            % Relative Displacement in Heave
            tol = 1e-10;
            org = testCase.regular.Rel.WEC_Sim_org.heave;
            new = testCase.regular.Rel.WEC_Sim_new.heave;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['Relative Displacement in Heave, Diff = '  ...
                     num2str(max(abs(org-new))) '\n']);
        end

        function body1_regCIC_disp_heave(testCase)
            % Body1 Displacement in Heave
            tol = 1e-10;
            org = testCase.regularCIC.B1.WEC_Sim_org.heave;
            new = testCase.regularCIC.B1.WEC_Sim_new.heave;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['Body1 Displacement in Heave, Diff = '     ...
                      num2str(max(abs(org-new))) '\n']);
        end

        function body2_regCIC_disp_heave(testCase)
            % Body2 Displacement in Heave
            tol = 1e-10;
            org = testCase.regularCIC.B2.WEC_Sim_org.heave;
            new = testCase.regularCIC.B2.WEC_Sim_new.heave;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['Body2 Displacement in Heave, Diff = '     ...
                      num2str(max(abs(org-new))) '\n']);
        end

        function bodyRel_regCIC_disp_heave(testCase)
            % Relative Displacement in Heave
            tol = 1e-10;
            org = testCase.regularCIC.Rel.WEC_Sim_org.heave;
            new = testCase.regularCIC.Rel.WEC_Sim_new.heave;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['Relative Displacement in Heave, Diff = '  ...
                     num2str(max(abs(org-new))) '\n']);
        end

        function body1_regSS_disp_heave(testCase)
            % Body1 Displacement in Heave
            tol = 1e-10;
            org = testCase.regularSS.B1.WEC_Sim_org.heave;
            new = testCase.regularSS.B1.WEC_Sim_new.heave;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['Body1 Displacement in Heave, Diff = '     ...
                     num2str(max(abs(org-new))) '\n']);
        end

        function body2_regSS_disp_heave(testCase)
            % Body2 Displacement in Heave
            tol = 1e-10;
            org = testCase.regularSS.B2.WEC_Sim_org.heave;
            new = testCase.regularSS.B2.WEC_Sim_new.heave;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['Body2 Displacement in Heave, Diff = '     ...
                     num2str(max(abs(org-new))) '\n']);
        end

        function bodyRel_regSS_disp_heave(testCase)
            % Relative Displacement in Heave
            tol = 1e-10;
            org = testCase.regularSS.Rel.WEC_Sim_org.heave;
            new = testCase.regularSS.Rel.WEC_Sim_new.heave;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['Relative Displacement in Heave, Diff = '  ...
                     num2str(max(abs(org-new))) '\n']);
        end

        function body1_irreg_disp_heave(testCase)
            % Body1 Displacement in Heave
            tol = 1e-10;
            org = testCase.irregularCIC.B1.WEC_Sim_org.heave;
            new = testCase.irregularCIC.B1.WEC_Sim_new.heave;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['Body1 Displacement in Heave, Diff = '     ...
                     num2str(max(abs(org-new))) '\n']);
        end

        function body2_irreg_disp_heave(testCase)
            % Body2 Displacement in Heave
            tol = 1e-10;
            org = testCase.irregularCIC.B2.WEC_Sim_org.heave;
            new = testCase.irregularCIC.B2.WEC_Sim_new.heave;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['Body2 Displacement in Heave, Diff = '     ...
                     num2str(max(abs(org-new))) '\n']);
        end

        function bodyRel_irreg_disp_heave(testCase)
            % Relative Displacement in Heave
            tol = 1e-10;
            org = testCase.irregularCIC.Rel.WEC_Sim_org.heave;
            new = testCase.irregularCIC.Rel.WEC_Sim_new.heave;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['Relative Displacement in Heave, Diff = '  ...
                      num2str(max(abs(org-new))) '\n']);
        end

        function irreg_0th_Spectral_Moment(testCase)
            % 0th Order Spectral Moment
            tol = 1e-10;
            org = testCase.irregularCIC.Sp.WEC_Sim_org.m0;
            new = testCase.irregularCIC.Sp.WEC_Sim_new.m0;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['0th Order Spectral Moment, Diff = '       ...
                      num2str(max(abs(org-new))) '\n']);
        end

        function irreg_2nd_Spectral_Moment(testCase)
            % 2nd Order Spectral Moment
            tol = 1e-10;
            org = testCase.irregularCIC.Sp.WEC_Sim_org.m2;
            new = testCase.irregularCIC.Sp.WEC_Sim_new.m2;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['2nd Order Spectral Moment, Diff = '       ...
                      num2str(max(abs(org-new))) '\n']);
        end

        function body1_irregSS_disp_heave(testCase)
            % Body1 Displacement in Heave
            tol = 1e-10;
            org = testCase.irregularSS.B1.WEC_Sim_org.heave;
            new = testCase.irregularSS.B1.WEC_Sim_new.heave;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['Body1 Displacement in Heave, Diff = '     ...
                      num2str(max(abs(org-new))) '\n']);
        end

        function body2_irregSS_disp_heave(testCase)
            % Body2 Displacement in Heave
            tol = 1e-10;
            org = testCase.irregularSS.B2.WEC_Sim_org.heave;
            new = testCase.irregularSS.B2.WEC_Sim_new.heave;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['Body2 Displacement in Heave, Diff = '     ...
                     num2str(max(abs(org-new))) '\n']);
        end

        function bodyRel_irregSS_disp_heave(testCase)
            % Relative Displacement in Heave
            tol = 1e-10;
            org = testCase.irregularSS.Rel.WEC_Sim_org.heave;
            new = testCase.irregularSS.Rel.WEC_Sim_new.heave;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['Relative Displacement in Heave, Diff = '  ...
                      num2str(max(abs(org-new))) '\n']);
        end

        function irregSS_0th_Spectral_Moment(testCase)
            % 0th Order Spectral Moment
            tol = 1e-10;
            org = testCase.irregularSS.Sp.WEC_Sim_org.m0;
            new = testCase.irregularSS.Sp.WEC_Sim_new.m0;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['0th Order Spectral Moment, Diff = '       ...
                     num2str(max(abs(org-new))) '\n']);
        end

        function irregSS_2nd_Spectral_Moment(testCase)
            % 2nd Order Spectral Moment
            tol = 1e-10;
            org = testCase.irregularSS.Sp.WEC_Sim_org.m2;
            new = testCase.irregularSS.Sp.WEC_Sim_new.m2;
            testCase.verifyEqual(new,org,'AbsTol',tol);
            fprintf(['2nd Order Spectral Moment, Diff = '       ...
                     num2str(max(abs(org-new))) '\n']);
        end

    end

end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function ykplus1 = sparseGalerkinControl_Discrete(t,y,u,p)

if isfield(p,'SelectVars')==0
    p.SelectVars = 1:length(y);
end
polyorder = p.polyorder;
usesine = p.usesine;
ahat = p.ahat;
yPool = poolData([y(p.SelectVars)' u],length(p.SelectVars)+1,polyorder,usesine);
ykplus1 = (yPool*ahat)';

author: eurika-kaiser

repo: SINDY-MPC


text:
function midi = extractPitch(configMfcc, configEsi, x, mfccHmm, esiHmm)
%EXTRACTPITCH   Function for A/U/V decision and pitch tracking.
%   Format: midi = extractPitch(configMfcc, configEsi, x, mfccHmm, esiHmm
%   Inputs:
%       configMfcc: Configuration for MFCC features.
%       configEsi:  Configuration for ESI features.
%       x:          Input signal.
%       mfccHmm:    HMM for MFCC features.
%       esiHmm:     HMM for ESI features.
%   Output:
%       midi:       Row vector containing the exxtracted pitch contour, in
%                   midi number scale, one number per frame. 0 stands for
%                   accompaniment, and -1 stands for unvoiced regions.

    mfcc = mfccFeature(configMfcc, x);
    sal = salience(configEsi, x);
    esi = esiFeature(configEsi, sal);

    midi = viterbi(mfcc, mfccHmm);
    segs = findSegs(midi == 3); % find voiced segs
    midi(midi == 1) = 0;        % accom
    midi(midi == 2) = -1;       % unvoiced
    for s = 1:size(segs,1)
        ind = segs(s,1) : segs(s,2);
        midi(ind) = refinePitch(configEsi, viterbi(esi(:,ind), esiHmm), sal(:,ind));
    end
end

function midi = refinePitch(config, state, sal)
% Refine pitch contour by picking the peaks in the salience map around the
%   coarse pitch contour.
    frames = length(state);
    midi = config.midisCoarse(state);
    for f = 1:frames
        ind = find((config.midisFine >= midi(f) - 0.5) & (config.midisFine <= midi(f) + 0.5));
        midis = config.midisFine(ind);
        [temp ind] = max(sal(ind, f));
        midi(f) = midis(ind);
    end
end

function segs = findSegs(bool)
% Find continuous segments of 1 in bool.
    bool = bool(:);
    start = bool & ~[0; bool(1:end-1)];
    finish = bool & ~[bool(2:end); 0];
    segs = [find(start), find(finish)];
end

author: MaigoAkisame


text:
function B = tocstruct(A)
    if iscstruct(A)
        B = A;
    elseif isnumeric(A)
        B.real = real(A);
        B.imag = imag(A);
    else
        error('Input does not have the expected format.');
    end
end

author: NicolasBoumal

repo: manopt


text:
function [cm] = m2cm(m)
% Convert length from meters to centimeters.
% Chad A. Greene 2012
cm = m*100;

author: Sable

repo: mcbench-benchmarks


text:
%% DNLS Example1
clc
% ODE System
 ode = @(t,z,p) [-p(1)*z(1) + 4;
                 2*z(1) - p(1)*z(2) + 5;
                 -4*z(1) - 2*z(2)*z(3) - p(2)];

% Initial Conditions
 z0 = [-1.5;1.25;1];

% True Parameter Values
 p = [2.345;1.1];

% Generate Fitting Data
 tm  = 0:0.1:2;                 %measurement times
 odeInt = @(t,z) ode(t,z,p);    %ODE function for ODE45
[~,zm] = ode45(odeInt,tm,z0);  %Solve ODEs

% Starting Guess
 theta0 = [1;0.5];

% Create OPTI Object
 Opt = opti('ode',ode,'data',tm,zm,'z0',z0,'theta0',theta0)

% Solve
[theta,fval,exitflag,info] = solve(Opt)

% Plot the Solution
plot(Opt)

%% DNLS Example2
clc
% Indicate Initial Conditions to Solve for Using NaNs
 z0 = [-1.5;NaN;NaN];

% Starting Guess [p;z0]
 theta0 = [1;0.5;0.5;0.5];

% Create OPTI Object
 Opt = opti('ode',ode,'data',tm,zm,'z0',z0,'theta0',theta0)

% Solve
[theta,fval,exitflag,info] = solve(Opt)

% Plot the Solution
plot(Opt)

%% Example 3
clc

% Add Bounds
 lb = [1.5;0;0.1;0.1];
 ub = [3.5;2;1.5;1.5];

% Create OPTI Object
 Opt = opti('ode',ode,'data',tm,zm,'bounds',lb,ub,'z0',z0,'theta0',theta0)

% Solve
[theta,fval,exitflag,info] = solve(Opt)

% Plot the Solution
plot(Opt)

%% DNLS Example4
clc
% Set Dynamic Options (specifying states of interest)
 dopts = optidynset('stateIndex',[2 3]);

% Add to General OPTI Options
 opts = optiset('display','iter','dynamicOpts',dopts);

% Index Measured Data
 zm23 = zm(:,[2 3]);

% Create OPTI Object
 Opt = opti('ode',ode,'data',tm,zm23,'z0',z0,'theta0',theta0,'options',opts)

% Solve
[theta,fval,exitflag,info] = solve(Opt)

% Plot the Solution
plot(Opt)

%% DNLS Example5
clc
% Initial Conditions
z0 = [-1.5;1.25;1];

% Measurement Times for each State
tm2  = 0:0.1:2;                %state 2
tm3  = 0.2:0.2:2;              %state 3

% Solve ODEs and index Measurement Data
[~,zm2] = ode45(odeInt,tm2,z0); zm2 = zm2(:,2);
% Note 0 is required below just to match initial condition in this example
[~,zm3] = ode45(odeInt,[0 tm3],z0); zm3 = zm3(2:end,3); %drop first point

% Group measurements and time stamps in cell arrays
tm_multi = {tm2;tm3};
zm_multi = {zm2;zm3};

% Indicate Initial Conditions to Solve for Using NaNs
 z0 = [-1.5;NaN;NaN];

% Create OPTI Object
Opt = opti('ode',ode,'data',tm_multi,zm_multi,'z0',z0,'theta0',theta0,'options',opts)

% Solve
[theta,fval,exitflag,info] = solve(Opt)

% Plot the Solution
plot(Opt)

%% DNLS Example6
clc
% Initial Conditions
z0 = [-1.5;1.25;1];

% Measurement Times for each State
tm1  = 0.5:0.1:2;   %state 1
tm2  = 1:0.1:2;    %state 2
tm3  = 0.2:0.2:2;   %state 3

% Solve ODEs and index Measurement Data
% Note 0 is required below just to match initial condition in this example
[~,zm1] = ode45(odeInt,[0 tm1],z0); zm1 = zm1((2:end),1);
[~,zm2] = ode45(odeInt,[0 tm2],z0); zm2 = zm2((2:end),2);
[~,zm3] = ode45(odeInt,[0 tm3],z0); zm3 = zm3((2:end),3);

% Group measurements and time stamps in cell arrays
tm_multi = {tm1;tm2;tm3};
zm_multi = {zm1;zm2;zm3};

% Indicate Initial Conditions to Solve for Using NaNs
 z0 = [NaN;NaN;NaN];

% New Initial Guess Vector [p;z0];
 theta0 = [1;0.5;0.5;0.5;0.5];

% Create OPTI Object
dopts = optidynset('initialT',0);
opts = optiset('display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm_multi,zm_multi,'z0',z0,'theta0',theta0,'options',opts)

% Solve
[theta,fval,exitflag,info] = solve(Opt)

% Plot the Solution
plot(Opt)

%% DNLS Example7
clc
% Initial Conditions
z0 = [-1.5;1.25;1];

% Measurement Times for each State
tm1  = 0.5:0.1:2;   %state 1
tm1b = 1:0.25:2;   %state 1
tm2  = 1:0.1:2;    %state 2
tm3  = 0.2:0.2:2;   %state 3

%ODE with modified parameters (representing a different run)
odeIntB = @(t,z) ode(t,z,p*0.85);

% Solve ODEs and index Measurement Data
% Note 0 is required below just to match initial condition in this example
[~,zm1]  = ode45(odeInt,[0 tm1],z0);  zm1  = zm1((2:end),1);
[~,zm1b] = ode45(odeIntB,[0 tm1b],z0); zm1b = zm1b((2:end),1);
[~,zm2]  = ode45(odeInt,[0 tm2],z0);  zm2  = zm2((2:end),2);
[~,zm3]  = ode45(odeInt,[0 tm3],z0);  zm3  = zm3((2:end),3);
[~,zm3b] = ode45(odeInt,[0 tm3],z0*0.75); zm3b = zm3b((2:end),3);

% Group measurements and time stamps in cell arrays
tm_multi = {[tm1 tm1b];tm2;[tm3 tm3]};
zm_multi = {[zm1;zm1b];zm2;[zm3;zm3b]};

% Indicate Initial Conditions to Solve for Using NaNs
 z0 = [NaN;NaN;NaN];

% New Initial Guess Vector [p;z0];
 theta0 = [1;0.5;0.5;0.5;0.5];

% Create OPTI Object
dopts = optidynset('initialT',0);
opts = optiset('display','iter','dynamicOpts',dopts);
Opt = opti('ode',ode,'data',tm_multi,zm_multi,'z0',z0,'theta0',theta0,'options',opts)

% Solve
[theta,fval,exitflag,info] = solve(Opt)

% Plot the Solution
plot(Opt)

%% DNLS Example8
clc
% Declare Symbolic Variables
syms p1 p2 z1 z2 z3
% Symbolic ODE
ODE = [-p1*z1 + 4;
      2*z1 - p1*z2 + 5;
      -4*z1 - 2*z2*z3 - p2];

% Evaluate Jacobians
dfdz_sym = jacobian(ODE,[z1 z2 z3])
dfdp_sym = jacobian(ODE,[p1 p2])

[dfdz,dfdp] = symDynJac(ode)

% Analytical Derivative Expressions
 dfdz = @(t,z,p) [-p(1) 0,       0;
                  2,    -p(1),   0;
                  -4,   -2*z(3), -2*z(2)];
 dfdp = @(t,z,p) [-z(1), 0;
                  -z(2), 0;
                  0,    -1];

% Set Dynamic Options
 dopts = optidynset('dfdz',dfdz,'dfdp',dfdp,'sensitivity','User','initialT',0);

% General OPTI Options
 opts = optiset('display','iter','derivCheck','on','dynamicOpts',dopts);

% Create OPTI Object
  Opt = opti('ode',ode,'data',tm_multi,zm_multi,'z0',z0,'theta0',theta0,...
             'options',opts)

% Solve
[theta,fval,exitflag,info] = solve(Opt)

% Plot the Solution
plot(Opt)

%% DNLS Ex9
clc
%ODE to fit
ode = @(t,z,p) p(1)*z(1)^2 - z(1)^3;
z0 = 0.001; %ic

%Generate measurement data
p = 1.5;                          %true parameter value
oi = @(t,z) ode(t,z,p);           %ode integrator function
tm = [0 65:70 80:1:90 100]*10;  %measurement times
[~,zm] = ode15s(oi,tm,z0);        %measurements

% Starting Guess
theta0 = 1;

%Build OPTI Object
dopts = optidynset('sensitivity','none');
opts1 = optiset('display','iter','iterfun',@optiplotlogfval,'dynamicOpts',dopts);
opts2 = optiset('display','iter','iterfun',@optiplotlogfval);

% Create OPTI Objects
OptNoDer = opti('ode',ode,'data',tm,zm,'z0',z0,'theta0',theta0,'options',opts1)
OptWDer = opti('ode',ode,'data',tm,zm,'z0',z0,'theta0',theta0,'options',opts2)

% Solve
[theta,fval,exitflag,info] = solve(OptNoDer)
[theta,fval,exitflag,info] = solve(OptWDer)

% Plot the Solution
subplot(121)
plot(OptNoDer); title('No Sensitivity');
subplot(122)
plot(OptWDer); title('With Sensitivity');

%% DNS Ex9b
clc
% Partial Derivatives
dfdz = @(t,z,p) 2*p(1)*z(1) - 3*z(1)^2;
dfdp = @(t,z,p) z(1)^2;

% Set Dynamic Options
dopts = optidynset('dfdz',dfdz,'dfdp',dfdp,'integrator','ode15s',...
                   'sensitivity','User');

% General OPTI Options
opts = optiset('display','iter','derivCheck','on','dynamicOpts',dopts);

% Create OPTI Object
Opt = opti('ode',ode,'data',tm,zm,'z0',z0,'theta0',theta0,'options',opts)

% Solve
[theta,fval,exitflag,info] = solve(Opt)

% Plot
plot(Opt)

%% DNS Ex9c
clc
% Set Dynamic Options
dopts = optidynset('integrator','ode15s','sensitivity','None');

% General OPTI Options
opts = optiset('solver','nomad','display','iter','dynamicOpts',dopts);

% Create OPTI Object
Opt = opti('ode',ode,'data',tm,zm,'bounds',0.5,2.5,'z0',z0,'theta0',theta0,'options',opts)

% Solve
[theta,fval,exitflag,info] = solve(Opt)

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [ success ] = write_lin_dyn_svr( location, means, w, b, cutoff)
%WRITE_LIN_SVR Summary of this function goes here
%   Detailed explanation goes here

    fileID = fopen(location, 'w');

    if(fileID ~= -1)

        % Write the regressor type 1 - linear dynamic SVR
        fwrite(fileID, 1, 'uint');

        fwrite(fileID, cutoff, 'float64');

        writeMatrixBin(fileID, means, 6);
        writeMatrixBin(fileID, w, 6);
        fwrite(fileID, b, 'float64');

        fclose(fileID);

        success = true;

    else
        success = false;
    end

end

author: TadasBaltrusaitis

repo: OpenFace


text:
function out = std(f)
%STD   Standard deviation of a CHEBFUN.
%   STD(F) is the standard deviation of the CHEBFUN F.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% The standard deviation is the square root SQRT() of the variance VAR():
out = sqrt(var(f));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
clear
path('./icon/',path);
Init;

%Constant value
RAD2DEG = 57.2957795;
DEG2RAD = 0.0174533;
%throttle when UAV is hovering
THR_HOVER = 0.609;
%% Initial condition
ModelInit_PosE = [0, 0, -100];
ModelInit_VelB = [0, 0, 0];
ModelInit_AngEuler = [0, 0, 0];
ModelInit_RateB = [0, 0, 0];
ModelInit_Rads = 557.142;

%max control angle,default 35deg
MAX_CONTROL_ANGLE_ROLL = 35;
MAX_CONTROL_ANGLE_PITCH  = 35;
%max control angle rate,rad/s
MAX_CONTROL_ANGLE_RATE_PITCH = 220;
MAX_CONTROL_ANGLE_RATE_ROLL = 220;
MAX_CONTROL_ANGLE_RATE_Y = 200;
%Maximum control speed, m/s
MAX_CONTROL_VELOCITY_XY = 5;
MAX_CONTROL_VELOCITY_Z = 3;
%Throttle amplitude
MAX_MAN_THR = 0.9;
MIN_MAN_THR = 0.05;
%% run simulink model
PosControl_tune

author: RflySim

repo: RflyExpCode

save_path: github-repos/MATLAB/RflySim-RflyExpCode


text:
function epsn = ReduceBoundary(eF, k, MaxK, cp)
% The shrink of the dynamic constraint boundary

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Ruwang Jiao

    z        = 1e-8;
    Nearzero = 1e-15;
    B        = MaxK./power(log((eF + z)./z), 1.0./cp);
    B(B==0)  = B(B==0) + Nearzero;
    f        = eF.* exp( -(k./B).^cp );
    tmp      = find(abs(f-z) < Nearzero);
    f(tmp)   = f(tmp).*0 + z;
    epsn     = f - z;
    epsn(epsn<=0) = 0;
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%% meshDistMarch
% Below is a demonstration of the features of the |meshDistMarch| function

%%
clear; close all; clc;

%% Syntax
% |[d,seedIndex]=meshDistMarch(F,V,indStart,optionStruct);|

%% Description
% The meshDistMarch function can be used to compute distances on meshes.
% The distances can be used for points sampling on the mesh or for
% remeshing.
% The function can operate on edge descriptions or face descriptions.
% Therefore for volumetric meshes (e.g. consisting of tetrahedra or
% hexahedra) appropriate face or edge data should be computed first to
% formulate the input.
%
% Input:
% E: the edges or faces for the mesh. E.g. an nx2 edge matrix or an nxm
% face matrix (n faces, m corners per face)
% V: the vertices for the mesh
% indStart: indices for one or more points to compute distances from
% optionStruct.toleranceLevel : The tolerance level for convergence. 0 is
% the default.
% optionStruct.numSeeds: Defines the number of seeds to generate on the mesh. Default is equal to startInd.
% optionStruct.waitBarOn=0; %Turn on/off waitbar
% optionStruct.unitEdgeOn=1; %Turn on/off the use of unit edge lengths
%
% Output:
% d: distances (from the start/seed points) on the mesh vertices
% seedIndex: nearest seed (or start) point indices for each vertex, forming
% a quasi-Voronoi tesselation.
% If the second output is not requested the performance is enhanced.

%% Examples

%%
% Plot settings
cMapDist=flipud(igviridis(250));
[cMapIndices,scrambleIndices]=scramble(viridis(250),1); %Colormap

faceAlpha1=1;
faceAlpha2=0.65;
fontSize=25;
markerSize=50;

%% Example: Edge or graph data

%%
% Create branching example. E defines edges and V is a vertex array

n=75;
x=linspace(0,0.9*pi,n);
y=sin(x);
V=[x(:) y(:) zeros(size(x(:)))];
V=evenlySampleCurve(V,n,'pchip',0);
x=V(:,1);
y=V(:,2);
dV=vecnormalize(diff(V,1,1));
VdV=V(1:1:end-1,:);
placePoints=1:3:n-1;
E=[(1:1:n-1)' (2:1:n)'];
numSteps=numel(placePoints);
l=linspace(0.75,0.25,numSteps);
a=linspace(0.25*pi,0.25*pi,numSteps);
for q=1:1:numSteps
    R=euler2DCM([0 0 a(q)]);
    vr1=dV(placePoints(q),:)*R.*l(q);
    vr2=dV(placePoints(q),:)*R'.*l(q);
    Vb=linspacen([x(placePoints(q)) y(placePoints(q)) 0],vr1+[x(placePoints(q)) y(placePoints(q)) 0],n)';
    Eb=[(1:1:n-1)' (2:1:n)'];
    E=[E; Eb+size(V,1)];
    V=[V; Vb];
    Vb=linspacen([x(placePoints(q)) y(placePoints(q)) 0],vr2+[x(placePoints(q)) y(placePoints(q)) 0],n)';
    Eb=[(1:1:n-1)' (2:1:n)'];
    E=[E; Eb+size(V,1)];
    V=[V; Vb];
end
E=[E;E+size(V,1);];
V2=V;
V2(:,2)=-V2(:,2);
V2(:,1)=-V2(:,1);
V=[V;V2];
E=[E;E+size(V,1);];
V2=V;
V2(:,2)=-V2(:,2);
V=[V;V2];
[E,V]=mergeVertices(E,V);

%%
% Compute distances on mesh

%Option set
indStart=1; %Index of the start point
optionStruct.toleranceLevel=0; %Tolerance for convergence
optionStruct.numSeeds=1; %Number of seeds
optionStruct.waitBarOn=0; %Turn on/off waitbar

%Compute distances on mesh description
d=meshDistMarch(E,V,indStart,optionStruct);

%%
% Visualization

cFigure; hold on;
title('Distances on an edge / graph model','fontSize',fontSize);
hp(1)=gpatch(E,V,'none',d,1,4);
hp(2)=plotV(V(indStart,:),'k.','MarkerSize',markerSize);
legend(hp,{'Mesh distances','Start point(s)'},'Location','SouthOutSide');
axisGeom;
view(2);
colorbar;
drawnow;

%% Example: Triangulated data

%%
% Get example triangulated mesh data
[F,V]=stanford_bunny;%graphicsModels(1);

%%
% Compute distances on mesh

%Option set
[~,indStart]=min(V(:,1)); %Index of the start point
optionStruct.toleranceLevel=0; %Tolerance for convergence
optionStruct.numSeeds=1; %Number of seeds
optionStruct.waitBarOn=0; %Turn on/off waitbar

%Compute distances on mesh description
d=meshDistMarch(F,V,indStart,optionStruct);

%%
% Visualization

cFigure;
subplot(1,2,1);hold on;
title('Distances on a triangulated surface model','fontSize',fontSize);
gpatch(F,V,'kw',d,1);
hp(1)=gpatch(F,V,'none',d,1,2);
hp(2)=plotV(V(indStart,:),'k.','MarkerSize',markerSize);
legend(hp,{'Mesh distances','Start point(s)'},'Location','SouthOutSide');
axisGeom;
colorbar;
camlight headlight;

subplot(1,2,2);hold on;
title('Distances on a triangulated surface model','fontSize',fontSize);
hp(1)=gpatch(F,V,d,'none',1);
hp(2)=plotV(V(indStart,:),'k.','MarkerSize',markerSize);
legend(hp,{'Mesh distances','Start point(s)'},'Location','SouthOutSide');
axisGeom;
camlight headlight;
colorbar;
drawnow;

%% Example: Quadrangulated data

%%
% Get example quadrangulated mesh data
n=4;
r=1;
[F,V]=quadSphere(n,r);

%%
% Compute distances on mesh

%Option set
[~,indStart]=min(V(:,1)); %Index of the start point
optionStruct.toleranceLevel=0; %Tolerance for convergence
optionStruct.numSeeds=1; %Number of seeds
optionStruct.waitBarOn=0; %Turn on/off waitbar

%Compute distances on mesh description
d=meshDistMarch(F,V,indStart,optionStruct);

%%
% Visualization

cFigure;
subplot(1,2,1);hold on;
title('Distances on a quad surface model','fontSize',fontSize);
gpatch(F,V,'kw',d,1);
hp(1)=gpatch(F,V,'none',d,1,2);
hp(2)=plotV(V(indStart,:),'k.','MarkerSize',markerSize);
legend(hp,{'Mesh distances','Start point(s)'},'Location','SouthOutSide');
axisGeom;
colorbar;
camlight headlight;

subplot(1,2,2);hold on;
title('Distances on a quad surface model','fontSize',fontSize);
hp(1)=gpatch(F,V,d,'none',1);
hp(2)=plotV(V(indStart,:),'k.','MarkerSize',markerSize);
legend(hp,{'Mesh distances','Start point(s)'},'Location','SouthOutSide');
axisGeom;
camlight headlight;
colorbar;
drawnow;

%% Example: Hexahedral mesh data

%%
% Get example hexahedral mesh data
cubeDimensions=[1 1 1];
cubeElementNumbers=[10 10 10];
outputStructType=2; %A structure compatible with mesh view
[meshStruct]=hexMeshBox(cubeDimensions,cubeElementNumbers,outputStructType);

%Access elements, nodes, and faces from the structure
HEX=meshStruct.elements; %The elements
F=meshStruct.faces; %The faces
V=meshStruct.nodes; %The nodes (vertices)

%Get edges for mesh marching
E=patchEdges(F); %Mesh edges

%%
% Compute distances on mesh

%Option set
[~,indStart]=min(V(:,1)); %Index of the start point
optionStruct.toleranceLevel=0; %Tolerance for convergence
optionStruct.numSeeds=1; %Number of seeds
optionStruct.waitBarOn=0; %Turn on/off waitbar

%Compute distances on mesh description
d=meshDistMarch(E,V,indStart,optionStruct);

%%
% Plotting model boundary surfaces and a cut view

cFigure; hold on
title('Distances on a hexahedral mesh','FontSize',fontSize);
hp(1)=gpatch(E,V,'none',d,1,3);
hp(2)=plotV(V(indStart,:),'k.','MarkerSize',markerSize);
legend(hp,{'Mesh distances','Start point(s)'},'Location','SouthOutSide');
colorbar;
axisGeom(gca,fontSize);
drawnow;

%% Example: Tetrahedral mesh data

%%
% Get example tetrahedral mesh data
[TET,V]=hex2tet(HEX,V,[],2);
F=element2patch(TET,V); %Element faces

%Get edges for mesh marching
E=patchEdges(F); %Mesh edges

%%
% Compute distances on mesh

%Option set
[~,indStart]=min(abs(V(:,1)-min(V(:,1)))+abs(V(:,2)-max(V(:,2)))+abs(V(:,3)-min(V(:,3)))); %Index of the start point
optionStruct.toleranceLevel=0; %Tolerance for convergence
optionStruct.numSeeds=1; %Number of seeds
optionStruct.waitBarOn=0; %Turn on/off waitbar

%Compute distances on mesh description
d=meshDistMarch(E,V,indStart,optionStruct);

%%
% Plotting model boundary surfaces and a cut view

cFigure; hold on
title('Distances on a tetrahedral mesh','FontSize',fontSize);
hp(1)=gpatch(E,V,'none',d,1,3);
hp(2)=plotV(V(indStart,:),'k.','MarkerSize',markerSize);
legend(hp,{'Mesh distances','Start point(s)'},'Location','SouthOutSide');
colorbar;
axisGeom(gca,fontSize);
drawnow;

%% Example: Mesh types effecting marching based distance computation
% In this example 3 mesh variations for a sphere are created. For each the
% start point is the point with the minimum x coordinate, i.e. a point on
% the far left of the plot. For each the true distance should be a smooth
% gradient with a maximum which equals pi.
% The marching algorithm uses the mesh to march and compute distances. One
% can imaging that straight paths (geodesically speaking) yeild the
% shortest distances. Hence zig-zag pattern (geodesically speaking) create
% false increased distances. This is not compensated for in this algorithm.
% Therefore the distance map depends on the mesh type and mesh
% connectivity. In the below example the regular quadrilateral mesh
% contains a straight path to the farthest point allong the equator. With
% mesh refinement the maximum distance would therefore converge on pi.
% However, other directions for the quadrilateral mesh produce zig-zag
% patterns causing distance to be altered. Hence the distance map departs
% from the smooth gradient allong the x-direction, one would expect. The
% triangulated mesh converted from the quadrilateral mesh contains both a
% straight path at the equator and also improved connectivity for other
% directions. The uniform geodesic triangulation contains the most smooth
% distance map.

%%
% Getting 3 different mesh types for a sphere.

n=3;
r=1;
[F1,V1]=quadSphere(n,r);
[F2,V2]=quad2tri(F1,V1);
[F3,V3]=geoSphere(n,r);

%%
% Compute distances on mesh

% -> Surface 1
%Option set
[~,indStart1]=min(V1(:,1)); %Index of the start point
optionStruct.toleranceLevel=0; %Tolerance for convergence
optionStruct.numSeeds=1; %Number of seeds
optionStruct.waitBarOn=0; %Turn on/off waitbar

%Compute distances on mesh description
[d1,i1]=meshDistMarch(F1,V1,indStart1,optionStruct);

% -> Surface 2
%Option set
[~,indStart2]=min(V2(:,1)); %Index of the start point
optionStruct.toleranceLevel=0; %Tolerance for convergence
optionStruct.numSeeds=1; %Number of seeds
optionStruct.waitBarOn=0; %Turn on/off waitbar

%Compute distances on mesh description
[d2,i2]=meshDistMarch(F2,V2,indStart2,optionStruct);

% -> Surface 3
%Option set
[~,indStart3]=min(V3(:,1)); %Index of the start point
optionStruct.toleranceLevel=0; %Tolerance for convergence
optionStruct.numSeeds=1; %Number of seeds
optionStruct.waitBarOn=0; %Turn on/off waitbar

%Compute distances on mesh description
[d3,i3]=meshDistMarch(F3,V3,indStart3,optionStruct);

%%
% Visualization

cFigure;
subplot(1,3,1);hold on;
title('Regular quad mesh','fontSize',fontSize);
hp(1)=gpatch(F1,V1,d1,d1);
hp(2)=plotV(V1(indStart1,:),'k.','MarkerSize',markerSize);
legend(hp,{'Mesh distances','Start point(s)'},'Location','SouthOutSide');
axisGeom;
colorbar; caxis([0 pi]);
camlight headlight;
view(2);

subplot(1,3,2);hold on;
title('Irregular quad mesh','fontSize',fontSize);
hp(1)=gpatch(F2,V2,d2,d2);
hp(2)=plotV(V2(indStart2,:),'k.','MarkerSize',markerSize);
legend(hp,{'Mesh distances','Start point(s)'},'Location','SouthOutSide');
axisGeom;
colorbar; caxis([0 pi]);
camlight headlight;
view(2);

subplot(1,3,3);hold on;
title('Geodesic triangulated mesh','fontSize',fontSize);
hp(1)=gpatch(F3,V3,d3,d3);
hp(2)=plotV(V3(indStart3,:),'k.','MarkerSize',markerSize);
legend(hp,{'Mesh distances','Start point(s)'},'Location','SouthOutSide');
axisGeom;
colorbar; caxis([0 pi]);
camlight headlight;
view(2);

drawnow;

%% Example: Using |meshDistMarch| for geodesic point sampling

%%
% Get example triangulated mesh data
[F,V]=graphicsModels(7);

%%
% Compute distances on mesh

numSeeds=250;

%Option set
[~,indStart]=min(V(:,1)); %Index of the start point
optionStruct.toleranceLevel=0; %Tolerance for convergence
optionStruct.numSeeds=numSeeds; %Number of seeds
optionStruct.waitBarOn=1; %Turn on/off waitbar

%Use weigths based on z-direction
% W=V(:,3);
% W=W-min(W);
% W=W./max(W);
% W=(W*9)+1;
% optionStruct.W=W;

%Compute distances on mesh description
[d,seedIndex]=meshDistMarch(F,V,indStart,optionStruct);
[indSeeds,~,ind2]=unique(seedIndex);

%%
% Visualization

cFigure;
subplot(1,2,1); hold on;
title('Distances on a triangulated surface model','fontSize',fontSize);
hp(1)=gpatch(F,V,d,'none',1); hp(1).FaceColor='Interp';
hp(2)=plotV(V(indSeeds,:),'k.','MarkerSize',markerSize);
legend(hp,{'Mesh distances','Seed point(s)'},'Location','SouthOutSide');
axisGeom;
camlight headlight;
colormap(gca,cMapDist); colorbar;

subplot(1,2,2); hold on;
title('Seed indices','fontSize',fontSize);
hp(1)=gpatch(F,V,'kw',ind2,1);
hp(2)=plotV(V(indSeeds,:),'k.','MarkerSize',markerSize);
legend(hp,{'Mesh distances','Seed point(s)'},'Location','SouthOutSide');
axisGeom;
camlight headlight;
colormap(gca,cMapIndices); %icolorbar;
drawnow;

%% Example: Using |meshDistMarch| for geodesic surface resampling
% See also: |remeshTriSurfDistMap|

[Fd,Vd,indSeed]=seedIndex2triangulation(F,V,seedIndex);

%%
% Visualization

cFigure;
subplot(1,2,1); hold on;
title('Distances on a triangulated surface model','fontSize',fontSize);
hp(1)=gpatch(F,V,d,'none',1); hp(1).FaceColor='Interp';
hp(2)=plotV(V(indSeeds,:),'k.','MarkerSize',markerSize);
legend(hp,{'Mesh distances','Seed point(s)'},'Location','SouthOutSide');
axisGeom;
camlight headlight;
colormap(gca,cMapDist); colorbar;

subplot(1,2,2); hold on;
title('Resampled surface model','fontSize',fontSize);
plotV(V(indSeed,:),'k.','MarkerSize',50);
hp(1)=gpatch(F,V,'kw','none',0.5);
hp(2)=gpatch(Fd,Vd,'gw','k',1,2);
legend(hp,{'Original mesh','Resampled mesh'},'Location','SouthOutSide');

axisGeom;
camlight headlight;
drawnow;

%% Example: Using unit edge lengths
% Forcing unit edge lenghts means each edge is considered equally long i.e.
% a length of 1. Hence when used in combination with resampling, this
% causes the algorithm to resample dense regions in a dense fashion and
% coarse regions in a coarse fashion. Therefore using
% optionStruct.unitEdgeOn=1; one can force the resampling to have similar
% degrees of relative density differences. In the example below a dinosaur
% mesh with a fine mesh at the limbs and a coarse mesh on the main body is
% resampled using unit edge lengths. The output can be seen to remain
% refined at the limbs.

%%
% Get example triangulated mesh data
[F,V]=graphicsModels(4);
[F,V]=subtri(F,V,2);

%%
% Compute distances on mesh

%Option set
[~,indStart]=min(V(:,1)); %Index of the start point
optionStruct.toleranceLevel=0; %Tolerance for convergence
optionStruct.numSeeds=1000; %Number of seeds
optionStruct.waitBarOn=1; %Turn on/off waitbar
optionStruct.unitEdgeOn=1;

%Compute distances on mesh description
[d,seedIndex]=meshDistMarch(F,V,indStart,optionStruct);
[indSeeds,~,ind2]=unique(seedIndex);

%Compute resampled surface
[Fd,Vd,indSeed]=seedIndex2triangulation(F,V,seedIndex);

%%
% Visualization

cFigure;
subplot(1,2,1); hold on;
title('Distances on a triangulated surface model','fontSize',fontSize);
hp(1)=gpatch(F,V,d,'none',1); hp(1).FaceColor='Interp';
hp(2)=plotV(V(indSeeds,:),'k.','MarkerSize',25);
legend(hp,{'Mesh distances','Seed point(s)'},'Location','SouthOutSide');
axisGeom;
camlight headlight;
colormap(gca,cMapDist); colorbar;

subplot(1,2,2); hold on;
title('Resampled surface model','fontSize',fontSize);
plotV(V(indSeed,:),'k.','MarkerSize',25);
hp(1)=gpatch(F,V,'kw','none',0.5);
hp(2)=gpatch(Fd,Vd,'gw','k',1,2);
legend(hp,{'Original mesh','Resampled mesh'},'Location','SouthOutSide');

axisGeom;
camlight headlight;
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.


text:
function f1 = p06_f1 ( x )

%*****************************************************************************80
%
%% P06_F1 evaluates the first derivative for problem 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the value of the variable.
%
%    Output, real F1, the first derivative of the
%    objective function.
%
  f1 = -1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [X, A, B, t, varargout] = VBMC(P, Y, options)
% Variational Bayesian Low Rank Matrix Completion
% Author (a.k.a. person to blame): S. Derin Babacan
%
% Last updated: January 31, 2012
% Usage:
% [X, A, B, t, varargout] = VSBL_MAT(P, Y, options)
% Solves for X = AB' in
% Y = P(X) + N = P(AB') + N  where X, A, B are low rank,
% N is dense Gaussian noise, and P is a binary sampling matrix.
%
% If you're using this code, please acknowledge
% Reference:
%    S. D. Babacan, M. Luessi, R. Molina, and A. K. Katsaggelos,
%    "Sparse Bayesian Methods for Low-Rank Matrix Estimation,"
%    IEEE Transactions on Signal Processing, 2012.
%
% This code is not optimized for speed. It implements the full variational
% Bayesian inference described in the paper, without any manipulation of
% the matrices. These are described in the paper above but not implemented
% in this version of the code. These manipulations might lead to significant
% decrease in running times depending on the size of matrix Y.
%
% -----------------------------------------------------------------------
%                                INPUTS
%   P               : binary (0-1) sampling matrix same size as Y
%   Y               : input matrix
%   options         : All options are *optional*. The default values are set
%                     automatically.
%
%   verbose       : output the progress? (0/1) default: 0.
%   init          : Initialization method.
%                   'rand': initialize A and B with random matrices
%                   'ml'  : Apply SVD to Y and initialize A and B using its factors.
%                   (default)
%   a_gamma0
%   b_gamma0      : hyperparameter values for gamma. default: 1e-6
%
%
%   a_beta0
%   b_beta0      : hyperparameter values for beta. default: 0
%
%   inf_flag      : Flag for inference of components of E. default: 1
%                    1: Standard
%                    2: fixed point MacKay
%   MAXITER       : max number of iterations. default: 100
%   UPDATE_BETA   : update noise variance? (0/1). default: 1
%   UPDATE_BETA_START : iteration number to start updating the noise variance
%   DIM_RED       : prune irrelevant dimensions during iterations? (0/1).  default: 1
%   DIMRED_THR    : threshold to remove columns from A and B. Used only if
%                   DIM_RED is 1.  default: 1e4
%   initial_rank  : The initial rank of A and B to start. Default is full
%                   rank, calculated from the SVD of Y. A smaller value is
%                   helpful if the data dimensions is large.
%   X_true,E_true : Original X and E matrices for simulations.
%                   If supplied and verbose=1,
%                   the errors are reported at each iteration
% -----------------------------------------------------------------------
%                                OUTPUTS
% X               : low rank component
% A, B            : low rank factors of X = AB'
% optional outputs (see the paper for the full definitions):
%   varargout{1} = gamma    --- hyperparameters for A and B
%   varargout{2} = beta     --- noise inverse variance
%   varargout{3} = it       --- number of iterations
%   varargout{4} = t        --- running time
%   varargout{5} = Sigma_A  --- covariance matrix of A
%   varargout{6} = Sigma_B  --- covariance matrix of B
% -----------------------------------------------------------------------
%
% Copyright (c): S. Derin Babacan, 2011
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License <http://www.gnu.org/licenses/> for more details.
%

%% Check input variables,
if ~exist('options','var')
    options = populate_vars([]);
else
    options = populate_vars(options);
end

verbose      = options.verbose;
a_gamma0     = options.a_gamma0;
b_gamma0     = options.b_gamma0;
a_beta0      = options.a_beta0;
b_beta0      = options.b_beta0;
MAXITER      = options.MAXITER;
inf_flag     = options.inf_flag;
UPDATE_BETA  = options.UPDATE_BETA;
DIMRED       = options.DIMRED;

if DIMRED,
    if isfield(options, 'DIMRED_THR')
        DIMRED_THR = options.DIMRED_THR;
    else
        DIMRED_THR = 1e4;
    end
end

UPDATE_BETA_START = options.UPDATE_BETA_START;

% For synthetic simulations
if isfield(options, 'X_true')
    X_true = options.X_true;
end

%% Initialization
[m n] = size(Y);
mn = m*n;
pmn = length(find(P == 1));
p = pmn/mn;

Y2sum = sum(Y(:).^2);
scale2 = Y2sum / (mn);
scale = sqrt(scale2);

% Initialize A and B
switch options.init,

    % Maximum likelihood
    case 'ml'
        [U, S, V] = svd(Y, 'econ');

        if strcmp(options.initial_rank, 'auto')
            r = min([m,n]);
        else
            r = options.initial_rank;
        end

        A = U(:,1:r)*(S(1:r,1:r)).^(0.5);
        B = (S(1:r,1:r)).^(0.5)*V(:,1:r)'; B = B';

        Sigma_A = repmat( scale*eye(r,r), [1 1 m] );
        Sigma_B = repmat( scale*eye(r,r), [1 1 n] );

        gammas = 1./diag(S.^2);
        gammas = gammas(1:r);
        gammas = ones(r,1)*scale;

        ssscale2 = 1;
        beta = ssscale2./scale2;

        if UPDATE_BETA == 0 & isfield(options, 'beta')
            beta = options.beta;
        end

    % Random initialization
    case 'rand'

        if strcmp(options.initial_rank, 'auto')
            r = min([m,n]);
        else
            r = options.initial_rank;
        end

        A = randn(m,r) * sqrt(scale);
        B = randn(n,r) * sqrt(scale);
        gammas = scale*ones(r,1);

        Sigma_A = repmat( scale*eye(r,r), [1 1 m] );
        Sigma_B = repmat( scale*eye(r,r), [1 1 n] );

        if UPDATE_BETA == 0 & isfield(options, 'beta')
            beta = options.beta;
        else
            beta = 1./scale2;
        end

end

X = A*B';

%% Iterations
tic
for it=1:MAXITER,

    old_X = X;

    Aw = diag(gammas);

    %% A step
    diagsA = 0;
    for i=1:m, %iterate over rows

        observed = find(P(i,:));
        Bi = B(observed,:);
        Sigma_A(:,:,i) = (beta*Bi'*Bi + beta*sum(Sigma_B(:,:,observed),3) + Aw)^(-1);
        A(i,:) = beta*Y(i,observed)*Bi*Sigma_A(:,:,i);
        diagsA = diagsA + diag( Sigma_A(:,:,i) );

    end

    %% B step
    diagsB = 0;
    for j=1:n, %Iterate over cols

        observed = find(P(:,j));
        Aj = A(observed,:);
        Sigma_B(:,:,j) = (beta*Aj'*Aj + beta*sum(Sigma_A(:,:,observed),3) + Aw)^(-1);
        B(j,:) = beta*Y(observed,j)'*Aj*Sigma_B(:,:,j);
        diagsB = diagsB + diag( Sigma_B(:,:,j) );

    end

    %% update X
    X = A*B';

    %% estimate gammas
    % Choose inference method (fixed point Mackay or standard)
    if inf_flag == 1, %Standard
        gammas = (m + n + a_gamma0)./( diag(B'*B) + diag(sum(Sigma_B,3)) + diag(A'*A)+ diag(sum(Sigma_A,3))+ b_gamma0);
    elseif inf_flag == 2, %Mackay
        gammas = (m + n + a_gamma0 - gammas.*( diag(sum(Sigma_A,3)) + diag(sum(Sigma_B,3)) ))./( diag(A'*A) + diag(B'*B)+b_gamma0);
    end

    %% estimate beta
    if  UPDATE_BETA & it>=UPDATE_BETA_START,

        err = sum(sum( abs(Y - P.*X).^2 ) );

        for l = 1: m
            observed = find(P(l, :));
            err = err + trace(B(observed,:)'*B(observed,:)*Sigma_A(:, :, l)) ...
                + trace(A(l, :)'*A(l, :) *sum(Sigma_B(:, :, observed), 3)) ...
                + trace( Sigma_A(:, :, l)*sum(Sigma_B(:, :, observed), 3));

        end

        beta = (pmn + a_beta0)/(err+b_beta0);
    end

    %% Prune irrelevant dimensions?
    if DIMRED,
        MAX_GAMMA = min(gammas) * DIMRED_THR;

        if sum(find(gammas > MAX_GAMMA)),

            indices = find(gammas <= MAX_GAMMA);

            A = A(:,indices);
            B = B(:,indices);
            gammas = gammas(indices);

            Sigma_A = Sigma_A(indices,indices,:);
            Sigma_B = Sigma_B(indices,indices,:);

            [m r] = size(A);
            [n r] = size(B);
        end

    end

    %% Check convergence and display progress
    conv = sum(sum( ( old_X - X ).^2 ) ) / mn;

    if verbose,

        if exist('X_true','var')
            % For synthetic simulations
            err = sum(sum( ( X - X_true ).^2 ) ) / sum(sum( ( X_true ).^2 ) );
            rms = sqrt( sum(sum( ( X - X_true ).^2 ) )/mn );
            relrec = norm(X_true-X,'fro')/norm(X_true,'fro');

            fprintf('it %d: Error = %g, beta = %g, conv = %g, r = %d\n', it, relrec, beta, conv, r);

        else

            fprintf('it %d: beta = %g, conv = %g, r = %d\n', it, beta, conv, r);

        end
    end

    % Check for convergence (do at least 5 iterations)
    if it > 5 & conv < options.thr
        break;
    end

end
t = toc;

varargout{1} = gammas;
varargout{2} = beta;
varargout{3} = it;
varargout{4} = t;
varargout{5} = Sigma_A;
varargout{6} = Sigma_B;

%% Function to populate options.
function [options] = populate_vars(options)

if isempty(options)
    options.init         = 'ml';
    options.verbose      = 1;
    options.a_gamma0     = 1e-6;
    options.b_gamma0     = 1e-6;
    options.a_beta0      = 0;
    options.b_beta0      = 0;
    options.MAXITER      = 100;
    options.DIMRED       = 1;
    options.UPDATE_BETA  = 1;
    options.mode         = 'VB';
    options.thr          = 1e-7;
    options.initial_rank = 'auto';
    options.inf_flag     = 1;
    options.UPDATE_BETA_START = 1;
else

    if ~isfield(options, 'init')
        options.init = 'ml';
    end

    if ~strcmp(options.init,'ml') | ~strcmp(options.init,'rand')
        options.init = 'ml';
    end

    if ~isfield(options, 'verbose')
        options.verbose   = 1;
    end

    if ~isfield(options, 'a_gamma0')
        options.a_gamma0   = 1e-6;
    end

    if ~isfield(options, 'b_gamma0')
        options.b_gamma0   = 1e-6;
    end

    if ~isfield(options, 'a_beta0')
        options.a_beta0   = 0;
    end

    if ~isfield(options, 'b_beta0')
        options.b_beta0   = 0;
    end

    if ~isfield(options, 'MAXITER')
        options.MAXITER   = 100;
    end

    if ~isfield(options, 'DIMRED')
        options.DIMRED   = 1;
    end

    if ~isfield(options, 'inf_flag')
        options.inf_flag   = 1;
    end

    if options.inf_flag ~= 1 & options.inf_flag ~= 2,
        options.inf_flag = 1; % Standard inference
    end

    if ~isfield(options, 'UPDATE_BETA')
        options.UPDATE_BETA = 1;
    end

    if ~isfield(options, 'UPDATE_BETA_START')
        options.UPDATE_BETA_START = 1;
    end

    if ~isfield(options, 'mode')
        options.mode = 'VB';
    end

    if ~isfield(options, 'thr')
        options.thr = 1e-7;
    end

    if ~isfield(options, 'initial_rank')
        options.initial_rank = 'auto';
    end

end

author: andrewssobral

repo: lrslibrary


text:
function grdwrite2(x,y,z,file);
%GRDWRITE2  Write a GMT grid file
%
% Uses built-in NetCDF capability (MATLAB R2008b or later) to
% write a COARDS-compliant netCDF grid file
% Duplicates (some) functionality of the program grdwrite (which requires
% compilation as a mexfile-based function on each architecture) using
% Matlab 2008b (and later) built-in NetCDF functionality
% instead of GMT libraries.
%
% GRDWRITE2(X,Y,Z,'filename') will create a grid file containing the
% data in the matrix Z.  X and Y should be either vectors with
% dimensions that match the size of Z or two-component vectors
% containing the max and min values for each.
%
% See also GRDREAD2, GRDINFO2

% For more information on GMT grid file formats, see:
% http://www.soest.hawaii.edu/gmt/gmt/doc/gmt/html/GMT_Docs/node70.html
% Details on Matlab's native netCDF capabilities are at:
% http://www.mathworks.com/access/helpdesk/help/techdoc/ref/netcdf.html

% GMT (Generic Mapping Tools, <http://gmt.soest.hawaii.edu>)
% was developed by Paul Wessel and Walter H. F. Smith

% Kelsey Jordahl
% Marymount Manhattan College
% http://marymount.mmm.edu/faculty/kjordahl/software.html

% Time-stamp: <Tue Jul 19 16:28:24 EDT 2011>

% Version 1.1.2, 19-Jul-2011
% Available at MATLAB Central
% <http://www.mathworks.com/matlabcentral/fileexchange/26290-grdwrite2>

if nargin < 4,
  help(mfilename);
  return,
end

% check for appropriate Matlab version (>=7.7)
V=regexp(version,'[ \.]','split');
if (str2num(V{1})<7) | (str2num(V{1})==7 & str2num(V{2})<7),
  ver
  error('grdread2: Requires Matlab R2008b or later!');
end

ncid = netcdf.create(file, 'NC_SHARE');
if isempty(ncid),
  return,
end

% set descriptive variables
conv='COARDS/CF-1.0';
title=file;
history='File written by MATLAB function grdwrite2.m';
desc=['Created ' datestr(now)];
vers='4.x';                             % is "x" OK?

% check X and Y
if (~isvector(x) | ~isvector(y)),
  error('X and Y must be vectors!');
end
if (length(x) ~= size(z,2)),    % number of columns don't match size of x
  minx=min(x); maxx=max(x);
  dx=(maxx-minx)/(size(z,2)-1);
  x=minx:dx:maxx;                       % write as a vector
end
if (length(y) ~= size(z,1)),    % number of rows don't match size of y
  miny=min(y); maxy=max(y);
  dy=(maxy-miny)/(size(z,1)-1);
  y=miny:dy:maxy;                       % write as a vector
end

% match Matlab class to NetCDF data type
switch class(z)
 case 'single'
  nctype='NC_FLOAT';
  nanfill=single(NaN);
 case 'double'
  nctype='NC_DOUBLE';
  nanfill=double(NaN);
 case 'int8'
  nctype='NC_BYTE';
  nanfill=intmin(class(z));
  disp(['Warning: ''No data'' fill value set to ' num2str(nanfill)])
 case 'int16'
  nctype='NC_SHORT';
  nanfill=intmin(class(z));
  disp(['Warning: ''No data'' fill value set to ' num2str(nanfill)])
 case 'int32'
  nctype='NC_INT';
  nanfill=intmin(class(z));
  disp(['Warning: ''No data'' fill value set to ' num2str(nanfill)])
 otherwise
  error(['Don''t know how to handle data of class ''' class(z) '''.  Try converting to a supported data type (int8, int16, int32, single or double).'])
end

% global
netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Conventions',conv);
netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'title',title);
netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'history',history);
netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'description',desc);
netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'GMT_version',vers);
% X
dimid = netcdf.defDim(ncid,'x',length(x));
varid = netcdf.defVar(ncid,'x','double',dimid);
netcdf.putAtt(ncid,varid,'long_name','x');
netcdf.putAtt(ncid,varid,'actual_range',[min(x) max(x)]);
netcdf.endDef(ncid);
netcdf.putVar(ncid,varid,x);
% Y
netcdf.reDef(ncid);
dimid = netcdf.defDim(ncid,'y',length(y));
varid = netcdf.defVar(ncid,'y','double',dimid);
netcdf.putAtt(ncid,varid,'long_name','y');
netcdf.putAtt(ncid,varid,'actual_range',[min(y) max(y)]);
netcdf.endDef(ncid);
netcdf.putVar(ncid,varid,y);
% Z
netcdf.reDef(ncid);
varid = netcdf.defVar(ncid,'z',nctype,[0 1]);
netcdf.putAtt(ncid,varid,'long_name','z');
netcdf.putAtt(ncid,varid,'_FillValue',nanfill);
netcdf.putAtt(ncid,varid,'actual_range',[min(z(:)) max(z(:))]);
netcdf.endDef(ncid);
netcdf.putVar(ncid,varid,z');
% close file
netcdf.close(ncid);

author: dbekaert

repo: TRAIN

save_path: github-repos/MATLAB/dbekaert-TRAIN


text:
function contourVisibility()
%Make contours whose visible flag is off hidden.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC stateS

indexS = planC{end};

for i=1:length(planC{indexS.structures})
    if ~isfield(planC{indexS.structures}, 'visible') || isempty(planC{indexS.structures}(i).visible)
        planC{indexS.structures}(i).visible = 1;
    end
end
vis = [planC{indexS.structures}.visible];

sGv = [];
for i = uint8(1:length(stateS.handle.CERRAxis))
    sG = getAxisInfo(i,'structureGroup');
    sGv = [sGv; sG.handles];
end

contours = findobj(sGv,'tag', 'structContour');

ud = get(contours, 'userdata');
if iscell(ud)
    ud = [ud{:}];
end
structNum = [ud.structNum];
tf = ismember(structNum, find(vis));
set(contours(~tf), 'visible', 'off');
set(contours(tf), 'visible', 'on');

if stateS.optS.structureDots
    contourDots = findobj(sGv,'tag', 'structContourDots');
    udCD = get(contourDots, 'userdata');
    udCD = [udCD{:}];
    structNumCD = udCD;
    tf = ismember(structNumCD, find(vis));
    set(contourDots(~tf), 'visible', 'off');
    set(contourDots(tf), 'visible', 'on');
end
return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [M, inliers] = RANSAC(x, fittingfn, distfn, degenfn, s, t)

    maxTrials = 1000;
    maxDataTrials = 100;

    [rows, npts] = size(x);

    % Desired probability of choosing at least one sample free from outliers (probably should be a parameter)
    p = 0.99;

    bestM = NaN;
    trialcount = 0;
    bestscore =  0;
    N = 1;

    while N > trialcount

        % Select at random s datapoints to form a trial model, M.
        degenerate = 1;
        count = 1;
        while degenerate
            % Generate s random indicies in the range 1..npts
            ind = randsample(npts, s);
            % Su Tan
%             while tooDense(x(:, ind), 10)
%                 ind = randsample(npts, s);
%             end
            % Test that these points are not a degenerate configuration.
            degenerate = feval(degenfn, x(:,ind));

            if ~degenerate
                M = feval(fittingfn, x(:,ind));
                if isempty(M)
                    degenerate = 1;
                end
            end

            % Safeguard against being stuck in this loop forever
            count = count + 1;
            if count > maxDataTrials
                disp('Unable to select a nondegenerate data set');
                break
            end
        end

%         if scoreX1 ~= 0
%             l = min([scoreX1(ind) scoreX2(ind)]);
%             t = (l + 0.6) * t;
%         end

        [inliers, score, M] = feval(distfn, M, x, t);

        % Find the number of inliers to this model.
        ninliers = length(inliers);

        if score > bestscore
            bestscore = score;
            bestinliers = inliers;
            bestM = M;

            % Update estimate of N
            fracinliers =  ninliers/npts;
            pNoOutliers = 1 -  fracinliers^s;
            pNoOutliers = max(eps, pNoOutliers);
            pNoOutliers = min(1-eps, pNoOutliers);
            N = log(1-p)/log(pNoOutliers);
        end

        trialcount = trialcount+1;

        % Safeguard against being stuck in this loop forever
        if trialcount > maxTrials
            fprintf('ransac reached the maximum number of %d trials\n', maxTrials);
            break
        end

        %fprintf('Iter: %d \n', trialcount);
    end

    if ~isnan(bestM)
        M = bestM;
        inliers = bestinliers;
    else
        M = [];
        inliers = [];
        disp('ransac was unable to find a useful solution');
    end
end

author: SuTanTank

repo: VideoStitchingViaShakinessRemoving


text:
function visualizeBoundary(X, y, model, varargin)
%VISUALIZEBOUNDARY plots a non-linear decision boundary learned by the SVM
%   VISUALIZEBOUNDARYLINEAR(X, y, model) plots a non-linear decision
%   boundary learned by the SVM and overlays the data on it

% Plot the training data on top of the boundary
plotData(X, y)

% Make classification predictions over a grid of values
x1plot = linspace(min(X(:,1)), max(X(:,1)), 100)';
x2plot = linspace(min(X(:,2)), max(X(:,2)), 100)';
[X1, X2] = meshgrid(x1plot, x2plot);
vals = zeros(size(X1));
for i = 1:size(X1, 2)
   this_X = [X1(:, i), X2(:, i)];
   vals(:, i) = svmPredict(model, this_X);
end

% Plot the SVM boundary
hold on
contour(X1, X2, vals, [1 1], 'b');
hold off;

end

author: atinesh-s

repo: Coursera-Machine-Learning-Stanford