text:
%% DEMO_febio_0035_blob_shear_contact_hex8
% Below is a demonstration for:
%
% * Building geometry for a hemi-spherical blob with hexahedral elements
% which is being sheared by a rigid wall.
% This demo consists off:
% * Defining the boundary conditions
% * Coding the febio structure
% * Running the model
% * Importing and visualizing the displacement results

%% Keywords
%
% * febio_spec version 3.0
% * febio, FEBio
% * indentation
% * contact, sliding, friction
% * rigid body constraints
% * hexahedral elements, hex8, hex20
% * quadrilaterial elements, quad4, quad8
% * shell elements
% * sphere
% * static, solid
% * hyperelastic, Ogden
% * hyperelastic, Fung orthotropic
% * displacement logfile

%%

clear; close all; clc;

%% Plot settings
fontSize=15;
faceAlpha1=0.8;
faceAlpha2=0.3;
markerSize=40;
markerSize2=15;
lineWidth=3;

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

% Hemi-sphere parameters
hemiSphereRadius=6;
numElementsMantel=4;
smoothEdge=1;
solidElementType='hex8';%'hex20';
membraneThickness=0.1;

% Ground plate parameters
plateRadius=1.3*hemiSphereRadius;

% Probe parameters
probeWidth=3*hemiSphereRadius;
filletProbe=hemiSphereRadius/2; %Fillet radius

% Define probe displacement
probeDisplacement=hemiSphereRadius*2;
probeOverlapFactor=0.3;
probeLength=hemiSphereRadius*2;

% Material parameter set
l=1e-3;
lt=1000*l;
m=1e-3;
mt=1000*m;

mu1=mt;
mu2=mt;
mu3=m;
lambda11=lt;
lambda22=lt;
lambda12=lt;
lambda33=l;
lambda23=l;
lambda13=l;
[E1,E2,E3,G12,G23,G31,v12,v23,v31]=lameInvertHookeOrthotropic(mu1,mu2,mu3,lambda11,lambda22,lambda33,lambda12,lambda23,lambda13);

materialPropertiesFung.E1=E1;
materialPropertiesFung.E2=E2;
materialPropertiesFung.E3=E3;
materialPropertiesFung.G12=G12;
materialPropertiesFung.G23=G23;
materialPropertiesFung.G31=G31;
materialPropertiesFung.v12=v12;
materialPropertiesFung.v23=v23;
materialPropertiesFung.v31=v31;
materialPropertiesFung.c=10;
materialPropertiesFung.k=1000*mean([G12 G23 G31]);

%Ogden parameters
materialPropertiesOgden.c1=1e-3; %Shear-modulus-like parameter
materialPropertiesOgden.m1=2; %Material parameter setting degree of non-linearity
materialPropertiesOgden.k=1000*materialPropertiesOgden.c1; %Bulk modulus

% FEA control settings
numTimeSteps=25;
max_refs=25; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=15; %Optimum number of iterations
max_retries=25; %Maximum number of retires
symmetric_stiffness=0;
min_residual=1e-20;
step_size=1/numTimeSteps;
dtmin=(1/numTimeSteps)/100; %Minimum time step size
dtmax=(1/numTimeSteps); %Maximum time step size

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

%% Creating model geometry and mesh
%

%Control settings
optionStruct.sphereRadius=hemiSphereRadius;
optionStruct.coreRadius=optionStruct.sphereRadius/2;
optionStruct.numElementsMantel=numElementsMantel;
optionStruct.numElementsCore=optionStruct.numElementsMantel*2;
optionStruct.outputStructType=2;
optionStruct.makeHollow=0;
optionStruct.cParSmooth.n=25;

% %Creating sphere
[meshStruct]=hexMeshHemiSphere(optionStruct);

% Access model element and patch data
Fb_blob=meshStruct.facesBoundary;
Cb_blob=meshStruct.boundaryMarker;
V_blob=meshStruct.nodes;
E_blob=meshStruct.elements;

F_blob=element2patch(E_blob);

%%

pointSpacingBlob=max(patchEdgeLengths(Fb_blob,V_blob));

%Smoothen edges
if smoothEdge==1
    %Get rigid region
    ind=1:1:size(V_blob,1); %Indices for all nodes
    indRigid1=find(ismember(ind,Fb_blob(Cb_blob==2,:)) & ~ismember(ind,Fb_blob(Cb_blob==1,:))); %Indices for new bottom surface nodes
    indRigid2=find(ismember(ind,Fb_blob(Cb_blob==1,:)) & ~ismember(ind,Fb_blob(Cb_blob==2,:))); %Indices for new bottom surface nodes
    indRigid=[indRigid1(:); indRigid2(:);];

    %Smoothing
    cPar.Method='HC';
    cPar.n=250;
    cPar.RigidConstraints=indRigid;
    [Vb_blob]=patchSmooth(F_blob,V_blob,[],cPar);
    indSmooth=unique(F_blob(:));
    V_blob(indSmooth,:)=Vb_blob(indSmooth,:);
    %Fix color data with new bottom surface
    Cb_blob=ones(size(Cb_blob));
    Cb_blob(all(ismember(Fb_blob,indRigid1),2))=2;

    meshStruct.nodes=V_blob;
end

if strcmp(solidElementType,'hex20')
    [E_blob,V_blob,~,~,Fb_blob]=hex8_hex20(E_blob,V_blob,{},Fb_blob);
    [Fb_blob_plot]=element2patch(Fb_blob,[],'quad8');

    meshStruct.elements=E_blob;
    meshStruct.nodes=V_blob;
    meshStruct.Fb=Fb_blob_plot;
    shellElementType='quad8';
else
    Fb_blob_plot=Fb_blob;
    shellElementType='quad4';
end

%%
% Visualize blob mesh

hFig=cFigure;
subplot(1,2,1); hold on;
hp=gpatch(Fb_blob_plot,V_blob,Cb_blob,'k',1);
hp.Marker='.';
hp.MarkerSize=markerSize2;

% patchNormPlot(Fb_blob,V_blob);
% plotV(V_blob(indRigid,:),'g.','MarkerSize',25);
axisGeom(gca,fontSize);
colormap(gjet); icolorbar;
camlight headlight;

hs=subplot(1,2,2); hold on;
title('Cut view of solid mesh','FontSize',fontSize);
optionStruct.hFig=[hFig hs];
gpatch(Fb_blob_plot,V_blob,'kw','none',0.25);
meshView(meshStruct,optionStruct);
axisGeom(gca,fontSize);
drawnow;

%% Creating rigid body shear surface

pointSpacingProbe=pointSpacingBlob/2;

%Sketching side profile
d=hemiSphereRadius*cos(asin(1-probeOverlapFactor));
x=[-probeLength-hemiSphereRadius -d -d];
y=[0 0 0];
z=[hemiSphereRadius*(1-probeOverlapFactor) hemiSphereRadius*(1-probeOverlapFactor) hemiSphereRadius*1.5];
V_probe_curve_sketch=[x(:) y(:) z(:)];

%Fillet sketch
np=100; %Number of points used to construct each fillet edge
[V_probe_curve]=filletCurve(V_probe_curve_sketch,filletProbe,np,0);
numPointsProbeCurve=ceil(max(pathLength(V_probe_curve))/pointSpacingProbe);
[V_probe_curve] = evenlySampleCurve(V_probe_curve,numPointsProbeCurve,'pchip',0);

% Extruding curve
% controlParametersExtrude.pointSpacing=pointSpacingProbe;
controlParametersExtrude.depth=hemiSphereRadius*2.5;
controlParametersExtrude.numSteps=ceil(controlParametersExtrude.depth/pointSpacingProbe);
controlParametersExtrude.numSteps=controlParametersExtrude.numSteps+iseven(controlParametersExtrude.numSteps); %Force uneven
controlParametersExtrude.patchType='quad';
controlParametersExtrude.dir=0;
controlParametersExtrude.n=[0 1 0];
controlParametersExtrude.closeLoopOpt=0;

[F_probe,V_probe]=polyExtrude(V_probe_curve,controlParametersExtrude);
F_probe=fliplr(F_probe); %Invert face orientation so normals point to blob

if strcmp(solidElementType,'hex20')
    [F_probe,V_probe]=quad4_quad8(F_probe,V_probe);
    [F_probe_plot]=element2patch(F_probe,[],'quad8');
else
    F_probe_plot=F_probe;
end

center_of_mass_probe=mean(V_probe,1);

%%
% Visualizing probe mesh

cFigure; hold on;
title('The probe surface mesh','fontSize',fontSize);
gpatch(Fb_blob_plot,V_blob,'kw','none',0.5);
hl(1)=plotV(V_probe_curve_sketch,'k.-.','lineWidth',3,'MarkerSize',25);
hl(2)=plotV(V_probe_curve,'r-','lineWidth',3,'MarkerSize',25);
hl(3)=gpatch(F_probe_plot,V_probe,'gw','k',1);
% hl(3).Marker='.';
% hl(3).MarkerSize=markerSize2;
legend(hl,{'Sketched probe curve','Rounded probe curve','Probe surface mesh'}); clear hl;
axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Join model node sets

V=[V_blob; V_probe];
F_probe=F_probe+size(V_blob,1);
F_probe_plot=F_probe_plot+size(V_blob,1);

%%
% Visualizing model

cFigure; hold on;
gtitle('Model components',fontSize);
hl(1)=gpatch(Fb_blob_plot,V,'rw','k',0.8);
hl(2)=gpatch(F_probe_plot,V,'gw','k',0.8);
legend(hl,{'Blob','Probe'}); clear hl;
axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Get local material axes

[N3,Vn]=patchNormal(Fb_blob_plot,V);
[N1,N2]=vectorOrthogonalPair(N3); %Get orthogonal vector pair

%%
% Visualizing axes

cFigure; hold on;
gtitle('Local material axes',fontSize);
hl(1)=gpatch(Fb_blob_plot,V,'w','k',1);
hl(2)=quiverVec(Vn,N1,pointSpacingBlob,'y');
hl(3)=quiverVec(Vn,N2,pointSpacingBlob,'g');
hl(4)=quiverVec(Vn,N3,pointSpacingBlob,'b');
legend(hl,{'Blob','1st axis','2nd axis','3rd axis'}); clear hl;
axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Visualize contact surfaces

cFigure;
title('Probe blob contact pair','fontsize',fontSize);
hl(1)=gpatch(F_probe_plot,V,'rw','k',1);
patchNormPlot(F_probe_plot,V);
hl(2)=gpatch(Fb_blob_plot,V,'gw','k',1);
patchNormPlot(Fb_blob_plot,V);
legend(hl,{'Secondary','Primary'}); clear hl;
axisGeom(gca,fontSize);
camlight headlight;

drawnow;

%% Setup boundary conditions
%

bcSupportList=unique(Fb_blob(Cb_blob==2,:));

%%
% Visualize boundary condition sets

cFigure; hold on;
title('Boundary conditions ','fontsize',fontSize);
gpatch(F_probe_plot,V,'kw','none',0.2);
gpatch(Fb_blob_plot,V,'kw','none',0.2);
hl(1)=plotV(V(bcSupportList,:),'k.','MarkerSize',markerSize);

legend(hl,{'Supported nodes'}); clear hl;
axisGeom(gca,fontSize);
camlight headlight;

drawnow;

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

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

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

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

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

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

