text:
clear; close all; clc;

%%
% Plot settings

markerSize=10;
markerSize2=10;
fontSize=15;
faceAlpha1=0.3;
faceAlpha2=0.5;
cMap=gjet(4);

%% 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=[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_sed=[febioFebFileNamePart,'_stress_out.txt']; %Log file name for exporting strain energy density

tireWidth = 120;
tireRadius = 150;
tireThickness = 8;
threadHeight = 8;
tireBulge=tireWidth/8;
pointSpacing = 8;
tireCircumference = 2*pi.*tireRadius;

numFeatures = 20;
groveWidth = 8;

distKeep=75; %Distance after which model is cropped, empty = full model
displacementVert=0.8*tireThickness;
displacementHor=-tireRadius/20;
contactInitialOffset=0.1;

%Material parameter set
c1=1e-3; %Shear-modulus-like parameter
m1=2; %Material parameter setting degree of non-linearity
k_factor=100; %Bulk modulus factor
k=c1*k_factor; %Bulk modulus

% FEA control settings
numTimeSteps=20; %Number of time steps desired
max_refs=40; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=25; %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
runMode='external';%'internal';
min_residual=1e-20;

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

%%

tireHalfWidth=tireWidth/2;
featureHeight = (tireCircumference-(numFeatures*groveWidth))./numFeatures;
featureWidth = tireHalfWidth-2*groveWidth;

% Coordinates of 1x1 square
v=[-1 -1 0;...
    1 -1 0;...
    1  1 0;...
    -1  1 0]/2;

% Deformation gradient tensor for shear
f=eye(3,3);
f(1,2)=1;

%Create boundary coordinates
V1=v; %Copy square
V1(:,1)=V1(:,1).*tireCircumference+tireCircumference/2; %Stretch to create circumference
V1(:,2)=V1(:,2).*tireHalfWidth+tireHalfWidth/2; %Stretch to get desired width
V1=evenlySpaceCurve(V1,pointSpacing,'linear',1,1:1:size(v,1)); %Evenly sample
V1=(f*V1')'; %Shear

%Create thread feature
Vf=v; %Copy square
Vf(:,1)=Vf(:,1).*featureHeight; %Stretch to create circumference
Vf(:,2)=Vf(:,2).*(featureWidth+groveWidth/2)+tireHalfWidth/2-groveWidth/4; %Stretch to get desired width
Vf=evenlySpaceCurve(Vf,pointSpacing,'linear',1,1:1:size(v,1)); %Evenly sample
Vf(:,1)=Vf(:,1)+featureHeight/2+groveWidth/2;

VF=cell(1,numFeatures);
for q=1:1:numFeatures
    Vff=Vf;
    Vff(:,1)=Vf(:,1)+(q-1).*(featureHeight+groveWidth); %Shift
    VF{q}=(f*Vff')'; %Shear and store
end

%%

%Defining 4 regions
regionSpec{1}={V1(:,[1 2])};
for q=1:1:numFeatures
    regionSpec{1}{q+1}=VF{q}(:,[1 2]);
end

for q=1:1:numFeatures
    regionSpec{q+1}{1}=VF{q}(:,[1 2]);
end

[Ft,Vt,Ct]=multiRegionTriMesh2D(regionSpec,pointSpacing,0,0);
Vt(:,3)=0;

Fb=fliplr(Ft);
Vb=Vt;
Vb(:,2)=-Vb(:,2);

Cb=Ct;
Cb(Cb>1)=Cb(Cb>1)-1+max(Ct);

[F,V,C]=joinElementSets({Ft,Fb},{Vt,Vb},{Ct,Cb});
[F,V]=mergeVertices(F,V);

%%

cFigure; hold on;
plotV(V1,'r.-','MarkerSize',markerSize);
plotV(VF,'b.-','MarkerSize',markerSize);
gpatch(F,V,C);

axisGeom(gca,fontSize); view(2);
colormap spectral; icolorbar;
gdrawnow;

%%

cMax=max(C);
numVert=size(V,1);
F2=F(C>1,:);
C2=C(C>1,:)+cMax;
Eb1=patchBoundary(F2);

Eb2=Eb1+numVert;
F=[F(C==1,:);F2+numVert];
V=[V;V];
C=[C(C==1,:);C2+max(C(:))];

%%
%

if tireBulge>eps(0)
    w=tireWidth;
    h=tireBulge;
    r=h/2+(w^2/(8*h));
    a=atan((w/2)./(r-h));

    TH = (V(:,2)./tireHalfWidth)*a;
    R = r.*ones(size(V,1),1);
    R(numVert+1:end)=R(numVert+1:end)+threadHeight;
    Z = V(:,1);
    [V(:,1),V(:,2),V(:,3)] = pol2cart(TH,R,Z);
    V=V(:,[3 2 1]);

    V(:,3)=V(:,3)-r+h;
else
    V(numVert+1:end,3)=V(numVert+1:end,3)+threadHeight;
end

%%
cFigure; hold on;
plotV(V1,'r.-','MarkerSize',markerSize);
plotV(VF,'b.-','MarkerSize',markerSize);
gpatch(F,V,C);

axisGeom(gca,fontSize); view(2);
colormap spectral; icolorbar;
gdrawnow;

%%

TH = V(:,1)./tireRadius;
R = tireRadius+V(:,3);

Z = V(:,2);
[V(:,1),V(:,2),V(:,3)] = pol2cart(TH,R,Z);
V=V(:,[1 3 2]); F=fliplr(F);

[Eb1,V1]=patchCleanUnused(Eb1,V);
[Eb2,V2]=patchCleanUnused(Eb2,V);

[F,V,~,indFix]=mergeVertices(F,V);

numSteps=ceil(threadHeight./pointSpacing)+1;

%Create quadrilateral faces
X=linspacen(V1(:,1),V2(:,1),numSteps);
Y=linspacen(V1(:,2),V2(:,2),numSteps);
Z=linspacen(V1(:,3),V2(:,3),numSteps);

Vn=[X(:) Y(:) Z(:)]; %Vertex set
Fn=repmat([fliplr(Eb1) Eb2],numSteps-1,1);
q=repmat(0:(numSteps-2),size(Eb1,1),1);
q=q(:);
Q=[q q q+1 q+1]*size(V1,1);
Fn=Fn+Q;

[Fn,Vn]=quad2tri(Fn,Vn,'a');

%%

[F,V]=patchCleanUnused(F,V);
C(C>1)=2; %Switch to use 2 for all thread profile parts
[F_out,V_out,C_out]=joinElementSets({F,Fn},{V,Vn},{C,max(C(:))+1*ones(size(Fn,1),1)});
[F_out,V_out]=patchCleanUnused(F_out,V_out);
[F_out,V_out]=mergeVertices(F_out,V_out);
Eb_out=patchBoundary(F_out);
%%

cFigure; hold on;
gpatch(F_out,V_out,C_out);
% patchNormPlot(Fq,Vq);
axisGeom(gca,fontSize); camlight headlight;
colormap spectral; icolorbar;
gdrawnow;

%%
% Create inner cylinder
inputStructCyl.cylRadius=tireRadius-tireThickness;
inputStructCyl.numRadial=ceil((2*pi*inputStructCyl.cylRadius)/pointSpacing);
inputStructCyl.cylHeight=tireWidth;
numheight=ceil(tireWidth./pointSpacing);
numheight=numheight+iseven(numheight); %Force uneven
inputStructCyl.numHeight=numheight;
inputStructCyl.meshType='tri';

%%
% Derive patch data for a cylinder
[F_in,V_in]=patchcylinder(inputStructCyl);
V_in=V_in(:,[1 3 2]);
Eb_in=patchBoundary(F_in);
C_in=(max(C_out)+1)*ones(size(F_in,1),1);

%%
inputStructGroup.outputType='label';
G_out=tesgroup(Eb_out,inputStructGroup);
G_in=tesgroup(Eb_in,inputStructGroup);

yMean_out=nan(2,1);
yMean_in=nan(2,1);
for q=1:1:2
    ind_out=Eb_out(G_out==q,:);
    yMean_out(q)=mean(V_out(ind_out(:),2));

    ind_in=Eb_in(G_in==q,:);
    yMean_in(q)=mean(V_in(ind_in(:),2));
end

[~,ind1_out]=min(yMean_out);
[~,ind2_out]=max(yMean_out);
indCurve1_out=edgeListToCurve(Eb_out(G_out==ind1_out,:));
indCurve1_out=indCurve1_out(1:end-1);
indCurve2_out=edgeListToCurve(Eb_out(G_out==ind2_out,:));
indCurve2_out=indCurve2_out(1:end-1);

[~,ind1_in]=min(yMean_in);
[~,ind2_in]=max(yMean_in);
indCurve1_in=edgeListToCurve(Eb_in(G_in==ind1_in,:));
indCurve1_in=indCurve1_in(1:end-1);
indCurve2_in=edgeListToCurve(Eb_in(G_in==ind2_in,:));
indCurve2_in=indCurve2_in(1:end-1);

[Fs1,Vs1]=regionTriMesh2D({V_out(indCurve1_out,[1 3]),V_in(indCurve1_in,[1 3])},pointSpacing,0);
Vs1(:,3)=-tireWidth/2;
Vs1=Vs1(:,[1 3 2]);
Cs1=(max(C_in)+1)*ones(size(Fs1,1),1);

[Fs2,Vs2]=regionTriMesh2D({V_out(indCurve2_out,[1 3]),V_in(indCurve2_in,[1 3])},pointSpacing,0);
Vs2(:,3)=tireWidth/2;
Vs2=Vs2(:,[1 3 2]);
Fs2=fliplr(Fs2);
Cs2=(max(Cs1)+1)*ones(size(Fs2,1),1);

%%

[FT,VT,CT]=joinElementSets({F_out,F_in,Fs1,Fs2},{V_out,V_in,Vs1,Vs2},{C_out,C_in,Cs1,Cs2});
[FT,VT]=patchCleanUnused(FT,VT);
[FT,VT]=mergeVertices(FT,VT);

%%

cFigure; hold on;

gpatch(FT,VT,CT,'none');
patchNormPlot(FT,VT);

axisGeom(gca,fontSize); camlight headlight;
colormap spectral; icolorbar;
gdrawnow;

%%

V_regions=getInnerPoint(FT,VT); %Define region points
V_holes=[]; %Define hole points
[regionTetVolumes]=tetVolMeanEst(FT,VT); %Volume estimate for regular tets
stringOpt='-pq1.2AaY'; %Options for tetgen

%%
% Mesh using TetGen

%Create tetgen input structure
inputStruct.stringOpt=stringOpt; %Tetgen options
inputStruct.Faces=FT; %Boundary faces
inputStruct.Nodes=VT; %Nodes of boundary
inputStruct.faceBoundaryMarker=CT;
inputStruct.regionPoints=V_regions; %Interior points for regions
inputStruct.holePoints=V_holes; %Interior points for holes
inputStruct.regionA=regionTetVolumes; %Desired tetrahedral volume for each region

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

%%
% Access mesh output structure

E=meshOutput.elements; %The elements
V=meshOutput.nodes; %The vertices or nodes
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(V_regions,'r.','MarkerSize',markerSize);
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);

axisGeom(gca,fontSize);
gdrawnow;

%%

if ~isempty(distKeep)
    Fb_full=Fb;
    V_full=V;
    [Fb_full,V_full]=patchCleanUnused(Fb_full,V_full);

    VE=patchCentre(E,V);
    logicKeepElements = VE(:,3) <= (min(VE(:,3))+distKeep);
    E=E(logicKeepElements,:);
    indKeepNodes=unique(E);
    logicKeep_Fb=all(ismember(Fb,indKeepNodes),2);
    Fb=Fb(logicKeep_Fb,:);
    Cb=Cb(logicKeep_Fb,:);

    logicKeep_F=all(ismember(F,indKeepNodes),2);
    F=F(logicKeep_F,:);

    [E,V,indFix]=patchCleanUnused(E,V);
    Fb=indFix(Fb);
    F=indFix(F);
    meshOutput.nodes=V;
    meshOutput.facesBoundary=Fb;
    meshOutput.boundaryMarker=Cb;
    meshOutput.faces=F;
    meshOutput.elements=E;
    meshOutput.elementMaterialID=ones(size(E,1),1);
else
    Fb_full=Fb;
    V_full=V;
end

%%

hf=cFigure;
subplot(1,2,1); hold on;
title('Input boundaries','FontSize',fontSize);
hp(1)=gpatch(Fb,V,Cb,'k',faceAlpha1);
hp(2)=plotV(V_regions,'r.','MarkerSize',markerSize);
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);

axisGeom(gca,fontSize);
gdrawnow;

%%

plateDim=[2*tireRadius tireWidth];
plateEl=ceil(plateDim./pointSpacing);
[Eg,Vg]=quadPlate(plateDim,plateEl);
Vg(:,3)=Vg(:,3)+min(VT(:,3))-contactInitialOffset;
center_of_mass=mean(Vg,1);

%% Joining node sets
Eg=Eg+size(V,1); %Fixed element indices
V=[V;Vg;]; %Combined node sets

%%

cFigure; hold on;
hp(1)=gpatch(Fb,V,Cb,'k',faceAlpha1);
hp(2)=gpatch(Eg,V,'kw','k',faceAlpha1);
% legend(hp,{'Input mesh','Interior point(s)'},'Location','NorthWestOutside');
axisGeom(gca,fontSize); camlight headlight;
colormap(cMap); icolorbar;
gdrawnow;

%% Define contact surfaces

% The rigid master surface of the sphere
F_contact_secondary=Eg;

% The deformable slave surface of the slab
logicContactSurf1=Cb==2;
F_contact_primary=fliplr(Fb(logicContactSurf1,:));

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

gpatch(Fb,V,'kw','none',faceAlpha2);
hl(1)=gpatch(F_contact_secondary,V,'g','k',1);
patchNormPlot(F_contact_secondary,V);
hl(2)=gpatch(F_contact_primary,V,'b','k',1);
patchNormPlot(F_contact_primary,V);

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

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

%% Define boundary conditions

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

%%
% Visualize BC's
hf=cFigure;
title('Boundary conditions model','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;

gpatch(Fb,V,'w','none',faceAlpha2);

hl2(1)=gpatch(Eg,V,'gw','k',1);
hl2(2)=plotV(V(bcSupportList,:),'k.','MarkerSize',markerSize);

legend(hl2,{'Rigid body plate','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';

%Create control structure for use by all steps
febio_spec.Control.analysis='STATIC';
stepStruct.Control.time_steps=numTimeSteps;
stepStruct.Control.step_size=1/numTimeSteps;
stepStruct.Control.solver.max_refs=max_refs;
stepStruct.Control.solver.max_ups=max_ups;
stepStruct.Control.solver.symmetric_stiffness=0;
febio_spec.Control.solver.min_residual=min_residual;
stepStruct.Control.time_stepper.dtmin=dtmin;
stepStruct.Control.time_stepper.dtmax=dtmax;
stepStruct.Control.time_stepper.max_retries=max_retries;
stepStruct.Control.time_stepper.opt_iter=opt_iter;

%Add template based default settings to proposed control section
[stepStruct.Control]=structComplete(stepStruct.Control,febio_spec.Control,1); %Complement provided with default if missing

%Remove control field (part of template) since step specific control sections are used
febio_spec=rmfield(febio_spec,'Control');

febio_spec.Step.step{1}.Control=stepStruct.Control;
febio_spec.Step.step{1}.ATTR.id=1;
febio_spec.Step.step{2}.Control=stepStruct.Control;
febio_spec.Step.step{2}.ATTR.id=2;

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

materialName2='Material2';
febio_spec.Material.material{2}.ATTR.name=materialName2;
febio_spec.Material.material{2}.ATTR.type='rigid body';
febio_spec.Material.material{2}.ATTR.id=2;
febio_spec.Material.material{2}.density=1;
febio_spec.Material.material{2}.center_of_mass=center_of_mass;

%Mesh section
% -> Nodes
febio_spec.Mesh.Nodes{1}.ATTR.name='nodeSet_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,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E; %The element matrix

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

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

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

febio_spec.MeshDomains.ShellDomain.ATTR.name=partName2;
febio_spec.MeshDomains.ShellDomain.ATTR.mat=materialName2;

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

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

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

%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';

%Rigid section
% -> Prescribed rigid body boundary conditions
febio_spec.Step.step{1}.Rigid.rigid_constraint{1}.ATTR.name='RigidFix_1';
febio_spec.Step.step{1}.Rigid.rigid_constraint{1}.ATTR.type='fix';
febio_spec.Step.step{1}.Rigid.rigid_constraint{1}.rb=2;
febio_spec.Step.step{1}.Rigid.rigid_constraint{1}.dofs='Rx,Ry';

febio_spec.Step.step{1}.Rigid.rigid_constraint{2}.ATTR.name='RigidPrescribe';
febio_spec.Step.step{1}.Rigid.rigid_constraint{2}.ATTR.type='prescribe';
febio_spec.Step.step{1}.Rigid.rigid_constraint{2}.rb=2;
febio_spec.Step.step{1}.Rigid.rigid_constraint{2}.dof='Rz';
febio_spec.Step.step{1}.Rigid.rigid_constraint{2}.value.ATTR.lc=1;
febio_spec.Step.step{1}.Rigid.rigid_constraint{2}.value.VAL=(displacementVert+contactInitialOffset);
febio_spec.Step.step{1}.Rigid.rigid_constraint{2}.relative=0;

febio_spec.Step.step{2}.Rigid.rigid_constraint{1}.ATTR.name='RigidFix_1';
febio_spec.Step.step{2}.Rigid.rigid_constraint{1}.ATTR.type='fix';
febio_spec.Step.step{2}.Rigid.rigid_constraint{1}.rb=2;
febio_spec.Step.step{2}.Rigid.rigid_constraint{1}.dofs='Ry,Rz,Ru,Rv,Rw';

febio_spec.Step.step{2}.Rigid.rigid_constraint{2}.ATTR.name='RigidPrescribe';
febio_spec.Step.step{2}.Rigid.rigid_constraint{2}.ATTR.type='prescribe';
febio_spec.Step.step{2}.Rigid.rigid_constraint{2}.rb=2;
febio_spec.Step.step{2}.Rigid.rigid_constraint{2}.dof='Rx';
febio_spec.Step.step{2}.Rigid.rigid_constraint{2}.value.ATTR.lc=2;
febio_spec.Step.step{2}.Rigid.rigid_constraint{2}.value.VAL=displacementHor;
% febio_spec.Step.step{2}.Rigid.rigid_constraint{2}.relative=1;

%Contact section
febio_spec.Contact.contact{1}.ATTR.type='sliding-elastic';
febio_spec.Contact.contact{1}.ATTR.surface_pair=febio_spec.Mesh.SurfacePair{1}.ATTR.name;
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; 2 1];

febio_spec.LoadData.load_controller{2}.ATTR.id=2;
febio_spec.LoadData.load_controller{2}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{2}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{2}.points.point.VAL=[0 0; 1 0; 2 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_sed;
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,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.runMode=runMode;

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

%% Import FEBio results

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

    %%
    % Importing nodal displacements from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_disp),1,1);

    %Access data
    N_disp_mat=dataStruct.data; %Displacement
    timeVec=dataStruct.time; %Time

    %Create deformed coordinate set
    V_DEF=N_disp_mat+repmat(V,[1 1 size(N_disp_mat,3)]);

    %%
    % Importing element stress from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_sed),1,1);

    %Access data
    E_sed_mat=dataStruct.data;
    E_sed_mat(isnan(E_sed_mat))=0;

        %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    [CV]=faceToVertexMeasure(E,V,E_sed_mat(:,:,end));

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    title('Strain energy density','Interpreter','Latex')
    gpatch(Fb_full,V_full,'w','none',0.25);
    hp=gpatch(Fb,V_DEF(:,:,end),CV,'none',1); %Add graphics object to animate
