text:
%% DEMO_febio_0065_clamp_tension_test
% Below is a demonstration for:
% 1) The creation of an FEBio model for clamped tensile testing
% 2) The use of multiple steps
% 4) Running an FEBio job with MATLAB
% 5) Importing FEBio results into MATLAB

%% Keywords
%
% * febio_spec version 3.0
% * febio, FEBio
% * tension, tensile
% * displacement control, displacement boundary condition
% * Hexahedral hex8
% * static, solid
% * hyperelastic, Ogden
% * displacement logfile
% * stress logfile

%%

clear; close all; clc;

%% Plot settings
fontSize=25;

%% Control parameters

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

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

%Specifying dimensions and number of elements
pointSpacing=1;
sampleWidth=10;
sampleThickness=2;
sampleClampedHeight=sampleWidth;
sampleGripGripHeight=sampleWidth.*2;
appliedLinearStrain=0.3;
clampCompressiveLinearStrain=0.3;

%Initial material parameter set
c1=1e-3;
m1=2;
k_factor=100;
k=c1*k_factor;

% FEA control settings
numTimeSteps=10; %Number of time steps desired
max_refs=25; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=12; %Optimum number of iterations
max_retries=25; %Maximum number of retires
dtmin=(1/numTimeSteps)/100; %Minimum time step size
dtmax=1/numTimeSteps; %Maximum time step size

min_residual=1e-20;
symmetric_stiffness=0;
runMode='external'; %'internal' or 'external'

%% Computing derived parameters
numElementsWidth=ceil(sampleWidth/pointSpacing);
numElementsWidth=numElementsWidth+iseven(numElementsWidth); %Force uneven so there is a middle element
numElementsThickness=ceil(sampleThickness/pointSpacing)+1;
numElementsGripGripHeight=ceil(sampleGripGripHeight/pointSpacing);
numElementsGripGripHeight=numElementsGripGripHeight+iseven(numElementsGripGripHeight); %Force uneven so there is a middle element
numElementsClampedHeight=ceil(sampleClampedHeight/pointSpacing);

clampCompressiveDisplacement=(sampleThickness.*clampCompressiveLinearStrain)/2;
clampTensionDisplacement=(sampleGripGripHeight.*appliedLinearStrain);

%% Creating strip region
% The region consists of three "boxes" which define the upper and lower
% clamped regions as well as the central region.

%Create box 1
boxDim=[sampleWidth sampleThickness sampleClampedHeight]; %Dimensions
boxEl=[numElementsWidth numElementsThickness numElementsClampedHeight]; %Number of elements
[box1]=hexMeshBox(boxDim,boxEl);
E1=box1.E;
V1=box1.V;
F1=box1.F;
Fb1=box1.Fb;
faceBoundaryMarker1=box1.faceBoundaryMarker;

%Create box 3 by copying the first
E3=E1;
V3=V1;
F3=F1;
Fb3=Fb1;
faceBoundaryMarker3=faceBoundaryMarker1;

%Shift first box up
V1(:,3)=V1(:,3)+sampleGripGripHeight/2+sampleClampedHeight/2;

%Shift third box down
V3(:,3)=V3(:,3)-sampleGripGripHeight/2-sampleClampedHeight/2;

%Create box 2
boxDim=[sampleWidth sampleThickness sampleGripGripHeight]; %Dimensions
boxEl=[numElementsWidth numElementsThickness numElementsGripGripHeight]; %Number of elements
[box2]=hexMeshBox(boxDim,boxEl);
E2=box2.E;
V2=box2.V;
F2=box2.F;
Fb2=box2.Fb;
faceBoundaryMarker2=box2.faceBoundaryMarker;

%% Merging box sets

%Join color data
faceBoundaryMarker_all=[faceBoundaryMarker1; faceBoundaryMarker2; faceBoundaryMarker3;];
faceBoundaryMarker_ind=[ones(size(Fb1,1),1);2*ones(size(Fb2,1),1); 3*ones(size(Fb3,1),1);];

%Join nodes, elements, and faces
V=[V1;V2;V3];
E=[E1;E2+size(V1,1);E3+size(V1,1)+size(V2,1)];
F=[F1;F2+size(V1,1);F3+size(V1,1)+size(V2,1)];
Fb=[Fb1;Fb2+size(V1,1);Fb3+size(V1,1)+size(V2,1)];

%Merge nodes
[F,V,ind1,ind2]=mergeVertices(F,V);
E=ind2(E);
Fb=ind2(Fb);

%%
% Plotting surface models
cFigure;
subplot(1,2,1); hold on;
title('Merged box sets','FontSize',fontSize);
gpatch(Fb,V,faceBoundaryMarker_all);
axisGeom(gca,fontSize);
colormap(gca,gjet(250)); icolorbar;

subplot(1,2,2); hold on;
title('Merged box sets','FontSize',fontSize);
gpatch(Fb,V,faceBoundaryMarker_ind);
axisGeom(gca,fontSize);
colormap(gca,gjet(250)); icolorbar;
drawnow;

%% Define clamping surfaces

logicContactSurf1=faceBoundaryMarker_all==3 & faceBoundaryMarker_ind==1;
Fc1=Fb(logicContactSurf1,:);

logicContactSurf2=faceBoundaryMarker_all==4 & faceBoundaryMarker_ind==1;
Fc2=Fb(logicContactSurf2,:);

logicContactSurf3=faceBoundaryMarker_all==4 & faceBoundaryMarker_ind==3;
Fc3=Fb(logicContactSurf3,:);

logicContactSurf4=faceBoundaryMarker_all==3 & faceBoundaryMarker_ind==3;
Fc4=Fb(logicContactSurf4,:);

%%
% Visualize clamping surfaces

cFigure; hold on;
title('Clamping surfaces','FontSize',fontSize);
gpatch(Fb,V,'kw','none',0.25);
hp(1)=gpatch(Fc1,V,'rw','k',1);
hp(2)=gpatch(Fc2,V,'gw','k',1);
hp(3)=gpatch(Fc3,V,'bw','k',1);
hp(4)=gpatch(Fc4,V,'yw','k',1);
legend(hp,{'Surf. 1','Surf. 2','Surf. 3','Surf. 4'})
axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Define BC's

bcPrescribeList1=unique(Fc1(:)); % Nodes of surface 1
bcPrescribeList2=unique(Fc2(:)); % Nodes of surface 2
bcPrescribeList3=unique(Fc3(:)); % Nodes of surface 3
bcPrescribeList4=unique(Fc4(:)); % Nodes of surface 4

%%
% Visualize boundary conditions

cFigure; hold on;
title('Complete model','FontSize',fontSize);

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

hp(1)=plotV(V(bcPrescribeList1,:),'r.','MarkerSize',25);
hp(2)=plotV(V(bcPrescribeList2,:),'g.','MarkerSize',25);
hp(3)=plotV(V(bcPrescribeList3,:),'b.','MarkerSize',25);
hp(4)=plotV(V(bcPrescribeList4,:),'y.','MarkerSize',25);
legend(hp,{'Node set 1','Node set 2','Node set 3','Node set 4'})
axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Get logic for middle elements (to help investigate strain in the middle)

searchTol=pointSpacing/1000;
VE=patchCentre(E,V);
logicMiddleElements= VE(:,3)<searchTol & VE(:,3)>-searchTol;

[FM]=element2patch(E(logicMiddleElements,:),[],'hex8');

%%
cFigure; hold on;
gpatch(Fb,V,'kw','none',0.25);

gpatch(FM,V,'rw','k',1);

% hp(1)=plotV(V(bcPrescribeList1,:),'r.','MarkerSize',25);
% legend(hp,{'Node set 1','Node set 2','Node set 3','Node set 4'})
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
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.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;

% Mesh section
% -> Nodes
febio_spec.Mesh.Nodes{1}.ATTR.name='Object1'; %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='hex8'; %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

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

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

nodeSetName3='bcPrescribeList3';
febio_spec.Mesh.NodeSet{3}.ATTR.name=nodeSetName3;
febio_spec.Mesh.NodeSet{3}.node.ATTR.id=bcPrescribeList3(:);

nodeSetName4='bcPrescribeList4';
febio_spec.Mesh.NodeSet{4}.ATTR.name=nodeSetName4;
febio_spec.Mesh.NodeSet{4}.node.ATTR.id=bcPrescribeList4(:);

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

%Boundary condition section
%STEP 1: Clamping compression
%Set 1
febio_spec.Step.step{1}.Boundary.bc{1}.ATTR.type='prescribe';
febio_spec.Step.step{1}.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
febio_spec.Step.step{1}.Boundary.bc{1}.dof='x';
febio_spec.Step.step{1}.Boundary.bc{1}.scale.ATTR.lc=1;
febio_spec.Step.step{1}.Boundary.bc{1}.scale.VAL=0;
febio_spec.Step.step{1}.Boundary.bc{1}.relative=1;

febio_spec.Step.step{1}.Boundary.bc{2}.ATTR.type='prescribe';
febio_spec.Step.step{1}.Boundary.bc{2}.ATTR.node_set=nodeSetName1;
febio_spec.Step.step{1}.Boundary.bc{2}.dof='y';
febio_spec.Step.step{1}.Boundary.bc{2}.scale.ATTR.lc=1;
febio_spec.Step.step{1}.Boundary.bc{2}.scale.VAL=clampCompressiveDisplacement;
febio_spec.Step.step{1}.Boundary.bc{2}.relative=1;

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

%Set 2
febio_spec.Step.step{1}.Boundary.bc{4}.ATTR.type='prescribe';
febio_spec.Step.step{1}.Boundary.bc{4}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{1}.Boundary.bc{4}.dof='x';
febio_spec.Step.step{1}.Boundary.bc{4}.scale.ATTR.lc=1;
febio_spec.Step.step{1}.Boundary.bc{4}.scale.VAL=0;
febio_spec.Step.step{1}.Boundary.bc{4}.relative=1;

febio_spec.Step.step{1}.Boundary.bc{5}.ATTR.type='prescribe';
febio_spec.Step.step{1}.Boundary.bc{5}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{1}.Boundary.bc{5}.dof='y';
febio_spec.Step.step{1}.Boundary.bc{5}.scale.ATTR.lc=1;
febio_spec.Step.step{1}.Boundary.bc{5}.scale.VAL=-clampCompressiveDisplacement;
febio_spec.Step.step{1}.Boundary.bc{5}.relative=1;

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

%Set 3
febio_spec.Step.step{1}.Boundary.bc{7}.ATTR.type='prescribe';
febio_spec.Step.step{1}.Boundary.bc{7}.ATTR.node_set=nodeSetName3;
febio_spec.Step.step{1}.Boundary.bc{7}.dof='x';
febio_spec.Step.step{1}.Boundary.bc{7}.scale.ATTR.lc=1;
febio_spec.Step.step{1}.Boundary.bc{7}.scale.VAL=0;
febio_spec.Step.step{1}.Boundary.bc{7}.relative=1;

febio_spec.Step.step{1}.Boundary.bc{8}.ATTR.type='prescribe';
febio_spec.Step.step{1}.Boundary.bc{8}.ATTR.node_set=nodeSetName3;
febio_spec.Step.step{1}.Boundary.bc{8}.dof='y';
febio_spec.Step.step{1}.Boundary.bc{8}.scale.ATTR.lc=1;
febio_spec.Step.step{1}.Boundary.bc{8}.scale.VAL=-clampCompressiveDisplacement;
febio_spec.Step.step{1}.Boundary.bc{8}.relative=1;