materialName2='Material2';
febio_spec.Material.material{2}.ATTR.name=materialName2;
febio_spec.Material.material{2}.ATTR.type='Fung orthotropic';
febio_spec.Material.material{2}.ATTR.id=2;
febio_spec.Material.material{2}.E1=materialPropertiesFung.E1;
febio_spec.Material.material{2}.E2=materialPropertiesFung.E2;
febio_spec.Material.material{2}.E3=materialPropertiesFung.E3;
febio_spec.Material.material{2}.G12=materialPropertiesFung.G12;
febio_spec.Material.material{2}.G23=materialPropertiesFung.G23;
febio_spec.Material.material{2}.G31=materialPropertiesFung.G31;
febio_spec.Material.material{2}.v12=materialPropertiesFung.v12;
febio_spec.Material.material{2}.v23=materialPropertiesFung.v23;
febio_spec.Material.material{2}.v31=materialPropertiesFung.v31;
febio_spec.Material.material{2}.c=materialPropertiesFung.c;
febio_spec.Material.material{2}.k=materialPropertiesFung.k;

materialName3='Material3';
febio_spec.Material.material{3}.ATTR.name=materialName3;
febio_spec.Material.material{3}.ATTR.type='rigid body';
febio_spec.Material.material{3}.ATTR.id=3;
febio_spec.Material.material{3}.density=1;
febio_spec.Material.material{3}.center_of_mass=center_of_mass_probe;

%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=solidElementType; %Element type
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E_blob,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E_blob; %The element matrix

partName2='Part2';
febio_spec.Mesh.Elements{2}.ATTR.name=partName2; %Name of this part
febio_spec.Mesh.Elements{2}.ATTR.type=shellElementType; %Element type
febio_spec.Mesh.Elements{2}.elem.ATTR.id=size(E_blob,1)+(1:1:size(Fb_blob,1))'; %Element id's
febio_spec.Mesh.Elements{2}.elem.VAL=Fb_blob;

partName3='Part3';
febio_spec.Mesh.Elements{3}.ATTR.name=partName3; %Name of this part
febio_spec.Mesh.Elements{3}.ATTR.type=shellElementType; %Element type
febio_spec.Mesh.Elements{3}.elem.ATTR.id=size(E_blob,1)+size(Fb_blob,1)+(1:1:size(F_probe,1))'; %Element id's
febio_spec.Mesh.Elements{3}.elem.VAL=F_probe;

% -> 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{1}.ATTR.name=partName2;
febio_spec.MeshDomains.ShellDomain{1}.ATTR.mat=materialName2;

febio_spec.MeshDomains.ShellDomain{2}.ATTR.name=partName3;
febio_spec.MeshDomains.ShellDomain{2}.ATTR.mat=materialName3;

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

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

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

%MeshData section
% -> ElementData
febio_spec.MeshData.ElementData{1}.ATTR.var='shell thickness';
febio_spec.MeshData.ElementData{1}.ATTR.elem_set=partName2;
febio_spec.MeshData.ElementData{1}.elem.ATTR.lid=(1:1:size(Fb_blob,1))';
febio_spec.MeshData.ElementData{1}.elem.VAL=membraneThickness.*ones(size(Fb_blob));

%MeshData section
% -> ElementData
febio_spec.MeshData.ElementData{2}.ATTR.elem_set=partName2;
febio_spec.MeshData.ElementData{2}.ATTR.var='mat_axis';

for q=1:1:size(N1,1)
    febio_spec.MeshData.ElementData{2}.elem{q}.ATTR.lid=q;
    febio_spec.MeshData.ElementData{2}.elem{q}.a=N1(q,:);
    febio_spec.MeshData.ElementData{2}.elem{q}.d=N2(q,:);
end

%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.Rigid.rigid_constraint{1}.ATTR.name='RigidFix_1';
febio_spec.Rigid.rigid_constraint{1}.ATTR.type='fix';
febio_spec.Rigid.rigid_constraint{1}.rb=3;
febio_spec.Rigid.rigid_constraint{1}.dofs='Ry,Rz,Ru,Rv,Rw';

febio_spec.Rigid.rigid_constraint{2}.ATTR.name='RigidPrescribe';
febio_spec.Rigid.rigid_constraint{2}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{2}.rb=3;
febio_spec.Rigid.rigid_constraint{2}.dof='Rx';
febio_spec.Rigid.rigid_constraint{2}.value.ATTR.lc=1;
febio_spec.Rigid.rigid_constraint{2}.value.VAL=probeDisplacement;
febio_spec.Rigid.rigid_constraint{2}.relative=0;

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

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

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

%% 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='external';%'internal';

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

%% Import FEBio results

if 1%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)]);

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

    DN_magnitude=sqrt(sum(N_disp_mat(:,:,end).^2,2)); %Current displacement magnitude

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; hold on;
    gtitle([febioFebFileNamePart,': Press play to animate']);
    hp1=gpatch(Fb_blob_plot,V_DEF(:,:,end),DN_magnitude,'k',1); %Add graphics object to animate
    hp1.FaceColor='interp';
    hp2=gpatch(F_probe_plot,V_DEF(:,:,end),'kw','none',0.5); %Add graphics object to animate