%     hp.Marker='.';
%     hp.MarkerSize=markerSize2;
    hp.FaceColor='interp';

    hp2=gpatch(Eg,V_DEF(:,:,end),'w','none',0.5); %Add graphics object to animate
    Ebb=patchBoundaryLabelEdges(Fb,V,Cb);
    hp3=gpatch(Ebb,V_DEF(:,:,end),'none','k',1,3);
    hp4=gpatch(patchBoundary(Eg),V_DEF(:,:,end),'none','kw',1,3);

    axisGeom(gca,fontSize);
    colormap(spectral(250)); colorbar;
    caxis([0 0.25*max(E_sed_mat(:))]);
    axis(axisLim(V_DEF,V_full)); %Set axis limits statically
    camlight headlight;

    % Set up animation features
    animStruct.Time=timeVec; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments

        [CV]=faceToVertexMeasure(E,V,E_sed_mat(:,:,qt));

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

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:
function c=context_menus(obj, tag, createmode, varargin)
    % manages context menus, which avoids repeatedly creating context menus
    % c=context_menus(obj, tag, createmode, varargin)
    % Contexts are attached to the figure, so deleting
    % objects that they are attached to will not delete the context menus.
    %
    % context menus can be reused.
    existing_contexts = findobj(obj.fig,'Type','uicontextmenu');
    c = findobj(existing_contexts, 'Tag',tag);

    switch createmode
        case 'overwrite'
            % delete existing context first, then recreate it and return handle
        case 'reuse'
            % if a context exists, just return a handle to it
    end
    switch tag

    end
end

author: CelsoReyes


text:
function [xi, em, gx, gy] = surface_interpolation(x, PARAM, INTERP, res, KLIM, nlev)
% surface_interpolation  Interpolate a surface from a scattered set of points
%
% [XI, EM, GX, GY] = surface_interpolation(X)
%
%   X is a 3-row matrix. Each column has the coordinates of a point that
%   belongs to the surface we want to interpolate.
%
%   XI is a 3-row matrix with the coordinates of the interpolated points.
%
%   EM is a 2-row matrix with the coordinates of the X points projected
%   onto the interpolation domain. Note that the interpolation domain may
%   change between methods, so don't expect results to be aligned if you do
%   e.g. plot3(gx(:), gy(:), y(:, 3), 'o').
%
%   GX, GY are the grid for the box that contains EM.
%
% ... = surface_interpolation(X, PARAM, INTERP, RES, KLIM, NLEV)
%
%   PARAM is a string with the method used to parametrize the surface and
%   X:
%
%     'xy' (default): No change, the X coordinates are kept the same.
%
%     'pca': X points are rotated according to their eigenvectors to make
%     the dominant plane of the points X as horizontal as possible before
%     interpolating.
%
%     'isomap': Use the Isomap method by [1] to "unfold" the curved surface
%     defined by X before interpolating. (This option requires function
%     IsomapII).
%
%   INTERP is a string with the interpolation method:
%
%      'tps' (default): Thin-plate spline. Global support.
%
%      'tsi': Matlab's TriScatteredInterp() function. Local support,
%      limited to the convex hull of the scattered points 2D projection on
%      the interpolation domain.
%
%      'gridfit': John D'Errico's gridfit() function [3] (note:
%      approximation, rather than interpolation). Local support with
%      extrapolation outside the convex hull.
%
%      'mba': Multilevel B-Spline Approximation Library by SINTEF ICT [4].
%      Local support, limited to a rectangle that tighly contains the
%      scattered points 2D projection on the interpolation domain.
%
%      'mbae': Like the 'mba' method, but first a thin-plate spline is used
%      to extrapolate values on the interpolation domain boundary. Then MBA
%      is used for the local support interpolation of X and the boundary
%      values set by the thin-plate spline.
%
%   RES is a 2-vector with the grid spacing in the x- and y-directions. By
%   default, RES=[1 1].
%
%   KLIM is a scalar factor for the extension of the interpolation domain.
%   By default, KLIM=1 and the interpolation domain is a rectangle that
%   tightly contains X. Sections of the interpolated surface that protude
%   from the image volume are removed.
%
%   NLEV is the number of levels in the hierarchical construction of 'mba'
%   and 'mbae'. For other INTERP options, it will be ignored.  By default,
%   NLEV = 7.
%
%
% [1] J.B. Tenenbaum, V. de Silva and J.C. Langford, "A Global Geometric
% Framework for Nonlinear Dimensionality Reduction", Science 290(5500):
% 2319-2323, 2000.
%
% [2] Isomap Homepage, http://isomap.stanford.edu/
%
% [3] Surface Fitting using gridfit by John D'Errico 11 Nov 2005 (Updated
% 29 Jul 2010). Code covered by the BSD License
% http://www.mathworks.com/matlabcentral/fileexchange/8998-surface-fitting-using-gridfit
%
% [4] MBA - Multilevel B-Spline Approximation Library
% http://www.sintef.no/Projectweb/Geometry-Toolkits/MBA/

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2010-2011 University of Oxford
% Version: 0.2.0
% $Rev: 782 $
% $Date: 2012-06-02 22:55:03 +0100 (Sat, 02 Jun 2012) $
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% 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. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
error(nargchk(2, 6, nargin, 'struct'));
error(nargoutchk(0, 4, nargout, 'struct'));