febio_spec.Step.step{1}.Boundary.bc{9}.ATTR.type='prescribe';
febio_spec.Step.step{1}.Boundary.bc{9}.ATTR.node_set=nodeSetName3;
febio_spec.Step.step{1}.Boundary.bc{9}.dof='z';
febio_spec.Step.step{1}.Boundary.bc{9}.scale.ATTR.lc=1;
febio_spec.Step.step{1}.Boundary.bc{9}.scale.VAL=0;
febio_spec.Step.step{1}.Boundary.bc{9}.relative=1;

%Set 4
febio_spec.Step.step{1}.Boundary.bc{10}.ATTR.type='prescribe';
febio_spec.Step.step{1}.Boundary.bc{10}.ATTR.node_set=nodeSetName4;
febio_spec.Step.step{1}.Boundary.bc{10}.dof='x';
febio_spec.Step.step{1}.Boundary.bc{10}.scale.ATTR.lc=1;
febio_spec.Step.step{1}.Boundary.bc{10}.scale.VAL=0;
febio_spec.Step.step{1}.Boundary.bc{10}.relative=1;

febio_spec.Step.step{1}.Boundary.bc{11}.ATTR.type='prescribe';
febio_spec.Step.step{1}.Boundary.bc{11}.ATTR.node_set=nodeSetName4;
febio_spec.Step.step{1}.Boundary.bc{11}.dof='y';
febio_spec.Step.step{1}.Boundary.bc{11}.scale.ATTR.lc=1;
febio_spec.Step.step{1}.Boundary.bc{11}.scale.VAL=clampCompressiveDisplacement;
febio_spec.Step.step{1}.Boundary.bc{11}.relative=1;

febio_spec.Step.step{1}.Boundary.bc{12}.ATTR.type='prescribe';
febio_spec.Step.step{1}.Boundary.bc{12}.ATTR.node_set=nodeSetName4;
febio_spec.Step.step{1}.Boundary.bc{12}.dof='z';
febio_spec.Step.step{1}.Boundary.bc{12}.scale.ATTR.lc=1;
febio_spec.Step.step{1}.Boundary.bc{12}.scale.VAL=0;
febio_spec.Step.step{1}.Boundary.bc{12}.relative=1;

%STEP 2 Tension
%Set 1
febio_spec.Step.step{2}.Boundary.bc{1}.ATTR.type='prescribe';
febio_spec.Step.step{2}.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
febio_spec.Step.step{2}.Boundary.bc{1}.dof='x';
febio_spec.Step.step{2}.Boundary.bc{1}.scale.ATTR.lc=2;
febio_spec.Step.step{2}.Boundary.bc{1}.scale.VAL=0;
febio_spec.Step.step{2}.Boundary.bc{1}.relative=1;

febio_spec.Step.step{2}.Boundary.bc{2}.ATTR.type='prescribe';
febio_spec.Step.step{2}.Boundary.bc{2}.ATTR.node_set=nodeSetName1;
febio_spec.Step.step{2}.Boundary.bc{2}.dof='y';
febio_spec.Step.step{2}.Boundary.bc{2}.scale.ATTR.lc=2;
febio_spec.Step.step{2}.Boundary.bc{2}.scale.VAL=0;
febio_spec.Step.step{2}.Boundary.bc{2}.relative=1;

febio_spec.Step.step{2}.Boundary.bc{3}.ATTR.type='prescribe';
febio_spec.Step.step{2}.Boundary.bc{3}.ATTR.node_set=nodeSetName1;
febio_spec.Step.step{2}.Boundary.bc{3}.dof='z';
febio_spec.Step.step{2}.Boundary.bc{3}.scale.ATTR.lc=2;
febio_spec.Step.step{2}.Boundary.bc{3}.scale.VAL=clampTensionDisplacement;
febio_spec.Step.step{2}.Boundary.bc{3}.relative=1;

%Set 2
febio_spec.Step.step{2}.Boundary.bc{4}.ATTR.type='prescribe';
febio_spec.Step.step{2}.Boundary.bc{4}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{2}.Boundary.bc{4}.dof='x';
febio_spec.Step.step{2}.Boundary.bc{4}.scale.ATTR.lc=2;
febio_spec.Step.step{2}.Boundary.bc{4}.scale.VAL=0;
febio_spec.Step.step{2}.Boundary.bc{4}.relative=1;

febio_spec.Step.step{2}.Boundary.bc{5}.ATTR.type='prescribe';
febio_spec.Step.step{2}.Boundary.bc{5}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{2}.Boundary.bc{5}.dof='y';
febio_spec.Step.step{2}.Boundary.bc{5}.scale.ATTR.lc=2;
febio_spec.Step.step{2}.Boundary.bc{5}.scale.VAL=0;
febio_spec.Step.step{2}.Boundary.bc{5}.relative=1;

febio_spec.Step.step{2}.Boundary.bc{6}.ATTR.type='prescribe';
febio_spec.Step.step{2}.Boundary.bc{6}.ATTR.node_set=nodeSetName2;
febio_spec.Step.step{2}.Boundary.bc{6}.dof='z';
febio_spec.Step.step{2}.Boundary.bc{6}.scale.ATTR.lc=2;
febio_spec.Step.step{2}.Boundary.bc{6}.scale.VAL=clampTensionDisplacement;
febio_spec.Step.step{2}.Boundary.bc{6}.relative=1;

%Set 3
febio_spec.Step.step{2}.Boundary.bc{7}.ATTR.type='prescribe';
febio_spec.Step.step{2}.Boundary.bc{7}.ATTR.node_set=nodeSetName3;
febio_spec.Step.step{2}.Boundary.bc{7}.dof='x';
febio_spec.Step.step{2}.Boundary.bc{7}.scale.ATTR.lc=2;
febio_spec.Step.step{2}.Boundary.bc{7}.scale.VAL=0;
febio_spec.Step.step{2}.Boundary.bc{7}.relative=1;

febio_spec.Step.step{2}.Boundary.bc{8}.ATTR.type='prescribe';
febio_spec.Step.step{2}.Boundary.bc{8}.ATTR.node_set=nodeSetName3;
febio_spec.Step.step{2}.Boundary.bc{8}.dof='y';
febio_spec.Step.step{2}.Boundary.bc{8}.scale.ATTR.lc=2;
febio_spec.Step.step{2}.Boundary.bc{8}.scale.VAL=0;
febio_spec.Step.step{2}.Boundary.bc{8}.relative=1;

febio_spec.Step.step{2}.Boundary.bc{9}.ATTR.type='prescribe';
febio_spec.Step.step{2}.Boundary.bc{9}.ATTR.node_set=nodeSetName3;
febio_spec.Step.step{2}.Boundary.bc{9}.dof='z';
febio_spec.Step.step{2}.Boundary.bc{9}.scale.ATTR.lc=2;
febio_spec.Step.step{2}.Boundary.bc{9}.scale.VAL=0;
febio_spec.Step.step{2}.Boundary.bc{9}.relative=1;

%Set 4
febio_spec.Step.step{2}.Boundary.bc{10}.ATTR.type='prescribe';
febio_spec.Step.step{2}.Boundary.bc{10}.ATTR.node_set=nodeSetName4;
febio_spec.Step.step{2}.Boundary.bc{10}.dof='x';
febio_spec.Step.step{2}.Boundary.bc{10}.scale.ATTR.lc=2;
febio_spec.Step.step{2}.Boundary.bc{10}.scale.VAL=0;
febio_spec.Step.step{2}.Boundary.bc{10}.relative=1;

febio_spec.Step.step{2}.Boundary.bc{11}.ATTR.type='prescribe';
febio_spec.Step.step{2}.Boundary.bc{11}.ATTR.node_set=nodeSetName4;
febio_spec.Step.step{2}.Boundary.bc{11}.dof='y';
febio_spec.Step.step{2}.Boundary.bc{11}.scale.ATTR.lc=2;
febio_spec.Step.step{2}.Boundary.bc{11}.scale.VAL=0;
febio_spec.Step.step{2}.Boundary.bc{11}.relative=1;

febio_spec.Step.step{2}.Boundary.bc{12}.ATTR.type='prescribe';
febio_spec.Step.step{2}.Boundary.bc{12}.ATTR.node_set=nodeSetName4;
febio_spec.Step.step{2}.Boundary.bc{12}.dof='z';
febio_spec.Step.step{2}.Boundary.bc{12}.scale.ATTR.lc=2;
febio_spec.Step.step{2}.Boundary.bc{12}.scale.VAL=0;
febio_spec.Step.step{2}.Boundary.bc{12}.relative=1;

%LoadData section
%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];

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.element_data{1}.ATTR.file=febioLogFileName_strain;
febio_spec.Output.logfile.element_data{1}.ATTR.data='Ez';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.element_data{1}.VAL=1: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;
febioAnalysis.maxLogCheckTime=120; %Max log file checking time

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