%     gpatch(Fb_all,V,0.5*ones(1,3),'none',0.25); %A static graphics object

    axisGeom(gca,fontSize);
    colormap(gjet(250)); colorbar;
    caxis([0 max(DN_magnitude(Fb_blob_plot(:)))]); caxis manual;
    axis(axisLim(V_DEF)); %Set axis limits statically
    camlight headlight;
    view(15,8);
    drawnow;
    % Set up animation features
    animStruct.Time=timeVec; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments
        DN_magnitude=sqrt(sum(N_disp_mat(:,:,qt).^2,2)); %Current displacement magnitude

        %Set entries in animation structure
        animStruct.Handles{qt}=[hp1 hp1 hp2]; %Handles of objects to animate
        animStruct.Props{qt}={'Vertices','CData','Vertices'}; %Properties of objects to animate
        animStruct.Set{qt}={V_DEF(:,:,qt),DN_magnitude,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 H = deriveBasis(A_equ)
[L,U,P] = lu(A_equ);
[L,U,P] = lu(A_equ');
r = colspaces(L');
AA = L';
H1 = AA(:,r);
H2 = AA(:,setdiff(1:size(AA,2),r));
H = P'*[-H1\H2;speye(size(H2,2))];
function  [indx]=colspaces(A)
indx = [];
for i = 1:size(A,2)
    s = max(find(A(:,i)));
    indx = [indx s];
end
indx = unique(indx);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function out = eval_chance( imdb, labelspath )
% RETRIEVAL_EVAL_CHANCE

% Read the files
labels = utls.readfile(labelspath);
poolSignatures = strsplit(labels{1}, ',');
% Compare the headers
numQueries = numel(labels) - 1;
imRetAps = zeros(1, numQueries); patchRetAps = zeros(1, numQueries);
numAllFeats = sum(cellfun(@(sign) imdb.getNumPatches(sign), poolSignatures));
for lni = 2:numel(labels)
  query_cluster_descs = strsplit(labels{lni}, ',');
  query_cluster_ims = remove_patchnum(query_cluster_descs);
  query_cluster_nfeats = sum(cellfun(@(sign) imdb.getNumPatches(sign), query_cluster_ims));

  imRetAps(lni) = (query_cluster_nfeats - 1) ./ (numAllFeats - 1);
  patchRetAps(lni) = (numel(query_cluster_descs) - 1) ./ (numAllFeats - 1);
end
out = struct('image_retr_ap', imRetAps, 'patch_retr_ap', patchRetAps);
end

function signs = remove_patchnum(signs)
% Get a sequence.imagename part of a signature
for si = 1:numel(signs)
  parts = strsplit(signs{si}, '.');
  assert(numel(parts) == 3, 'Invalid patch signature.');
  signs{si} = strjoin(parts(1:2), '.');
end
end

author: hpatches

repo: hpatches-benchmark


text:
function i4vec_index_test ( )

%*****************************************************************************80
%
%% I4VEC_INDEX_TEST tests I4VEC_INDEX;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_INDEX_TEST\n' );
  fprintf ( 1, '  For an integer vector:\n' );
  fprintf ( 1, '  I4VEC_INDEX:              first index of given value;\n' );

  seed = 123456789;
  b = -n;
  c = n;

  [ a, seed ] = i4vec_uniform_ab ( n, b, c, seed );

  i4vec_print ( n, a, '  Input vector:' );

  fprintf ( 1, '\n' );

  aval = a(floor(n/2));
  fprintf ( 1, '\n' );
  j = i4vec_index ( n, a, aval );
  fprintf ( 1, '  Index of first occurrence of %d is %d\n', aval, j );

  aval = aval + 1;
  j = i4vec_index ( n, a, aval );
  fprintf ( 1, '  Index of first occurrence of %d is %d\n', aval, j );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
      function [A, phi, f] = amplitude_spectrum(td)
         %amplitude_spectrum   Simple method to compute amplitude
         %   spectrum for a trace. Uses the MATLAB fft function.
         %   [A, phi, f] = amplitude_spectrum(td)
         %
         %   Inputs:
         %       td - a single TraceData
         %
         %   Outputs:
         %       A - the amplitude coefficients
         %       phi - the phase coefficients
         %       f - the frequency values corresponding to elements of A and phi
         %
         %   Example:
         %       [A phi, f] = amplitude_spectrum(td)
         %       plot(f,A);
         %       xlabel('Frequency (Hz)')
         %       ylabel('Amplitude');
         %
         %   See also fft
         %   Glenn Thompson, November 21, 2014

         %TODO: Rename this to amplitudespectrum or move it elsewhere

         N = length(td.data);
         NFFT = 2 ^ nextpow2(N); % Next power of 2 from length of y
         Y = fft(td.data, NFFT); % X will have same length as signal, and will be complex with a magnitude and phase
         A = 2 * abs(Y(1:NFFT/2+1)) / N;
         phi = angle(Y);
         f = td.samplerate / 2 * linspace(0,1,NFFT/2+1);
      end

author: geoscience-community-codes

repo: GISMO


text:
function mr = mrLoadAnalyze(pth);
% Load an ANALYZE 7.5-format file as a mrVista mr struct.
%
%  mr = mrLoadAnalyze(pth);
%
%
% ras, 01/06.
mr = mrCreateEmpty;

[mr.data, mr.voxelSize, mr.hdr] = mrReadAnalyzeData(pth);
mr.dims = size(mr.data);

%% check for file with extra fields
[p f ext] = fileparts(pth);
extraFieldsFile = fullfile(p, [f '.mat']);
if exist(extraFieldsFile, 'file')
	tmp = load(extraFieldsFile);
	fields = setdiff(fieldnames(tmp), 'M');
	for i = 1:length(fields)
		field = fields{i};
		if isfield(mr, field)
			% combine the existing and additional info
			mr.(field) = mergeStructures(mr.(field), tmp.(field));
		else
			mr.(field) = tmp.(field);
		end
	end
end

%% add spaces
mr.spaces = mrStandardSpaces(mr);

mr.spaces(end+1).name = 'R|A|S';
mr.spaces(end).xform = mr.hdr.mat;
mr.spaces(end).dirLabels =  {'R <--> L' 'A <--> P' 'S <--> I'};
mr.spaces(end).sliceLabels = {'Sagittal' 'Coronal' 'Axial'};

mr.spaces(end+1) = mr.spaces(end);
mr.spaces(end).name = 'I|P|R';
mr.spaces(end).xform = ras2ipr(mr.spaces(end-1).xform, mr.dims);
mr.spaces(end).dirLabels =  {'S <--> I' 'A <--> P' 'L <--> R'};
mr.spaces(end).sliceLabels = {'Axial' 'Coronal' 'Sagittal'};

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
classdef figurePrinter < handle

    properties (Access = protected)
        x0 = 1920;
        y0 = 500;
        width = 600;
        height = 600;
        %labelFontSize = 30;
        labelFontSize = 200;
        %fontSize = 20;
        fontSize = 30;
        outFormat = '-dpng';
        figID
    end

    methods (Access = public)

        function print(obj,figName)
            obj.printFigure(figName)
        end

    end

    methods (Access = protected)

        function setFontSizes(obj)
            obj.setCurrentAxisFontSize();
            obj.setXlabelFontSize();
            obj.setYlabelFontSize();
        end

        function setPositionAndSize(obj)
            set(obj.figID,'units','points','position',[obj.x0,obj.y0,obj.width,obj.height])
        end

        function printFigure(obj,figName)
            print(obj.figID,figName,obj.outFormat);
        end

    end

    methods (Access = private)

        function setCurrentAxisFontSize(obj)
            set(gca,'fontsize',obj.fontSize)
        end

        function setXlabelFontSize(obj)
            ca = get(obj.figID,'CurrentAxes');
            xl = get(ca,'xlabel');
            set(xl,'FontSize',obj.labelFontSize)
        end

        function setYlabelFontSize(obj)
            ca = get(obj.figID,'CurrentAxes');
            yl = get(ca,'ylabel');
            set(yl,'FontSize',obj.labelFontSize)
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function ax = setAxis(fh, axis_id)
% Enabled a certain axis in the figure, update MATLAB curfigure status
% useful to reset the current axis in case of new plotting
%
% INPUT
%   fh       figurehandle
%   axis_id  number of the axis (1 is the first inserted)
%
% SYNTAX
%   setAxis(fh, axis_id)

%--- * --. --- --. .--. ... * ---------------------------------------------
%               ___ ___ ___
%     __ _ ___ / __| _ | __|
%    / _` / _ \ (_ |  _|__ \
%    \__, \___/\___|_| |___/
%    |___/                    v 1.0RC1
%
%--------------------------------------------------------------------------
%  Copyright (C) 2021 Geomatics Research & Development srl (GReD)
%  Written by: Giulio Tagliaferro
%  Contributors:     ...
%  A list of all the historical goGPS contributors is in CREDITS.nfo
%--------------------------------------------------------------------------
%
%   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/>.
%
%--------------------------------------------------------------------------
% 01100111 01101111 01000111 01010000 01010011
%--------------------------------------------------------------------------

    if not(isa(fh, 'matlab.ui.Figure'))
        axis_id = fh;
        fh = gcf;
    end
    set(0, 'CurrentFigure', fh);
    try
        ax = fh.Children(max(1, end + 1 - axis_id));
        subplot(ax);
    catch
        ax = axes();
    end
end

author: goGPS-Project


text:
% Dasslc test problem
% requires files: dydt1.m

t0  = 0.0;        % initial value for independent variable
tf  = 1.0;        % final value for independent variable
y0  = [3 6 5]';   % initial state variables

[t,y]=dasslc('dydt1',[t0 tf],y0);
plot(t,y);

author: Sable

repo: mcbench-benchmarks


text:
function varargout = subsref(grains,s)
% implements grains(1:3)
%
% Syntax
%   grains(1:10)            % the 10 first grains
%   grains('Fe')            % only Fe grains
%   grains('id',5)          % give the grain with id 5
%   grains(5)               % give the 5th grain in the list
%   grains( ~grains('fe') ) % all grains but Fe
%   grains(cond)
%
% Input
%  grains - @grain2d
%  cond   - logical array with same size as grains
%

% some special cases to speed things up
if strcmp(s(1).type,'()') && ...
    length(s)>1 && strcmp(s(2).type,'.') && strcmp(s(2).subs,'meanOrientation')

  if strcmp(s(1).type,'{}')
    ind = grains.id2ind(s(1).subs{1});
  else
    ind = subsind(grains,s(1).subs);
  end

  phId = unique(grains.phaseId(ind));

  if length(phId) > 1
    error('MTEX:MultiplePhases',['\n' ...
      '----------------------------------------------------------------\n'...
      'Your variable contains the phases: ' ...
      grains.mineralList{phId(1)} ', ' grains.mineralList{phId(2)} '\n\n' ...
      'However, you are executing a command that is only permitted for a single phase!\n\n' ...
      'Please read the chapter ' doclink('EBSDSelect','"select EBSD data"')  ...
      ' for how to restrict grains to a single phase.\n' ...
      '----------------------------------------------------------------\n']);
  end
  if isempty(ind)
    ori = orientation;
  else
    ori = orientation(grains.prop.meanRotation(ind),grains.CSList{phId});
  end

  if numel(s)>2
    [varargout{1:nargout}] = builtin('subsref',ori,s(3:end));
  else
    varargout{1} = ori;
  end
  return
end

if strcmp(s(1).type,'()') || strcmp(s(1).type,'{}')

  if strcmp(s(1).type,'{}')
    ind = grains.id2ind(s(1).subs{1});
  else
    ind = subsind(grains,s(1).subs);
  end

  grains = subSet(grains,ind);

  % is there something more to do?
  if numel(s)>1
    s = s(2:end);
  else
    varargout{1} = grains;
    return
  end
end

% maybe reference to a dynamic property
if isProperty(grains,s(1).subs)

  [varargout{1:nargout}] = subsref@dynProp(grains,s);

else

  [varargout{1:nargout}] = builtin('subsref',grains,s);

end

end

author: mtex-toolbox

repo: mtex

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


text:
function test_issue1585

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_read_header ft_read_data ft_read_event

%%

cd(dccnpath('/home/common/matlab/fieldtrip/data/test/issue1585'))

% cfg = [];
% cfg.dataset      = 'eeglab_data.mat';
% cfg.datatype     = 'eeg';

% cfg.headerformat = 'eeglab_set';
% cfg.dataformat   = 'eeglab_set';
% cfg.eventformat  = 'eeglab_set';
%
% cfg.method    = 'convert';
% cfg.bidsroot  = './bids';
% cfg.sub       = '001';
% cfg.task      = 'whatever';
% data2bids(cfg);

% after creating the BIDS dataset, I used a text editor to remove a row from the channels.tsv

%%

try
  % this should give an error
  hdr = ft_read_header('bids/sub-001/eeg/sub-001_task-whatever_eeg.eeg', 'readbids', 'yes');
  dat = ft_read_data('bids/sub-001/eeg/sub-001_task-whatever_eeg.eeg', 'readbids', 'yes');
  evt = ft_read_event('bids/sub-001/eeg/sub-001_task-whatever_eeg.eeg', 'readbids', 'yes');
  failed = false;
catch
  failed = true;
end
assert(failed, 'it should have failed');

% this should give a warning
hdr = ft_read_header('bids/sub-001/eeg/sub-001_task-whatever_eeg.eeg', 'readbids', 'ifmakessense');
dat = ft_read_data('bids/sub-001/eeg/sub-001_task-whatever_eeg.eeg', 'readbids', 'ifmakessense');
evt = ft_read_event('bids/sub-001/eeg/sub-001_task-whatever_eeg.eeg', 'readbids', 'ifmakessense');

% this should work without any complaints
hdr = ft_read_header('bids/sub-001/eeg/sub-001_task-whatever_eeg.eeg', 'readbids', 'no');
dat = ft_read_data('bids/sub-001/eeg/sub-001_task-whatever_eeg.eeg', 'readbids', 'no');
evt = ft_read_event('bids/sub-001/eeg/sub-001_task-whatever_eeg.eeg', 'readbids', 'no');

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%  Figure 10.21      Feedback Control of Dynamic Systems, 6e
%                        Franklin, Powell, Emami
%
%   fig10_21.m is a script to generate Fig 10.21, the
%   frequency response of the LQR symmetric rootlocus compensator for the
%   satellite position control, non-colocated case WITH ESTIMATOR

m=[1, 0.1]; k=[0, 0.091] ; d=[0, 0.0036]; k1=[0, 0.4];
[f,g,h,j] = twomass(m,k,d);
s=[f, g;h, 0];r=[0*g;1]; n=s\r;nx=n(1:4);nu=n(5);
[f1,g,h,j] = twomass(m,k1,d);
a=[f, 0*f;
-h'*h, -f'];
b=[g;0*g];
c=[0*h, g'];
d=[0];

hold off; clf

P=eig(a-b*c*0.1621);
pc=P(real(P<0)==1);
K=place(f,g,pc);
nbar=nu+K*nx;
% eig(f-g*K)
P=eig(a-b*c*3.056e7);
pe=P(real(P<0)==1);
L=place(f',h',pe)';
ac=f-g*K-L*h ;bc=L;cc=K;dc=0;
ac=f-g*K-L*h ;bc=L;cc=K;dc=0;
% [numc denc]=ss2tf(ac,bc,cc,dc)

w=logspace(-1,1);
w(26) = 1; w(25) = .94;w(25)=.9;
[mag ph]=bode(ac,bc,cc,dc,1,w);
subplot(211); loglog(w,mag); grid;
xlabel('\omega (rad/sec)');
ylabel('Magnitude, |D_4(j\omega)|');
title('Fig. 10.21: Bode plot of the optimal compensator D_4(s)');
subplot(212); semilogx(w,ph);grid;
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');

%Bode grid
bodegrid


text:
%  Figure 3.24      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%script to generate Fig. 3.24

clf;
zeta=.5;
k=1/zeta;
den=[1 1 1];

a=10;
num=[k/a 1];
t=0:.1:10;
y1=step(num,den,t);

a=4;
num=[k/a 1];
t=0:.1:10;
y2=step(num,den,t);

a=2;
num=[k/a 1];
t=0:.1:10;
y3=step(num,den,t);

a=1;
num=[k/a 1];
t=0:.1:10;
y4=step(num,den,t);

plot(t,y1,'-',t,y2,'-',t,y3,'-',t,y4,'-'),grid
title('Fig. 3.24 Step response with \xi = 0.5')
xlabel('\omega_n t')
ylabel('Step response of H(s)')
text(.1,.9,'\alpha=')
text(.6,.9,'1')
text(1.1,.9,'2')
text(1.5,.85,'4')
text(1.8,.8,'100')


text:
function demmet1(plot_wait)
%DEMMET1 Demonstrate Markov Chain Monte Carlo sampling on a Gaussian.
%
%	Description
%	The problem consists of generating data from a Gaussian in two
%	dimensions using a Markov Chain Monte Carlo algorithm. The points are
%	plotted one after another to show the path taken by the chain.
%
%	DEMMET1(PLOTWAIT) allows the user to set the time (in a whole number
%	of seconds) between the plotting of points.  This is passed to PAUSE
%
%	See also
%	DEMHMC1, METROP, GMM, DEMPOT
%

%	Copyright (c) Ian T Nabney (1996-2001)

if nargin == 0 | plot_wait < 0
  plot_wait = 0; % No wait if not specified or incorrect
end
dim = 2;            	% Data dimension
ncentres = 1;		% Number of centres in mixture model

seed = 42;              % Seed for random weight initialization.
randn('state', seed);
rand('state', seed);

clc
disp('This demonstration illustrates the use of the Markov chain Monte Carlo')
disp('algorithm to sample from a Gaussian distribution.')
disp('The mean is at [0 0].')
disp(' ')
disp('First we set up the parameters of the mixture model we are sampling')
disp('from.')
disp(' ')
disp('Press any key to continue.')
pause

% Set up mixture model to sample from
mix = gmm(dim, ncentres, 'spherical');
mix.centres(1, :) = [0 0];
x = [0 4];  % Start vector

% Set up vector of options for hybrid Monte Carlo.

nsamples = 150;		% Number of retained samples.

options = foptions;     % Default options vector.
options(1) = 0;		% Switch off diagnostics.
options(14) = nsamples;	% Number of Monte Carlo samples returned.
options(18) = 0.1;

clc
disp('Next we take 150 samples from the distribution.')
disp('Sampling starts at the point [0 4].')
disp('The new state is accepted if the threshold value is greater than')
disp('a random number between 0 and 1.')
disp(' ')
disp('Press any key to continue.')
pause

[samples, energies] = metrop('dempot', x, options, '', mix);

clc
disp('The plot shows the samples generated by the MCMC function in order')
disp('as an animation to show the path taken by the Markov chain.')
disp('The different colours are used to show that the first few samples')
disp('should be discarded as they lie too far from the mean.')
disp(' ')
disp('Press any key to continue.')
pause
probs = exp(-energies);
fh1 = figure;
g1end = floor(nsamples/4);

for n = 1:nsamples

  if n < g1end
    Marker = 'k.';
    p1 = plot(samples(n,1), samples(n,2), Marker, ...
      'EraseMode', 'none', 'MarkerSize', 12);
    if n == 1
      axis([-3 5 -2 5])
    end
  else
    Marker = 'r.';
    p2 = plot(samples(n,1), samples(n,2), Marker, ...
      'EraseMode', 'none', 'MarkerSize', 12);
  end
  hold on
  drawnow;  % Force drawing immediately
  pause(plot_wait);
end
lstrings = char(['Samples 1-' int2str(g1end)], ...
  ['Samples ' int2str(g1end+1) '-' int2str(nsamples)]);
legend([p1 p2], lstrings, 1);

disp(' ')
disp('Press any key to exit.')
pause
close(fh1);
clear all;

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
% Script: view_coulomb.m
% Display Coulomb stress change map together with seismicity
% Allow to analyze seismicity rate changes using the Ztools
%
% J. Woessner, jochen.woessner@sed.ethz.ch, 05.07.2004

report_this_filefun(mfilename('fullpath'));

% Default value for map view
bMap =1;

% Load CFS file
[sFilename, sPathname] = uigetfile('*.mat', 'Pick CFS change MAT-file');
sHelp = [sPathname sFilename];
sFile = [sFilename(1:length(sFilename)-4)];
load(sHelp)
rCfs = eval(sFile);

% Load rate change file
[sFilename1, sPathname1] = uigetfile('*.mat', 'Pick Rate Change MAT-file');
sHelp1 = [sPathname1 sFilename1];
sFile1 = [sFilename1(1:length(sFilename1)-4)];
load(sHelp1)

% Set up the Seismicity Map window Enviroment
% Find out of figure already exists
[existFlag,figNumber]=figure_exists('Coulomb-map',1);
newcfsmapWindowFlag=~existFlag;

if newcfsmapWindowFlag
    oldfig_button = 0;
end

if oldfig_button == 0
    cfsmap = figure_w_normalized_uicontrolunits( ...
        'Name','Coulomb-map',...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'NextPlot','add', ...
        'backingstore','on',...
        'Visible','off', ...
        'Position',[ (fipo(3:4) - [600 400]) ZmapGlobal.Data.map_len]);
    % make menu bar
    matdraw

    % Display
    add_menu_divider();
    add_symbol_menu('eq_plot');

    % Menu: Ztools
    op1 = uimenu('Label','Ztools');
    uimenu(op1,'Label','Rate change, p-,c-,k-value map in aftershock sequence (MLE) ',...
        'Callback','sel= ''in'';,rcvalgrid_a2');
    uimenu(op1,'Label','Load existing  Rate change, p-,c-,k-value map (MLE)',...
        'Callback','sel= ''lo'';rcvalgrid_a2');

    % Menu: ZAnalyze
    op2 = uimenu('Label',' ZAnalyze ');
    uimenu(op2,'Label','Refresh ', 'Callback','view_coulombmap')
    uimenu(op2,'Label','Select EQ in Circle - Constant R',...
         'Callback','h1 = gca;met = ''ra''; ho=false;plot_circbootfit_a2;watchoff(cfsmap)')
    uimenu(op2,'Label','Select EQ with const. number',...
         'Callback','h1 = gca;ho2=true;ho=true;plot_constnrbootfit_a2;watchoff(cfsmap)')

    % Set default colormap
    colormap(jet)
end  % This is the end of the figure setup.

% Load the data from the Coulomb data file
% Define colormap
% mColormap = gui_Colormap_Rastafari(256);
colormap(jet);

% Get the gridding
vY = linspace(rCfs.fMinLat,rCfs.fMaxLat,rCfs.nNy);
vY = fliplr(vY);
vX = linspace(rCfs.fMinLon,rCfs.fMaxLon,rCfs.nNx);

% Plot Coulomb stress map with seismicity
figure_w_normalized_uicontrolunits(cfsmap)
% Fix color scale for imagesc
vClims = [-2 2];
hCoulomb = imagesc(vX,vY,rCfs.mCfs,vClims);
shading interp;
set(gca,'Ydir','normal');
% Colorbar
hColor = colorbar;
chl = get(hColor,'Ylabel');
set(chl,'String','\Delta CFS [bar]','FontS',10,'Rot',270);

% Labeling
ylabel('Latitude [deg]');
xlabel('Longitude [deg]');

% Add rate change map
hold on;
normlap2=ones(length(ll),1)*nan;
% Relative rate change
normlap2(ll)= mRcGrid(:,15);
mRelchange = reshape(normlap2,length(yvect),length(xvect));
hRate = pcolor(xvect,yvect,mRelchange)
set(hRate, 'AlphaData', 0.6, 'AlphaDataMapping', 'none');
shading(gca,'flat')
% Plot seismicity
ploeq=plot(a.Longitude,a.Latitude,'Markersize',3,'Marker','o','Linestyle','none','Color',[0 0 0],'Tag','eq_plot');
hold off;

author: CelsoReyes


text:
    function [ ] = showNoisyStatistics(noisyStats)
    %% Plot statistics for a dataset collection based on single noisy structure

    %% Extract the fields
    collectionTitle = noisyStats.collectionTitle;
    statisticsTitles = noisyStats.statisticsTitles;
    s = noisyStats.statisticsIndex;
    statistics = noisyStats.statistics;
    divideColor = [0.85, 0.85, 0.85];

    %% Plot median versus SDR overall channel deviation
    minDev = min(statistics(:, s.medDev));
    maxDev = max(statistics(:, s.medDev));
    minrSD = min(statistics(:, s.rSDDev));
    maxrSD = max(statistics(:, s.rSDDev));
    baseTitle = 'Median versus SDR overall channel deviation';
    figure ('Name', baseTitle, 'Color', [1, 1, 1]);
    hold on
    plot(statistics(:, s.medDev), statistics(:, s.rSDDev), 'ok')
    xlabel(statisticsTitles{s.medDev});
    ylabel(statisticsTitles{s.rSDDev});
    title(collectionTitle)
    medDev = median(statistics(:, s.medDev));
    rSDDev = median(statistics(:, s.rSDDev));
    plot(medDev, rSDDev, '+r', 'MarkerSize', 14, 'LineWidth', 3);
    set(gca, 'XLim', [minDev, maxDev], 'XLimMode', 'manual', ...
        'YLim', [minrSD, maxrSD], 'YLimMode', 'manual');
    title({collectionTitle; baseTitle; ...
        ['[Median of dataset medians: ' num2str(medDev) '(med dev) ' ...
        num2str(rSDDev) '(rSD dev)]']}, 'interpreter', 'none');
    legend('Dataset', 'Collection median', 'Location', 'NorthWest')
    box on
    hold off

    %% Plot median versus SDR window channel deviation
    minWinDev = min(statistics(:, s.medWinDev));
    maxWinDev = max(statistics(:, s.medWinDev));
    minWinrSD = min(statistics(:, s.rSDWinDev));
    maxWinrSD = max(statistics(:, s.rSDWinDev));
    baseTitle = 'Median versus SDR window channel deviation';
    figure ('Name', baseTitle, 'Color', [1, 1, 1]);
    hold on
    plot(statistics(:, s.medWinDev), statistics(:, s.rSDWinDev), 'ok')
    xlabel(statisticsTitles{s.medWinDev});
    ylabel(statisticsTitles{s.rSDWinDev});
    title(collectionTitle)
    medWinDev = median(statistics(:, s.medWinDev));
    rSDWinDev = median(statistics(:, s.rSDWinDev));
    plot(medWinDev, rSDWinDev, '+r', 'MarkerSize', 14, 'LineWidth', 3);
    line([minWinDev; maxWinDev], [minWinrSD; maxWinrSD], 'LineWidth', 3, ...
        'Color', divideColor);
    set(gca, 'XLim', [minWinDev, maxWinDev], 'XLimMode', 'manual', ...
        'YLim', [minWinrSD, maxWinrSD], 'YLimMode', 'manual');
    title({collectionTitle; baseTitle; ...
        ['[Median of dataset medians: ' num2str(medWinDev) '(med dev) ' ...
        num2str(rSDWinDev) '(rSD dev)]']}, 'interpreter', 'none');
    legend('Dataset', 'Collection median', 'Location', 'NorthWest')
    box on
    hold off

%% Plot median maximum correlation versus median window channel deviation
    minWinDev = min(statistics(:, s.medWinDev));
    maxWinDev = max(statistics(:, s.medWinDev));
    minCorr = min(statistics(:, s.medCor));
    baseTitle = 'Median maximum correlation versus median window channel deviation';
    figure ('Name', baseTitle, 'Color', [1, 1, 1]);
    hold on
    plot(statistics(:, s.medCor), statistics(:, s.medWinDev), 'ok')

    xlabel(statisticsTitles{s.medCor});
    ylabel(statisticsTitles{s.medWinDev});
    medWinDev = median(statistics(:, s.medWinDev));
    medCorr = median(statistics(:, s.medCor));
    plot(medCorr, medWinDev, '+r', 'MarkerSize', 14, 'LineWidth', 3);
    set(gca, 'XLim', [minCorr, 1], 'XLimMode', 'manual', ...
        'YLim', [minWinDev, maxWinDev], 'YLimMode', 'manual');
    title({collectionTitle; baseTitle; ...
        ['[Median of dataset medians: ' num2str(medCorr) '(corr) ' ...
        num2str(medWinDev) '(med dev)]']}, 'interpreter', 'none');
    legend('Dataset', 'Collection median', 'Location', 'NorthWest')
    box on
    hold off

    %% Plot median maximum correlation versus window channel deviation ratio
    devRatio = statistics(:, s.rSDWinDev)./statistics(:, s.medWinDev);
    minWinRatio = min(devRatio);
    maxWinRatio = max(devRatio);
    yLabelString = 'rSDR/median window channel deviation';
    minCorr = min(statistics(:, s.medCor));
    baseTitle = 'Median maximum correlation versus window channel deviation ratio';
    figure ('Name', baseTitle, 'Color', [1, 1, 1]);
    hold on
    plot(statistics(:, s.medCor), devRatio, 'ok')
    xlabel(statisticsTitles{s.medCor});
    ylabel(yLabelString);
    medRatio = median(devRatio);
    medCorr = median(statistics(:, s.medCor));
    plot(medCorr, medRatio, '+r', 'MarkerSize', 14, 'LineWidth', 3);
    set(gca, 'XLim', [minCorr, 1], 'XLimMode', 'manual', ...
        'YLim', [minWinRatio, maxWinRatio], 'YLimMode', 'manual');
    title({collectionTitle; baseTitle; ...
        ['[Median of dataset medians: ' num2str(medCorr) '(corr) ' ...
        num2str(medRatio) '(rSD/med dev)]']}, 'interpreter', 'none');
    legend('Dataset', 'Collection median', 'Location', 'NorthWest')
    box on
    hold off

%% Plot average maximum correlation versus median window channel deviation
    minWinDev = min(statistics(:, s.medWinDev));
    maxWinDev = max(statistics(:, s.medWinDev));
    minCorr = min(statistics(:, s.medCor));
    baseTitle = 'Average maximum correlation versus median window channel deviation';
    figure ('Name', baseTitle, 'Color', [1, 1, 1]);
    hold on
    plot(statistics(:, s.aveCor), statistics(:, s.medWinDev), 'ok')

    xlabel(statisticsTitles{s.aveCor});
    ylabel(statisticsTitles{s.medWinDev});
    medWinDev = median(statistics(:, s.medWinDev));
    aveCorr = median(statistics(:, s.aveCor));
    plot(aveCorr, medWinDev, '+r', 'MarkerSize', 14, 'LineWidth', 3);
    set(gca, 'XLim', [minCorr, 1], 'XLimMode', 'manual', ...
        'YLim', [minWinDev, maxWinDev], 'YLimMode', 'manual');
    title({collectionTitle; baseTitle; ...
        ['[Median of dataset medians: ' num2str(aveCorr) '(corr) ' ...
        num2str(medWinDev) '(med dev)]']}, 'interpreter', 'none');
    legend('Dataset', 'Collection median', 'Location', 'NorthWest')
    box on
    hold off

    %% Plot average maximum correlation versus window channel deviation ratio
    devRatio = statistics(:, s.rSDWinDev)./statistics(:, s.medWinDev);
    minWinRatio = min(devRatio);
    maxWinRatio = max(devRatio);
    yLabelString = 'rSDR/median window channel deviation';
    minCorr = min(statistics(:, s.aveCor));
    baseTitle = 'Average maximum correlation versus window channel deviation ratio';
    figure ('Name', baseTitle, 'Color', [1, 1, 1]);
    hold on
    plot(statistics(:, s.aveCor), devRatio, 'ok')
    xlabel(statisticsTitles{s.aveCor});
    ylabel(yLabelString);
    medRatio = median(devRatio);
    aveCorr = median(statistics(:, s.aveCor));
    plot(aveCorr, medRatio, '+r', 'MarkerSize', 14, 'LineWidth', 3);
    set(gca, 'XLim', [minCorr, 1], 'XLimMode', 'manual', ...
        'YLim', [minWinRatio, maxWinRatio], 'YLimMode', 'manual');
    title({collectionTitle; baseTitle; ...
        ['[Median of dataset medians: ' num2str(aveCorr) '(corr) ' ...
        num2str(medRatio) '(rSD/med dev)]']}, 'interpreter', 'none');
    legend('Dataset', 'Collection median', 'Location', 'NorthWest')
    box on
    hold off

    %% Plot median versus SDR overall HF noise
    minHF = min(statistics(:, s.medHF));
    maxHF = max(statistics(:, s.medHF));
    minrSD = min(statistics(:, s.rSDHF));
    maxrSD = max(statistics(:, s.rSDHF));
    baseTitle = 'Median versus SDR overall channel HF noise';
    figure ('Name', baseTitle, 'Color', [1, 1, 1]);
    hold on
    plot(statistics(:, s.medHF), statistics(:, s.rSDHF), 'ok')
    xlabel(statisticsTitles{s.medHF});
    ylabel(statisticsTitles{s.rSDHF});
    title(collectionTitle)
    medHF = median(statistics(:, s.medHF));
    rSDHF = median(statistics(:, s.rSDHF));
    plot(medHF, rSDHF, '+r', 'MarkerSize', 14, 'LineWidth', 3);
    set(gca, 'XLim', [minHF, maxHF], 'XLimMode', 'manual', ...
        'YLim', [minrSD, maxrSD], 'YLimMode', 'manual');
    title({collectionTitle; baseTitle; ...
        ['[Median of dataset medians: ' num2str(medHF) '(med HF) ' ...
        num2str(rSDDev) '(rSD HF)]']}, 'interpreter', 'none');
    legend('Dataset', 'Collection median', 'Location', 'NorthWest')
    box on
    hold off

    %% Plot median versus SDR window channel HF noise
    minWinHF = min(statistics(:, s.medWinHF));
    maxWinHF = max(statistics(:, s.medWinHF));
    minWinrSD = min(statistics(:, s.rSDWinHF));
    maxWinrSD = max(statistics(:, s.rSDWinHF));
    baseTitle = 'Median versus SDR window channel HF noise';
    figure ('Name', baseTitle, 'Color', [1, 1, 1]);
    hold on
    plot(statistics(:, s.medWinHF), statistics(:, s.rSDWinHF), 'ok')
    xlabel(statisticsTitles{s.medWinHF});
    ylabel(statisticsTitles{s.rSDWinHF});
    title(collectionTitle)
    medWinHF = median(statistics(:, s.medWinHF));
    rSDWinHF = median(statistics(:, s.rSDWinHF));
    plot(medWinHF, rSDWinHF, '+r', 'MarkerSize', 14, 'LineWidth', 3);
    set(gca, 'XLim', [minWinHF, maxWinHF], 'XLimMode', 'manual', ...
        'YLim', [minWinrSD, maxWinrSD], 'YLimMode', 'manual');
    title({collectionTitle; baseTitle; ...
        ['[Median of dataset medians: ' num2str(medWinHF) '(med HF) ' ...
        num2str(rSDWinHF) '(rSD HF)]']}, 'interpreter', 'none');
    legend('Dataset', 'Collection median', 'Location', 'NorthWest')
    box on
    hold off

%% Plot median maximum correlation versus median window channel HF noise
    minWinHF = min(statistics(:, s.medWinHF));
    maxWinHF = max(statistics(:, s.medWinHF));
    minCorr = min(statistics(:, s.medCor));
    baseTitle = 'Median maximum correlation versus median window channel HF noise';
    figure ('Name', baseTitle, 'Color', [1, 1, 1]);
    hold on
    plot(statistics(:, s.medCor), statistics(:, s.medWinHF), 'ok')

    xlabel(statisticsTitles{s.medCor});
    ylabel(statisticsTitles{s.medWinHF});
    medWinHF = median(statistics(:, s.medWinHF));
    medCorr = median(statistics(:, s.medCor));
    plot(medCorr, medWinHF, '+r', 'MarkerSize', 14, 'LineWidth', 3);
    set(gca, 'XLim', [minCorr, 1], 'XLimMode', 'manual', ...
        'YLim', [minWinHF, maxWinHF], 'YLimMode', 'manual');
    title({collectionTitle; baseTitle; ...
        ['[Median of dataset medians: ' num2str(medCorr) '(corr) ' ...
        num2str(medWinHF) '(med HF)]']}, 'interpreter', 'none');
    legend('Dataset', 'Collection median', 'Location', 'NorthWest')
    box on
    hold off

    %% Plot median maximum correlation versus window channel HF ratio
    HFRatio = statistics(:, s.rSDWinHF)./statistics(:, s.medWinHF);
    minWinRatio = min(HFRatio);
    maxWinRatio = max(HFRatio);
    yLabelString = 'rSDR/median window channel HF noise';
    minCorr = min(statistics(:, s.medCor));
    baseTitle = 'Median maximum correlation versus window channel HF noise ratio';
    figure ('Name', baseTitle, 'Color', [1, 1, 1]);
    hold on
    plot(statistics(:, s.medCor), HFRatio, 'ok')
    xlabel(statisticsTitles{s.medCor});
    ylabel(yLabelString);
    medRatio = median(HFRatio);
    medCorr = median(statistics(:, s.medCor));
    plot(medCorr, medRatio, '+r', 'MarkerSize', 14, 'LineWidth', 3);
    set(gca, 'XLim', [minCorr, 1], 'XLimMode', 'manual', ...
        'YLim', [minWinRatio, maxWinRatio], 'YLimMode', 'manual');
    title({collectionTitle; baseTitle; ...
        ['[Median of dataset medians: ' num2str(medCorr) '(corr) ' ...
        num2str(medRatio) '(rSD/med HF)]']}, 'interpreter', 'none');
    legend('Dataset', 'Collection median', 'Location', 'NorthWest')
    box on
    hold off

%% Plot average maximum correlation versus median window channel deviation
    minWinHF = min(statistics(:, s.medWinHF));
    maxWinHF = max(statistics(:, s.medWinHF));
    minCorr = min(statistics(:, s.medCor));
    baseTitle = 'Average maximum correlation versus median window HF noise';
    figure ('Name', baseTitle, 'Color', [1, 1, 1]);
    hold on
    plot(statistics(:, s.aveCor), statistics(:, s.medWinHF), 'ok')

    xlabel(statisticsTitles{s.aveCor});
    ylabel(statisticsTitles{s.medWinHF});
    medWinHF = median(statistics(:, s.medWinHF));
    aveCorr = median(statistics(:, s.aveCor));
    plot(aveCorr, medWinHF, '+r', 'MarkerSize', 14, 'LineWidth', 3);
    set(gca, 'XLim', [minCorr, 1], 'XLimMode', 'manual', ...
        'YLim', [minWinHF, maxWinHF], 'YLimMode', 'manual');
    title({collectionTitle; baseTitle; ...
        ['[Median of dataset medians: ' num2str(aveCorr) '(corr) ' ...
        num2str(medWinHF) '(med HF)]']}, 'interpreter', 'none');
    legend('Dataset', 'Collection median', 'Location', 'NorthWest')
    box on
    hold off

    %% Plot average maximum correlation versus window channel HF noise ratio
    devRatio = statistics(:, s.rSDWinHF)./statistics(:, s.medWinHF);
    minWinRatio = min(devRatio);
    maxWinRatio = max(devRatio);
    yLabelString = 'rSDR/median window channel HF noise';
    minCorr = min(statistics(:, s.aveCor));
    baseTitle = 'Average maximum correlation versus window channel HF noise ratio';
    figure ('Name', baseTitle, 'Color', [1, 1, 1]);
    hold on
    plot(statistics(:, s.aveCor), devRatio, 'ok')
    xlabel(statisticsTitles{s.aveCor});
    ylabel(yLabelString);
    medRatio = median(devRatio);
    aveCorr = median(statistics(:, s.aveCor));
    plot(aveCorr, medRatio, '+r', 'MarkerSize', 14, 'LineWidth', 3);
    set(gca, 'XLim', [minCorr, 1], 'XLimMode', 'manual', ...
        'YLim', [minWinRatio, maxWinRatio], 'YLimMode', 'manual');
    title({collectionTitle; baseTitle; ...
        ['[Median of dataset medians: ' num2str(aveCorr) '(corr) ' ...
        num2str(medRatio) '(rSD/med HF)]']}, 'interpreter', 'none');
    legend('Dataset', 'Collection median', 'Location', 'NorthWest')
    box on
    hold off

author: VisLab

repo: EEG-Clean-Tools


text:
function G = Kernel_Integration_Approx(dt, para)

G = zeros(length(dt(:)), size(para.g,2));

M = size(para.g,1);
Nums = ceil(dt./para.dt);
for i = 1:length(dt(:))
    if Nums(i)<=M
        G(i,:) = sum(para.g(1:Nums(i),:)).*para.dt;
    else
        G(i,:) = sum(para.g).*para.dt;
    end
end

author: HongtengXu

repo: Hawkes-Process-Toolkit


text:
function opt_option_out = fcnChooseOption(opt_option,tol_opt, u0)
% specify and configure optimization method
%refer to goo.gl/aTQFwr
%(http://www.mathworks.com/help/optim/ug/fmincon.html#inputarg_options)
    warning off all;
    if ( opt_option == 0 )
        opt_option_out = optimset('Display','off',...
            'TolFun', tol_opt,...
            'MaxIter', 10000,...
            'Algorithm', 'active-set',...
            'FinDiffType', 'forward',...
            'RelLineSrchBnd', [],...
            'RelLineSrchBndDuration', 1,...
            'TolConSQP', 1e-6);
    elseif ( opt_option == 1 )
        opt_option_out = optimset('Display','off',...
            'TolFun', tol_opt,...
            'MaxIter', 2000,...
            'Algorithm', 'interior-point',...
            'AlwaysHonorConstraints', 'bounds',...
            'FinDiffType', 'forward',...
            'HessFcn', [],...
            'Hessian', 'bfgs',...
            'HessMult', [],...
            'InitBarrierParam', 0.1,...
            'InitTrustRegionRadius', sqrt(size(u0,1)*size(u0,2)),...
            'MaxProjCGIter', 2*size(u0,1)*size(u0,2),...
            'ObjectiveLimit', -1e20,...
            'ScaleProblem', 'obj-and-constr',...
            'SubproblemAlgorithm', 'cg',...
            'TolProjCG', 1e-2,...
            'TolProjCGAbs', 1e-10);
    %                       'UseParallel','always',...
    elseif ( opt_option == 2 )
        opt_option_out = optimset('Display','off',...
            'TolFun', tol_opt,...
            'MaxIter', 2000,...
            'Algorithm', 'trust-region-reflective',...
            'Hessian', 'off',...
            'MaxPCGIter', max(1,floor(size(u0,1)*size(u0,2)/2)),...
            'PrecondBandWidth', 0,...
            'TolPCG', 1e-1);
    end

end

author: juchengquan


text:
classdef MimImageDatabaseSeries < handle
    % MimImageDatabaseSeries.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties (Constant)
        CurrentVersionNumber = 3
    end

    properties (SetAccess = private)
        Name
        StudyName
        ImageMap
        Modality
        Date
        Time
        SeriesUid
        PatientId

        Version
        FirstImagePath
        ImageFileFormat
        StudyUid

    end

    methods
        function obj = MimImageDatabaseSeries(series_id, single_image_metainfo)
            if nargin > 0
                obj.Name = single_image_metainfo.SeriesDescription;
                obj.StudyName = single_image_metainfo.StudyDescription;
                obj.Modality = single_image_metainfo.Modality;
                obj.Date = single_image_metainfo.Date;
                obj.Time = single_image_metainfo.Time;
                obj.SeriesUid = series_id;
                obj.PatientId = single_image_metainfo.PatientId;

                obj.FirstImagePath = single_image_metainfo.ImagePath;
                obj.ImageFileFormat = single_image_metainfo.ImageFileFormat;
                obj.StudyUid = single_image_metainfo.StudyUid;

                obj.ImageMap = containers.Map;

                obj.Version = obj.CurrentVersionNumber;
            end
        end

        function visible_path = GetVisiblePath(obj)
            if obj.ImageMap.Count == 1
                filenames = obj.ImageMap.values;
                first_filename = filenames{1};
                visible_path = first_filename.FullFile;
            else
                visible_path = obj.FirstImagePath;
            end
        end

        function AddImage(obj, single_image_metainfo)
            obj.ImageMap(single_image_metainfo.ImageUid) = CoreFilename(single_image_metainfo.ImagePath, single_image_metainfo.ImageFilename);
        end

        function num_images = NumberOfImages(obj)
            num_images = obj.ImageMap.Count;
        end

        function image_info = GetImageInfo(obj)
            filenames_cells = [];
            filenames_objects = obj.ImageMap.values;
            for image_index = 1 : numel(filenames_objects)
                filenames_cells{end + 1} = filenames_objects{image_index};
            end
            path = obj.FirstImagePath;
            image_type = obj.ImageFileFormat;
            uid = obj.SeriesUid;
            study_uid = obj.StudyUid;
            modality = obj.Modality;
            image_info = PTKImageInfo(path, filenames_cells, image_type, uid, study_uid, modality);
        end

    end

    methods (Static)
        function obj = loadobj(a)
            % This method is called when the object is loaded from disk.

            if isa(a, 'MimImageDatabaseSeries')
                obj = a;
            else
                % In the case of a load error, loadobj() gives a struct
                obj = MimImageDatabaseSeries;
                for field = fieldnames(a)'
                    if isprop(obj, field{1})
                        mp = findprop(obj, (field{1}));
                        if (~mp.Constant) && (~mp.Dependent) && (~mp.Abstract)
                            obj.(field{1}) = a.(field{1});
                        end
                    end
                end
            end

        end
    end

end

author: tomdoel

repo: pulmonarytoolkit


text:
function [boxesOut uniqueIdx] = BoxRemoveDuplicates(boxesIn)
% function boxOut = BoxRemoveDuplicates(boxIn)
%
% Removes duplicate boxes. Leaves the boxes in the same order
% Keeps the first box of each kind.
%
% boxesIn:          N x 4 array containing boxes
%
% boxexOut:         M x 4 array of boxes witout duplicates. M <= N
% uniqueIdx:        Indices of retained boxes from boxesIn
%
%     Jasper Uijlings - 2013

[dummy uniqueIdx] = unique(boxesIn, 'rows', 'first');
uniqueIdx = sort(uniqueIdx);
boxesOut = boxesIn(uniqueIdx,:);

author: zouchuhang

repo: LayoutNet


text:
function HD = VARhd(VAR)
% =======================================================================
% Compute the historical decomposition of the time series in a VAR
% estimated with VARmodel and identified with VARir/VARfevd
% =======================================================================
% HD = VARhd(VAR)
% -----------------------------------------------------------------------
% INPUTS
%   - VAR: structure, result of VARmodel -> VARir/VARfevd function
% -----------------------------------------------------------------------
% OUTPUT
%   - HD: structure including the historical decomposition
% =======================================================================
% Ambrogio Cesa Bianchi, March 2015
% ambrogio.cesabianchi@gmail.com

% I thank Andrey Zubarev for finding a bug in the contribution of the
% exogenous variables when nvar_ex~=0 and nlag_ex>0.

%% Check inputs
%===============================================
if ~exist('VAR','var')
    error('You need to provide VAR structure, result of VARmodel');
end
% If there is VAR check that the inverse of the A matrix is not empty
invA = VAR.invA;
if isempty(invA)
    error('You need to identify the VAR before running VARhd. Run VARir/VARfevd first.');
end

%% Retrieve and initialize variables
%===============================================
Fcomp   = VAR.Fcomp;                      % Companion matrix
const   = VAR.const;                      % constant and/or trends
F       = VAR.Ft';                        % make comparable to notes
eps     = invA\transpose(VAR.residuals);  % structural errors
nvar    = VAR.nvar;                       % number of endogenous variables
nvar_ex = VAR.nvar_ex;                    % number of exogenous (excluding constant and trend)
nvarXeq = VAR.nvar * VAR.nlag;            % number of lagged endogenous per equation
nlag    = VAR.nlag;                       % number of lags
nlag_ex = VAR.nlag_ex;                    % number of lags of the exogenous
Y       = VAR.Y;                          % left-hand side
X       = VAR.X(:,1+const:nvarXeq+const); % right-hand side (no exogenous)
nobs    = size(Y,1);                      % number of observations

%% Compute historical decompositions
%===============================================

% Contribution of each shock
    invA_big = zeros(nvarXeq,nvar);
    invA_big(1:nvar,:) = invA;
    Icomp = [eye(nvar) zeros(nvar,(nlag-1)*nvar)];
    HDshock_big = zeros(nlag*nvar,nobs+1,nvar);
    HDshock = zeros(nvar,nobs+1,nvar);
    for j=1:nvar; % for each variable
        eps_big = zeros(nvar,nobs+1); % matrix of shocks conformable with companion
        eps_big(j,2:end) = eps(j,:);
        for i = 2:nobs+1
            HDshock_big(:,i,j) = invA_big*eps_big(:,i) + Fcomp*HDshock_big(:,i-1,j);
            HDshock(:,i,j) =  Icomp*HDshock_big(:,i,j);
        end
    end

% Initial value
    HDinit_big = zeros(nlag*nvar,nobs+1);
    HDinit = zeros(nvar, nobs+1);
    HDinit_big(:,1) = X(1,:)';
    HDinit(:,1) = Icomp*HDinit_big(:,1);
    for i = 2:nobs+1
        HDinit_big(:,i) = Fcomp*HDinit_big(:,i-1);
        HDinit(:,i) = Icomp *HDinit_big(:,i);
    end

% Constant
    HDconst_big = zeros(nlag*nvar,nobs+1);
    HDconst = zeros(nvar, nobs+1);
    CC = zeros(nlag*nvar,1);
    if const>0
        CC(1:nvar,:) = F(:,1);
        for i = 2:nobs+1
            HDconst_big(:,i) = CC + Fcomp*HDconst_big(:,i-1);
            HDconst(:,i) = Icomp * HDconst_big(:,i);
        end
    end

% Linear trend
    HDtrend_big = zeros(nlag*nvar,nobs+1);
    HDtrend = zeros(nvar, nobs+1);
    TT = zeros(nlag*nvar,1);
    if const>1;
        TT(1:nvar,:) = F(:,2);
        for i = 2:nobs+1
            HDtrend_big(:,i) = TT*(i-1) + Fcomp*HDtrend_big(:,i-1);
            HDtrend(:,i) = Icomp * HDtrend_big(:,i);
        end
    end

% Quadratic trend
    HDtrend2_big = zeros(nlag*nvar, nobs+1);
    HDtrend2 = zeros(nvar, nobs+1);
    TT2 = zeros(nlag*nvar,1);
    if const>2;
        TT2(1:nvar,:) = F(:,3);
        for i = 2:nobs+1
            HDtrend2_big(:,i) = TT2*((i-1)^2) + Fcomp*HDtrend2_big(:,i-1);
            HDtrend2(:,i) = Icomp * HDtrend2_big(:,i);
        end
    end

% Exogenous
    HDexo_big = zeros(nlag*nvar,nobs+1);
    HDexo = zeros(nvar,nobs+1);
    EXO = zeros(nlag*nvar,nvar_ex*(nlag_ex+1));
    if nvar_ex>0;
        VARexo = VAR.X_EX;
        EXO(1:nvar,:) = F(:,nvar*nlag+const+1:end); % this is c in my notes
        for i = 2:nobs+1
            HDexo_big(:,i) = EXO*VARexo(i-1,:)' + Fcomp*HDexo_big(:,i-1);
            HDexo(:,i) = Icomp * HDexo_big(:,i);
        end
    end

% All decompositions must add up to the original data
HDendo = HDinit + HDconst + HDtrend + HDtrend2 + HDexo + sum(HDshock,3);

%% Save and reshape all HDs
%===============================================
HD.shock = zeros(nobs+nlag,nvar,nvar);  % [nobs x shock x var]
    for i=1:nvar
        for j=1:nvar
            HD.shock(:,j,i) = [nan(nlag,1); HDshock(i,2:end,j)'];
        end
    end
HD.init   = [nan(nlag-1,nvar); HDinit(:,1:end)'];    % [nobs x var]
HD.const  = [nan(nlag,nvar);   HDconst(:,2:end)'];   % [nobs x var]
HD.trend  = [nan(nlag,nvar);   HDtrend(:,2:end)'];   % [nobs x var]
HD.trend2 = [nan(nlag,nvar);   HDtrend2(:,2:end)'];  % [nobs x var]
HD.exo    = [nan(nlag,nvar);   HDexo(:,2:end)'];     % [nobs x var]
HD.endo   = [nan(nlag,nvar);   HDendo(:,2:end)'];    % [nobs x var]

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function rules(nx,ny,nz,alpha)

ssz=get(0,'ScreenSize');

hf=figure;
x0=ssz(3)/2;
y0=ssz(4)/2;
kfsz=0.9;
fpos=[x0-kfsz*ssz(3)/2 y0-kfsz*ssz(4)/2 kfsz*ssz(3) kfsz*ssz(4)];
set(hf,'MenuBar','none','Color',[1 1 1],'name','rules','NumberTitle','off','units','pixels','position',fpos);

ha = axes('parent',hf);
set(ha,'Visible','off');

set(ha,'UserData',[nx,ny,nz,alpha]);

% menu:
 hm=uimenu(hf,'Label','Save Image / Print');
    uimenu(hm,'Label','Save Image As...','Callback',['save_im(' num2str(hf) ')']);
    uimenu(hm,'Label','Page Setup...','Callback',['pgs(' num2str(hf) ')'],'Separator','on');
    uimenu(hm,'Label','Print Preview...','Callback',['pp(' num2str(hf) ')']);
    uimenu(hm,'Label','Print...','Callback',['pd('  num2str(hf) ')']);

euler_rule(ha);

h1 = uicontrol('Style', 'pushbutton', 'String', 'Euler to axis angle',...
    'Position', [0.25*fpos(3) 0.03*fpos(4) 0.2*fpos(4) 0.04*fpos(4)], 'Callback', ['euler_rule(' num2str(ha,'%20.20e') ')']);

h2 = uicontrol('Style', 'pushbutton', 'String', 'axis angle to Euler',...
    'Position', [0.45*fpos(3) 0.03*fpos(4) 0.2*fpos(4) 0.04*fpos(4)], 'Callback', ['axan_rule(' num2str(ha,'%20.20e') ')']);

h3 = uicontrol('Style', 'pushbutton', 'String', 'quatenion form',...
    'Position', [0.65*fpos(3) 0.03*fpos(4) 0.2*fpos(4) 0.04*fpos(4)], 'Callback', ['quat_rule(' num2str(ha,'%20.20e') ')']);


text:
function [im_x, im_y] = project( camera, world_X, world_Y, world_Z )
%PROJECT: project a 3D point into an image
%
%   [IM_X,IM_Y] = PROJECT(CAMERA,WORLD_X,WORLD_Y,WORLD_Z) projects one or
%   more 3D point in the world coordinate frame into image coordinates

%   Copyright 2005-2009 The MathWorks, Inc.
%   $Revision: 1.0 $    $Date: 2006/06/30 00:00:00 $

z = camera.rawP(3,1) * world_X ...
    + camera.rawP(3,2) * world_Y ...
    + camera.rawP(3,3) * world_Z ...
    + camera.rawP(3,4);
im_y = round( (camera.rawP(2,1) * world_X ...
    + camera.rawP(2,2) * world_Y ...
    + camera.rawP(2,3) * world_Z ...
    + camera.rawP(2,4)) ./ z);
im_x = round( (camera.rawP(1,1) * world_X ...
    + camera.rawP(1,2) * world_Y ...
    + camera.rawP(1,3) * world_Z ...
    + camera.rawP(1,4)) ./ z);


text:
function [d, xd, bw] = kernelDensity(x, bins, h, kernel)
%KERNELDENSITY Calculate the kernel density of a data set
%
%   [D, XD] = IOSR.STATISTICS.KERNELDENSITY(X) calculate the kernel density
%   D of a dataset X for query points XD. The kernel density is calculated
%   for 100 query points equally spaced between the minimum and maximum of
%   the data X. The density is estimated using a gaussian kernel with a
%   width that is optimal for normal data. X may be a vector, matrix, or
%   multi-dimensional array; the entire array is treated as the sample. D
%   and XD are 100-point column vectors. NaN are excluded from the
%   calculations.
%
%   ... = IOSR.STATISTICS.KERNELDENSITY(X, BINS) calculates the density for
%   the query points specified by BINS. If BINS is a scalar, then BINS
%   points are queried between MIN(X(:)) and MAX(X(:)); if BINS is a
%   vector, then the values are used as the query points directly. D and XD
%   are column vectors.
%
%   ... = IOSR.STATISTICS.KERNELDENSITY(X, BINS, H) uses the bandwidth H to
%   calculate the kernel density. H must be a scalar. BINS may be an empty
%   array in order to use the default described above.
%
%   ... = IOSR.STATISTICS.KERNELDENSITY(X, BINS, [], KERNEL) uses the
%   kernel function specified by KERNEL to calculate the density. The
%   kernel may be:
%       - 'normal' (default),
%       - 'uniform',
%       - 'triangular',
%       - 'epanechnikov',
%       - 'quartic',
%       - 'triweight',
%       - 'tricube',
%       - 'cosine',
%       - 'logistic',
%       - 'sigmoid', or
%       - 'silverman'.
%   For the uniform case the bandwidth is set to 15% of the range of the
%   data [1]. Otherwise the bandwidth is chosen to be optimal for normal
%   data assuming a gaussian kernel.
%
%   ... = IOSR.STATISTICS.KERNELDENSITY(X, BINS, H, KERNEL) allows the
%   bins BINS and bandwidth H to be specified directly.
%
%   [D, XD, BW] = IOSR.STATISTICS.KERNELDENSITY(...) returns the badwidth
%   BW.
%
%   Examples
%
%     Example 1: Plot the kernel density of gaussian data
%       figure
%       % gaussian random numbers
%       y = randn(100000, 1);
%       % density
%       [d, xd] = iosr.statistics.kernelDensity(y);
%       % plot
%       plot(xd, d);
%
%     Example 2: Density trace with 200 bins of width of 10% of data range
%       figure
%       % random numbers
%       y = randn(100000, 1);
%       % y range
%       range = max(y(:)) - min(y(:));
%       % density trace
%       [d, xd] = iosr.statistics.kernelDensity(y, 200, 0.1*range, 'uniform');
%       % plot
%       plot(xd, d);
%
%   References
%
%   [1] Hintze, Jerry L.; Nelson, Ray D. (1998). "Violin Plots: A Box
%       Plot-Density Trace Synergism". The American Statistician. 52 (2):
%       181?4.

    %% input check

    x = x(:);
    x = x(~isnan(x));
    assert(numel(x) > 1, 'X must be a vector, matrix, or array.')

    % x bins
    if nargin < 2
        bins = [];
    end
    if isempty(bins)
        bins = 100;
    end
    if isscalar(bins)
        bins = linspace(min(x), max(x), round(bins));
    end
    bins = bins(:);

    % bin width
    if nargin < 3
        h = [];
    end

    % kernel
    if nargin < 4
        kernel = [];
    end
    if isempty(kernel)
        kernel = 'normal';
    end
    % return kernel function
    switch lower(kernel)
        case 'uniform'
            K = @(u) 0.5*(abs(u) <= 1);
            if isempty(h)
                h = 0.15 * (max(x) - min(x));
            end
        case 'normal'
            K = @(u) ((1/sqrt(2*pi)) * exp(-0.5*(u.^2)));
        case 'triangular'
            K = @(u) ((1-abs(u)) .* (abs(u) <= 1));
        case 'epanechnikov'
            K = @(u) ((0.75*(1-(u.^2))) .* (abs(u) <= 1));
        case 'quartic'
            K = @(u) (((15/16)*(1-(u.^2)).^2) .* (abs(u) <= 1));
        case 'triweight'
            K = @(u) (((35/32)*(1-(u.^2)).^3) .* (abs(u) <= 1));
        case 'tricube'
            K = @(u) (((70/81)*(1-(abs(u).^3)).^3) .* (abs(u) <= 1));
        case 'cosine'
            K = @(u) (((pi/4)*cos((pi/2)*u)) .* (abs(u) <= 1));
        case 'logistic'
            K = @(u) (1 / (exp(u) + 2 + exp(-u)));
        case 'sigmoid'
            K = @(u) ((2/pi) * (1 / (exp(u) + exp(-u))));
        case 'silverman'
            K = @(u) (0.5 * exp((-abs(u))/(sqrt(2))) .* sin((abs(u))/(sqrt(2)) + (pi/4)));
        otherwise
            error('Unknown kernel specified');
    end
    if isempty(h)
        h = ((4*(std(x).^5))/(3*numel(x))).^(1/5);
    end

    assert(isscalar(h), 'h must be a scalar')

    %% calculate kernel density

    xd = sort(bins);
    d = zeros(size(xd));

    for i = 1:numel(xd)
        d(i) = sum(K((x-xd(i))/h))./(numel(x)*h);
    end
    d(isnan(d)) = 0;
    bw = h;

end

author: IoSR-Surrey

repo: MatlabToolbox


text:
d = 500;
A = randn(d);
b = randn(d,1);
niter = 100;
fprintf('cols:')
tic; for i = 1:niter scale_cols(A,b); end; toc
fprintf('rows:')
tic; for i = 1:niter scale_rows(A,b); end; toc

author: andrewssobral

repo: lrslibrary


text:
addpath('testImage_Video/')
addpath('model/')
%for pic
load('cars_meta.mat');
load('cifar10NetRCNN.mat') %for detect
load('AlexNet_New.mat');%for recognition

%input image:
frame=imread('6.jpg');
%crop:
frame=imresize(frame,[480 640]);
outputImage=frame;
   % DetectRCNN(frame,cifar10NetRCNN);
   [bboxes, ~, ~] = detect(cifar10NetRCNN, frame);

   if  ~isempty(bboxes)
   size_array=size(bboxes);
   length=size_array(1);
   for i=1:length
       box=bboxes(i,:);
       frame_=imcrop(frame,box);
       frame_=imresize(frame_,[227 227]);
       type_num=classify(AlexNet_New,frame_);
       outputImage=insertObjectAnnotation(outputImage, 'rectangle', box, class_names{type_num},'LineWidth',3);
   end
   end
   imshow(outputImage);

author: chenjoya


text:
function list = dirext(path, varargin)
% DIREXT Extended version of standard dir command.
% ------------------------------------------------------------------------------
% SYNTAX
% 1 list = dirext(path)
%     Lists all files and folders contained in path.
%
% 2 list = dirext(path, Filter)
%     With Filter a filtering of the files and folders can be performed.
% ------------------------------------------------------------------------------
% INPUT
% 1 path
%     Search path.
%
% 2 [Filter]
%     String or cell array of strings containing filter informations. All
%     options offered by the original dir command are possible. Among them
%     possible choices are:
%       * 'isdir'    Output list contains just folders.
% ------------------------------------------------------------------------------
% OUTPUT
% 1 list
%     List of files and/or folders, found in path. Result is cell array.
% ------------------------------------------------------------------------------
% EXAMPLES
% 1 List all files and folders in current folder.
%     dirext('*')
%
% 2 List all files and folders in Y:\Temp\.
%     dirext('Y:\Temp')
%
% 3 List just folders in Y:\Temp\.
%     dirext('Y:\Temp', 'isdir')
% ------------------------------------------------------------------------------
% philipp.glira@gmail.com
% ------------------------------------------------------------------------------

% Input parsing ----------------------------------------------------------------

p = inputParser;
p.addRequired('path',     @ischar);
p.addOptional('Filter', [], @(x) ischar(x) || iscell(x));
p.parse(path, varargin{:});
p = p.Results;
% Clear required inputs to avoid confusion
clear path

% Create raw list of files/folders ---------------------------------------------

listStruct = dir(p.path);

% Error if no files or folders are found
if isempty(listStruct), error('No files or folders found.'); end

% Delete two first entries ('.' and '..') if present
idxLogPoint  = strcmp({listStruct.name}, '.');  % where is point?
idxLogDPoint = strcmp({listStruct.name}, '..'); % where is double point?
idxLog = idxLogPoint | idxLogDPoint;
% Delete if found
if any(idxLog), listStruct(idxLog) = []; end

% Filtering of results? --------------------------------------------------------

if ~isempty(p.Filter)

    % Just folders
    if strcmpi(p.Filter, 'isdir')
        LogIdx = ([listStruct.isdir] == 1);
        listStruct = listStruct(LogIdx);
    end

end

list = {listStruct.name};

author: pglira


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

%%
clear; close all; clc;

%% Syntax
% |P=rayleighpdf(M,S);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<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


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
function Xloci1 = get_X_cell_efficient(hull,Xs,ind_Xs,d_Xs,p,n);

temps1 = [];
for j=1:p
    temps1(j) = d_Xs{j}(hull(j));
end

Xloci1 = ones(n,prod(temps1));

for j=1:prod(temps1)
    temp = ind2subv(temps1,j);
    for k=1:p
        Xloci1(:,j) = Xloci1(:,j) .* Xs(:,ind_Xs{k}{hull(k)}(temp(k)));
    end
end

repo: REST


text:
%kckhueshift 'shifts the hues on a Hue-only (1 band) image '
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros ckhueshift.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'First Input data object'
% Double: shift 'Shift by', default: 0: 'shift float value'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kckhueshift(i, {'i','';'shift',0;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% ckhueshift - shifts the hues on a Hue-only (1 band) image
%
%  DESCRIPTION
% Adds a constant value to each value on the object, correcting for circularity. This will make a "hue shift" on the object.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
% Expects a 1-band object.
%
%  REFERENCES
% All references for the COLOR toolbox are on the $COLOR/Readme.html file.
%
%  COPYRIGHT
% Copyright (C) 1997 Rafael Santos. Khoros (C) Khoral Research, Inc.
%

function varargout = kckhueshift(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kckhueshift(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'shift', 0;'o', '__output'};
maxval={0,360,0};
minval={0,-360,0};
istoggle=[0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','Double','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'ckhueshift"  '],Inputs,narglist);


text:
classdef TopologyMonitoring < handle

    properties (Access = public)
        type
        designVariable
        cost
        constraint
        dualVariable
        historyPrinter
        convergenceVars
        monitorDocker
        postProcess
        incrementalScheme
        monitor
    end

    properties (Access = private)
        nIter
        hasFinished
        lineSearch
        lineSearchTrials
    end

    methods (Access = public)

        function obj = TopologyMonitoring(cParams)
            obj.init(cParams)
        end

        function compute(obj,cParams)
            switch obj.type
                case 'fmincon'
                    obj.plotFmincon(cParams);
                case 'NullSpace'
                    obj.plotNullSpace(cParams);
                case 'AlternatingPrimalDual'
                    obj.plotAugmentedLagrangian(cParams);
                case 'DualNestedInPrimal'
                    obj.plotBisection(cParams);
                case 'IPOPT'
                    obj.plotIPOPT(cParams);
                case 'MMA'
                    obj.plotMMA(cParams);
                otherwise
                    error('Optimizer not implemented')
            end
        end

        function create(obj,cParams)
            obj.createHistoryPrinter(cParams.historyPrinterSettings);
            obj.createConvergenceVariables(cParams);
            obj.createMonitorDocker(cParams.monitoringDockerSettings);
            obj.createPostProcess(cParams.postProcessSettings);
        end

    end

    methods (Access = private)

        function init(obj,cParams)
            obj.type     = cParams.type;
            obj.cost              = cParams.cost;
            obj.constraint        = cParams.constraint;
            obj.designVariable    = cParams.designVar;
            obj.dualVariable      = cParams.dualVariable;
            obj.incrementalScheme = cParams.incrementalScheme;
        end

        function plotFmincon(obj,cParams)
            obj.dualVariable.value = zeros(obj.constraint.nSF,1);
            obj.nIter              = cParams.nIter;
            obj.hasFinished        = cParams.hasFinished;
            foOpt                  = cParams.firstorderopt;
            normXsquare            = obj.designVariable.computeL2normIncrement();
            incX                   = sqrt(normXsquare);
            obj.designVariable.updateOld();
            switch cParams.algorithm
                case 'sqp'
                    stepL = cParams.stepsize;
                case 'interior-point'
                    stepL = cParams.trustregionradius;
                otherwise
            end
            obj.printOptimizerVariable();
            obj.convergenceVars.reset();
            obj.convergenceVars.append(incX);
            obj.convergenceVars.append(foOpt);
            obj.convergenceVars.append(stepL);
            obj.refreshMonitoring();
            obj.printHistory();
        end

        function plotNullSpace(obj,cParams)
            deltaCost              = obj.cost.value - cParams.oldCost;
            obj.nIter              = cParams.nIter;
            normXsquare            = obj.designVariable.computeL2normIncrement();
            obj.lineSearch         = cParams.tau;
            obj.lineSearchTrials   = cParams.lineSearchTrials;
            obj.hasFinished        = cParams.hasFinished;
            incX                   = sqrt(normXsquare);
            obj.designVariable.updateOld();
            obj.printOptimizerVariable();
            obj.convergenceVars.reset();
            obj.convergenceVars.append(deltaCost);
            obj.convergenceVars.append(incX);
            obj.convergenceVars.append(obj.lineSearch);
            obj.convergenceVars.append(obj.lineSearchTrials);
            obj.convergenceVars.append(cParams.meritNew);
            obj.refreshMonitoring();
            obj.printHistory();
        end

        function plotAugmentedLagrangian(obj,cParams)
            deltaCost              = obj.cost.value - cParams.oldCost;
            obj.nIter              = cParams.nIter;
            normXsquare            = obj.designVariable.computeL2normIncrement();
            obj.lineSearch         = cParams.tau;
            obj.lineSearchTrials   = cParams.lineSearchTrials;
            obj.hasFinished        = cParams.hasFinished;
            incX                   = sqrt(normXsquare);
            obj.designVariable.updateOld();
            obj.printOptimizerVariable();
            obj.convergenceVars.reset();
            obj.convergenceVars.append(deltaCost);
            obj.convergenceVars.append(incX);
            obj.convergenceVars.append(obj.lineSearch);
            obj.convergenceVars.append(obj.lineSearchTrials);
            obj.convergenceVars.append(cParams.meritNew);
            obj.refreshMonitoring();
            obj.printHistory();
        end

        function plotBisection(obj,cParams)
            deltaCost              = obj.cost.value - cParams.oldCost;
            obj.nIter              = cParams.nIter;
            normXsquare            = obj.designVariable.computeL2normIncrement();
            obj.lineSearch         = cParams.tau;
            obj.lineSearchTrials   = cParams.lineSearchTrials;
            obj.hasFinished        = cParams.hasFinished;
            incX                   = sqrt(normXsquare);
            obj.designVariable.updateOld();
            obj.printOptimizerVariable();
            obj.convergenceVars.reset();
            obj.convergenceVars.append(deltaCost);
            obj.convergenceVars.append(incX);
            obj.convergenceVars.append(obj.lineSearch);
            obj.convergenceVars.append(obj.lineSearchTrials);
            obj.refreshMonitoring();
            obj.printHistory();
        end

        function plotIPOPT(obj,cParams)
            obj.hasFinished = cParams.hasFinished;
            obj.nIter       = cParams.nIter;
            normXsquare     = obj.designVariable.computeL2normIncrement();
            obj.printOptimizerVariable();
            obj.convergenceVars.reset();
            obj.convergenceVars.append(cParams.inf_pr);
            obj.convergenceVars.append(cParams.inf_du);
            obj.convergenceVars.append(sqrt(normXsquare));
            obj.refreshMonitoring();
            obj.printHistory();
        end

        function plotMMA(obj,cParams)
            obj.hasFinished = cParams.hasFinished;
            obj.nIter       = cParams.nIter;
            obj.printOptimizerVariable();
            obj.convergenceVars.reset();
            obj.convergenceVars.append(cParams.KKTnorm);
            obj.convergenceVars.append(cParams.outitFrac);
            obj.refreshMonitoring();
            obj.printHistory();
        end

        function createHistoryPrinter(obj,cParams)
            cParams.optimizer  = obj;
            cParams.cost       = obj.cost;
            cParams.constraint = obj.constraint;
            obj.historyPrinter = OptimizationMetricsPrinterFactory.create(cParams);
        end

        function createConvergenceVariables(obj,cParams)
            s                   = cParams.optimizerNames;
            cVarD               = ConvergenceVarsDispatcher.dispatchNames(s);
            n                   = numel(cVarD);
            cVar                = ConvergenceVariables(n);
            obj.convergenceVars = cVar;
        end

        function createMonitorDocker(obj,s)
            s.designVariable  = obj.designVariable;
            s.dualVariable    = obj.dualVariable;
            s.cost            = obj.cost;
            s.constraint      = obj.constraint;
            s.convergenceVars = obj.convergenceVars;
            obj.monitorDocker = MonitoringDocker(s);
        end

        function createPostProcess(obj,cParams)
            if cParams.shallPrint
                d                  = obj.createPostProcessDataBase(cParams);
                d.printMode        = cParams.printMode;
                d.nDesignVariables = obj.designVariable.nVariables;
                obj.postProcess    = Postprocess('TopOptProblem',d);
            end
        end

        function d = createPostProcessDataBase(obj,cParams)
            d.mesh    = obj.designVariable.mesh;
            d.outName = cParams.femFileName;
            d.ptype   = cParams.ptype;
            ps = PostProcessDataBaseCreator(d);
            d  = ps.create();
            d.pdim    = cParams.pdim;
            d.ndim    = obj.designVariable.mesh.ndim;
            d.optimizer  = obj.type;
            d.cost       = obj.cost;
            d.constraint = obj.constraint;
            d.designVar  = obj.designVariable.type;
        end

        function printOptimizerVariable(obj)
            if ~isempty(obj.postProcess)
            d.fields     = obj.designVariable.getVariablesToPlot();
            d.cost       = obj.cost;
            d.constraint = obj.constraint;
%             obj.postProcess.print(obj.nIter,d);
            end
        end

        function printHistory(obj)
            iStep = obj.incrementalScheme.iStep;
            obj.historyPrinter.print(obj.nIter,iStep);
        end

        function printHistoryFinalValues(obj)
            obj.historyPrinter.printFinal();
        end

        function refreshMonitoring(obj)
            iStep = obj.incrementalScheme.iStep;
            nStep = obj.incrementalScheme.nSteps;
            obj.monitorDocker.refresh(obj.nIter,obj.hasFinished,iStep,nStep);
        end

    end

end

author: SwanLab


text:
function p = addMultipliedEqualityCuts(p)

newRows = [];
for i =  p.linears
    for j = 1:p.K.f
        row = p.F_struc(j,:);
        used = find(row(2:end));
        if ~any(p.variabletype(used))
            % Linear constraint
            ok = 1;
            for k = 1:length(used)
                find1 = findrows(p.bilinears(:,2:end),[i used(k)]);
                find2 = findrows(p.bilinears(:,2:end),[used(k) i]);
                if (isempty(find1) & isempty(find2))
                    ok = 0;
                    break
                else
                    monom(k) = [find1 find2];
                end
            end
            if ok
                monom = p.bilinears(monom,1);
                row(1+monom) = row(1+used);
                row(1+used) = 0;
                row(i) = row(1);
                row(1)=0;
                newRows = [newRows;row];
            end
        end
    end
end
if ~isempty(newRows)
    p.F_struc = [newRows;p.F_struc];
    p.K.f = p.K.f + size(newRows,1);
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function R = strsimilarity(a,b)
%STRSIMILARITY Similarity measure between two character vectors.
%   R = STRSIMILARITY(A,B) computes the similarity measure, R, defined
%   in Eq. (14-22) between character vectors A and B. The vectors do not
%   have to be of the same length, but they should be otherwise
%   registered for this method to make sense. All blanks in both vectors
%   are deleted, so they should not be used as valid characters in
%   defining A and B. Inputs can be either character vectors or scalar
%   strings, in the sense defined by MATLAB. If the inputs are strings
%   they are converted to character vectors.
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

% If inputs are scalar strings, covert them to character vectors. If
% they already are character vectors, they are not affected by the
% following function.
if (isstring(a) && isstring(b)) || (ischar(a) && ischar(b))
   a = convertStringsToChars(a);
   b = convertStringsToChars(b);
else
   error('a and b must both character vectors or both be strings')
end

% Make sure the inputs are character vectors and not higher-dimensional
% character arrays.
if ~(size(a,1) == 1 || size(a,2) == 1) && ~(size(b,1) == 1 ||...
                                                         size(b,2) == 1)
   error('a and b must be character vectors')
end

% Make sure they are rows for dimensional consistency later on.
a = a(:)';
b = b(:)';

% Remove all blanks.
a = a(~isspace(a));
b = b(~isspace(b));

% Pad the end of the shorter string.
La = length(a);
Lb = length(b);
if La > Lb
   b = [b,blanks(La - Lb)];
else
   a = [a,blanks(Lb - La)];
end

% Compute the similarity measure.
I = find(a == b);
alpha = numel(I);
den = max(La,Lb) - alpha;
R = alpha/den;

author: dipum

repo: dipum-toolbox

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