% defaults
if (nargin < 2 || isempty(PARAM))
    PARAM = 'xy';
end
if (nargin < 3 || isempty(INTERP))
    INTERP = 'tps';
end
if (nargin < 4 || isempty(res))
    res = [1 1];
end
if (nargin < 5 || isempty(KLIM))
    KLIM = 1;
end
if (nargin > 5 && ~isempty(nlev) ...
        && ~(strcmp(INTERP, 'mba') || strcmp(INTERP, 'mbae')))
    warning('NLEV input argument ignored for this INTERP option')
end
if (nargin < 6 || isempty(nlev))
    nlev = 7;
end

%% map the 3D points (x,y,z) to a 2D domain (u,v)

% this is analogous to computing the knot vector for a curve interpolation.
% The idea is that (x,y)->z is not necessarily a function, due to the valve
% folding over. However, we hope that (u,v)->(x,y,z) is a function
switch PARAM

    case 'xy'

        % (u,v) is simply (x,y)
        em = x(1:2, :);

    case 'pca'

        % rotate valve points to make the valve surface as horizontal as
        % possible
        m = mean(x, 2);
        em = x - m(:, ones(1, size(x, 2)));
        eigv = pts_pca(em);
        if any(isnan(eigv(:)))
            error('Cannot interpolate surface');
        end
        em = eigv' * em;
        em = em(1:2, :);

    case 'isomap'

        % compute distance matrix
        d = dmatrix(x, x, 'euclidean');

        % compute 2-d projection of the 3-d data
        options.dims = 2;
        options.display = 0;
        options.overlay = 0;
        options.verbose = 0;
        em = IsomapII(d, 'k', round(size(x, 2)/3), options);
        em = em.coords{1};

    otherwise
        error('Parametrization method not implemented')
end

%% compute interpolation domain

% find box that contains embedded coordinates
emmin = min(em, [], 2);
emmax = max(em, [], 2);

% box size and centroid
delta = emmax - emmin;
boxm = mean([emmax emmin], 2);

% extend the box
emmin = boxm - delta/2*KLIM;
emmax = boxm + delta/2*KLIM;

% generate grid for the embedding box
[gy, gx] = ndgrid(emmin(2):res(1):emmax(2), emmin(1):res(2):emmax(1));

%% compute interpolating surface

% source and target points that will define the warp
%s = em; % don't duplicate data in memory
%t = x; % don't duplicate data in memory