%% Import FEBio results

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

    %%
    % Importing nodal displacements from a log file

    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 strain data from a log file
    dataStruct=importFEBio_logfile(fullfile(savePath,febioLogFileName_strain),1,1);

    %Access data
    E_strain=dataStruct.data;

    [F,CF]=element2patch(E,E_strain(:,:,end));
    CF_V=faceToVertexMeasure(F,V,CF);

    E_strain_middle_mean=squeeze(mean(E_strain(logicMiddleElements,:,:),1));

    %% Compute grimp implied strain
    gripStrainLinear=double(timeVec>1).*(timeVec-1).*appliedLinearStrain;
    gripStrain=1/2*((gripStrainLinear+1).^2-1); %Green-Lagrange strain
    maxGridStrain=max(abs(gripStrain));
      %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    axLim=[min(min(V_DEF,[],3),[],1); max(max(V_DEF,[],3),[],1)];
    indBc=[bcPrescribeList1;bcPrescribeList2;bcPrescribeList3;bcPrescribeList4;];

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; hold on;%Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    title('$E_{zz}$','Interpreter','Latex')
    hp1=gpatch(Fb,V_DEF(:,:,end),CF_V,'k',1);
    hp1.FaceColor='interp';
    hp2=plotV(V(indBc,:),'k.','MarkerSize',25);

    axisGeom(gca,fontSize);
    colormap(warmcold(250)); hc=colorbar;
    caxis([-maxGridStrain maxGridStrain]);
    hc.Ticks=linspace(-maxGridStrain,maxGridStrain,7);
    axis(axLim(:)'); %Set axis limits statically

    camlight headlight; axis off;

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

        [~,CF]=element2patch(E,E_strain(:,:,qt));
        CF_V=faceToVertexMeasure(F,V,CF);

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp1 hp1 hp2 hp2 hp2]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData','XData','YData','ZData'}; %Properties of objects to animate
        animStruct.Set{qt}={V_def,CF_V,V_def(indBc,1),V_def(indBc,2),V_def(indBc,3)}; %Property values for to set in order to animate
    end
    anim8(hf,animStruct); %Initiate animation feature
    gdrawnow;

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

    indBc=[bcPrescribeList1;bcPrescribeList2;bcPrescribeList3;bcPrescribeList4;];

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure;
    subplot(1,2,1); hold on;
    title('$E_{zz}$','Interpreter','Latex');
    gtitle([febioFebFileNamePart,': Press play to animate']);
    hp1=gpatch(Fb,V_DEF(:,:,end),CF_V,'none',1,0.5); hp1.FaceColor='interp';
    hp2=plotV(V_DEF(indBc,:,end),'k.','MarkerSize',1);
    hp3=gpatch(FM,V_DEF(:,:,end),CF_V,'k',1,2); hp3.FaceColor='interp';

    colormap(warmcold(250)); hc=colorbar;
    caxis([-maxGridStrain maxGridStrain]);
    hc.Ticks=linspace(-maxGridStrain,maxGridStrain,7);
    axisGeom(gca,fontSize); camlight headlight;
    axis(axisLim(V_DEF)); %Set axis limits statically
    view(0,0); axis off;

    subplot(1,2,2); hold on;
    xlabel('Time (s)'); ylabel('$E_{zz}$','Interpreter','Latex');
    hpl4=plot(timeVec,E_strain_middle_mean,'g.-','LineWidth',3);
    hpl5=plot(timeVec,gripStrain,'r.-','LineWidth',3);
    hp4=plot(timeVec(end),E_strain_middle_mean(end),'g.','MarkerSize',50);
    hp5=plot(timeVec(end),gripStrain(end),'r.','MarkerSize',50);
    legend([hpl4 hpl5],{'True strain $E_{zz}$','"Intended applied" strain $E_{zz}$'},'Location','NorthOutside','Interpreter','Latex');
    axis tight; box on; grid on; set(gca,'FontSize',fontSize);

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

        [~,CF]=element2patch(E,E_strain(:,:,qt));
        CF_V=faceToVertexMeasure(F,V,CF);

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp1 hp1 hp2 hp2 hp2 hp3 hp3 hp4 hp4 hp5 hp5]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData','XData','YData','ZData','Vertices','CData','XData','YData','XData','YData'}; %Properties of objects to animate
        animStruct.Set{qt}={V_def,CF_V,V_def(indBc,1),V_def(indBc,2),V_def(indBc,3),V_def,CF_V,timeVec(qt),E_strain_middle_mean(qt),timeVec(qt),gripStrain(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


function [tracker, qscore, f] = track_frame(tracker, frame_id, frame_image, bboxes_det, index_det, seq_name, opt)

    global client_tcp;

    % tracked, decide to tracked or occluded
    if tracker.state == opt.STATE_TRACKED || tracker.state == opt.STATE_ACTIVATED
        tracker = ECO_tracking(frame_id, frame_image, bboxes_det, tracker, opt);
        is_tracking_score_good = tracker.eco.is_confident;
        is_overlap_ratio_good = mean(tracker.eco.bb_overlaps);
        tracker.bb = tracker.eco.bb';
        qscore = 0;
        if is_tracking_score_good == 1
            if seq_name == 'MOT16-03'
                overlap_threshold = 0.1;
            else
                overlap_threshold = 0.8;
            end
            if tracker.state == opt.STATE_TRACKED && is_overlap_ratio_good > overlap_threshold
                label = 1;
            elseif tracker.state == opt.STATE_ACTIVATED %&& is_overlap_ratio_good > 0.8
                label = 1;
            else
                label = -1;
            end
        else
            label = -1;
        end

        % make a decision
        if label > 0
            if tracker.num_tracked > (0.25 * tracker.fps)
                tracker.state = opt.STATE_TRACKED;
            else
                tracker.state = opt.STATE_ACTIVATED;
            end
        else
            if tracker.state == opt.STATE_ACTIVATED
                tracker.state = opt.STATE_STOP;
                fprintf('Target %d end.\n', tracker.target_id);
            else
                tracker.state = opt.STATE_LOST;
            end
        end

        % build the bboxes structure new
        if tracker.state == opt.STATE_TRACKED || tracker.state == opt.STATE_ACTIVATED
            bbox.x = tracker.bb(1);
            bbox.y = tracker.bb(2);
            bbox.w = tracker.bb(3) - tracker.bb(1);
            bbox.h = tracker.bb(4) - tracker.bb(2);
            bbox.r = 1;
            bbox.eco_x = tracker.eco.bb(1);
            bbox.eco_y = tracker.eco.bb(2);
            bbox.eco_w = tracker.eco.bb(3) - tracker.eco.bb(1);
            bbox.eco_h = tracker.eco.bb(4) - tracker.eco.bb(2);

            if ~isempty(bboxes_det.fr)
                o = calc_overlap(bbox, 1, bboxes_det, 1:numel(bboxes_det.fr));
                tracker.eco.bb_overlaps(1:end-1) = tracker.eco.bb_overlaps(2:end);
                if max(o) > 0.5
                    tracker.eco.bb_overlaps(end) = 1;
                else
                    tracker.eco.bb_overlaps(end) = 0;
                end
            else
                tracker.eco.bb_overlaps(1:end-1) = tracker.eco.bb_overlaps(2:end);
                tracker.eco.bb_overlaps(end) = 0;
            end
        else
            bbox = sub_bboxes(tracker.bboxes, numel(tracker.bboxes.fr));
        end
        if tracker.state == opt.STATE_ACTIVATED && tracker.eco.bb_overlaps(end) == 0
            tracker.state = opt.STATE_STOP;
        end
        bbox.fr = frame_id;
        bbox.id = tracker.target_id;
        bbox.state = tracker.state;
        tracker.bboxes = concatenate_bboxes(tracker.bboxes, bbox);

    % occluded, decide to tracked or occluded
    elseif tracker.state == opt.STATE_LOST
        % association
        if isempty(index_det)
            qscore = 0;
            label = -1;
            f = [];
        else
            bboxes = sub_bboxes(bboxes_det, index_det);
            ctrack = motion_predict(frame_id, tracker, opt);
            bboxes_predict.w = tracker.bboxes.w(end);
            bboxes_predict.h = tracker.bboxes.h(end);
            bboxes_predict.x = ctrack(1) - (bboxes_predict.w - 1) / 2;
            bboxes_predict.y = ctrack(2) - (bboxes_predict.h - 1) / 2;
            o_predict = calc_overlap(bboxes_predict, 1, bboxes, 1:numel(bboxes.fr));
            [motion_score, motion_ind] = max(o_predict);
            bboxes_tracker = tracker.bboxes;
            index_state = find(bboxes_tracker.state == opt.STATE_TRACKED | bboxes_tracker.state == opt.STATE_ACTIVATED);
            bboxes_tracker = sub_bboxes(bboxes_tracker, index_state);
            if ~isempty(bboxes_tracker.fr)
                fr_tracker = bboxes_tracker.fr(end);
            else % for target which is suppressed since the first frame it appears
                fr_tracker = 1;
            end
            traj_dir = ['img_traj/' seq_name '/' num2str(tracker.target_id) '/'];
            frame_id_double = double(frame_id);
            save('mot_py.mat', 'traj_dir', 'bboxes', 'frame_id_double', 'seq_name');

            fwrite(client_tcp, 'client ok');
            fread(client_tcp, 9); % size is 9 for 'server ok' message
            load('similarity.mat');

            [ass_score, ind] = max(similarity);
            index_tracked = find(tracker.bboxes.state == opt.STATE_TRACKED | tracker.bboxes.state == opt.STATE_ACTIVATED);
            bboxes_tracked = sub_bboxes(tracker.bboxes, index_tracked);
            fr_tracked = bboxes_tracked.fr(end);
            if frame_id - fr_tracker == 1
                prev_c = [bboxes_tracked.x(end) + bboxes_tracked.w(end) / 2, bboxes_tracked.y(end) + bboxes_tracked.h(end) / 2];
                prev_w = bboxes_tracked.w(end);
                prev_h = bboxes_tracked.h(end);
                bboxes_det_one = sub_bboxes(bboxes_det, index_det(ind));
                cur_c = [bboxes_det_one.x + bboxes_det_one.w/2, bboxes_det_one.y + bboxes_det_one.h/2];
                dis = norm(cur_c - prev_c) / prev_w;
                ratio = bboxes_det_one.h / prev_h;
                ratio = min(ratio, 1/ratio);
            end

            if ass_score > opt.association_score_thre
                label = 1;
                bboxes_det_one = sub_bboxes(bboxes_det, index_det(ind));
                fprintf('Target %d associated by appearance.\n', tracker.target_id);

            elseif frame_id - fr_tracker < 5 && ass_score > 0.4 && motion_score > 0.5
                label = 1;
                bboxes_det_one = sub_bboxes(bboxes, motion_ind);
                fprintf('Target %d associated by motion.\n', tracker.target_id);
            else
                label = -1;
            end

            if label == 1
                tracker.eco.rect_position = double(round([bboxes_det_one.x, bboxes_det_one.y, bboxes_det_one.w, bboxes_det_one.h]));
                tracker.eco.pos = single([bboxes_det_one.y + (bboxes_det_one.h - 1)/2, bboxes_det_one.x + (bboxes_det_one.w - 1)/2]);
                tracker.eco.target_sz = double([bboxes_det_one.h, bboxes_det_one.w]);
                tracker.eco.base_target_sz = double(tracker.eco.target_sz / tracker.eco.currentScaleFactor);
                tracker.eco.base_target_sz(2) = tracker.eco.base_target_sz(2)* 0.5;
                tracker = ECO_tracking(frame_id, frame_image, bboxes_det_one, tracker, opt);
            end
        end

        if label > 0
            tracker.bb = tracker.eco.bb';

            % association
            if tracker.num_tracked > 0.25 * tracker.fps
                tracker.state = opt.STATE_TRACKED;
            else
                tracker.state = opt.STATE_ACTIVATED;
            end

            % build the bboxes structure
            bbox = [];
            bbox.fr = frame_id;
            bbox.id = tracker.target_id;
            bbox.x = tracker.bb(1);
            bbox.y = tracker.bb(2);
            bbox.w = tracker.bb(3) - tracker.bb(1);
            bbox.h = tracker.bb(4) - tracker.bb(2);
            bbox.r = 1;
            bbox.state = tracker.state;
            bbox.eco_x = tracker.eco.bb(1);
            bbox.eco_y = tracker.eco.bb(2);
            bbox.eco_w = tracker.eco.bb(3) - tracker.eco.bb(1);
            bbox.eco_h = tracker.eco.bb(4) - tracker.eco.bb(2);

            if tracker.bboxes.fr(end) == frame_id
                bboxes = tracker.bboxes;
                index = 1:numel(bboxes.fr)-1;
                tracker.bboxes = sub_bboxes(bboxes, index);
            end

            tracker.bboxes = interpolate_traj(tracker.bboxes, bbox, tracker.fps, opt);

            if isempty(bboxes_det.fr) == 0
                o = calc_overlap(bbox, 1, bboxes_det, 1:numel(bboxes_det.fr));
                tracker.eco.bb_overlaps(1:end-1) = tracker.eco.bb_overlaps(2:end);
                if max(o) > 0.5
                    tracker.eco.bb_overlaps(end) = 1;
                else
                    tracker.eco.bb_overlaps(end) = 0;
                end
            else
                tracker.eco.bb_overlaps(1:end-1) = tracker.eco.bb_overlaps(2:end);
                tracker.eco.bb_overlaps(end) = 0;
            end
        else
            % no association
            tracker.state = opt.STATE_LOST;
            bbox = sub_bboxes(tracker.bboxes, numel(tracker.bboxes.fr));
            bbox.fr = frame_id;
            bbox.id = tracker.target_id;
            bbox.state = tracker.state;

            if tracker.bboxes.fr(end) == frame_id
                bboxes = tracker.bboxes;
                index = 1:numel(bboxes.fr)-1;
                tracker.bboxes = sub_bboxes(bboxes, index);
            end
            tracker.bboxes = concatenate_bboxes(tracker.bboxes, bbox);
        end
    end


text:
function ft_realtime_asaproxy(cfg)

% FT_REALTIME_ASAPROXY reads continuous data from the ASA acquisition system and
% writes it to a FieldTrip buffer. This function uses the NeuroSDK software, which
% can be obtained from ANT.
%
% The FieldTrip buffer is a network transparent server that allows the acquisition
% client to stream data to it. An analysis client can connect to read the data upon
% request. Multiple clients can connect simultaneously, each analyzing a specific
% aspect of the data concurrently.
%
% Use as
%   ft_realtime_asaproxy(cfg)
%
% The configuration should contain
%   cfg.channel              = cell-array, see FT_CHANNELSELECTION (default = 'all')
%
% The target to write the data to is configured as
%   cfg.target.datafile      = string, target destination for the data (default = 'buffer://localhost:1972')
%   cfg.target.dataformat    = string, default is determined automatic
%
% To stop this realtime function, you have to press Ctrl-C
%
% See also FT_REALTIME_SIGNALPROXY, FT_REALTIME_SIGNALVIEWER

% Copyright (C) 2009, 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/>.
%
% $Id$

% set the defaults
if ~isfield(cfg, 'target'),               cfg.target = [];                                  end
if ~isfield(cfg.target, 'headerformat'),  cfg.target.headerformat = [];                     end % default is detected automatically
if ~isfield(cfg.target, 'dataformat'),    cfg.target.dataformat = [];                       end % default is detected automatically
if ~isfield(cfg.target, 'datafile'),      cfg.target.datafile = 'buffer://localhost:1972';  end
if ~isfield(cfg, 'blocksize'),            cfg.blocksize = 1;                                end % in seconds
% if ~isfield(cfg, 'minblocksize'),         cfg.minblocksize = 0;                             end % in seconds
% if ~isfield(cfg, 'maxblocksize'),         cfg.maxblocksize = 1;                             end % in seconds
% if ~isfield(cfg, 'channel'),              cfg.channel = 'all';                              end

sHost               = 'asa';       % change this to 'asa', 'eemagine', 'tmsi' or 'test' other work with other hosts
bDirectAmp          = true;

% connect to amplifier server
if (strcmpi(sHost, 'asa') || strcmpi(sHost, 'eemagine'))
  bDirectAmp   = false;
  oAcquisition =  actxserver('Eemagine.Acquisition');
end

oDevice = device(sHost);
oDevice = connect(oDevice);

try
  pause(1);
  count = 0;

  while true

    % get new EEG data
    dat = getEEG(oDevice, cfg.blocksize);
    dat = dat';

    % FIXME this data is the last N seconds, which means that there will be overlap
    % between subsequent blocks that are processed by this function

    count = count + 1;
    nchans = size(dat,1);
    nsamples = size(dat,2);
    fprintf('processing segment %d, %d nchans, %d nsamples\n', count, nchans, nsamples);

    hdr = [];
    hdr.Fs = (nsamples-1)/cfg.blocksize;
    hdr.nChans = nchans;
    hdr.nSamples = nsamples;
    % hdr.nSamplesPre = 0;
    % hdr.label = {};
    % for i=1:nchans
    %   hdr.label{i} = sprintf('%d', i);
    % end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % from here onward it is specific to writing the data to another stream
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if count==1
      % flush the file, write the header and subsequently write the data segment
      ft_write_data(cfg.target.datafile, dat, 'header', hdr, 'dataformat', cfg.target.dataformat, 'append', false);
    else
      % write the data segment
      ft_write_data(cfg.target.datafile, dat, 'header', hdr, 'dataformat', cfg.target.dataformat, 'append', true);
    end

    % pause for some time
    pause(cfg.blocksize);
  end
catch
  % disconnect the device
  oDevice = disconnect(oDevice);
  disp(lasterr);
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Copyright (C) 2014-2016, 2019, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @deftypemethod  @@sym {} latex (@var{x})
%% @deftypemethodx @@sym {@var{s} =} latex (@var{x})
%% Display or return LaTeX typesetting code for symbolic expression.
%%
%% Example:
%% @example
%% @group
%% syms x
%% latex(sin(x/2))
%%   @print{} \sin@{\left(\frac@{x@}@{2@} \right)@}
%% @end group
%%
%% @group
%% A = [sym(1) 2; sym(3) 4];
%% s = latex(A)
%%   @result{} s = \left[\begin@{matrix@}1 & 2\\3 & 4\end@{matrix@}\right]
%% @end group
%% @end example
%%
%% @seealso{@@sym/disp, @@sym/pretty}
%% @end deftypemethod

function varargout = latex(x)

  if (nargin ~= 1)
    print_usage ();
  end

  cmd = { 'return sp.latex(*_ins),' };

  s = pycall_sympy__ (cmd, x);

  if (nargout == 0)
    disp(s)
  else
    varargout = {s};
  end

end

%!test
%! syms x
%! y = sin(x);
%! assert (strcmp (latex (y), '\sin{\left(x \right)}'))

%!assert (strcmp (latex (exp (sym('x'))), 'e^{x}'))

repo: octsympy


text:
function [ C , coord , MN ] = ccv ( A , N , F )
% Generates a crystal consists of lattice sites and bases, Calculates
% connecting vectors from point to point, Finds the number of nearest
% neighbors of each atomic site.
%
% function [ C , coord , MN ] = ccv ( A , N , F )
%
% arguments: ( input )
%
%  A - ( class - double ) a ( 3 * 3 ) matrix that each row is a vector that
%  shows primary generator vectors of a lattice in three dimensions.
%
%  N - ( class - double ) a matrix with three possitive integers that shows
%  the number of atoms in each dimensions.
%
%  F - ( class - double ) a ( m * 3 ) matrix that show the bases
%  coordinates. ( m is the number of bases )
%
% arguments: ( output )
%
%  C - ( class - double ) A ( N ( 1 ) * N ( 2 ) * N ( 3 ) ) * ( N ( 1 ) *
%  N ( 2 ) * N ( 3 ) ) * 3 matrix connecting vector between each two atomic
%  sites
%
%  coord - ( class - double ) Coordination number ( Number of nearest
%  neighbors ) of each atomic site
%
%  MN - ( class - double ) Mean distance between the nearest neighbor of
%  each atomic site
%
% Example:
%  A = [ 2 3 3 ; 6 5 3 ; 2 5 3 ] ;
%  N = [ 4 3 5 ] ;
%  F = [ 0 0 0 ; 0 0 1 ; 0 -2 0 ] ;
%  [ CV , COORD , MN ] = ccv ( A , N , F )
%
% See also rcv.
%
% Copyright 2009
%
% Release Date: 2009-10-12

% check for simple errors

if nargin < 3
    F = [ 0 0 0 ] ;
end % end of if loop

if nargin < 2
    N = [ 3 3 3 ] ;
end % end of if loop

if nargin < 1
    A = [ 1 0 0 ; 0 1 0 ; 0 0 1 ] ;
end % end of if loop

if ( N ( 1 ) <= 0 ) || ( N ( 2 ) <= 0 ) || ( N ( 3 ) <= 0 ) % condition of negative numbers
    error ' Enterd demensions must be positive integers. ' % error message
end % end of if loop

V = dot ( A ( 1 , : ) , cross ( A ( 2 , : ) , A ( 3 , : ) ) ) ; % Primitive Cell Volume

if V == 0 % condition of same plane vectors
    error ' Vectors must not be in same plane. ' % error message
end % end of if loop

% end of error checking

n = ( N ( 1 ) * N ( 2 ) ) * N ( 3 ) ; % calculates the total number of atomic sites
xy = zeros ( n , 3 ) ; % Preallocating

for i = 0 : N ( 1 ) - 1 % i is the numerator of for loop
    for j = 0 : N ( 2 ) - 1 % j is the numerator of for loop
        for k = 0 : N ( 3 ) - 1 % k is the numerator of for loop
            xy ( ( ( N ( 1 ) * j ) + ( i + 1 ) ) + ( N ( 1 ) * N ( 2 ) ) * k , : ) = A ( 1 , : ) * i + A ( 2 , : ) * j + A( 3 , : ) * k ;
        end % end of for loop
    end % end of for loop
end % end of for loop

xyold = xy ;
f = size ( F ) ;
xyo = [ ] ;

for q = 1 : f ( 1 ) % q is the numerator of for loop
    for p = 1 : 3 % p is the numerator of for loop
        xy ( : , p ) = xyold ( : , p ) + F ( q , p ) ;
    end % end of for loop
    xyo = [ xyo ; xy ] ;
end % end of for loop

n = length ( xyo ) ;
xy = xyo ;

C = zeros ( n , n , 3 ) ; % Preallocating
D = zeros ( n , n , 3 ) ; % Preallocating
E = zeros ( n , 3 ) ; % Preallocating
coord = zeros ( 1 , n ) ; % Preallocating
MIN = zeros ( n-1 , n ); % Preallocating
k = 1 ; % Preallocating

for P = 1 : n % P is the numerator of for loop
    for j = 1 : 3 % j is the numerator of for loop
        C ( P , 1 : n , j ) = xy ( 1 : n , j ) - xy ( P , j ) ; % Connecting vector elements between to lattice atomic sites
        D ( P , 1 : n , j ) = ( xy ( 1 : n , j ) - xy ( P , j ) ) .^ 2 ;
    end % end of for loop
    E ( P , 1 : n ) = sqrt ( D ( P , 1 : n , 1 ) + D ( P , 1 : n , 2 ) + D ( P , 1 : n , 3 ) ) ;
end % end of for loop

for P = 1 : n % P is the numerator of for loop
    for i = 1 : n % i is the numerator of for loop
        if E ( P , i ) ~= 0
            MIN ( k ) = E ( P , i ) ; % Eliminates zero value elements
            k = k + 1 ;
        end % end of if loop
    end % end of for loop
end % end of for loop

MN = min ( MIN ) ; % Nearest neighbors distance

for P = 1 : n % P is the numerator of for loop
    t = 0 ; % Preallocating
    for i = 1 : n % i is the numerator of for loop
        if E ( P , i ) == MN ( P )
            t = t + 1 ; % Counts number of nearest neighbors
        end % end of if loop
    end % end of for loop
    coord ( P ) = t ;
end % end of for loop

% With special thanks to John D'Errico
% By Ali Mohammad Razeghi
% My Email ( am_razeghi@yahoo.com ) is also ready to get full detailed
% commentation.

author: Sable

repo: mcbench-benchmarks


text:
function [ n_data_new, n, c ] = phi_values ( n_data )

%*****************************************************************************80
%
%% PHI_VALUES returns some values of the PHI function.
%
%  Definition:
%
%    PHI(N) is the number of integers between 1 and N which are
%    relatively prime to N.  I and J are relatively prime if they
%    have no common factors.  The function PHI(N) is known as
%    "Euler's totient function".
%
%    By convention, 1 and N are relatively prime.
%
%  First values:
%
%     N  PHI(N)
%
%     1    1
%     2    1
%     3    2
%     4    2
%     5    4
%     6    2
%     7    6
%     8    4
%     9    6
%    10    4
%    11   10
%    12    4
%    13   12
%    14    6
%    15    8
%    16    8
%    17   16
%    18    6
%    19   18
%    20    8
%
%  Formula:
%
%    PHI(U*V) = PHI(U) * PHI(V) if U and V are relatively prime.
%
%    PHI(P**K) = P**(K-1) * ( P - 1 ) if P is prime.
%
%    PHI(N) = N * Product ( P divides N ) ( 1 - 1 / P )
%
%    N = Sum ( D divides N ) PHI(D).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%  Parameters:
%
%    Input, integer N_DATA, indicates the index of the previous test data
%    returned, or is 0 if this is the first call.  For repeated calls,
%    set the input value of N_DATA to the output value of N_DATA_NEW
%    from the previous call.
%
%    Output, integer N_DATA_NEW, the index of the test data.
%
%    Output, integer N, the argument of the PHI function.
%
%    Output, integer C, the value of the PHI function.
%
  n_max = 20;
  c_vec = [ ...
      1,   1,   2,   2,   4,   2,   6,   4,   6,   4, ...
      8,   8,  16,  20,  16,  40, 148, 200, 200, 648 ];
  n_vec = [ ...
      1,   2,   3,   4,   5,   6,   7,   8,   9,  10, ...
     20,  30,  40,  50,  60, 100, 149, 500, 750, 999 ];

  n_data_new = n_data;

  if ( n_data_new < 0 )
    n_data_new = 0;
  end

  n_data_new = n_data_new + 1;

  if ( n_max < n_data_new )
    n_data_new = 0;
    n = 0;
    c = 0;
  else
    n = n_vec(n_data_new);
    c = c_vec(n_data_new);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function calpak_test337 ( )

%*****************************************************************************80
%
%% CALPAK_TEST337 tests MONTH_EG_LUNAR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 March 2013
%
%  Author:
%
%    John Burkardt
%
  n_test = 2;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST337\n' );
  fprintf ( 1, '  For the Egyptian Lunar calendar,\n' );
  fprintf ( 1, '  MONTH_LENGTH_EG_LUNAR returns month lengths.\n' );

  y_test(1) = 1;
  y_test(2) = 2;

  for i_test = 1 : n_test

    y = y_test(i_test);
    sy = y_to_s_eg_lunar ( y );
    months = year_length_months_eg_lunar ( y );
    days = year_length_eg_lunar ( y );

    fprintf ( 1, '\n' );
    fprintf ( 1, '  %d\n', y );
    fprintf ( 1, '  %s\n', sy );
    fprintf ( 1, '  Year length in months = %d\n', months );
    fprintf ( 1, '  Year length in days = %d\n', days );
    fprintf ( 1, '\n' );

    for m = 1 : months
      month_name = month_to_month_name_eg_lunar ( m );
      fprintf ( 1, '  %10s  %2d\n', month_name, month_length_eg_lunar ( y, m ) );
    end

  end

  return
end


text:
function [ hdls ] = tracks_viewer( dates, track_data, background_data, ...
    pixel_dir, track_indexes, eddy_attribute_data, eddy_attribute_dir, eddy_plotting_data, chelton_eddies )
%TRACKS_VIEWER Start the eddy viewers with the following components
%   Current date box: display current date that the data will be plotted
%   Eddy type box: what type of eddies that will be plotted
%   Track type: what type of tracks will be plotted
%   Track length: minimum and maximum length of the tracks which will be plotted
%   Backgrond panel:
%       Dropdown box: select what background will be plotted
%       Min value, max value edit box: min and max values of the background data which will be plotted
%       Auto panel: auto adjusting background values when zoomed in/out
%   Show contours: whether or not showing eddies' contour
%   Show hurricanes: whether or not showing hurricanes
%   Select eddy button: select an eddy on the map to see more information of that eddy

hdls = ViewerHandle();
hdls.dates = dates;
hdls.track_data = track_data;
hdls.pixel_dir = pixel_dir;
hdls.track_indexes = track_indexes;
hdls.attribute_dir = eddy_attribute_dir;
hdls.eddy_attribute_data = eddy_attribute_data;
hdls.chelton_eddies = chelton_eddies;

hdls.constants.eddy_types = eddy_plotting_data.eddy_types;
hdls.constants.eddy_markers = eddy_plotting_data.eddy_markers;
hdls.constants.eddy_colors = eddy_plotting_data.eddy_colors;
hdls.constants.eddy_type_names = eddy_plotting_data.eddy_type_names;

hdls.constants.PAST_TRACK_COLOR = [0 1 0];
hdls.constants.CYC_FUTURE_TRACK_COLOR = [0 0 1];
hdls.constants.ANT_FUTURE_TRACK_COLOR = [1 0 0];
hdls.constants.CHELTON_CYC_FUTURE_TRACK_COLOR = 'k';
hdls.constants.CHELTON_ANT_FUTURE_TRACK_COLOR = 'm';
hdls.constants.CHELTON_CURRENT_TRACK_COLOR = 'k';
hdls.constants.CHELTON_MARKER = 'd';

hdls.constants.CURRENT_EDDY_SIZE = 9;
hdls.constants.OTHER_EDDY_SIZE = 6;

hdls.constants.COAST_COLOR = [205 133 63] / 255;

horizontalBarFontSize = 12;

%% Initialize GUI Components
hdls.fig = figure('Name', 'Tracks Viewer', ...
    'NumberTitle', 'off', ...
    'Units', 'normalized', ...
    'OuterPosition', [0 .05 1 .9], ...
    'NextPlot', 'replace');

set(hdls.fig,'toolbar','figure');
% Initialize map axes, set properties.
hdls.mapAx = axesm('pcarree', ...
    'Grid', 'on', ...
    'parallellabel', 'on', ...
    'meridianlabel', 'on',...
    'maplatlimit', [-90 90], ...
    'maplonlimit', [20 380]);
set(gca, ...
    'Units', 'normalized', ...
    'Position', [0 0 .85 .8]);

zoom_handle = zoom(hdls.fig);
set(zoom_handle,'ActionPostCallback',@postZoomCallback);
set(zoom_handle,'ActionPreCallback',@preZoomCallback);
mapLatLimit = getm(hdls.mapAx, 'maplatlimit');
mapLonLimit = getm(hdls.mapAx, 'maplonlimit');
[hdls.mapXLimit(1), hdls.mapYLimit(1)] = mfwdtran(mapLatLimit(1), mapLonLimit(1)); % The axis coordinates at which the map starts
[hdls.mapXLimit(2), hdls.mapYLimit(2)] = mfwdtran(mapLatLimit(2), mapLonLimit(2));
hdls.axisXLimit = [-4 4]; %Initial axis limits, will be used to figure out zoomed in region
hdls.axisYLimit = [-2 2];

colorbar();

hdls.coast_plot = draw_coast_line(hdls.mapAx, hdls.constants.COAST_COLOR);

hdls.currentDateIndex = 1;
initializeHorizontalBar();

initializeDateListBox(dates, hdls);

display_axis_legend();

%% GUI initialization functions
    function initializeHorizontalBar()
        width = .85;
        height = .075;
        separator = .005; % The distance between horizontal panels, corresponding to horizontal bar width
        horizontalBar = uipanel(hdls.fig, ...
            'Units', 'normalized', ...
            'Position', [0 .875 width height],...
            'BorderType', 'none', ...
            'BackgroundColor', [.8 .8 .8]);

        panel_width = 0.1;

        currentSize = [panel_width 1];
        currentLeftBottom = [0 0];
        initializeDateBoxPanel(horizontalBar, currentLeftBottom, currentSize, hdls, horizontalBarFontSize);

        lastSize = currentSize;
        lastLeftBottom = currentLeftBottom;
        currentSize = [panel_width 1];
        currentLeftBottom = ...
            [(lastLeftBottom(1) + lastSize(1) + separator) 0];
        initializeEddyTypePanel(horizontalBar, currentLeftBottom, currentSize, hdls, horizontalBarFontSize);

        lastSize = currentSize;
        lastLeftBottom = currentLeftBottom;
        currentSize = [panel_width 1];
        currentLeftBottom = ...
            [(lastLeftBottom(1) + lastSize(1) + separator) 0];
        initializeTrackTypePanel(horizontalBar, currentLeftBottom, currentSize, hdls, horizontalBarFontSize);

        lastSize = currentSize;
        lastLeftBottom = currentLeftBottom;
        currentSize = [panel_width 1];
        currentLeftBottom = ...
            [(lastLeftBottom(1) + lastSize(1) + separator) 0];
        initializeLengthFilterPanel(horizontalBar, currentLeftBottom, currentSize, hdls);

        lastSize = currentSize;
        lastLeftBottom = currentLeftBottom;
        currentSize = [panel_width 1];
        currentLeftBottom = ...
            [(lastLeftBottom(1) + lastSize(1) + separator) 0];
        initializeBackgroundPanel(horizontalBar, currentLeftBottom, currentSize, hdls, horizontalBarFontSize, background_data);

        lastSize = currentSize;
        lastLeftBottom = currentLeftBottom;
        currentSize = [panel_width 1];
        currentLeftBottom = ...
            [(lastLeftBottom(1) + lastSize(1) + separator) 0];
        initializeContourPanel(horizontalBar, currentLeftBottom, currentSize, hdls, horizontalBarFontSize);

        lastSize = currentSize;
        lastLeftBottom = currentLeftBottom;
        currentSize = [panel_width 1];
        currentLeftBottom = ...
            [(lastLeftBottom(1) + lastSize(1) + separator) 0];
        initializeCheltonPanel(horizontalBar, currentLeftBottom, currentSize, hdls, horizontalBarFontSize);

        lastSize = currentSize;
        lastLeftBottom = currentLeftBottom;
        currentSize = [panel_width 1];
        currentLeftBottom = ...
            [(lastLeftBottom(1) + lastSize(1) + separator) 0];
        initializeSelectEddyPanel(horizontalBar, currentLeftBottom, currentSize, hdls, horizontalBarFontSize);
    end

    function preZoomCallback(~, ~)
    end

    function postZoomCallback(~, evd)
        hdls.axisXLimit = get(evd.Axes, 'XLim');
        hdls.axisYLimit = get(evd.Axes, 'YLim');
        if hdls.is_autoadjusting_background
            hdls.update_background();
        end
    end

    function display_axis_legend()
        % Create dumb plot objects to legend plots on the map

        cyc_past_track_plot = plotm([NaN NaN], ...
            'Color', hdls.constants.PAST_TRACK_COLOR, ...
            'Parent', hdls.mapAx);
        cyc_future_track_plot = plotm([NaN NaN], ...
            'Color', hdls.constants.CYC_FUTURE_TRACK_COLOR, ...
            'Parent', hdls.mapAx);
        ant_future_track_plot = plotm([NaN NaN], ...
            'Color', hdls.constants.ANT_FUTURE_TRACK_COLOR, ...
            'Parent', hdls.mapAx);
        chelton_cyc_future_track_plot = plotm([NaN NaN], ...
            'Color', hdls.constants.CHELTON_CYC_FUTURE_TRACK_COLOR, ...
            'Parent', hdls.mapAx);
        chelton_ant_future_track_plot = plotm([NaN NaN], ...
            'Color', hdls.constants.CHELTON_ANT_FUTURE_TRACK_COLOR, ...
            'Parent', hdls.mapAx);

        eddy_plots = zeros(1, length(hdls.constants.eddy_types));
        for i = 1:length(hdls.constants.eddy_type_names)
            eddy_plots(i) = plotm([NaN NaN], ...
                'Color', hdls.constants.eddy_colors{i}, ...
                'marker', hdls.constants.eddy_markers{i}, ...
                'LineStyle', 'None', ...
                'Parent', hdls.mapAx);
        end

        chelton_eddy_plot = plotm([NaN NaN], ...
            'Color', hdls.constants.CHELTON_CURRENT_TRACK_COLOR, ...
            'Marker', hdls.constants.CHELTON_MARKER, ...
            'LineStyle', 'None');

        legend_names{1} = 'Coast';
        legend_names{2} = 'Tracks in the past';
        legend_names{3} = 'Future cyclonic tracks';
        legend_names{4} = 'Future anticyclonic tracks';
        legend_names{5} = 'Chelton future cyclonic tracks';
        legend_names{6} = 'Chelton future anticyclonic tracks';
        for i = 1:length(eddy_plots)
            legend_names{6 + i} = [hdls.constants.eddy_type_names{i} ' features'];
        end
        legend_names{end + 1} = 'Chelton features';

        legend( [hdls.coast_plot cyc_past_track_plot cyc_future_track_plot ant_future_track_plot ...
            chelton_cyc_future_track_plot chelton_ant_future_track_plot eddy_plots chelton_eddy_plot], ...
            legend_names);
    end
end

author: jfaghm

repo: OceanEddies

save_path: github-repos/MATLAB/jfaghm-OceanEddies


text:
function [out,Xt,str,ts] = rfnn_miso_grid(t,Xt,u,flag,itaVector,alphaVector, NumInVars,NumInTerms,x0,T)

% This program is an implementation of the on line RFNN (MISO) system.
% The structure of the network is determined by the user.
% The input space is partitioned using the grid-type method.
% All parameters of the network are estimated by Gradient Descent (GD)
% through error backpropagation.

    ninp = NumInVars;
    nout = 1;
   ninps = ninp+nout+1;  % number of inputs to sfunction [ x y LE ]
   NumRules = NumInTerms^NumInVars;  % Grid-Type Input Space Partitioning.
       ns = 4*NumInVars*NumInTerms + NumRules;
     nds = 3*NumInVars*NumInTerms + NumRules;
     % Learning Rates
     ita1 = itaVector(1); ita2 = itaVector(2);
     ita3 = itaVector(3); ita4 = itaVector(4);
     % Momentum Constants.
     alpha1 = alphaVector(1); alpha2 = alphaVector(2);
     alpha3 = alphaVector(3); alpha4 = alphaVector(4);
%  ----------------------- % initial informations --------------
if abs(flag)==0

    out = [0,ns+nds,1+ns+nds,ninps,0,1,1];    % states, outputs, inputs, ?, df, #ts
    str = [];                                 % API block consistency
     ts = T;                                  % sample time
    Xt = x0;
%  ----------------------- % state derivatives -----------------
elseif abs(flag) == 2

          x = u(1:ninp);
          e = u(ninp+1:ninp+nout);
   learning = u(ninp+nout+1);

if learning == 1

   % Unroll the states:
   off=1;
   off_end=NumInVars*NumInTerms;
   mean2=reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off + NumInVars*NumInTerms-1;
   sigma2=reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off+NumInVars*NumInTerms-1;
   Theta2=reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off + NumRules-1;
   W = reshape(Xt(off:off_end),NumRules,1);

   off=off_end+1;
   off_end=off+NumInVars*NumInTerms-1;
   Out2 = reshape(Xt(off:off_end),NumInVars,NumInTerms);

   % Unroll the differential states:
   off=off_end+1;
   off_end=off + NumInVars*NumInTerms-1;
   dmean2=reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off + NumInVars*NumInTerms-1;
   dsigma2=reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off + NumInVars*NumInTerms-1;
   dTheta2=reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off + NumRules-1;
   dW = reshape(Xt(off:off_end),NumRules,1);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                                                    FEEDFORWARD OPERATION                                                      %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % LAYER 2 - INPUT TERM NODES
  Out2_pr = Out2;

 In2 = x*ones(1,NumInTerms) + Out2_pr.*Theta2;
 Out2 = exp(-((In2-mean2)./sigma2).^2);

% LAYER 3 - RULE (PRODUCT) NODES
precond = comb(Out2);
 Out3 = prod(precond,2);

%%%%%%%%%%%% END OF NETWORK FUNCTIONALITY SECTION %%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%			 					                 PARAMETER LEARNING SECTION	                       			                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BACKWARD PASS. Error Backpropagation
% LAYERS 4-3
delta3 = e*W;

% LAYER 2
Q = zeros(NumInVars,NumInTerms,NumRules);
for i=1:NumInVars
     for j=1:NumInTerms
          for k=1:NumRules
       	       if Out2(i,j)== precond(k,i) && Out2(i,j)~=0
      			  Q(i,j,k) = (Out3(k)/Out2(i,j))*delta3(k);
               end
          end
   	  end
 end

 ThetaE = sum(Q,3);

% LAYER 1 PARAMETER ADJUSTMENT BY GRADIENT DESCENT.
deltamean2  =  2*ThetaE.*Out2.*(In2-mean2)./((sigma2).^2);
deltasigma2 =  2*ThetaE.*Out2.*((In2-mean2).^2)./(sigma2.^3);
deltaTheta2 = -2*ThetaE.*Out2.*(In2-mean2).*Out2_pr./sigma2.^2;

dmean2 = ita2*deltamean2 + alpha2*dmean2;
 mean2  = mean2 + dmean2;

dsigma2 = ita3*deltasigma2 + alpha3*dsigma2;
  sigma2 = sigma2 + dsigma2;

dTheta2 = ita4*deltaTheta2 + alpha4*dTheta2;
 Theta2  = Theta2 + dTheta2;

% LAYER 4 PARAMETER ADJUSTMENT
  deltaW = e*Out3;
  dW = ita1*deltaW + alpha1*dW;
   W = W + dW;

%%%%%%%%%%%   END OF PARAMETER LEARNING PROCESS %%%%%%%%%
% State Vector Storage.
% Xt = [mean2 sigma2 Theta2 W Out2 dmean2 dsigma2 dTheta2 dW];

Xt = [ reshape(mean2,NumInVars*NumInTerms,1);
          reshape(sigma2,NumInVars*NumInTerms,1);
          reshape(Theta2,NumInVars*NumInTerms,1);
          W;
          reshape(Out2,NumInVars*NumInTerms,1);
          reshape(dmean2,NumInVars*NumInTerms,1);
          reshape(dsigma2,NumInVars*NumInTerms,1);
          reshape(dTheta2,NumInVars*NumInTerms,1);
          dW;];
end

out=Xt;

%  ----------------------- % outputs -------------------------
elseif flag == 3

  % Unpack the network's parameters first...
   off=1;
   off_end=NumInVars*NumInTerms;
   mean2=reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off + NumInVars*NumInTerms-1;
   sigma2=reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off+NumInVars*NumInTerms-1;
   Theta2 =reshape(Xt(off:off_end),NumInVars,NumInTerms);

   off=off_end+1;
   off_end=off + NumRules - 1;
   W = Xt(off:off_end);

   off=off_end+1;
   off_end=off+NumInVars*NumInTerms-1;
   Out2 = reshape(Xt(off:off_end),NumInVars,NumInTerms);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %                                                       FEEDFORWARD OPERATION                                                       %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % LAYER 1 - INPUT TERM NODES
  Out2_pr = Out2;

  x = u(1:ninp);
 In2 = x*ones(1,NumInTerms) + Out2_pr.*Theta2;
 Out2 = exp(-((In2-mean2)./sigma2).^2);

% LAYER 3 - RULE (PRODUCT) NODES
precond = comb(Out2);
 Out3 = prod(precond,2);

 % LAYER 4
  outact = W.'*Out3;

  % Block Outputs Vector Formation.
   out=[outact;Xt];

else
   out=[];
end


text:
classdef prtClusterSphericalKmeans < prtCluster
    % prtClusterSphericalKmeans Spherical K-Means Clustering with Cosine
    %     Metric
    %
    % sk = prtClusterSphericalKmeans generates a spherical K-means
    %   clustering object based on the algorithm description in [1] below.
    %   Spherical K-means operates under the assumption that the input data
    %   rows have zero-mean and unit standard deviation.  Otherwise the
    %   resulting clusters may not work very well.
    %
    % sk = prtClusterSphericalKmeans(varargin) enables the inclusion of
    %   various parameter/value pairs.
    %
    %
    %  A prtClusterSphericalKmeans object has the following properites:
    %
    %   nClusters - 3 - The number of clusters to learn
    %   nIters - 10 - the number of clustering iterations to use
    %
    %  A prtClusterSphericalKmeans object also inherits all properties
    %   and functions from the prtCluster class
    %
    % [1] Learning Feature Representations with K-means, Adam Coates and
    % Andrew Y. Ng
    %
    % Example usage:
    %    See the entry in ]blogs\torrione_2013.03.15_CoatesNg_Kmeans for
    %    example usage

    properties (SetAccess=private)
        name = 'Spherical K-Means'
        nameAbbreviation = 'SKM'
    end

    properties
        nClusters = 3;
        nIters = 10;
    end

    properties (SetAccess = protected)
        clusterCenters = [];   % The cluster centers
    end

    methods

        function self = set.nClusters(self,value)
            if ~prtUtilIsPositiveScalarInteger(value)
                error('prt:prtClusterKmeans:nClusters','value (%s) must be a positive scalar integer',mat2str(value));
            end
            self.nClusters = value;
        end

        function self = prtClusterSphericalKmeans(varargin)
            % Allow for string, value pairs
            self = prtUtilAssignStringValuePairs(self,varargin{:});
        end
    end

    methods (Access=protected, Hidden = true)

        function self = trainAction(self,dataSet)

            X = dataSet.X;
            d1 = dataSet.bootstrap(self.nClusters); %random inits
            clusters = d1.X';

            clusters = bsxfun(@rdivide,clusters,sqrt(sum(clusters.^2,1)));
            clusters(isnan(clusters)) = 0;

            inner = X*clusters;
            [val,indJ] = max(inner,[],2);

            matInd = sub2ind(size(inner),(1:length(indJ))',indJ);
            boolMat = false(size(inner));
            boolMat(matInd) = true;
            inner(~boolMat) = 0;

            for i = 1:self.nIters
                clusters = clusters + X'*inner;

                clusters = bsxfun(@rdivide,clusters,sqrt(sum(clusters.^2,1)));
                clusters(isnan(clusters)) = 0;

                inner = X*clusters;
                [val,indJ] = max(inner,[],2);
                matInd = sub2ind(size(inner),(1:length(indJ))',indJ);
                boolMat = false(size(inner));
                boolMat(matInd) = true;
                inner(~boolMat) = 0;
            end
            self.clusterCenters = clusters;
        end

        function dataSet = runAction(self,dataSet)

            X = dataSet.X;
            inner = X*self.clusterCenters;
            dataSet.X = inner;
            %             [val,indJ] = max(inner,[],2);
            %             matInd = sub2ind(size(inner),(1:length(indJ))',indJ);
            %             boolMat = false(size(inner));
            %             boolMat(matInd) = true;
            %
            %             dataSet.X = double(boolMat);

        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function whl = Find_RB_LED(file)

% This code tests an approximate median filter for separating out the red
% and blue leds mounted on top of the subjects head
% I was using the DataMax system, which was synced using a red LED. Hence,
% I am tracking the x,y position of the animal as well as when the sync
% light comes on.

%file       = 'Mouse12-120808-01.mpg';
%[~,name,~] = fileparts(file);
% Create Text file for writing out LED locations and sync trigger
%fid        = fopen(sprintf('%s.whl',name),'w+');
% Creater readerobj of file
readerobj  = VideoReader(file);
width      = readerobj.Width;
height     = readerobj.Height;
threshF    = 150;    % threshold on foreground pixel intensity
% Initial frame
Fint       = 1;

% Initialize grid for locating centroid of LED
[X,Y] = meshgrid(1:width,1:height);

% Initialize background as a grayscale image of the first frame
bg_bw     = rgb2gray(read(readerobj,Fint));

% Initialize foreground image and difference image
fg          = zeros(size(bg_bw));
fr_diff     = zeros(size(bg_bw));

% Initialize color mask
rmask  = zeros(height,width,3,'uint8');
bmask  = zeros(height,width,3,'uint8');

% Initialize fr in case first frame reading returns error
fr  = zeros(height,width,3,'uint8');

% Initialize whl matrix
whl = zeros(readerobj.NumberOfFrames,4);

for i = Fint:readerobj.NumberOfFrames
%     % Initialize red and blue LEDS as missing (=> -1)
    Red = [-1 -1];
    Blue = [-1 -1];

    %fprintf('%i',i);

    % Access frame of interest - if error (mostly for the last frames,
    % don't know why), will analyze previous frame...
    try
        fr    = read(readerobj,i);
    end

    % Spatial filter
    filt = fspecial('average',50);
    frs = imfilter(fr,filt);
    fr2 = fr-frs;

    filt = fspecial('average',10);
    fr3 = imfilter(fr2, filt);

    % convert frame to thresholdable
%     fr_bw = rgb2gray(fr);
    fr_r = double(fr2(:,:,1));
    fr_b = double(fr2(:,:,3));
    template = zthreshimage(rgb2gray(fr2),0.1);%pixels that aren't too dim after filtering... ie worth paying attention to
    fr_r_ok = fr_r;
    fr_b_ok = fr_b;
    fr_b_ok(find(~template)) = 0;
    fr_r_ok(find(~template)) = 0;

    rbdiff = fr_r-fr_b;
    brdiff = fr_b-fr_r;
    %would have to brightness threshold first
    rbratio = log(fr_r_ok./fr_b_ok);
    brratio = log(fr_b_ok./fr_r_ok);

    % bright enough in each color
    zbright_r = zthreshimage(fr_r,1);
    zbright_b = zthreshimage(fr_b,1);
    absbright_r = threshimage(fr_r,60);
    absbright_b = threshimage(fr_b,60);

    % color selective enough by diff
    zdiffsel_r = zthreshimage(rbdiff,3);% diff above z score
    zdiffsel_b = zthreshimage(brdiff,3);
    absdiffsel_r = threshimage(rbdiff,20);% diff above absolute level
    absdiffsel_b = threshimage(brdiff,20);

    % color selective enough by ratio
    zratsel_r = zthreshimage(rbratio,1);% diff above z score
    zratsel_b = zthreshimage(brratio,1);

%     ratsel_r = im2bw(rbratio,.51);% above 55st %ile ratio
%     ratsel_b = im2bw(brratio,.51);% above 51st %ile ratio

    % bright and selective enough
%     combo_r = diffsel_r.*ratsel_r.*bright_r;
%     combo_b = diffsel_b.*ratsel_r.*bright_b;
    combo_r = zbright_r.*absbright_r.*zdiffsel_r.*absdiffsel_r.*zratsel_r;
    combo_b = zbright_b.*absbright_b.*zdiffsel_b.*absdiffsel_b.*zratsel_b;

% keep only the largest object in each thresholded image
    % blue
    [l,numobjs] = bwlabel(combo_b);
    if numobjs>0
        numeach = zeros(1,numobjs);
        for a = 1:numobjs
            numeach(a) = sum(l(:)==a);
        end
        [~,bgroup] = max(numeach);
        pl = regionprops(l,'PixelList');
        Blue = round(mean(pl(bgroup).PixelList,1));
    end
    % red
    [l,numobjs] = bwlabel(combo_r);
    if numobjs>0
        numeach = zeros(1,numobjs);
        for a = 1:numobjs
            numeach(a) = sum(l(:)==a);
        end
        [~,rgroup] = max(numeach);
        pl = regionprops(l,'PixelList');
        Red = round(mean(pl(rgroup).PixelList,1));
    end

    whl(i,:) = [Red(1),Red(2),Blue(1),Blue(2)];

    % End processing time, now outputting/plotting
    if mod(i,100)==0 & i<1000
        ixr = whl(i-99:i,1);
        ixr = ixr>-1;
        ixb = whl(i-99:i,3);
        ixb = ixb>-1;
        ok = ixr.*ixb;
        figure(1),clf,
        subplot(3,1,1);imshow(uint8(fr));title('If looks wrong, Ctrl+C and change threshold \!!')

        subplot(3,3,5);imagesc(cat(3,combo_r,zeros(size(combo_r)),combo_b));
        xlim([0 size(fr,1)])
        ylim([0 size(fr,2)])
        axis tight

        subplot(3,3,8)
            plot(whl(i-99:i,1),whl(i-99:i,2),'r.')
            hold on
            plot(whl(i-99:i,3),whl(i-99:i,4),'b.')
            mx = mean([whl(i-100+find(ok),1) whl(i-100+find(ok),3)],2);
            my = mean([whl(i-100+find(ok),2) whl(i-100+find(ok),4)],2);
            plot(mx,my,'Marker','.','color',[0.5 0.5 0.5],'LineStyle','none')
            xlim([0 size(fr,1)])
            ylim([0 size(fr,2)])
            axis ij
        fprintf('Detection of Red LED failed %i times (%i times for the blue LED) in the last 100 frames\n',sum(~ixr),sum(~ixb));
    end
    if i == 1;
        h = waitbar(i/readerobj.NumberOfFrames);
    else
        waitbar(i/readerobj.NumberOfFrames,h);
    end

end
close(h)

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function idx = getDealiasingIndexes(~, N, nVars)
%GETGDEALIASINGINDEXES  % Returns indexes for dealiasing procedure (2/3-rule).
%
% See also SPINOP2.

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

toOne = floor(N/2) + 1 - ceil(N/6):floor(N/2) + ceil(N/6);
idx = false(N, N);
idx(toOne, toOne) = 1;
idx = repmat(idx, nVars, 1);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function gp = gp_unpak(gp, w, param)
%GP_UNPAK  Set GP parameters from vector to structure
%
%  Description
%    GP = GP_UNPAK(GP, W, PARAM) takes an Gaussian Process data
%    structure GP and a parameter vector W, and returns a Gaussian
%    Process structure identical to the input, except that the
%    parameters has been set to the ones in W. PARAM defines which
%    parameters are present in the W vector. If PARAM is not given
%    the function unpacks all parameters.
%
%    Each of the following strings in PARAM defines one group of
%    parameters to unpack:
%      covariance  - unpack parameters of covariance function
%      likelihood  - unpack parameters of likelihood
%      inducing    - unpack inducing inputs (in sparse approximations):
%                    W = gp.X_u(:)
%
%    By combining the strings one can unpack more than one group of
%    parameters. For example:
%      covariance+inducing   - unpack covariance function parameters
%                              and inducing inputs
%      covariance+likelihood - unpack covariance function parameters
%                              and likelihood parameters
%
%    Inside each group (such as covariance functions) the
%    parameters to be unpacked is defined by the existence of a
%    prior structure. For example, if GP has two covariance
%    functions but only the first one has prior for its parameters
%    then only the parameters of the first one are unpacked. Thus,
%    also inducing inputs require prior if they are to be
%    optimized.
%
%    GP_PAK and GP_UNPAK functions are used, e.g., when GP
%    parameters are optimized with GP_OPTIM or sampled with GP_MC.
%    See GP_SET and option 'infer_params'.
%
%  See also
%    GP_PAK, GP_SET
%
% Copyright (c) 2007-2010 Jarno Vanhatalo

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if isempty(w)
  return
end

if nargin < 3
  param = gp.infer_params;
end

if size(w,1) > 1
  error(' The vector to be packed has to be row vector! \n')
end

% Unpack SVI variational parameters
if isfield(gp, 'latent_method') && isequal(gp.latent_method, 'SVI') ...
    && isfield(gp, 't1')
  w=w(:);
  gp.t1=w(1:length(gp.t1));
  w(1:length(gp.t1))=[];
  t2=reshape(w(1:numel(gp.t2)), length(gp.t1), length(gp.t1));
  [V,D]=eig(t2);
  gp.t2=t2;
  w(1:numel(gp.t2))=[];
  dtmp=-0.5./diag(D);
%   if any(dtmp<eps)
%     % Force positive definite covariance matrix
%     fprintf(['Variational covariance not positive definite. ' ...
%       'Forcing positive definitiness.\n']);
%     dtmp(dtmp<eps)=eps;
%     %gp.t2=-0.5.*(V*diag(1./dtmp))/V;
%   end
  dtmp=diag(dtmp);
  gp.S=(V*dtmp)/V;
  gp.m=gp.S*gp.t1;
  w=w';
end

% Unpack the parameters of covariance functions
if ~isempty(strfind(param, 'covariance'))
  ncf = length(gp.cf);

  for i=1:ncf
    gpcf = gp.cf{i};
    [gpcf, w] = gpcf.fh.unpak(gpcf, w);
    gp.cf{i} = gpcf;
  end

end

% Unpack the parameters of likelihood function
if ~isempty(strfind(param, 'likelihood'))
  [gp.lik, w] = gp.lik.fh.unpak(gp.lik, w);
  if isfield(gp, 'latent_method') && isequal(gp.latent_method, 'SVI') ...
      && ~isequal(gp.lik.type, 'Gaussian') && ~isempty(gp.lik.p.sigma2)
    gp.lik.sigma2=exp(w(1));
    w(1)=[];
  end
  % Unpack the parameters of the second likelihood function (monotonicity)
  if isfield(gp, 'lik_mono')
    [gp.lik_mono, w] = gp.lik_mono.fh.unpak(gp.lik_mono, w);
    if isfield(gp, 'latent_method') && isequal(gp.latent_method, 'SVI') ...
         && ~isempty(gp.lik_mono.p.sigma2)
      gp.lik_mono.sigma2=exp(w(1));
      w(1)=[];
    end
  end
end

% Unpack the inducing inputs
if ~isempty(strfind(param, 'inducing'))
  if isfield(gp,'p') && isfield(gp.p, 'X_u') && ~isempty(gp.p.X_u)
    if ~iscell(gp.p.X_u)
      % One prior for all inducing inputs
      lu = length(gp.X_u(:));
      gp.X_u = reshape(w(1:lu), size(gp.X_u'))';
      if lu < length(w)
        w = w(lu+1:end);
      end
      [gp.p.X_u, w]=gp.p.X_u.fh.unpak(gp.p.X_u, w);
    else
      % Own prior for each inducing input
      d=numel(gp.X_u(1,:));
      for i=1:size(gp.X_u)
        gp.X_u(i,:) = reshape(w(1:d), size(gp.X_u(1,:)));
        w(1:d)=[];
        [gp.p.X_u{i}, w] = gp.p.X_u{i}.fh.unpak(gp.p.X_u{i}, w);
      end
    end
  end
end

% Unpack the covariance function parameters
if ~isempty(strfind(param, 'mean'))
  nmf = length(gp.meanf);

  for i=1:nmf
    gpmf = gp.meanf{i};
    [gpmf, w] = gpmf.fh.unpak(gpmf, w);
    gp.meanf{i} = gpmf;
  end

end

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
% demo SIFT using LabelMe toolbox

img = imread('demo1.jpg');
img = imresize(img, .5, 'bilinear');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SIFT parameters:
SIFTparam.grid_spacing = 1; % distance between grid centers
SIFTparam.patch_size = 16; % size of patch from which to compute SIFT descriptor (it has to be a factor of 4)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% CONSTANTS (you can not change this)
w = SIFTparam.patch_size/2; % boundary

% COMPUTE SIFT: the output is a matrix [nrows x ncols x 128]
tic;
SIFT = LMdenseSift(img, '', SIFTparam);
toc

figure
subplot(121)
imshow(img(w:end-w+1,w:end-w+1,:))
title('cropped image')
subplot(122)
showColorSIFT(SIFT)
title('SIFT color coded')

author: CSAILVision

repo: LabelMeToolbox


text:
function datatype = GetDatatype(handles)
global maingui

datatype=0;

if nargin==0
    return
end
if isempty(handles)
    return
end

plotRawVal  = get(handles.radiobuttonPlotRaw, 'value');
plotODVal   = get(handles.radiobuttonPlotOD, 'value');
plotConcVal = get(handles.radiobuttonPlotConc, 'value');
plotHRFVal  = get(handles.checkboxPlotHRF, 'value');

c = maingui.buttonVals;

if plotRawVal && ~plotHRFVal
    datatype      = c.RAW;
elseif plotRawVal && plotHRFVal
    datatype      = c.RAW_HRF;
elseif plotODVal && ~plotHRFVal
    datatype      = c.OD;
elseif plotODVal && plotHRFVal
    datatype      = c.OD_HRF;
elseif plotConcVal && ~plotHRFVal
    datatype      = c.CONC;
elseif plotConcVal && plotHRFVal
    datatype      = c.CONC_HRF;
end
maingui.datatype = datatype;

author: BUNPC


text:
function bbci_trigger_print(value)
%BBCI_TRIGGER_PRINT - Proxi for trigger function, just printing

v= clock;
time_str= sprintf('%02d:%02d:%06.3f', v(4), v(5), v(6));
fprintf('[TRIGGER PRINT] trig #%03d at %s\n', value, time_str);

author: bbci


text:
function [x,y,z]=ell2xyz(lat,lon,h,a,e2)
% ELL2XYZ  Converts ellipsoidal coordinates to cartesian.
%   Vectorized.
% Version: 2011-02-19
% Useage:  [x,y,z]=ell2xyz(lat,lon,h,a,e2)
%          [x,y,z]=ell2xyz(lat,lon,h)
% Input:   lat - vector of ellipsoidal latitudes (radians)
%          lon - vector of ellipsoidal E longitudes (radians)
%          h   - vector of ellipsoidal heights (m)
%          a   - ref. ellipsoid major semi-axis (m); default GRS80
%          e2  - ref. ellipsoid eccentricity squared; default GRS80
% Output:  x \
%          y  > vectors of cartesian coordinates in CT system (m)
%          z /

% Copyright (c) 2011, Michael R. Craymer
% All rights reserved.
% Email: mike@craymer.com

if nargin ~= 3 & nargin ~= 5
  warning('Incorrect number of input arguments');
  return
end
if nargin == 3
  [a,b,e2]=refell('grs80');
end

v=a./sqrt(1-e2*sin(lat).*sin(lat));
x=(v+h).*cos(lat).*cos(lon);
y=(v+h).*cos(lat).*sin(lon);
z=(v.*(1-e2)+h).*sin(lat);


text:
% Fig. 5.12   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami

clf
np=[1 1];
dp=[1 4 0 0];
rlocus(np,dp)
axis([-6 2 -3 3])
grid on
title('Root locus for Figure 5.12')
z=0:.1:.9;
 wn= .5:.5:6;
 sgrid(z, wn)
 hold off


text:
%--- help for load_varlist ---
%
%  INTERNAL FUNCTION: retrieves the list of all the variables
%  entering the object or equation
%
%

author: jmaih


text:
function [ ] = save_callback_denoising_general( theta, info, state, eI, varargin)
% Copyright (c) 2014-present University of Illinois at Urbana-Champaign
% All rights reserved.
%
% Developed by:     Po-Sen Huang, Paris Smaragdis
%                   Department of Electrical and Computer Engineering
%                   Department of Computer Science
%
% save model and run evaluation while minfunc is running

if mod(info.iteration, 100) == 0
    if isfield(eI, 'iterStart')
      info.iteration = info.iteration+eI.iterStart;
    end
    if ~strcmp(state,'init')
        if isfield(eI,'ioffset'),
            saveName = sprintf('%smodel_off%d_%d.mat',eI.saveDir,eI.ioffset, info.iteration);
        else
            saveName = sprintf('%smodel_%d.mat',eI.saveDir,info.iteration);
        end
        save(saveName, 'theta', 'eI',  'info');
    end
    % run evaluation
%       if mod(info.iteration, 20) == 0 && info.iteration>=0
%          test_mir1k_general_kl_bss3(eI.modelname, theta, eI, 'iter', info.iteration);
%       end
end;

end

author: posenhuang

repo: deeplearningsourceseparation


text:
function omega = angle(f,ori,varargin)
% angle fibre to orientation or fibre to fibre
%
% Syntax
%
%   omega = angle(f,ori) % angle orientation to fibre
%   omega = angle(f1,f2) % angle fibre to fibre
%
% Input
%  f, f1, f2 - @fibre
%  ori - @orientation
%
% Output
%  omega - double
%
% See also
% orientation/angle

if isa(ori,'orientation')
  omega = angle(ori .\ f.r,f.h,varargin{:});

else
  omega = max(angle(f,orientation(ori),varargin{:}));

  % in the non symmetric case we have also
  %omega = min(angle(f.h,ori.h) + angle(f.r,ori.r), angle(f.h,-ori.h) + angle(f.r,-ori.r));

end

end

author: mtex-toolbox

repo: mtex

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


text:
function out = proj_psd(X)
%PROJ_PSD computes the orthogonal projection onto the cone of positive semidefinite matrics
%                                     (psd cone)   {X: X psd }
%
%  Usage:
%  out = PROJ_PSD(X)
%  ===========================================
%  INPUT:
%  X - matrix to be projected
%  ===========================================
%  Assumptions:
%  X is symmetric
%  ===========================================
%  Output:
%  out - projection matrix

% This file is part of the FOM package - a collection of first order methods for solving convex optimization problems
% Copyright (C) 2017 Amir and Nili Beck
%
% 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/>.

if (nargin < 1)
    error ('usage: proj_psd(X)') ;
end

eps = 1e-10 ;   % defalut value for eps : 1e-10
if ((size(X,1) ~= size(X,2)) || (norm( X - X') > eps))
    error('usage: proj_psd(X) - X should be a symmetric matrix') ;
end

X = 0.5 * (X + X');

[V,D] = eig(X) ;

out = V * max(D,0) *V' ;

author: hiroyuki-kasai

repo: SGDLibrary


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

%%
clear; close all; clc;

%% Syntax
% |h=axis_inc(s);|

%% Description
% This function widens the axis limits for the axis with the handle |hs|
% using the scale factor |s| (e.g. if s=2 the axis width is doubled). The
% scaling parameter |s| can be a single scalar or a vector such that each
% axis direction is scaled differently.
%
%% Examples
%

%%
% Plot settings
fontSize=25;
lineWidth=3;

%% Scaling limits on a 2D plot

scaleFactor=2;

cFigure;
h1=subplot(1,2,1); hold on;
title('Original axis limits');
plot([-1 1 1 -1],[-1 -1 1 1],'b-','LineWidth',lineWidth);
axis tight; axis equal; view(2); box on; grid on;
set(h1,'FontSize',fontSize);

h2=subplot(1,2,2); hold on;
title(['Axis limits expanded using scale factor of ',sprintf('%.2f',scaleFactor)]);
plot([-1 1 1 -1],[-1 -1 1 1],'b-','LineWidth',lineWidth);
axis tight; axis equal; view(2); box on; grid on;
axis_inc(scaleFactor,h2);
set(h2,'FontSize',fontSize);
drawnow;

%% Scaling limits on a 3D plot

[F,V]=graphicsModels(8);

scaleFactor=[2 4 2];

cFigure;
h1=subplot(1,2,1); hold on;
title('Original axis limits');
gpatch(F,V,'bw','none');
axisGeom(h1,fontSize);
set(h1,'FontSize',fontSize);
camlight headlight;

h2=subplot(1,2,2); hold on;
title(['Axis limits expanded using scale factor of ',sprintf('%.2f, %.2f, %.2f',scaleFactor)]);
gpatch(F,V,'bw','none');
axisGeom(h1,fontSize);
axis_inc(scaleFactor,h2);
set(h2,'FontSize',fontSize);
camlight headlight;
drawnow;

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


text:
function SR

clc;
clear;
close all;
%% the example fuction
datasize=200;
var_noise=0.1;
x=linspace(-10,10,datasize)';
    f=sin(x)./x;
noise=var_noise*randn(datasize,1);
y=f+noise;
x_raw=x;
plot(x_raw,f,'r-','LineWidth',2),hold on;
%%
Inset=30;
dataindex=randsample(datasize,Inset);
xI=x(dataindex);
yI=y(dataindex);
%%
%%test
testnum=500;
x_test=linspace(-15,15,testnum)';

K_mn=Kernel_func(xI,x,0.7,1.2);
K_mm=Kernel_func(xI,xI,0.7,1.2);
 rank(K_mm)
K_mat=K_mn*1/var_noise*eye(datasize)*K_mn'+K_mm;
K_sm=Kernel_func(x_test,xI,0.7,1.2);
L=chol(K_mat,'lower');
b=K_mn*1/var_noise*eye(datasize)*y;
f_mean=K_sm*(L'\(L\b));

for i=1:testnum
K_sm=Kernel_func(x_test(i),xI,0.7,1.2);
v=L\(K_sm');
  f_var(i)=v'*v;
end
f_var=f_var';

plot(x_test,f_mean,'b-','LineWidth',2),hold on;

x_area=[x_test;flipdim(x_test,1)];
ybound_std=[(f_mean+1.96*sqrt(f_var));flipdim(f_mean-1.96*sqrt(f_var),1)];
fill(x_area,ybound_std,[187 255 0]/255,'EdgeColor','none');
plot(x_raw,f,'r-','LineWidth',2),hold on;
plot(xI,yI,'b+'),hold on;
plot(x_test,f_mean,'b-','LineWidth',2),hold on;
%%
legend1=legend('target function','training points','preditive distribution','95% confidence interval');
 set(legend1,'Box','off','Color','none',...
     'Location','NorthEast');
xlabel('(f)SR with 30 points ')
xlim([-15,15]);
ylim([-2,2]);
set(gca,'XTick',[-15:3:15])
matlab2tikz( 'SubReg30.tex' )

author: epfl-lasa


text:
function obj=MOxUnitTestNode(name)
% Initialize MOxUnitTestNode instance
%
% obj=MOxUnitTestNode()
%
% Input:
%   name            name of this TestNode
%
% Output:
%   obj             MOxUnitTestNode instance
%
% Notes:
%   - This class is intended as an 'abstract' superclass, and should not
%     be initialized in typical use cases. Instead, use its subclasses,
%     MoxUnitTestSuite and MoxUnitTestCase.
%
% NNO 2015

    s=struct();
    s.name=name;
    obj=class(s,'MOxUnitTestNode');

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function mean = truncated_normal_ab_mean ( mu, sigma, a, b )

%*****************************************************************************80
%
%% TRUNCATED_NORMAL_AB_MEAN returns the mean of the truncated Normal PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 August 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real MU, SIGMA, the mean and standard deviation of the
%    parent Normal distribution.
%
%    Input, real A, B, the lower and upper truncation limits.
%
%    Output, real MEAN, the mean of the PDF.
%
  alpha = ( a - mu ) / sigma;
  beta = ( b - mu ) / sigma;

  alpha_cdf = normal_01_cdf ( alpha );
  beta_cdf = normal_01_cdf ( beta );

  alpha_pdf = normal_01_pdf ( alpha );
  beta_pdf = normal_01_pdf ( beta );

  mean = mu + sigma * ( alpha_pdf - beta_pdf ) / ( beta_cdf - alpha_cdf );

  return
end


text:
function r = trirnd(varargin)
%TRIRND Pseudorandom numbers drawn from the triangular distribution
%
%   R = IOSR.STATISTICS.TRIRND(N) returns an N-by-N matrix containing
%   pseudorandom values drawn from the triangular distribution constrained
%   to (-1,1) and mode = 0.
%
%   IOSR.STATISTICS.TRIRND(M,N) or IOSR.STATISTICS.TRIRND([M,N]) returns an
%   M-by-N matrix.
%
%   IOSR.STATISTICS.TRIRND(M,N,P,...) or
%   IOSR.STATISTICS.TRIRND([M,N,P,...]) returns an M-by-N-by-P-by-...
%   array.
%
%   IOSR.STATISTICS.TRIRND returns a scalar.
%
%   TRIRND(SIZE(A)) returns an array the same size as A.
%
%   Note: The size inputs M, N, P, ... should be nonnegative integers.
%   Negative integers are treated as 0.
%
%   The sequence of numbers produced by TRIRND is determined by the
%   settings of the uniform random number generator that underlies RAND,
%   RANDI, and RANDN. Control that shared random number generator using
%   RNG.
%
%   See also IOSR.STATISTICS.LAPRND, RAND, RANDN, RANDI, RNG.

%   Based on code (Matlab FE File ID: #13705) written by Elvis Chen, 2007.

%   Copyright 2016 University of Surrey.

    % Generate traingular noise
    u1 = rand(varargin{:})-0.5;
    u2 = rand(varargin{:})-0.5;
    r = u1+u2;

end

author: IoSR-Surrey

repo: MatlabToolbox