% interpolate
switch INTERP
    case 'tps' % thin-plate spline
        xi = pts_tps_map(em', x', [gx(:) gy(:)]);
    case 'tsi' % Matlab's TriScatteredInterp
        fx = TriScatteredInterp(em(1, :)', em(2, :)', x(1, :)', 'natural');
        fy = TriScatteredInterp(em(1, :)', em(2, :)', x(2, :)', 'natural');
        fz = TriScatteredInterp(em(1, :)', em(2, :)', x(3, :)', 'natural');
        fx = fx(gx, gy);
        fy = fy(gx, gy);
        fz = fz(gx, gy);
        xi = [fx(:) fy(:) fz(:)];
    case 'gridfit'
        fx = gridfit(em(1, :)', em(2, :)', x(1, :)', ...
            emmin(1):res(2):emmax(1), emmin(2):res(1):emmax(2), ...
            'tilesize', 150);
        fy = gridfit(em(1, :)', em(2, :)', x(2, :)', ...
            emmin(1):res(2):emmax(1), emmin(2):res(1):emmax(2), ...
            'tilesize', 150);
        fz = gridfit(em(1, :)', em(2, :)', x(3, :)', ...
            emmin(1):res(2):emmax(1), emmin(2):res(1):emmax(2), ...
            'tilesize', 150);
        xi = [fx(:) fy(:) fz(:)];
    case 'mba' % Multilevel B-Spline Approximation Library
        xi = [...
            mba_surface_interpolation(em(1, :)', em(2, :)', x(1, :)', gx(:), gy(:), nlev) ...
            mba_surface_interpolation(em(1, :)', em(2, :)', x(2, :)', gx(:), gy(:), nlev) ...
            mba_surface_interpolation(em(1, :)', em(2, :)', x(3, :)', gx(:), gy(:), nlev)];
    case 'mbae' % MBA extrapolated using a TPS
        % use TPS to interpolate the points, but we are only interested in
        % the edges and corners of the interpolation, where the TPS is
        % extrapolating linearly; also, we decimate the points in the edges
        em2 = [gx(1, 1:end)' gy(1, 1:end)' ; ... % top edge
            gx(2:end, end) gy(2:end, end) ; ... % right edge
            gx(end, 1:end-1)' gy(end, 1:end-1)' ; ... % bottom edge
            gx(2:end-1, 1) gy(2:end-1, 1) ; ... % left edge
            ]';
        xi = pts_tps_map(em', x', em2')';
        x = [x xi(:, 1:10:end)];
        % local support interpolation with boundary conditions provided by
        % the TPS
        em = [em em2(:, 1:10:end)];
        xi = [...
            mba_surface_interpolation(em(1, :)', em(2, :)', x(1, :)', gx(:), gy(:), nlev) ...
            mba_surface_interpolation(em(1, :)', em(2, :)', x(2, :)', gx(:), gy(:), nlev) ...
            mba_surface_interpolation(em(1, :)', em(2, :)', x(3, :)', gx(:), gy(:), nlev)];
    otherwise
        error('Interpolation method not implemented')
end

author: tomdoel

repo: pulmonarytoolkit


text:
%  internal function

%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)

function [img,hdr] = load_nii_img(hdr,filetype,fileprefix,machine,img_idx,dim5_idx,dim6_idx,dim7_idx,old_RGB)

   if ~exist('hdr','var') | ~exist('filetype','var') | ~exist('fileprefix','var') | ~exist('machine','var')
      error('Usage: [img,hdr] = load_nii_img(hdr,filetype,fileprefix,machine,[img_idx],[dim5_idx],[dim6_idx],[dim7_idx],[old_RGB]);');
   end

   if ~exist('img_idx','var') | isempty(img_idx) | hdr.dime.dim(5)<1
      img_idx = [];
   end

   if ~exist('dim5_idx','var') | isempty(dim5_idx) | hdr.dime.dim(6)<1
      dim5_idx = [];
   end

   if ~exist('dim6_idx','var') | isempty(dim6_idx) | hdr.dime.dim(7)<1
      dim6_idx = [];
   end

   if ~exist('dim7_idx','var') | isempty(dim7_idx) | hdr.dime.dim(8)<1
      dim7_idx = [];
   end

   if ~exist('old_RGB','var') | isempty(old_RGB)
      old_RGB = 0;
   end

   %  check img_idx
   %
   if ~isempty(img_idx) & ~isnumeric(img_idx)
      error('"img_idx" should be a numerical array.');
   end

   if length(unique(img_idx)) ~= length(img_idx)
      error('Duplicate image index in "img_idx"');
   end

   if ~isempty(img_idx) & (min(img_idx) < 1 | max(img_idx) > hdr.dime.dim(5))
      max_range = hdr.dime.dim(5);

      if max_range == 1
         error(['"img_idx" should be 1.']);
      else
         range = ['1 ' num2str(max_range)];
         error(['"img_idx" should be an integer within the range of [' range '].']);
      end
   end

   %  check dim5_idx
   %
   if ~isempty(dim5_idx) & ~isnumeric(dim5_idx)
      error('"dim5_idx" should be a numerical array.');
   end

   if length(unique(dim5_idx)) ~= length(dim5_idx)
      error('Duplicate index in "dim5_idx"');
   end

   if ~isempty(dim5_idx) & (min(dim5_idx) < 1 | max(dim5_idx) > hdr.dime.dim(6))
      max_range = hdr.dime.dim(6);

      if max_range == 1
         error(['"dim5_idx" should be 1.']);
      else
         range = ['1 ' num2str(max_range)];
         error(['"dim5_idx" should be an integer within the range of [' range '].']);
      end
   end

   %  check dim6_idx
   %
   if ~isempty(dim6_idx) & ~isnumeric(dim6_idx)
      error('"dim6_idx" should be a numerical array.');
   end

   if length(unique(dim6_idx)) ~= length(dim6_idx)
      error('Duplicate index in "dim6_idx"');
   end

   if ~isempty(dim6_idx) & (min(dim6_idx) < 1 | max(dim6_idx) > hdr.dime.dim(7))
      max_range = hdr.dime.dim(7);

      if max_range == 1
         error(['"dim6_idx" should be 1.']);
      else
         range = ['1 ' num2str(max_range)];
         error(['"dim6_idx" should be an integer within the range of [' range '].']);
      end
   end

   %  check dim7_idx
   %
   if ~isempty(dim7_idx) & ~isnumeric(dim7_idx)
      error('"dim7_idx" should be a numerical array.');
   end

   if length(unique(dim7_idx)) ~= length(dim7_idx)
      error('Duplicate index in "dim7_idx"');
   end

   if ~isempty(dim7_idx) & (min(dim7_idx) < 1 | max(dim7_idx) > hdr.dime.dim(8))
      max_range = hdr.dime.dim(8);

      if max_range == 1
         error(['"dim7_idx" should be 1.']);
      else
         range = ['1 ' num2str(max_range)];
         error(['"dim7_idx" should be an integer within the range of [' range '].']);
      end
   end

   [img,hdr] = read_image(hdr,filetype,fileprefix,machine,img_idx,dim5_idx,dim6_idx,dim7_idx,old_RGB);

   return					% load_nii_img

%---------------------------------------------------------------------
function [img,hdr] = read_image(hdr,filetype,fileprefix,machine,img_idx,dim5_idx,dim6_idx,dim7_idx,old_RGB)

   switch filetype
   case {0, 1}
      fn = [fileprefix '.img'];
   case 2
      fn = [fileprefix '.nii'];
   end

   fid = fopen(fn,'r',machine);

   if fid < 0,
      msg = sprintf('Cannot open file %s.',fn);
      error(msg);
   end

   %  Set bitpix according to datatype
   %
   %  /*Acceptable values for datatype are*/
   %
   %     0 None                     (Unknown bit per voxel) % DT_NONE, DT_UNKNOWN
   %     1 Binary                         (ubit1, bitpix=1) % DT_BINARY
   %     2 Unsigned char         (uchar or uint8, bitpix=8) % DT_UINT8, NIFTI_TYPE_UINT8
   %     4 Signed short                  (int16, bitpix=16) % DT_INT16, NIFTI_TYPE_INT16
   %     8 Signed integer                (int32, bitpix=32) % DT_INT32, NIFTI_TYPE_INT32
   %    16 Floating point    (single or float32, bitpix=32) % DT_FLOAT32, NIFTI_TYPE_FLOAT32
   %    32 Complex, 2 float32      (Use float32, bitpix=64) % DT_COMPLEX64, NIFTI_TYPE_COMPLEX64
   %    64 Double precision  (double or float64, bitpix=64) % DT_FLOAT64, NIFTI_TYPE_FLOAT64
   %   128 uint8 RGB                 (Use uint8, bitpix=24) % DT_RGB24, NIFTI_TYPE_RGB24
   %   256 Signed char            (schar or int8, bitpix=8) % DT_INT8, NIFTI_TYPE_INT8
   %   511 Single RGB              (Use float32, bitpix=96) % DT_RGB96, NIFTI_TYPE_RGB96
   %   512 Unsigned short               (uint16, bitpix=16) % DT_UNINT16, NIFTI_TYPE_UNINT16
   %   768 Unsigned integer             (uint32, bitpix=32) % DT_UNINT32, NIFTI_TYPE_UNINT32
   %  1024 Signed long long              (int64, bitpix=64) % DT_INT64, NIFTI_TYPE_INT64
   %  1280 Unsigned long long           (uint64, bitpix=64) % DT_UINT64, NIFTI_TYPE_UINT64
   %  1536 Long double, float128  (Unsupported, bitpix=128) % DT_FLOAT128, NIFTI_TYPE_FLOAT128
   %  1792 Complex128, 2 float64  (Use float64, bitpix=128) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128
   %  2048 Complex256, 2 float128 (Unsupported, bitpix=256) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128
   %
   switch hdr.dime.datatype
   case   1,
      hdr.dime.bitpix = 1;  precision = 'ubit1';
   case   2,
      hdr.dime.bitpix = 8;  precision = 'uint8';
   case   4,
      hdr.dime.bitpix = 16; precision = 'int16';
   case   8,
      hdr.dime.bitpix = 32; precision = 'int32';
   case  16,
      hdr.dime.bitpix = 32; precision = 'float32';
   case  32,
      hdr.dime.bitpix = 64; precision = 'float32';
   case  64,
      hdr.dime.bitpix = 64; precision = 'float64';
   case 128,
      hdr.dime.bitpix = 24; precision = 'uint8';
   case 256
      hdr.dime.bitpix = 8;  precision = 'int8';
   case 511
      hdr.dime.bitpix = 96; precision = 'float32';
   case 512
      hdr.dime.bitpix = 16; precision = 'uint16';
   case 768
      hdr.dime.bitpix = 32; precision = 'uint32';
   case 1024
      hdr.dime.bitpix = 64; precision = 'int64';
   case 1280
      hdr.dime.bitpix = 64; precision = 'uint64';
   case 1792,
      hdr.dime.bitpix = 128; precision = 'float64';
   otherwise
      error('This datatype is not supported');
   end

   hdr.dime.dim(find(hdr.dime.dim < 1)) = 1;

   %  move pointer to the start of image block
   %
   switch filetype
   case {0, 1}
      fseek(fid, 0, 'bof');
   case 2
      fseek(fid, hdr.dime.vox_offset, 'bof');
   end

   %  Load whole image block for old Analyze format or binary image;
   %  otherwise, load images that are specified in img_idx, dim5_idx,
   %  dim6_idx, and dim7_idx
   %
   %  For binary image, we have to read all because pos can not be
   %  seeked in bit and can not be calculated the way below.
   %
   if hdr.dime.datatype == 1 | isequal(hdr.dime.dim(5:8),ones(1,4)) | ...
	(isempty(img_idx) & isempty(dim5_idx) & isempty(dim6_idx) & isempty(dim7_idx))

      %  For each frame, precision of value will be read
      %  in img_siz times, where img_siz is only the
      %  dimension size of an image, not the byte storage
      %  size of an image.
      %
      img_siz = prod(hdr.dime.dim(2:8));

      %  For complex float32 or complex float64, voxel values
      %  include [real, imag]
      %
      if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792
         img_siz = img_siz * 2;
      end

      %MPH: For RGB24, voxel values include 3 separate color planes
      %
      if hdr.dime.datatype == 128 | hdr.dime.datatype == 511
	 img_siz = img_siz * 3;
      end

      img = fread(fid, img_siz, sprintf('*%s',precision));

      d1 = hdr.dime.dim(2);
      d2 = hdr.dime.dim(3);
      d3 = hdr.dime.dim(4);
      d4 = hdr.dime.dim(5);
      d5 = hdr.dime.dim(6);
      d6 = hdr.dime.dim(7);
      d7 = hdr.dime.dim(8);

      if isempty(img_idx)
         img_idx = 1:d4;
      end

      if isempty(dim5_idx)
         dim5_idx = 1:d5;
      end

      if isempty(dim6_idx)
         dim6_idx = 1:d6;
      end

      if isempty(dim7_idx)
         dim7_idx = 1:d7;
      end
   else

      d1 = hdr.dime.dim(2);
      d2 = hdr.dime.dim(3);
      d3 = hdr.dime.dim(4);
      d4 = hdr.dime.dim(5);
      d5 = hdr.dime.dim(6);
      d6 = hdr.dime.dim(7);
      d7 = hdr.dime.dim(8);

      if isempty(img_idx)
         img_idx = 1:d4;
      end

      if isempty(dim5_idx)
         dim5_idx = 1:d5;
      end

      if isempty(dim6_idx)
         dim6_idx = 1:d6;
      end

      if isempty(dim7_idx)
         dim7_idx = 1:d7;
      end

      %  compute size of one image
      %
      img_siz = prod(hdr.dime.dim(2:4));

      %  For complex float32 or complex float64, voxel values
      %  include [real, imag]
      %
      if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792
         img_siz = img_siz * 2;
      end

      %MPH: For RGB24, voxel values include 3 separate color planes
      %
      if hdr.dime.datatype == 128 | hdr.dime.datatype == 511
         img_siz = img_siz * 3;
      end

      % preallocate img
      img = zeros(img_siz, length(img_idx)*length(dim5_idx)*length(dim6_idx)*length(dim7_idx) );
      currentIndex = 1;

      for i7=1:length(dim7_idx)
         for i6=1:length(dim6_idx)
            for i5=1:length(dim5_idx)
               for t=1:length(img_idx)

                  %  Position is seeked in bytes. To convert dimension size
                  %  to byte storage size, hdr.dime.bitpix/8 will be
                  %  applied.
                  %
                  pos = sub2ind([d1 d2 d3 d4 d5 d6 d7], 1, 1, 1, ...
			img_idx(t), dim5_idx(i5),dim6_idx(i6),dim7_idx(i7)) -1;
                  pos = pos * hdr.dime.bitpix/8;

                  if filetype == 2
                     fseek(fid, pos + hdr.dime.vox_offset, 'bof');
                  else
                     fseek(fid, pos, 'bof');
                  end

                  %  For each frame, fread will read precision of value
                  %  in img_siz times
                  %
                  img(:,currentIndex) = fread(fid, img_siz, sprintf('*%s',precision));
                  currentIndex = currentIndex +1;

               end
            end
         end
      end
   end

   %  For complex float32 or complex float64, voxel values
   %  include [real, imag]
   %
   if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792
      img = reshape(img, [2, length(img)/2]);
      img = complex(img(1,:)', img(2,:)');
   end

   fclose(fid);

   %  Update the global min and max values
   %
   hdr.dime.glmax = double(max(img(:)));
   hdr.dime.glmin = double(min(img(:)));

   %  old_RGB treat RGB slice by slice, now it is treated voxel by voxel
   %
   if old_RGB & hdr.dime.datatype == 128 & hdr.dime.bitpix == 24
      % remove squeeze
      img = (reshape(img, [hdr.dime.dim(2:3) 3 hdr.dime.dim(4) length(img_idx) length(dim5_idx) length(dim6_idx) length(dim7_idx)]));
      img = permute(img, [1 2 4 3 5 6 7 8]);
   elseif hdr.dime.datatype == 128 & hdr.dime.bitpix == 24
      % remove squeeze
      img = (reshape(img, [3 hdr.dime.dim(2:4) length(img_idx) length(dim5_idx) length(dim6_idx) length(dim7_idx)]));
      img = permute(img, [2 3 4 1 5 6 7 8]);
   elseif hdr.dime.datatype == 511 & hdr.dime.bitpix == 96
      img = double(img(:));
      img = single((img - min(img))/(max(img) - min(img)));
      % remove squeeze
      img = (reshape(img, [3 hdr.dime.dim(2:4) length(img_idx) length(dim5_idx) length(dim6_idx) length(dim7_idx)]));
      img = permute(img, [2 3 4 1 5 6 7 8]);
   else
      % remove squeeze
      img = (reshape(img, [hdr.dime.dim(2:4) length(img_idx) length(dim5_idx) length(dim6_idx) length(dim7_idx)]));
   end

   if ~isempty(img_idx)
      hdr.dime.dim(5) = length(img_idx);
   end

   if ~isempty(dim5_idx)
      hdr.dime.dim(6) = length(dim5_idx);
   end

   if ~isempty(dim6_idx)
      hdr.dime.dim(7) = length(dim6_idx);
   end

   if ~isempty(dim7_idx)
      hdr.dime.dim(8) = length(dim7_idx);
   end

   return						% read_image

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function isdir_or_mkdir(p)

% ISDIR_OR_MKDIR Checks that a directory exists, or if not, creates the directory and
% all its parent directories.
%
% See also FOPEN_OR_ERROR

tok = tokenize(p, filesep);
if isempty(tok{1})
  tok{1} = filesep;
end
for i=1:numel(tok)
  d = fullfile(tok{1:i});
  if ~isfolder(d)
    ft_notice('creating directory %s', d);
    mkdir(d);
  end
end

author: spm


text:
function [ x ] = psit_fp( X,h,n,image)
% function [ x ] = psit_fp( X,h,n,image)
% PSI_FP maps pixels to coefficients
% Input:
%       X       : pixel values
%       h   : wavelet transform filter coefficients
%       n   : signal size
%       image  : 1 or 0 based on whether or not x is a square image
%Output:
%       X   : the coefficient values associated with the pixel values X
if image==1
    X=reshape(X,[sqrt(n),sqrt(n)]);
end
x=mdwt(X,h,log2(sqrt(n)));
x=x(:);
end

author: ricedsp


text:
% DEMSILHOUETTEAVERAGE Shows the average of the poses.

% FORMAT
% DESC show the average of the training poses from Agarwal and Triggs.
%
% SEEALSO : demSilhouetteGp1, demSilhouetteGp2
%
% COPYRIGHT : Neil D. Lawrence, 2008

% GP

randn('seed', 1e7)
rand('seed', 1e7)

dataSetName = 'silhouette';

capName = dataSetName;;
capName(1) = upper(capName(1));

% load data
[X, y, XTest, yTest] = mapLoadData(dataSetName);

yPred = repmat(mean(y), size(yTest, 1), 1);
xyzankurAnimCompare(yPred, yTest);

yDiff = (yPred - yTest);
rmsError = sqrt(sum(sum(yDiff.*yDiff))/prod(size(yDiff)));

counter = 0;
if printDiagram
  fileBase = ['dem' capName 'Average'];
  figure
  handle = xyzankurVisualise(yPred(1, :), 1);
  printPlot([fileBase], '../tex/diagrams', '../html')
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function c=cuttype(t,j)
%CUTTYPE Type of cut used at branches in decision tree.
%   C=CUTTYPE(T) returns an N-element cell array C indicating the type of
%   cut at each node in the tree T, where N is the number of nodes in the
%   tree.  For each node K, C{K} is 'continuous' if the cut is defined in
%   the form Z<V for a variable Z and cutpoint V, 'categorical' if the cut
%   is defined by whether Z takes a value in a set of categories, or '' if
%   K is not a branch node.  The CUTVAR method returns the cutpoints for
%   'continuous' cuts, and the CUTCATEGORIES method returns the set of
%   categories for 'categorical' cuts.
%
%   V=CUTTYPE(T,J) takes an array J of node numbers and returns the cut
%   types for the specified nodes.
%
%   See also CLASSREGTREE, CLASSREGTREE/NUMNODES, CLASSREGTREE/CUTVAR.

%   Copyright 2006-2007 The MathWorks, Inc.
%   $Revision: 1.1.6.2 $  $Date: 2007/02/15 21:48:06 $

if nargin>=2 && ~validatenodes(t,j)
    error('stats:classregtree:parent:InvalidNode',...
          'J must be an array of node numbers or a logical array of the proper size.');
end

% Get variable defining each cut
if nargin<2
    v = t.var;
else
    v = t.var(j,:);
end
v = v(:);

% Create array of cut types
c = repmat({''},numel(v),1);
c(v>0) = {'continuous'};
c(v<0) = {'categorical'};

author: Cloud-CV

repo: object-proposals


text:
function [im] = im_pad(im)

[w,h,~] = size(im);

if mod(w,2)==1
    im = cat(1,im, im(end,:,:)) ;
end
if mod(h,2)==1
    im = cat(2,im, im(:,end,:)) ;
end

end

author: cszn

repo: SRMD

save_path: github-repos/MATLAB/cszn-SRMD


text:
function [theta, J_history] = gradientDescent(X, y, theta, alpha, num_iters)
%GRADIENTDESCENT Performs gradient descent to learn theta
%   theta = GRADIENTDESENT(X, y, theta, alpha, num_iters) updates theta by
%   taking num_iters gradient steps with learning rate alpha

% Initialize some useful values
m = length(y); % number of training examples
J_history = zeros(num_iters, 1);

for iter = 1:num_iters

    % ====================== YOUR CODE HERE ======================
    % Instructions: Perform a single gradient step on the parameter vector
    %               theta.
    %
    % Hint: While debugging, it can be useful to print out the values
    %       of the cost function (computeCost) and gradient here.
    %

    delta = (1/m)*sum(X.*repmat((X*theta - y), 1, size(X,2)));

    theta = (theta' - (alpha * delta))';

    % ============================================================

    % Save the cost J in every iteration
    J_history(iter) = computeCost(X, y, theta);

end

end

author: Borye

path: github-repos/MATLAB/Borye-machine-learning-coursera-1/machine-learning-coursera-1-033fdc2e6da393eeb1179a09aafe92362021effb/Week 2 Assignments/Linear Regression with Multiple Variables/mlclass-ex1/gradientDescent.m


text:
function b = inCircle(point, circle)
%INCIRCLE test if a point is located inside a given circle.
%
%   B = inCircle(POINT, CIRCLE)
%   return true if point is located inside the circle
%
%   Example:
%   inCircle([1 0], [0 0 1])
%   inCircle([0 0], [0 0 1])
%   returns true, whereas
%   inCircle([1 1], [0 0 1])
%   return false
%
%   See also:
%   circles2d, onCircle

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2004-04-07
% Copyright 2004 INRA - TPV URPOI - BIA IMASTE

% deprecation warning
warning('geom2d:deprecated', ...
    '''inCircle'' is deprecated, use ''isPointInCircle'' instead');

d = sqrt(sum(power(point - circle(:,1:2), 2), 2));
b = d-circle(:,3)<=1e-12;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function test_bug3192

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_plot_mesh ft_plot_box ft_plot_headmodel ft_plot_dipole ft_plot_headshape

% one color for all vertex
cfg = [];
cfg.sourcemodel.xgrid  = -20:5:20;
cfg.sourcemodel.ygrid  = -20:5:20;
cfg.sourcemodel.zgrid  = -20:5:20;
grid  = ft_prepare_sourcemodel(cfg);
figure, ft_plot_mesh(grid, 'vertexcolor', 'blue', 'facecolor', 'brain', 'edgecolor', 'skull')

% different colors for each vertex
c = [];
temp = {'r','b'};
for iPos=1:length(grid.pos), c = [c temp{mod(iPos,2)+1}]; end
figure, ft_plot_mesh(grid, 'vertexcolor', c)

figure, ft_plot_box([-1 1 2 3], 'facecolor', 'brain')

% ft_plot_sens, _vol, _headshape, and _dipole just forward to ft_plot_mesh
elecs = [];
elecs.elecpos = [23 42 -31; 69 52 1; 61 67 26; 52 65 45];
elecs.label = {'1' '2' '3' '4'};
figure, ft_plot_sens(elecs, 'edgecolor', 'red')

load(dccnpath('/home/common/matlab/fieldtrip/template/headmodel/standard_singleshell.mat'))
figure, ft_plot_headmodel(vol, 'edgecolor', 'blue', 'vertexcolor', 'red', 'facecolor', 'brain')

vol.pos = [23 42 -31; 69 52 1; 61 67 26; 52 65 45];
figure, ft_plot_headshape(vol, 'edgecolor', 'blue', 'vertexcolor', 'red', 'facecolor', 'brain')

figure, ft_plot_dipole([1 2 3], [1 2 3], 'color', 'brain')

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Set axes limit.

clear all;
addpath('../lib');

%% lets plot 3 cycles of 50Hz AC voltage
f = 50;
Vm = 10;
phi = 0;

% generate the signal
t = [0:0.0001:3/f];
th = 2*pi*f*t;
v = Vm*sin(th+phi);

figure;
plot(t*1E3, v);

%% plot now
opt.XLabel = 'Time, t (ms)'; % xlabel
opt.YLabel = 'Voltage, V (V)'; %ylabel
opt.XLim = [0, 40]; % set x axis limit
opt.YLim = [-11, 11]; % set y axis limit

% Save? comment the following line if you do not want to save
opt.FileName = 'plotAxisLimit.jpg';

% create the plot
setPlotProp(opt);

author: masumhabib

repo: PlotPub

save_path: github-repos/MATLAB/masumhabib-PlotPub


text:
function [I,src] = rayMeasure(ray,Xmes,rad,rMax)
%+========================================================================+
%|                                                                        |
%|           OPENRAY - LIBRARY FOR TRI-DIMENSIONAL RAY TRACING            |
%|           openRay is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : rayMeasure.m                                  |
%|    #    |   VERSION    : 0.41                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 01.04.2018                                    |
%| ( === ) |   SYNOPSIS   : Spherical measurement                         |
%|  `---'  |                                                              |
%+========================================================================+

% Initialization
I   = cell(1,length(ray.pos)-1);
src = cell(length(ray.pos)-1,1);
dst = zeros(length(ray),1);

% Loop
for i = 1:size(I,2)
    % Ray positions
    Pi   = ray.pos{i};
    Pip1 = ray.pos{i+1};

    % Ray direction and length
    U   = (Pip1 - Pi);
    lgt = sqrt(sum(U.^2,2));
    U   = U ./ (lgt * [1 1 1]);

    % Initial position -> micro
    PM = ones(size(Pi,1),1)*Xmes - Pi;

    % Measure triangle (pythagore)
    hyp2 = sum(PM.^2,2);
    adj  = sum(PM.*U,2);
    opp2 = hyp2 - adj.^2;

    % Measured ray
    I{i} = find( (hyp2 < lgt.^2) & (adj >= 0) & (opp2 <= rad^2) & (dst + sqrt(hyp2) < rMax) );

    % Sources (focusing)
    if ~isempty(I{i})
        src{i} = Pi(I{i},:) - (dst(I{i})*[1 1 1]) .* U(I{i},:);
    else
        src{i} = zeros(0,3);
    end

    % Total distance from initial position
    dst = dst + lgt;
end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
classdef ParadigmFBCSP < ParadigmDataflowSimplified
    % Paradigm for complex oscillatory processes using the filter-bank CSP algorithm.
    % Result = para_multiband_csp(Input-Data, Operation-Mode, Options...)
    %
    % Filter-bank CSP [1,2] is a simple extension of the basic CSP method (see ParadigmCSP), in which for
    % each of several time/frequency filters a set of CSP filters is learned, followed by log-variance
    % feature extraction, concatenation of all features (over all chosen spectral filters) and
    % subsequent machine learning. It is not a general replacement for CSP due to the problem of
    % overfitting, but is very useful whenever oscillatory processes in different frequency bands (and
    % with different spatial topographies) are jointly active, and their concerted behavior must be
    % taken into account for a given prediction task. Filter-bank CSP can also be used to capture
    % oscillations in multiple time windows, instead of frequency windows (for example for the detection
    % of complex event-related dynamics).
    %
    % Since the dimensionality of the feature space is larger than in CSP, and since complex
    % interactions may be present, a more complex classifier than the default LDA may be necessary to
    % learn an appropriate model. On the other hand, more flexibility amplifies the risk of overfitting
    % (especially with only little calibration data), so that the performance should always be compared
    % to standard CSP (and Spec-CSP). Another reason is that complex (relevant) interactions between
    % different frequency bands are seemingly rarely observed in practice. The most important
    % user-configurable parameters are the selection regions in time and frequency and the learner
    % component.
    %
    % Typical applications would be those in which either complex event-related oscillatory dynamics
    % happen (for example when reacting to a particular stimulus) and/or where non-trivial interactions
    % between frequency bands (e.g. alpha/theta) are relevant, such as, for example, in workload
    % measurements.
    %
    % Example: Consider a calibration data set in which a subject is maintaining and updating
    % different number of items in his/her working memory at different times, e.g. while performing
    % the n-Back task [2]. Events with types 'n1','n2','n3' indicate challenge stimuli in which the
    % respective number of items is being processed by the person. The goal is to be able to predict
    % the working-memory load of the person following the presentation of such a memory-related
    % challenge. An epoch of 3 seconds relative to each challenge is selected, and three different
    % regions are chosen, two of them over the entire interval, covering the theta and alpha ryhthm,
    % respectively, and one region that is restricted to a window around the time of heaviest
    % cognitive processing. The three regions are specified as a cell array of flt_select
    % parameters.
    %
    %   data = io_loadset('data sets/mary/nback.eeg')
    %   myapproach = {'FBCSP' 'SignalProcessing',{'EpochExtraction',[-0.5 2.5]}, ...
    %       'Prediction', {'FeatureExtraction',{'FreqWindows',[4 6; 7 15; 7 15],'TimeWindows',[-0.5 2.5; -0.5 2.5; 0.25 1.25]}, ...
    %                      'MachineLearning',{'Learner','logreg'}}}
    %   [loss,model,stats] = bci_train('Data',data, 'Approach','ParadigmFBCSP, 'TargetMarkers',{'n1','n2','n3'})
    %
    % References;
    %   [1] Quadrianto Novi, Cuntai Guan, Tran Huy Dat, and Ping Xue, "Sub-band Common Spatial Pattern (SBCSP) for Brain-Computer Interface"
    %       Proceedings of the 3rd International IEEE EMBS Conference on Neural Engineering Kohala Coast, Hawaii, USA, May 2-5, 2007
    %   [2] Kai K. Ang, Zhang Y. Chin, Haihong Zhang, Cuntai Guan, "Filter Bank Common Spatial Pattern (FBCSP) in Brain-Computer Interface"
    %       In 2008 IEEE International Joint Conference on Neural Networks (IEEE World Congress on Computational Intelligence) (June 2008), pp. 2390-2397.
    %   [3] Owen, A. M., McMillan, K. M., Laird,A. R. & Bullmore, E. "N-back working memory paradigm: A meta-analysis of normative functional neuroimaging studies."
    %       Human Brain Mapping, 25, 46-59, 2005
    %
    % Name:
    %   Filter-Bank CSP
    %
    %                           Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
    %                           2010-04-29

    methods

        function defaults = preprocessing_defaults(self)
            % define the default pre-processing parameters of this paradigm
            defaults = {'EpochExtraction',[0.5 3.5],'Resampling',200};
        end

        function model = feature_adapt(self,varargin)
            % adapt a feature representation using the CSP algorithm
            args = arg_define(varargin, ...
                arg_norep('signal'), ...
                arg({'patterns','PatternPairs'},3,uint32([1 1 64 10000]),'CSP patterns per band (times two).','cat','Feature Extraction'), ...
                arg({'shrinkage_cov','ShrinkageCovariance'},false,[],'Shrinkage covariance estimator. Whether to use shrinkage to estimate the covariance matrices.'), ...
                arg({'robust_cov','RobustCovariance'},false,[],'Robust covariance estimation. Whether to use robust cov estimation.'), ...
                arg({'freqwnds','FreqWindows'},[0.5 3; 4 7; 8 12; 13 30; 31 42],[0 0.5 200 1000],'Frequency bands of interest. Matrix containing one row for the start and end of each frequency band from which CSP patterns shall be computed. Values in Hz.','cat','Feature Extraction'), ...
                arg({'timewnds','TimeWindows'},[],[],'Time windows of interest. Matrix containing one row for the start and end of each time window from which CSP patterns shall be computed. Values in seconds. If both this and the freqwnds parameter are non-empty, they should have the same number of rows.','cat','Feature Extraction'), ...
                arg({'winfunc','WindowFunction'},'rect',{'barthann','bartlett','blackman','blackmanharris','bohman','cheb','flattop','gauss','hamming','hann','kaiser','nuttall','parzen','rect','taylor','triang','tukey'},'Type of window function. Typical choices are rect (rectangular), hann, gauss, blackman and kaiser.'),...
                arg({'winparam','WindowParameter','param'},[],[],'Parameter of the window function. This is mandatory for cheb, kaiser and tukey and optional for some others.','shape','scalar'));

            if args.signal.nbchan == 1
                error('Multi-band CSP does intrinsically not support single-channel data (it is a spatial filter).'); end
            if args.signal.nbchan < args.patterns
                error('Multi-band CSP prefers to work on at least as many channels as you request output patterns. Please reduce the number of pattern pairs.'); end
            if ~isempty(args.freqwnds) && ~isempty(args.timewnds) && size(args.freqwnds,1) ~= size(args.timewnds,1)
                error('If both time and frequency windows are specified, both arrays must have the same number of rows (together they define the windows in time and frequency).'); end
            if isempty(args.timewnds)
                args.timewnds = zeros(size(args.freqwnds,1),0); end
            if isempty(args.freqwnds)
                args.freqwnds = zeros(size(args.timewnds,1),0); end

            filters = [];
            patterns = [];
            % for each window
            for w = 1:max(size(args.freqwnds,1),size(args.timewnds,1))
                % pre-parse arguments for flt_window and flt_spectrum (for fast subsequent online use)
                time_args{w} = arg_report('vals',@flt_window,{'time',{args.timewnds(w,:),args.winfunc,args.winparam}});
                freq_args{w} = arg_report('vals',@flt_spectrum,{'freq',args.freqwnds(w,:)});
                for k=1:2
                    % filter trial subrange in time and frequency
                    data = exp_eval_optimized(flt_spectrum('signal',flt_window('signal',set_picktrials(args.signal,'rank',k),time_args{w}),freq_args{w}));
                    if args.robust_cov
                        covar{k} = hlp_diskcache('featuremodels',@cov_blockgeom,reshape(data.data,size(data.data,1),[])',max([data.nbchan*2,data.srate*2,size(data,3)]));
                    else
                        if args.shrinkage_cov
                            covar{k} = hlp_diskcache('featuremodels',@cov_shrink,reshape(data.data,size(data.data,1),[])');
                        else
                            covar{k} = cov(reshape(data.data,size(data.data,1),[])');
                        end
                    end
                    covar{k}(~isfinite(covar{k})) = 0;
                end
                [V,D] = eig(covar{1},covar{1}+covar{2}); %#ok<NASGU>
                P = inv(V);
                filters = [filters V(:,[1:args.patterns end-args.patterns+1:end])];
                patterns = [patterns P([1:args.patterns end-args.patterns+1:end],:)'];
            end
            model = struct('filters',{filters},'patterns',{patterns},'time_args',{time_args},'freq_args',{freq_args},'chanlocs',{args.signal.chanlocs});
        end

        function features = feature_extract(self,signal,featuremodel)
            W = length(featuremodel.freq_args);
            F = size(featuremodel.filters,2);
            T = size(signal.data,3);
            features = zeros(T,F);
            for w = 1:W
                % filter data in time & frequency
                data = exp_eval_optimized(flt_spectrum('signal',flt_window('signal',signal,featuremodel.time_args{w}),featuremodel.freq_args{w}));
                inds = (w-1)*(F/W)+(1:(F/W));
                for t=1:T
                    features(t,inds) = log(var(data.data(:,:,t)' * featuremodel.filters(:,inds))); end
            end
        end

        function visualize_model(self,varargin) %#ok<*INUSD>
            args = arg_define([0 3],varargin, ...
                arg_norep({'myparent','Parent'},[],[],'Parent figure.'), ...
                arg_norep({'featuremodel','FeatureModel'},[],[],'Feature model. This is the part of the model that describes the feature extraction.'), ...
                arg_norep({'predictivemodel','PredictiveModel'},[],[],'Predictive model. This is the part of the model that describes the predictive mapping.'), ...
                arg({'patterns','PlotPatterns'},true,[],'Plot patterns instead of filters. Whether to plot spatial patterns (forward projections) rather than spatial filters.'), ...
                arg({'weight_scaled','WeightScaled'},false,[],'Scaled by weight. Whether to scale the patterns by weight.'));
            arg_toworkspace(args);

            % find the relevant components
            scores = predictivemodel.model.w;
            scores = sqrt(abs(scores));
            % optionally remove the bias if included in w
            if length(scores) == size(featuremodel.patterns,2)+1
                scores = scores(1:end-1); end
            % frequency labels
            % titles = repmat({'delta','theta','alpha','beta','gamma'},8,1); titles = titles(:);
            % extract relevant patterns
            patterns = featuremodel.patterns(:,find(scores)); %#ok<FNDSB>
            filters = featuremodel.filters(:,find(scores)); %#ok<FNDSB>
            % plot them
            if args.weight_scaled
                if args.patterns
                    topoplot_grid(patterns,featuremodel.chanlocs,'scales',scores(find(scores))/max(scores)*1);
                else
                    topoplot_grid(filters,featuremodel.chanlocs,'scales',scores(find(scores))/max(scores)*1);
                end
            else
                if args.patterns
                    topoplot_grid(patterns,featuremodel.chanlocs);
                else
                    topoplot_grid(filters,featuremodel.chanlocs);
                end
            end
            % figure;
        end

        function layout = dialog_layout_defaults(self)
            % define the default configuration dialog layout
            layout = {'SignalProcessing.Resampling.SamplingRate', 'SignalProcessing.EpochExtraction', '', ...
                'Prediction.FeatureExtraction.FreqWindows', 'Prediction.FeatureExtraction.TimeWindows', ...
                'Prediction.FeatureExtraction.WindowFunction', '', 'Prediction.FeatureExtraction.PatternPairs', '', ...
                'Prediction.MachineLearning.Learner'};
        end

        function tf = needs_voting(self)
            tf = true;
        end
    end
end

repo: REST


text:
function [pnt, tri] = mesh_cube()

% MESH_CUBE creates a triangulated cube
%
% Use as
%   [pos, tri] = mesh_cube()
%
% See also MESH_TETRAHEDRON, MESH_OCTAHEDRON, MESH_ICOSAHEDRON, MESH_SPHERE, MESH_CONE

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

pnt = [
  -1 -1 -1
  -1  1 -1
   1  1 -1
   1 -1 -1
  -1 -1  1
  -1  1  1
   1  1  1
   1 -1  1
  ];

tri = [
  1 2 4
  2 3 4
  1 5 6
  1 6 2
  2 6 7
  2 7 3
  3 7 8
  3 8 4
  4 8 5
  4 5 1
  8 6 5
  8 7 6
  ];


text:
function out = mtimes(a,b,varargin)
% outer quaternion multiplication
%
% Syntax
%   out = a * b
%   out = a * v
%
% Input
%  a - @SO3Grid
%  b - @quaternion
%  v - @vector3d
%
% Output
%  out - @SO3Grid / @vector3d

if isa(a,'SO3Grid') % right multiplication

  if isa(b,'SO3Grid'), b = orientation(b); end

  out = orientation(a) * b;

elseif isa(a,'quaternion')

  if length(a) == 1 % rotate center only

    r = mtimes@orientation(a,b);
    out = b;
    out.a = r.a;
    out.b = r.b;
    out.c = r.c;
    out.d = r.d;
    out.i = r.i;
    if isempty(b.center)
      out.center = a;
    else
      out.center = a * out.center;
    end

  else

    out = a * orientation(b);

  end

else
  error('type mismatch!')
end

author: mtex-toolbox

repo: mtex

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


text:
function exact = p19_exact ( )

%*****************************************************************************80
%
%% P19_EXACT returns the exact integral for problem 19.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real EXACT, the value of the integral.
%
  exact = - 1.0;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
 function pre = newfft_exact_p(st, om)
%function pre = newfft_exact_p(st, om)
%|
%| precomputed (large) matrix for exact forward NUFFT

if ~isvar('om') || isempty(om)
	om = st.om; % default
end
if isempty(om), fail 'om or st.om required', end

dd = st.dd;
Nd = st.Nd;
n_shift = st.n_shift;

if dd > 3
	fail 'only up to 3D is done'
else
	Nd = [Nd(:); ones(3-length(Nd),1)];
	n_shift = [n_shift(:); zeros(3-length(n_shift),1)];
end

for id=1:3 % fix: dd
	nn{id} = [0:(Nd(id)-1)]-n_shift(id);
end

[nn{1}, nn{2}, nn{3}] = ndgrid(nn{1}, nn{2}, nn{3});

pre = 0;
for id=1:dd
	pre = pre + om(:,id) * col(nn{id})';
end
pre = exp(-1i*pre);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% OP_V_GRADP: assemble the matrix B = [b(i,j)], b(i,j) = (epsilon grad p_i, v_j).
%
%   mat = op_v_gradp (spv, spp, msh, epsilon);
%   [rows, cols, values] = op_v_gradp (spv, spp, msh, epsilon);
%
% INPUT:
%
%   spv:     structure representing the space of vectorial trial functions  (see sp_vector/sp_evaluate_col)
%   spp:     structure representing the space of scalar test functions (see sp_scalar/sp_evaluate_col)
%   msh:     structure containing the domain partition and the quadrature rule (see msh_cartesian/msh_evaluate_col)
%   epsilon: physical parameter
%
% OUTPUT:
%
%   mat:    assembled matrix
%   rows:   row indices of the nonzero entries
%   cols:   column indices of the nonzero entries
%   values: values of the nonzero entries
%
% Copyright (C) 2009, 2010 Carlo de Falco
% Copyright (C) 2011, 2017 Rafael Vazquez
%
%    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/>.

function varargout = op_v_gradp (spv, spp, msh, coeff)

  ndir = size (spp.shape_function_gradients, 1);

  if (ndir ~= spv.ncomp)
    error ('Inconsistent dimensions between the number of components of v and the gradient of p')
  end

  rows = zeros (msh.nel * spp.nsh_max * spv.nsh_max, 1);
  cols = zeros (msh.nel * spp.nsh_max * spv.nsh_max, 1);
  values = zeros (msh.nel * spp.nsh_max * spv.nsh_max, 1);

  jacdet_weights = msh.jacdet .* msh.quad_weights .* coeff;

  ncounter = 0;
  for iel = 1:msh.nel
    if (all (msh.jacdet(:, iel)))
      shpv_iel = reshape (spv.shape_functions(:, :, :, iel), spv.ncomp, msh.nqn, 1, spv.nsh_max);
      gradp_iel = reshape (spp.shape_function_gradients(:, :, :, iel), ndir, msh.nqn, spp.nsh_max, 1);

      jacdet_iel = reshape (jacdet_weights(:,iel), [1,msh.nqn,1,1]);

      jacdet_shpv = bsxfun (@times, jacdet_iel, shpv_iel);
      tmp1 = sum (bsxfun (@times, jacdet_shpv, gradp_iel), 1);
      elementary_values = reshape (sum (tmp1, 2), spp.nsh_max, spv.nsh_max);

      [rows_loc, cols_loc] = ndgrid (spp.connectivity(:,iel), spv.connectivity(:,iel));
      indices = rows_loc & cols_loc;
      rows(ncounter+(1:spv.nsh(iel)*spp.nsh(iel))) = rows_loc(indices);
      cols(ncounter+(1:spv.nsh(iel)*spp.nsh(iel))) = cols_loc(indices);
      values(ncounter+(1:spv.nsh(iel)*spp.nsh(iel))) = elementary_values(indices);
      ncounter = ncounter + spv.nsh(iel)*spp.nsh(iel);
    else
      warning ('geopdes:jacdet_zero_at_quad_node', 'op_v_gradp: singular map in element number %d', iel)
    end
  end

  if (nargout == 1 || nargout == 0)
    varargout{1} = sparse (rows(1:ncounter), cols(1:ncounter), ...
                           values(1:ncounter), spp.ndof, spv.ndof);
  elseif (nargout == 3)
    varargout{1} = rows(1:ncounter);
    varargout{2} = cols(1:ncounter);
    varargout{3} = values(1:ncounter);
  else
    error ('op_v_gradp: wrong number of output arguments')
  end

end

%% COPY OF THE FIRST VERSION OF THE FUNCTION (MORE UNDERSTANDABLE)
% function mat = op_v_gradp (spv, spp, msh, coeff)
%
%   mat = spalloc(spp.ndof, spv.ndof, 1);
%   ndir = size (spp.shape_function_gradients, 1);
%   for iel = 1:msh.nel
%     if (all (msh.jacdet(:,iel)))
%       mat_loc = zeros (spp.nsh(iel), spv.nsh(iel));
%       for idof = 1:spp.nsh(iel)
%         ishg = reshape(spp.shape_function_gradients(:,:,idof,iel),ndir,[]);
%         for jdof = 1:spv.nsh(iel)
%           jshg = reshape(spv.shape_functions(:,:,jdof,iel),ndir,[]);
% % The cycle on the quadrature points is vectorized
%           %for inode = 1:msh.nqn
%             mat_loc(idof, jdof) = mat_loc(idof, jdof) + ...
%              sum (msh.jacdet(:, iel) .* msh.quad_weights(:, iel) .* ...
%              sum (ishg .* jshg, 1).' .* coeff(:, iel));
%           %end
%         end
%       end
%       mat(spp.connectivity(:, iel), spv.connectivity(:, iel)) = ...
%         mat(spp.connectivity(:, iel), spv.connectivity(:, iel)) + mat_loc;
%     else
%       warning ('geopdes:jacdet_zero_at_quad_node', 'op_v_gradp: singular map in element number %d', iel)
%     end
%   end
%
% end

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
function [ Z, H, dnorm] = cnmf ( X, k, y, varargin )
% Matrix sizes
% X: m x n
% Z: m x num_of_components
% H: num_of_components x num_of_components

% Process optional arguments
pnames = {'z0' 'h0' 'bUpdateH' 'maxiter' 'nonlinearity_function', 'TolFun'};

% Do SVD initialisation of the init components

if 1
    [z0, h0] = NNDSVD(abs(X), k, 0);
else
    z0 = rand(size(X, 1), k);
    h0 = rand(k, size(X,2));
end

dflts  = {z0, h0, 1, 300, @(x) x, 1e-5};

[z0, h0, bUpdateH, max_iter, g, tolfun] = ...
        internal.stats.parseArgs(pnames,dflts,varargin{:});

% X =>  % p x n
Z = z0; % p x k

A = ind2vec([y; [length(y)+1:size(X, 2)]']')';
H = max(abs(A' * pinv(h0)), eps); % c x k

for i = 1:max_iter
    if bUpdateH
        numer = A' * X' * Z;
        H = H .* (numer ./ (((A' * A) * H * (Z' * Z)) + eps(numer)));
    end

    numer = X * A * H;
    Z = Z .* (numer ./  (Z * (H' * A') * (A * H) + eps(numer)));

    if mod(i, 10) == 0 || mod(i+1, 10) == 0
        s = X - Z * H' * A';
        dnorm = sqrt(sum(s(:).^2));

        if mod(i+1, 10) == 0
            dnorm0 = dnorm;
            continue
        end

%         if mod(i, 100) == 0
            display(sprintf('...CNMF iteration #%d out of %d, error: %f\n', i, max_iter, dnorm));
%         end

%         if exist('dnorm0')
%             assert(dnorm <= dnorm0, sprintf('Rec. error increasing! From %f to %f. (%d)', dnorm0, dnorm, k));
%         end

        % Check for convergence
        if exist('dnorm0') && dnorm0-dnorm <= tolfun*max(1,dnorm0)
            display(sprintf('Stopped at %d: dnorm: %f, dnorm0: %f', i, dnorm, dnorm0));
            break;
        end

    end
end

H = A * H;
H = H(length(y) +1 : end, :)';

author: jwyang


text:
% VL_PEGASOS PEGASOS linear SVM solver
%   W = VL_PEGASOS(X, Y, LAMBDA) learns a linear SVM W given training
%   vectors X, their labels Y, and the regularization parameter LAMBDA
%   using the PEGASOS [1] solver. The algorithm finds a minimizer W of
%   the objective function
%
%     LAMBDA/2 |W|^2 + 1/N SUM_i LOSS(W, X(:,i), Y(i))
%
%   where LOSS(W,X,Y) = MAX(0, 1 - Y W'X) is the hinge loss and N is
%   the number of training vectors in X.
%
%   ALGORITHM. PEGASOS is an implementation of stochastic subgradient
%   descent. At each iteration a data point is selected at random, the
%   subgradient of the cost function relative to that data point is
%   computed, and a step is taken in that direction. The step size is
%   inversely proportional to the iteration number. See [1] for
%   details.
%
%   NumIterations:: [10 / LAMBDA]
%     Sets the maximum number of iterations.
%
%   BiasMultiplier:: [0]
%     Appends to the data X the specified scalar value B. This
%     approximates the training of a linear SVM with bias.  The bias
%     can be recovered from the optimal weight vector W as W(end) *
%     B. It is often useful to precondition the computation of the
%     gradient with respect to this element by a small value,
%     e.g. 0.1/B (see PRECONDITIONER).
%
%   StartingModel:: [null vector]
%     Specify the initial value for the weight vector W.
%
%   StartingIteration:: [1]
%     Specify the iteration number to start from. The only effect
%     is to change the step size, as this is inversely proportional
%     to the iteration number.
%
%   Permutation:: [empty]
%     Specify a permutation PERM to be used to sample the data (this
%     disables random sampling). Specifically, at the T-th iteration
%     the algorithm takes a step w.r.t. the PERM[T']-th data point,
%     where T' is T modulo the number of data samples
%     (i.e. MOD(T'-1,NUMSAMPLES)+1). PERM needs not to be
%     bijective. This allows specifying certain data points more or
%     less frequently, implicitly increasing their relative weight in
%     the error term. A common application is to balance an unbalanced
%     dataset.
%
%   Preconditioner:: [empty]
%     Specify a diagonal preconditioner PREC. The elements of this
%     vector are multiplied to the function subgradient before adding
%     the latter to the current model estimate. The dimension of the
%     vector PREC is the same of the model, plus one if the SVM has
%     bias (see BIASMULTIPLIER).
%
%   Verbose::
%     Be verbose.
%
%   Example::
%     The options StartingModel and StartingIteration can be used
%     to continue training. I.e., the command
%
%       vl_twister('state',0) ;
%       w = vl_pegasos(x,y,lambda,'NumIterations',1000) ;
%
%     produces the same result as the sequence
%
%       vl_twister('state',0) ;
%       w = vl_pegasos(x,y,lambda,'NumIterations',500) ;
%       w = vl_pegasos(x,y,lambda,'NumIterations',500, ...
%                      'StartingIteration', 501, ...
%                      'StartingModel', w) ;
%
%   REFERENCES::
%     [1] S. Shalev-Shwartz, Y. Singer, N. Srebro, and
%     A. Cotter. Pegasos: Primal Estimated sub-GrAdient SOlver for
%     SVM. MBP, 2010.
%
%   See also: VL_HOMKERMAP(), VL_HELP().

% AUTHORIGHTS

author: jianxiongxiao

repo: ProfXkit


text:
RMS1=RMS_multifractal;
RMS2=RMS_monofractal;
RMS3=RMS_ordinary;
YMatrix1=[multifractal,mean(multifractal).*ones(size(multifractal)),(mean(multifractal)+RMS1).*ones(size(multifractal)),(mean(multifractal)-RMS1).*ones(size(multifractal))];
YMatrix2=[monofractal,mean(monofractal).*ones(size(monofractal)),(mean(monofractal)+RMS2).*ones(size(monofractal)),(mean(monofractal)-RMS2).*ones(size(monofractal))];
YMatrix3=[whitenoise,mean(whitenoise).*ones(size(whitenoise)),(mean(whitenoise)+RMS3).*ones(size(whitenoise)),(mean(whitenoise)-RMS3).*ones(size(whitenoise))];

% Create figure
figure1 = figure('PaperSize',[20.98 29.68],'Color',[1 1 1]);

% Create axes
axes1 = axes('Parent',figure1,'XTick',zeros(1,0),...
    'Position',[0.13 0.6555 0.775 0.2695],...
    'LineWidth',2,...
    'FontSize',14);
% Uncomment the following line to preserve the Y-limits of the axes
ylim(axes1,[-7 20]);
hold(axes1,'all');

% Create multiple lines using matrix input to plot
pplot1 = plot(YMatrix1,'Parent',axes1,'Color',[1 0 0],'LineWidth',2);
set(pplot1(1),'Color',[0 0 1],'DisplayName','Noise like time series',...
    'LineWidth',0.5);
set(pplot1(2),'LineStyle','--','DisplayName','Mean');
set(pplot1(3),'DisplayName','+/- 1 RMS');

% Create axes
axes2 = axes('Parent',figure1,'YTick',[-4 -2 0 2 4],'XTick',zeros(1,0),...
    'Position',[0.13 0.3876 0.775 0.2687],...
    'LineWidth',2,...
    'FontSize',14);
% Uncomment the following line to preserve the Y-limits of the axes
ylim(axes2,[-6 6]);
hold(axes2,'all');

% Create multiple lines using matrix input to plot
pplot2 = plot(YMatrix2,'Parent',axes2,'Color',[1 0 0],'LineWidth',2);
set(pplot2(1),'Color',[0 0 1],'LineWidth',0.5);
set(pplot2(2),'LineStyle','--');

% Create axes
axes3 = axes('Parent',figure1,'YTick',[-4 -2 0 2 4],...
    'Position',[0.13 0.1198 0.775 0.2692],...
    'LineWidth',2,...
    'FontSize',14);
% Uncomment the following line to preserve the Y-limits of the axes
ylim(axes3,[-6 6]);
hold(axes3,'all');

% Create multiple lines using matrix input to plot
pplot3 = plot(YMatrix3,'Parent',axes3,'Color',[1 0 0],'LineWidth',2);
set(pplot3(1),'Color',[0 0 1],'LineWidth',0.5);
set(pplot3(2),'LineStyle','--');

% Create xlabel
xlabel('Time (sample number)','FontSize',14);

% Create ylabel
ylabel('Amplitude (measurement units)','FontSize',14);

% Create legend
legend1 = legend(axes1,'Noise like time series','Mean','+/- 1 RMS');
set(legend1,'Position',[0.6322 0.8219 0.2066 0.1059]);

% Create textbox
annotation(figure1,'textbox',[0.264 0.8856 0.1917 0.05384],...
    'String',{'Multifractal time series'},...
    'FontSize',14,...
    'FitBoxToText','off',...
    'LineStyle','none');

% Create textbox
annotation(figure1,'textbox',[0.2684 0.3234 0.1917 0.05384],...
    'String',{'White noise'},...
    'FontSize',14,...
    'FitBoxToText','off',...
    'LineStyle','none');

% Create textbox
annotation(figure1,'textbox',[0.2641 0.5937 0.1917 0.05384],...
    'String',{'Monofractal time series'},...
    'FontSize',14,...
    'FitBoxToText','off',...
    'LineStyle','none');

clear pplot1 pplot2 pplot3 legend1 axes1 figure1 axes2 figure2 axes3 figure3 ans YMatrix1 YMatrix2 YMatrix3 RMS1 RMS2 RMS3

author: Sable

repo: mcbench-benchmarks


text:
function weights=realized_kernel_weights(options)
% Computes weights for Realized Kernels
%
% USAGE:
%   [WEIGHTS] = realized_kernel_weights(OPTIONS)
%
% INPUTS:
%   OPTIONS - A realized kernel options structure.  See help realized_options for details
%
% OUTPUTS:
%   WEIGHTS - H by 1 vector of kernel weights where H depends on the bandwidth and kernel
%
% COMMENTS:
%   This is a helper function for REALIZED_KERNEL. See Barndorf-Nielsen,
%   Hansen, Lunde and Shephard (2008a, 2008b) for details about realized kernels and
%   their properties.
%
%  See also REALIZED_KERNEL, REALIZED_KERNEL_BANDWIDTH,

% Copyright: Kevin Sheppard
% kevin.sheppard@economics.ox.ac.uk
% Revision: 1    Date: 5/1/2008

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin~=1
    error('One input required.')
end

% List of flat top kernels
flatTopKernelList = {'bartlett','twoscale','2ndorder','epanechnikov',...
    'cubic','multiscale','5thorder','6thorder','7thorder','8thorder','parzen',...
    'th1','th2','th5','th16'};

% List of non flat top kernels
nonFlatTopKernelList = {'nonflatparzen','qs','fejer','thinf','bnhls'};

% Combined kernel list
kernelList = [flatTopKernelList nonFlatTopKernelList];

if ~isfield(options,'kernel') || ~ismember(options.kernel,kernelList)
    error('KERNEL must be a field of OPTIONS and one of the listed types.')
end

if ~isfield(options,'bandwidth') || options.bandwidth<0
    error('BANDWIDTH must be a field of OPTIONS and a non-negative scalar.')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input Checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Take relevant fields
kernel = options.kernel;
bandwidth = options.bandwidth;

% A big lookup table
if bandwidth>0
    if ismember(options.kernel,flatTopKernelList)
        % If it is flat top then the first lag is actually 0 which produces a weights of 1
        H = round(bandwidth);
        x = (1:H)';
        x = (x-1)/H;
    else
        H = bandwidth;
    end

    switch lower(kernel)
        case {'bartlett','twoscale'}
            weights = 1-x;
        case {'2ndorder'}
            weights = 1-2*x+x.^2;
        case {'epanechnikov'}
            weights = 1-x.^2;
        case {'cubic','multiscale'}
            weights = 1-2*x.^2+2*x.^3;
        case '5thorder'
            weights = 1-10*x.^3+15*x.^4-6*x.^5;
        case '6thorder'
            weights = 1-15*x.^4+24*x.^5-10*x.^6;
        case '7thorder'
            weights = 1-21*x.^5+35*x.^6-15*x.^7;
        case '8thorder'
            weights = 1-28*x.^6+48*x.^7-21*x.^8;
        case 'parzen'
            weights = (1-6*x.^2+6*x.^3).*(x>=0 & x<=1/2) + 2*(1-x).^3.*(x>1/2 & x<1);
        case 'th1'
            weights = sin(pi/2*(1-x)).^2;
        case 'th2'
            weights = sin(pi/2*(1-x).^2).^2;
        case 'th5'
            weights = sin(pi/2*(1-x).^5).^2;
        case 'th16'
            weights = sin(pi/2*(1-x).^16).^2;
        case 'nonflatparzen'
            x = (1:H)';
            x = x./(H+1);
            weights = (1-6*x.^2+6*x.^3).*(x>=0 & x<=1/2) + 2*(1-x).^3.*(x>1/2 & x<1);
        case 'qs'
            % Truncate at 30 * H where H is BW
            x = unique(round(1:30*H))';
            x = x./(H+1);
            weights = 3./x.^2 .* (sin(x)./x - cos(x));
        case 'fejer'
            % Truncate at 30 * H where H is BW
            x = unique(round(1:30*H))';
            x = x./(H+1);
            weights = (sin(x)./x).^2;
        case 'thinf'
            % Truncate at 4 * H where H is BW
            x = unique(round(1:4*H))';
            x = x./(H+1);
            weights = sin((pi/2).*exp(-x)).^2;
        case 'bnhls'
            % Truncate at 10 * H where H is BW
            x = unique(round(1:10*H))';
            x = x./(H+1);
            weights = (1+x).*exp(-x);
        otherwise
            error('OPTIONS.KERNEL must be one of the listed types.')
    end
else
    weights=[];
end

author: bashtage

repo: mfe-toolbox


text:
function f2 = p08_f2 ( x )

%*****************************************************************************80
%
%% P08_F2 evaluates the second derivative for problem 8.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the values of the variables.
%
%    Output, real F2, the second derivative.
%
  if ( x == pi )
    f2 = 1.0;
  else
    f2 = 6.0 + ( - 2.0 / ( x - pi ) / ( x - pi ) ) / pi^4;
  end

  return
end


text:
function [t]=tt_heaviside(n,L,ind)
% Heaviside vector in the TT format
% function [t]=tt_heaviside(n,L,ind)
% Returns the tt_tensor with elements
%   theta(i-ind), i=1:N, ind=1,...,N, N=prod(n), theta(0)=1.
%
% See also: tt_unit, tt_shift
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et. al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

if (numel(n)==1)
    n = n*ones(1,L);
end;

if (numel(ind)==1)
    ind = tt_ind2sub(reshape(n, 1, []), ind);
end;

t = cell(L,1);

t{1} = zeros(1,n(1),2);
t{1}(1,:,1) = heaviside((1:n(1))-ind(1));
t{1}(1,:,2) = 1;

for i=2:L-1
    t{i} = zeros(2,n(i),2);
    t{i}(1,:,1)=[zeros(ind(i)-1,1); 1; zeros(n(i)-ind(i), 1)];
    t{i}(2,:,1)=heaviside((1:n(i))-ind(i)-1);
    t{i}(2,:,2)=1;
end;

t{L} = zeros(2,n(L),1);
t{L}(1,:,1)=[zeros(ind(L)-1,1); 1; zeros(n(L)-ind(L), 1)];
t{L}(2,:,1)=heaviside((1:n(L))-ind(L)-1);

if (L==1)
    t{1} = sum(t{1}, 1);
end;

t = cell2core(tt_tensor,t);
end

function [y]=heaviside(x)
% Heaviside function of a _real_ array
% function [y]=heaviside(x)
% Zeroes strictly negative values in x, all others are cast to 1

y = sign(x)+1;
y(y==1)=2;
y = y*0.5;
end

author: oseledets

repo: TT-Toolbox


text:
function epoch = spm_cfg_eeg_epochs
% Configuration file for M/EEG epoching
%__________________________________________________________________________
% Copyright (C) 2008-2016 Wellcome Trust Centre for Neuroimaging

% Stefan Kiebel
% $Id: spm_cfg_eeg_epochs.m 7095 2017-06-07 10:15:59Z vladimir $

D        = cfg_files;
D.tag    = 'D';
D.name   = 'File Name';
D.filter = 'mat';
D.num    = [1 1];
D.help   = {'Select the M/EEG mat file.'};

bc        = cfg_menu;
bc.tag    = 'bc';
bc.name   = 'Baseline correction';
bc.help   = {'Perform baseline correction when epoching, or not.'};
bc.labels = {'Yes','No'};
bc.values = {1 0};
bc.val    = {1};

% input via trl file
trlfile        = cfg_files;
trlfile.tag    = 'trlfile';
trlfile.name   = 'Trial definition file';
trlfile.filter = 'mat';
trlfile.num    = [1 1];
trlfile.help   = {'Select the trialfile mat file.'};

eventpadding         = cfg_entry;
eventpadding.tag     = 'eventpadding';
eventpadding.name    = 'Event padding';
eventpadding.strtype = 'r';
eventpadding.val     = {0};
eventpadding.num     = [1 1];
eventpadding.help    = {['In seconds: the additional time period around each trial',...
    'for which the events are saved with the trial (to let the ',...
    'user keep and use for analysis events which are outside ',...
    'trial borders). Default is 0 s.']};

% input via trialdef
timewin         = cfg_entry;
timewin.tag     = 'timewin';
timewin.name    = 'Time window';
timewin.strtype = 'r';
timewin.num     = [1 2];
timewin.help    = {'Start and end of epoch [ms].'};

conditionlabel         = cfg_entry;
conditionlabel.tag     = 'conditionlabel';
conditionlabel.name    = 'Condition label';
conditionlabel.strtype = 's';
conditionlabel.help    = {''};

eventtype         = cfg_entry;
eventtype.tag     = 'eventtype';
eventtype.name    = 'Event type';
eventtype.strtype = 's';
eventtype.help    = {''};

eventvalue         = cfg_entry;
eventvalue.tag     = 'eventvalue';
eventvalue.name    = 'Event value';
eventvalue.strtype = 'e';
eventvalue.help    = {''};

trlshift         = cfg_entry;
trlshift.tag     = 'trlshift';
trlshift.name    = 'Shift';
trlshift.strtype = 'r';
trlshift.num     = [1 1];
trlshift.val     = {0};
trlshift.help    = {'shift the triggers by a fixed amount (ms)',...
                   '(e.g. projector delay).'};

trialdef      = cfg_branch;
trialdef.tag  = 'trialdef';
trialdef.name = 'Trial';
trialdef.val  = {conditionlabel eventtype eventvalue, trlshift};
trialdef.help = {''};

define1        = cfg_repeat;
define1.tag    = 'unused';
define1.name   = 'Trial definitions';
define1.values = {trialdef};
define1.help   = {''};

define      = cfg_branch;
define.tag  = 'define';
define.name = 'Define trial';
define.val  = {timewin define1};
define.help = {''};

% input via trialdef
trialength         = cfg_entry;
trialength.tag     = 'trialength';
trialength.name    = 'Trial length';
trialength.strtype = 'r';
trialength.num     = [1 1];
trialength.help    = {'Arbitary trial length [ms].'};

arbitrary      = cfg_branch;
arbitrary.tag  = 'arbitrary';
arbitrary.name = 'Arbitrary trials';
arbitrary.val  = {trialength, conditionlabel};
arbitrary.help = {'Epoch the data in arbitray fixed length trials (e.g. for spectral analysis).'};

trlchoice        = cfg_choice;
trlchoice.tag    = 'trialchoice';
trlchoice.name   = 'How to define trials';
trlchoice.help   = {'Choose one of the two options how to define trials.'}';
trlchoice.values = {trlfile define arbitrary};

prefix         = cfg_entry;
prefix.tag     = 'prefix';
prefix.name    = 'Filename Prefix';
prefix.help    = {'Specify the string to be prepended to the filenames of the epoched dataset. Default prefix is ''e''.'};
prefix.strtype = 's';
prefix.num     = [1 Inf];
prefix.val     = {'e'};

epoch          = cfg_exbranch;
epoch.tag      = 'epoch';
epoch.name     = 'Epoching';
epoch.val      = {D, trlchoice, bc, eventpadding, prefix};
epoch.help     = {'Epoch continuous EEG/MEG data.'};
epoch.prog     = @eeg_epochs;
epoch.vout     = @vout_eeg_epochs;
epoch.modality = {'EEG'};

%==========================================================================
function out = eeg_epochs(job)
% construct the S struct
S.D = job.D{1};

if isfield(job.trialchoice, 'define')
    S.timewin  = job.trialchoice.define.timewin;
    S.trialdef = job.trialchoice.define.trialdef;
elseif isfield(job.trialchoice, 'arbitrary')
    S.trialength = job.trialchoice.arbitrary.trialength;
    S.conditionlabels = job.trialchoice.arbitrary.conditionlabel;
else
    trlfile = load(char(job.trialchoice.trlfile));
    usetrl  = 0;

    % In the new code trl file contains both trl matrix and trial definition
    % struct. trl is usually only applicable to the file on which it was
    % defined and there it must have priority because it can be manually
    % adjusted in the GUI. Otherwise trialdef has priority if present.
    if isfield(trlfile, 'trl')
        if ~all(isfield(trlfile, {'trialdef', 'timewin'}))
             usetrl = 1;
        elseif isfield(trlfile, 'source')
                D = spm_eeg_load(S.D);
                if isequal(D.fname, trlfile.source)
                    usetrl = 1;
                end
        end
    elseif ~all(isfield(trlfile, {'trialdef', 'timewin'}))
        error('The trial definition file could not be inetrpreted');
    end

    if usetrl
        S.trl = trlfile.trl;
        if isfield(trlfile, 'conditionlabels')
            S.conditionlabels = trlfile.conditionlabels;
        end
    else
        S.trialdef = trlfile.trialdef;
        S.timewin = trlfile.timewin;
    end
end

S.bc = job.bc;
S.prefix = job.prefix;
S.eventpadding = job.eventpadding;

out.D      = spm_eeg_epochs(S);
out.Dfname = {fullfile(out.D.path, out.D.fname)};

%==========================================================================
function dep = vout_eeg_epochs(job)
% Output is always in field "D", no matter how job is structured
dep = cfg_dep;
dep.sname = 'Epoched Data';
% reference field "D" from output
dep.src_output = substruct('.','D');
% this can be entered into any evaluated input
dep.tgt_spec   = cfg_findspec({{'strtype','e'}});

dep(2) = cfg_dep;
dep(2).sname = 'Epoched Datafile';
% reference field "Dfname" from output
dep(2).src_output = substruct('.','Dfname');
% this can be entered into any file selector
dep(2).tgt_spec   = cfg_findspec({{'filter','mat'}});

author: spm