text:
%% STUDY_02_wrist_scaffold_single_fiber_iFEA_viscoelastic
% Below is a demonstration for:
%
% * Inverse FEA based material parameter optimisation
% * febio_spec version 3.0
% * febio, FEBio
% * hexahedral elements, hex8
% * static, solid
% * hyperelastic, Ogden
% * displacement logfile
% * stress logfile

clear; close all; clc;

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

%% Control parameters
% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'2021_Nataliya_Perevoshchikova_wrist_project','data','temp');
loadPathMatProp=fullfile(defaultFolder,'2021_Nataliya_Perevoshchikova_wrist_project','data','mat_prop','PCL');

%%
% Building a quadrilateral circular mesh
pointSpacingHeight=0.1;
d=0.35;
h=10;
xSpacing=0.5;
ySpacing=0.4;
nCopies_x=8;
nCopies_y=8;

r=d/2;
ne=2; %Elements in radius
f=0.6; %Fraction (with respect to outer radius) where central square appears
NumStrands=4*7;
%Create the mesh
[Fq,Vq]=discQuadMesh(ne,r,f);

%%
% Lofting to a fiber
Vq(:,3)=0;
Vq1=Vq;
Vq1(:,3)=Vq1(:,3)-h/2;
Vq2=Vq;
Vq2(:,3)=Vq2(:,3)+h/2;

numStepsSweep=ceil(h/pointSpacingHeight);
xc=zeros(numStepsSweep,1);
yc=zeros(numStepsSweep,1);
zc=linspace(-h/2,h/2,numStepsSweep)';
Vc=[xc yc zc];

[~,~,~,S]=sweepLoft(Vq1,Vq2,[0 0 1],[0 0 1],Vc,numStepsSweep,0,0);

%Compose hexahedral elements
X=S.X'; Y=S.Y'; Z=S.Z'; %Coordinate matrices
V=[X(:) Y(:) Z(:)]; %Create node list

I=size(Vq,1)*((1:1:numStepsSweep)-1);
I=I(ones(size(Fq,1),1),:);
I=I(:);

FQ=repmat(Fq,numStepsSweep,1)+I(:,ones(size(Fq,2),1));
E=[FQ(1:end-size(Fq,1),:) FQ(size(Fq,1)+1:end,:)]; %The hexahedral elements
F_start=FQ(1:size(Fq,1),:);
F_end=FQ(end-size(Fq,1)+1:end,:);

[E,V,ind1,ind2]=mergeVertices(E,V); %Merge nodes (start and end are not shared yet)

[F]=element2patch(E); %Element faces
F_start=ind2(F_start); %Start faces for BC's
F_end=ind2(F_end); %End faces for BC's
ind=tesBoundary(F,V); %Indices of boundary faces
Fb=F(ind,:); %Boundary faces (for plotting)

%%
% Visualizing mesh

cFigure; hold on;
title('Hexahedral mesh one fiber','fontSize',fontSize);
gpatch(Fb,V,'gw','k',0.5);
gpatch(F_start,V,'rw','k',1);
gpatch(F_end,V,'bw','k',1);
axisGeom;
camlight headlight;
set(gca,'FontSize',fontSize);
drawnow;

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

C1=57.21894783089972;
C2=57.23071414361289;
K=11444.966197451262;
M1=2.00012552045284;
M2=-1.9997730268725347;
%Initial Viscoelastic parameters
g1_ini=1.25;%Viscoelastic QLV proportional coefficient
t1_ini=25;%Viscoelastic QLV time coefficient
P=[t1_ini g1_ini C1 M1 C2 M2];
displacementMagnitude=0.6030;

% Experimental viscoelastic properties
matname = fullfile(loadPathMatProp,'Viscoelastic_properties.mat');
load(matname);

cFigure;
title('Viscoelastic region','fontSize',fontSize);
hold on;
xlabel('Time, s','FontSize',fontSize);
ylabel('Load rate','FontSize',fontSize);
plot(V_time,lc,'k.-','LineWidth',3);
set(gca,'FontSize',fontSize);
camlight headlight;

[M,I]=max(lc);
time_load=V_time(1:I);
time_unload=V_time(I+1:end);
lc_load=lc(1:I);
lc_unload=lc(I+1:end);

% FEA control settings
max_refs=50; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=6; %Optimum number of iterations
max_retries=10; %Maximum number of retires
k_factor=100;
NumStrands=4*7;

t_load=max(time_load); %Time from start to max load
t_step_ini1=t_load/50; %Initial desired step size
numTimeSteps1=round(t_load/t_step_ini1); %Number of time steps desired
t_step1=t_load/numTimeSteps1; %Step size
dtmin1=t_step1/100; %Smallest allowed step size
dtmax1=t_step1; %Largest allowed step size

t_unload=V_time(end)-max(time_load);
t_step_ini2=t_step_ini1; %Initial desired step size
numTimeSteps2=round(t_unload/t_step_ini2); %Number of time steps desired
t_step2=t_unload/numTimeSteps2; %Step size
dtmin2=t_step2/100; %Smallest allowed step size
dtmax2=t_step2; %Largest allowed step size

cFigure;
title('Viscoelastic region','fontSize',fontSize);
hold on;
xlabel('Time, s','FontSize',fontSize);
ylabel('Load rate','FontSize',fontSize);
plot(time_load,lc_load,'r.-','LineWidth',3);
plot(time_unload,lc_unload,'k.-','LineWidth',3);
set(gca,'FontSize',fontSize);
camlight headlight;

%% DEFINE BC's
%Supported nodes
bcRigidList=unique(F_start(:));

%Prescribed force nodes
bcPrescribeList=unique(F_end(:));
bcPrescribeMagnitudes=displacementMagnitude(ones(1,numel(bcPrescribeList)),:);

%Visualize BC's
cFigure; hold on;
title('Boundary conditions','FontSize',fontSize);
gpatch(F,V,'y','k',1);

hl(1)=plotV(V(bcRigidList,:),'k.','MarkerSize',markerSize);
hl(2)=plotV(V(bcPrescribeList,:),'r.','MarkerSize',markerSize);

legend(hl,{'BC full support','BC prescribed Z displacement'})
axisGeom;
camlight headlight;
set(gca,'FontSize',fontSize);
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';

analysisType='STATIC';
%Control sections for each step
febio_spec.Step.step{1}.Control=febio_spec.Control; %Copy from template
febio_spec.Step.step{1}.ATTR.id=1;
febio_spec.Step.step{1}.Control.analysis=analysisType;
febio_spec.Step.step{1}.Control.time_steps=numTimeSteps1;
febio_spec.Step.step{1}.Control.step_size=t_step1;
febio_spec.Step.step{1}.Control.solver.max_refs=max_refs;
febio_spec.Step.step{1}.Control.solver.max_ups=max_ups;
febio_spec.Step.step{1}.Control.time_stepper.dtmin=dtmin1;
febio_spec.Step.step{1}.Control.time_stepper.dtmax=dtmax1;
febio_spec.Step.step{1}.Control.time_stepper.max_retries=max_retries;
febio_spec.Step.step{1}.Control.time_stepper.opt_iter=opt_iter;

febio_spec.Step.step{2}.Control=febio_spec.Control; %Copy from template
febio_spec.Step.step{2}.ATTR.id=2;
febio_spec.Step.step{2}.Control.analysis=analysisType;
febio_spec.Step.step{2}.Control.time_steps=numTimeSteps2;
febio_spec.Step.step{2}.Control.step_size=t_step2;
febio_spec.Step.step{2}.Control.solver.max_refs=max_refs;
febio_spec.Step.step{2}.Control.solver.max_ups=max_ups;
febio_spec.Step.step{2}.Control.time_stepper.dtmin=dtmin2;
febio_spec.Step.step{2}.Control.time_stepper.dtmax=dtmax2;
febio_spec.Step.step{2}.Control.time_stepper.max_retries=max_retries;
febio_spec.Step.step{2}.Control.time_stepper.opt_iter=opt_iter;

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

materialName1='Material1';
febio_spec.Material.material{1}.ATTR.name=materialName1;

%Viscoelastic part
febio_spec.Material.material{1}.ATTR.type='uncoupled viscoelastic';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.g1=g1_ini;
febio_spec.Material.material{1}.t1=t1_ini;

%Elastic part
febio_spec.Material.material{1}.elastic{1}.ATTR.type='Ogden';
febio_spec.Material.material{1}.elastic{1}.c1=C1;
febio_spec.Material.material{1}.elastic{1}.m1=M1;
febio_spec.Material.material{1}.elastic{1}.c2=C2;
febio_spec.Material.material{1}.elastic{1}.m2=M2;
febio_spec.Material.material{1}.elastic{1}.k=K;

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

% % -> Elements
partName1='Part1';
febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of the element set
febio_spec.Mesh.Elements{1}.ATTR.mat=1; %material index for this set
febio_spec.Mesh.Elements{1}.ATTR.type='hex8'; %Element type of this set
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E;
%

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

% % -> NodeSets
% -> NodeSets
nodeSetName1='bcSupportList';
nodeSetName2='bcPrescribeList';

febio_spec.Mesh.NodeSet{1}.ATTR.name=nodeSetName1;
febio_spec.Mesh.NodeSet{1}.node.ATTR.id=bcRigidList(:);
%
febio_spec.Mesh.NodeSet{2}.ATTR.name=nodeSetName2;
febio_spec.Mesh.NodeSet{2}.node.ATTR.id=bcPrescribeList(:);
%
% %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';

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

febio_spec.Boundary.bc{3}.ATTR.type='fix';
febio_spec.Boundary.bc{3}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{3}.dofs='z';

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

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

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

%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;time_load lc_load;time_unload lc_unload];

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

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

febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_stress;
febio_spec.Output.logfile.element_data{1}.ATTR.data='sx';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.element_data{1}.VAL=1:size(E,1);

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

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

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

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

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

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.disp_log_on=1; %Display convergence information in the command window
febioAnalysis.runMode='internal';%'external';

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

%% Import FEBio results

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

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

    N_disp_mat=N_disp_mat(:,2:end,:);
    sizImport=size(N_disp_mat);
    sizImport(3)=sizImport(3)+1;
    N_disp_mat_n=zeros(sizImport);
    N_disp_mat_n(:,:,2:end)=N_disp_mat;
    N_disp_mat=N_disp_mat_n;
    DN=N_disp_mat(:,:,end);
    DN_magnitude=sqrt(sum(DN(:,1).^2,2)); %X
    V_def=V+DN;
    [CF]=vertexToFaceMeasure(F,DN_magnitude);

    %Importaing force from log file
    [time_mat, Force_mat,~]=importFEBio_logfile(fullfile(savePath,febioLogFileName_force));
    time_mat=[0; time_mat(:)]; %Time
    Force_mat=Force_mat(:,2:end,:);
    sizImport=size(Force_mat);
    sizImport(3)=sizImport(3)+1;
    Force_mat_n=zeros(sizImport);
    Force_mat_n(:,:,2:end)=Force_mat;
    Force_mat=Force_mat_n;
    time_sim=time_mat;
    %%
    % Plotting the simulated results using |anim8| to visualize and animate
    % deformations

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    hp=gpatch(F,V_def,CF,'k',1); %Add graphics object to animate
    gpatch(F,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)]);
    axis([min(V_def(:,1)) max(V_def(:,1)) min(V_def(:,2)) max(V_def(:,2)) min(V_def(:,3)) max(V_def(:,3))]); %Set axis limits statically
    view(130,25); %Set view direction
    camlight headlight;

    % Set up animation features
    animStruct.Time=time_mat; %The time vector
    for qt=1:1:size(N_disp_mat,3) %Loop over time increments
        DN=N_disp_mat(:,:,qt); %Current displacement
        DN_magnitude=sqrt(sum(DN.^2,2)); %Current displacement magnitude
        V_def=V+DN; %Current nodal coordinates
        [CF]=vertexToFaceMeasure(F,DN_magnitude); %Current color data to use

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

    %%
    % Calculate the simulated applied uniaxial stretch
    DZ_set=N_disp_mat(bcPrescribeList,3,:); %Z displacements of the prescribed set
    DZ_set=mean(DZ_set,1); %Calculate mean Z displacements across nodes
    DZ_set=squeeze(DZ_set);

    %Calculated the simulated force
    FZ_set=Force_mat(bcPrescribeList,3,:);
    FZ_set=sum(FZ_set,1); %Calculate mean X force across nodes
    FZ_set=squeeze(FZ_set);

    disp_sim=DZ_set;
    load_sim=FZ_set*NumStrands;

        %%
    % Visualize curves

    hf=cFigure;
    hold on;
    title('Cyclic load-displacement curve','FontSize',fontSize);
    xlabel('Time','FontSize',fontSize); ylabel('Load','FontSize',fontSize);

    hl(1)=plot(time_mat,load_sim,'r.-','lineWidth',lineWidth*1.1,'markerSize',markerSize*1.3);
    hl(2)=plot(V_time,V_load,'k-','lineWidth',lineWidth);

    legend(hl,{'Simulation','Experiment'},'Location','northwest');

    view(2); axis tight;  grid on; axis square; box on;
    set(gca,'FontSize',fontSize);
    drawnow;

end

%% Create structures for optimization
% Material structure
mat_struct.id=1; %Material id

mat_struct.par_names={'t1','g1','c1','m1','c2','m2','k'}; %Parameter names
mat_struct.par_values={t1_ini g1_ini C1 M1 C2 M2 K}; %Parameter values

objectiveStruct.parNormFactors=abs(P); %This will normalize the parameters to ones(size(P))
objectiveStruct.Pb_struct.xx_c=P; %Parameter constraining centre
objectiveStruct.Pb_struct.xxlim=[[P(1)/100 0.01 P(3)/100 -30 P(5)/100 -30]' [P(1)*1e3 50 P(3)*100 30 P(5)*100 30]']; %Parameter bounds

febioAnalysis.time_on=0;
febioAnalysis.time_log_on=0;

%What should be known to the objective function:
objectiveStruct.h=hl(1);
objectiveStruct.bcPrescribeList=bcPrescribeList;
objectiveStruct.time_exp=V_time;
objectiveStruct.load_exp=V_load;
objectiveStruct.febioAnalysis=febioAnalysis;
objectiveStruct.febio_spec=febio_spec;
objectiveStruct.febioFebFileName=febioFebFileName;
objectiveStruct.mat_struct=mat_struct;
objectiveStruct.k_factor=k_factor;
objectiveStruct.NumStrands=NumStrands;

%Optimisation settings
maxNumberIterations=100; %Maximum number of optimization iterations
maxNumberFunctionEvaluations=maxNumberIterations*10; %Maximum number of function evaluations, N.B. multiple evaluations are used per iteration
functionTolerance=1e-25; %Tolerance on objective function value
parameterTolerance=0.001; %Tolerance on parameter variation
displayTypeIterations='iter';

objectiveStruct.method=2;

%File names of output files
output_names.load=fullfile(savePath,febioLogFileName_force);

objectiveStruct.run_output_names=output_names;

%% start optimization

Pn=P./objectiveStruct.parNormFactors;

switch objectiveStruct.method
    case 1 %fminsearch and Nelder-Mead
        OPT_options=optimset('fminsearch'); % 'Nelder-Mead simplex direct search'
        OPT_options = optimset(OPT_options,'MaxFunEvals',maxNumberFunctionEvaluations,...
            'MaxIter',maxNumberIterations,...
            'TolFun',functionTolerance,...
            'TolX',parameterTolerance,...
            'Display',displayTypeIterations,...
            'FinDiffRelStep',1e-2,...
            'DiffMaxChange',0.5);
        [Pn_opt,OPT_out.fval,OPT_out.exitflag,OPT_out.output]= fminsearch(@(Pn) objectiveFunctionIFEA(Pn,objectiveStruct,savePath),Pn,OPT_options);
    case 2 %lsqnonlin and Levenberg-Marquardt
        OPT_options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt');
        OPT_options = optimoptions(OPT_options,'MaxFunEvals',maxNumberFunctionEvaluations,...
            'MaxIter',maxNumberIterations,...
            'TolFun',functionTolerance,...
            'TolX',parameterTolerance,...
            'Display',displayTypeIterations,...
            'FinDiffRelStep',1e-2,...
            'DiffMaxChange',0.5);
        [Pn_opt,OPT_out.resnorm,OPT_out.residual]= lsqnonlin(@(Pn) objectiveFunctionIFEA(Pn,objectiveStruct,savePath),Pn,[],[],OPT_options);
end

%%
[Fopt,OPT_stats_out]=objectiveFunctionIFEA(Pn_opt,objectiveStruct,savePath);

%% Unnormalize and constrain parameters

P_opt=Pn_opt.*objectiveStruct.parNormFactors; %Scale back, undo normalization

%Constraining parameters
for q=1:1:numel(P)
    [P(q)]=boxconstrain(P(q),objectiveStruct.Pb_struct.xxlim(q,1),objectiveStruct.Pb_struct.xxlim(q,2),objectiveStruct.Pb_struct.xx_c(q));
end

disp_text=sprintf('%6.16e,',P_opt); disp_text=disp_text(1:end-1);
disp(['P_opt=',disp_text]);

%%
hf1=cFigure;
%title('displacement load curves optimised','FontSize',fontSize);
xlabel('Time (s)','FontSize',fontSize); ylabel('Force (N)','FontSize',fontSize); zlabel('Z','FontSize',fontSize); hold on;

Hn(1)=plot(OPT_stats_out.time_sim,OPT_stats_out.load_sim,'r.-','lineWidth',lineWidth*1.1,'markerSize',markerSize*1.3);
%Hn(1)=plot(OPT_stats_out.time_sim,OPT_stats_out.load_sim,'r-','lineWidth',lineWidth,'markerSize',markerSize);
%Hn(2)=plot(V_time,V_load,'g-','lineWidth',lineWidth);
Hn(2)=plot(V_time,V_load,'k-','lineWidth',lineWidth);

legend(Hn,{'Simulation','Experiment'},'Location','northwest');
view(2); axis tight;  grid on;
set(gca,'FontSize',fontSize);
drawnow;

function [Fopt,OPT_stats_out]=objectiveFunctionIFEA(Pn,objectiveStruct,savePath)

%%

febioFebFileName=objectiveStruct.febioFebFileName;
febio_spec=objectiveStruct.febio_spec;

%% Unnormalize and constrain parameters

P=Pn.*objectiveStruct.parNormFactors; %Scale back, undo normalization
P_in=P; %Proposed P

%Constraining parameters
for q=1:1:numel(P)
    [P(q)]=boxconstrain(P(q),objectiveStruct.Pb_struct.xxlim(q,1),objectiveStruct.Pb_struct.xxlim(q,2),objectiveStruct.Pb_struct.xx_c(q));
end

%% Setting material parameters

%Acces material parameters

mat_struct=objectiveStruct.mat_struct;
mat_struct.par_values={P(1) P(2) P(3) P(4) P(5) P(6) (P(3)+P(5))*objectiveStruct.k_factor};

disp('SETTING MATERIAL PARAMETERS...');
disp(['Proposed (norm.): ',sprintf(repmat('%6.6e ',[1,numel(Pn)]),Pn)]);
disp(['Proposed        : ',sprintf(repmat('%6.6e ',[1,numel(P_in)]),P_in)]);
disp(['Set (constr.)   : ',sprintf(repmat('%6.6e ',[1,numel(P)]),P)]);

%Assign material parameters
matId=mat_struct.id;
for q=1:1:2
    parNameNow=mat_struct.par_names{q};
    parValuesNow=mat_struct.par_values{q};
    febio_spec.Material.material{matId}.(parNameNow)=mat2strIntDouble(parValuesNow);
end

for q=3:1:numel(mat_struct.par_names)
    parNameNow=mat_struct.par_names{q};
    parValuesNow=mat_struct.par_values{q};
    febio_spec.Material.material{matId}.elastic{1}.(parNameNow)=mat2strIntDouble(parValuesNow);
end

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

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

%% START FEBio

[runFlag]=runMonitorFEBio(objectiveStruct.febioAnalysis);

pause(0.1);

bcPrescribeList=objectiveStruct.bcPrescribeList;
time_exp=objectiveStruct.time_exp;
load_exp=objectiveStruct.load_exp;

if runFlag==1

    %Importing force from log file
    [time_mat, Force_mat,~]=importFEBio_logfile(objectiveStruct.run_output_names.load);
    time_mat=[0; time_mat(:)]; %Time
    Force_mat=Force_mat(:,2:end,:); %Remove index
    sizImport=size(Force_mat);
    sizImport(3)=sizImport(3)+1;
    Force_mat_n=zeros(sizImport);
    Force_mat_n(:,:,2:end)=Force_mat;
    Force_mat=Force_mat_n; %With added zeros for t=0

    %Calculated the simulated force
    FZ_set=Force_mat(bcPrescribeList,3,:);
    FZ_set=sum(FZ_set,1); %Calculate mean X force across nodes
    FZ_set=squeeze(FZ_set);

    time_sim=time_mat;
    load_sim=FZ_set*objectiveStruct.NumStrands; %Multiply force by number of strands

    if ~isempty(objectiveStruct.h)
        objectiveStruct.h.XData=time_sim;
        objectiveStruct.h.YData=load_sim;
        drawnow;
    end

    load_sim_exp = interp1(time_sim,load_sim,time_exp,'pchip');

    %Derive Fopt
    loadDev=load_exp-load_sim_exp;

    switch objectiveStruct.method
        case 1
            Fopt=sum((loadDev).^2)/length(loadDev); %Sum of squared differences
        case 2
            Fopt=sum((loadDev).^2)/length(loadDev);%loadDev(:); %Squared differences
    end

    OPT_stats_out.load_sim=load_sim;
    OPT_stats_out.time_sim=time_sim;
    OPT_stats_out.loadDev=loadDev;
    OPT_stats_out.Fopt=Fopt;

else %Output NaN
    switch objectiveStruct.method
        case 1
            Fopt=NaN;
        case 2
            Fopt=NaN(size(load_exp)); %Squared differences
    end
    OPT_stats_out=[];
end

end

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [] = generateSubsystemLayouts( minerva, cobra_model, color )
% Generates subsystem layouts
%
% USAGE:
%
%    generateSubsystemLayouts( minerva, cobra_model, color )
%
% INPUTS:
%    minerva:           Struct with the information of minerva instance:
%                       address, login, password and model (map)
%    cobra_model:       COBRA model structure
%    color:             Color

    subsystems = getModelSubSystems(cobra_model);

    for i= 1:length(subsystems)

        if nargin < 3
            response = generateSubsytemsLayout(cobra_model, subsystems(i));
        else
            response = generateSubsytemsLayout(cobra_model, subsystems(i), color);
        end

        if ~isempty(regexp(response, strcat('"creator":"', minerva.login, '"')))
            result = [subsystems(i), ' successfully sent to ReconMap.'];
            disp(result)
        else
            result = [subsystems(i), ' failed.'];
            disp(result)
        end

    end
end

author: opencobra

repo: cobratoolbox


text:
function D = generate_2d_3x3_sparse_operator(mat, sz)
% D = generate_first_difference_matrix_for_image(sz, {dif_order, spacing})
% sz -- image size 2x1 or 3x1
% dif_order -- 1 or 2. 1st or 2nd derivative
% spacing -- size of pixel in physical units. 2x1 or 3x1
%
% Creates derivative matrix, if IMG is rect. image or volume, then
% G = D * IMG(:);
% G is [prod(size(IMG))*ndims]x[1] vector, containing drivatives in
% 1, 2, {3} directions.
%
    P = prod(sz);
    s1 = sz(1);
    s2 = sz(2);

    [m1, m2] = ndgrid([1 : s1-1, 1], 1 : s2);
    tp = m1(:) + (m2(:) - 1) * s1;
    tm = (m1(:)+1) + (m2(:) - 1) * s1;
    K = numel(tp);
    A1 = (sparse(1:K, tm, 1, K, P) - sparse(1:K, tp, 1, K, P)) / spacing(1);

    [m1, m2] = ndgrid(1 : s1, [1 : s2-1, 1]);
    tp = m1(:) + (m2(:) - 1) * s1;
    tm = m1(:) + (m2(:)) * s1;
    K = numel(tp);
    A2 = (sparse(1:K, tm, 1, K, P) - sparse(1:K, tp, 1, K, P)) / spacing(2);
    A = [A1; A2];
end

repo: pTVreg


text:
function T = slreadtext(filename)
%SLREADTEXT Reads in a text file into a cell array
%
% $ Syntax $
%   T = slreadtext(filename)
%
% $ Description $
%   - T = slreadtext(filename) Reads in a text file to a cell array with
%     each cell storing a line string.
%
% $ Remarks $
%   1. No extra operations will be performed on the strings besides deleting
%      the trailing spaces in end of lines.
%
% $ History $
%   Created by Dahua Lin, on 2005-06-02
%

fid = fopen(filename, 'r');
if (fid <= 0)
    error(['Fail to open file ', filename]);
end

T = {};
icount = 0;
strline = fgets(fid);
while ~isequal(strline, -1)
    icount = icount + 1;
    T{icount, 1} = deblank(strline);
    strline = fgets(fid);
end

fclose(fid);

author: lmthang


text:
function analyzeAirwayTreeSegments(ptDir,baseName,followupName,segmentsFile)
% function analyzeAirwayTreeSegments(ptDir,doseFile,baseTreeFile,followupTreeFile,segmentsFile)
%
% Example:
% ptDir = 'path/to/patient/dir';
% baseName = 'Baseline';
% followupName = 'Stenosis';
% segmentsFile = 'optional/path/to/segments/file.mat'; % .mat file containing segments
% segmentsFile = ''; % leave empty if no segments need to be loaded
% analyzeAirwayTreeSegments(ptDir,baseName,followupName,segmentsFile)
%
% APA, 6/9/2021

global stateS
%mergedFile = [strtok(baseTreeFile,'_'),'_',strtok(followupTreeFile,'_')];
% mergedFile = [strtok(followupTreeFile,'_')];
mergedFileName = fullfile(ptDir,'merged_files',[followupName,'.mat']);
folllowScanNum = 2;
sliceCallBack('INIT');
sliceCallBack('OPENNEWPLANC',mergedFileName);
sliceCallBack('layout', 3)
hAxis = stateS.handle.CERRAxis(2);
sliceCallBack('selectaxisview', hAxis, 'transverse');
setAxisInfo(hAxis, 'scanSelectMode', 'manual', 'scanSets', folllowScanNum,...
    'xRange',[],'yRange',[]);
setAxisInfo(hAxis, 'structSelectMode', 'manual',...
    'doseSelectMode', 'manual',...
    'structureSets', [],...
    'doseSets', []);
sliceCallBack('planeLocatorToggle')
sliceCallBack('refresh');
showPatientOrientation

%% Plot trees
% radInd = strfind(followupTreeFile,'_radius');
%vfFile = fullfile(ptDir,'registered',[followupTreeFile(1:radInd-1),'_vf.mat']);
% baseRadiusInd = strfind(baseTreeFile,'_radius');
vfFile = fullfile(ptDir,'registered',[baseName,'_',followupName,'_vf.mat']);
%doseFile = fullfile(ptDir,'registered',doseFile);
baseTreeFile = fullfile(ptDir,'AirwayTree',[baseName,'_radius.mat']);
followupTreeFile = fullfile(ptDir,'AirwayTree',[followupName,'_radius.mat']);

% Longitudinal selection
if ~exist('segmentsFile','var')
    segmentsFile = '';
end
select_segments_on_base_and_followup_trees(mergedFileName,baseTreeFile,...
followupTreeFile,vfFile,segmentsFile)

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
% The COBRAToolbox: testVerifyModel.m
%
% Purpose:
%     - Tests the Model Verification function
%

model = getDistributedModel('ecoli_core_model.mat');

% test whether the correct "invalid" subSystem is found
modelSub = model;
modelSub.subSystems(20) = {'blubb'};

res = verifyModel(modelSub,'silentCheck',true);

assert(~isempty(res))
if isfield(res,'Errors')
    assert(isequal(res.Errors.propertiesNotMatched.subSystems, sprintf('Field does not match the required properties at the following positions: \n            20')));
end

% test whether rules are checked crrectly
modelRule = model;
modelRule.rules(3) = {'x(17) )'};
res = verifyModel(modelRule,'silentCheck',true);

assert(~isempty(res))
assert(isequal(res.Errors.propertiesNotMatched.rules, sprintf('Field does not match the required properties at the following positions: \n       3')));


text:
function [ readerobj ] = open_tsx_reader( filename )
%OPEN_TSX_READER Intiates a reader object for TSX XML description file
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

%% Read in XML metadata
xp=javax.xml.xpath.XPathFactory.newInstance.newXPath();
domnode=xmlread(filename); % Read main TSX XML descriptor file
filepath=fileparts(filename);
geoxmlfile=dir(fullfile(filepath,'ANNOTATION','GEOREF.xml'));
if ~isempty(geoxmlfile)
    geo_dom_node = xmlread(fullfile(filepath,'ANNOTATION','GEOREF.xml'));
else
    geo_dom_node = [];
end
meta=meta2sicd_tsxxml(domnode, geo_dom_node);
for i=1:length(meta)
    meta{i}.native.tsxxml=domnode;
end
symmetry=[(meta{1}.SCPCOA.SideOfTrack=='L') 0 1]; % COSAR written in range lines

%% Setup readers for each component
num_files=str2double(xp.evaluate('count(level1Product/productComponents/imageData/file/location)',domnode));
filelist=cell(1,num_files); readerobj=filelist; % Create empty cell arrays
for i=1:num_files
    basepathname=fileparts(filename);
    xpath_str_base=['level1Product/productComponents/imageData[' num2str(i) ']/file/'];
    host=char(xp.evaluate([xpath_str_base 'location/host'],domnode));
    pathname=char(xp.evaluate([xpath_str_base 'location/path'],domnode));
    datafilename=char(xp.evaluate([xpath_str_base 'location/filename'],domnode));
    filelist{i}=fullfile(basepathname,host,pathname,datafilename);
    format = char(xp.evaluate('level1Product/productInfo/imageDataInfo/imageDataFormat',domnode));
    if strcmp(format,'COSAR')
        readerobj{i}=open_cos_reader_noxml(filelist{i},symmetry);
        cos_meta=readerobj{i}.get_meta();
        meta{i}.native.cos=cos_meta.native.cos;
        if (meta{i}.ImageData.NumCols~=cos_meta.native.cos.az)||...
            (meta{i}.ImageData.NumRows~=cos_meta.native.cos.rs)
            error('OPEN_TSX_READER:INCONSISTENT_SIZE','COSAR and XML metadata do not match.');
        end
        % Reading ValidData region from COSAR is slow, so we comment it out
        % meta{i}.ImageData.ValidData = cosar_valid_data(filelist{i},symmetry(1));
        % valid_latlons=point_slant_to_ground(...
        %     [[meta{i}.ImageData.ValidData.Vertex.Row]; ...
        %     [meta{i}.ImageData.ValidData.Vertex.Col]], meta{i});
        % for j = 1:size(valid_latlons,2)
        %     meta{i}.GeoData.ValidData.Vertex(j).Lat=valid_latlons(1,j);
        %     meta{i}.GeoData.ValidData.Vertex(j).Lon=valid_latlons(2,j);
        % end
    elseif strcmp(format,'GEOTIFF')
        readerobj{i}=open_tiff_reader(filelist{i}); % Not necessarily oriented properly
        meta{i}.native.tiff=readerobj{i}.get_meta();
    end
    readerobj{i}.get_meta = @() meta{i};
end
end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function [Jsl] = singleSL(model, cutoff, eliList, atpm)
%% [Jsl] = singleSL(model,cutoff,eliList,atpm)
% INPUT
% model (the following fields are required - others can be supplied)
%   S            Stoichiometric matrix
%   b            Right hand side = dx/dt
%   c            Objective coefficients
%   lb           Lower bounds
%   ub           Upper bounds
%   rxns         Reaction Names
% OPTIONAL
% cutoff         cutoff percentage value for lethality.Default is 0.01.
% eliList        List of reactions to be ignored for lethality
% analysis:Exchange Reactions, ATPM etc.
% atpm           ATPM Reaction Id in model.rxns if other than 'ATPM'
% OUTPUT
% Jsl            Single lethal reactions identified
% Aditya Pratapa       6/26/14.

if exist('cutoff', 'var')
    if isempty(cutoff)
        cutoff = 0.01;
    end
else
    cutoff = 0.01;
end

if exist('atpm', 'var')
    if isempty(atpm)
        atpm = 'ATPM'; % Reaction Id of ATP maintenance reaction- by default it takes 'ATPM'
    end
else
    atpm = 'ATPM';
end

if exist('eliList', 'var')
    if isempty(eliList)
        eliList = model.rxns(ismember(model.rxns, atpm)); % To eliminate ATPM.
    end
else
        eliList = model.rxns(ismember(model.rxns, atpm));
end

Jsl = [];

% Step1 Identify Single Lethal Reactions...
% Identify minNorm flux distribution
solWT = optimizeCbModel(model, 'max', 'one');
grWT = solWT.f;
Jnz = find(~eq(solWT.x, 0));
if (~isempty(eliList))
    eliIdx = find(ismember(model.rxns, eliList)); % Index of reactions not considered for lethality analysis
    Jnz = Jnz(~ismember(Jnz, eliIdx)); % Jnz
end
h = waitbar(0, '0.00', 'Name', 'Identifying Jsl...');

% Identify Single Lethal Reaction Deletions...
modeldel = model;

for iRxn = 1:length(Jnz)
    delIdx_i = Jnz(iRxn);
    modeldel.lb(delIdx_i) = 0;
    modeldel.ub(delIdx_i) = 0;
    solKO_i = optimizeCbModel(modeldel);
    if (solKO_i.f < cutoff*grWT || isnan(solKO_i.f))
        Jsl = [Jsl; delIdx_i];
    end
   % Reset bounds on idx reaction
    modeldel.lb(delIdx_i) = model.lb(delIdx_i);
    modeldel.ub(delIdx_i) = model.ub(delIdx_i);
    waitbar(iRxn / length(Jnz), h, [num2str(round(iRxn * 100 / length(Jnz))) '% completed...']);
end

close(h);
Jsl = model.rxns(Jsl);
end


text:
%SerialLink.animate   Update a robot animation
%
% R.animate(q) updates an existing animation for the robot R.  This will have
% been created using R.plot().
%
% Updates graphical instances of this robot in all figures.
%
% Notes::
% - Not a general purpose method, used for Simulink robot animation.
%
% See also SerialLink.plot.

function animate(robot, q, handles)

    if nargin < 3
        handles = findobj('Tag', robot.name);
    end

    for handle=handles'     % for each graphical robot instance
        h = get(handle, 'UserData');
        animate2( h, q);
    end

function animate2(h, q)

    robot = h.robot;
    n = robot.n;
    L = robot.links;

    mag = h.mag;

    % compute the link transforms, and record the origin of each frame
    % for the animation.

    % points on the multiline, skeleton of the arm
    x = zeros(1,n);
    y = zeros(1,n);
    z = zeros(1,n);

    % points on the multiline, skeleton of the shadow
    xs = zeros(1,n);
    ys = zeros(1,n);
    zs = zeros(1,n);

    % add first point to the link/shadow line, the base

    t = robot.base;

    x(1) = t(1,4);
    y(1) = t(2,4);
    z(1) = t(3,4);
    xs(1) = t(1,4);
    ys(1) = t(2,4);
    zs(1) = h.zmin;

    Tn = t;  % robot base

    if robot.mdh == 0
        % standard DH parameterization

        % Tn(:,:,j) is pose of joint j
        for j=1:n
            Tn(:,:,j) = t;

            t = t * L(j).A(q(j));

            % add point to arm and shadow skeleton
            %   first point is the base
            x(j+1) = t(1,4);
            y(j+1) = t(2,4);
            z(j+1) = t(3,4);
            xs(j+1) = t(1,4);
            ys(j+1) = t(2,4);
            zs(j+1) = h.zmin;
        end
        t = t *robot.tool;
    else
        % modified DH parameterization
        for j=1:n
                        x(j+1) = t(1,4);
            y(j+1) = t(2,4);
            z(j+1) = t(3,4);
            xs(j+1) = t(1,4);
            ys(j+1) = t(2,4);
            zs(j+1) = h.zmin;

            t = t * L(j).A(q(j));
                        Tn(:,:,j) = t;

            % add point to arm and shadow skeleton
            %   first point is the base

        end
        t = t *robot.tool;

    end

    %
    % draw the robot stick figure and the shadow
    %
    set(h.links,'xdata', x, 'ydata', y, 'zdata', z);
    if isfield(h, 'shadow')
        set(h.shadow,'xdata', xs, 'ydata', ys, 'zdata', zs);
    end

    %
    % display the joints as cylinders with rotation axes
    %
    if isfield(h, 'joint')
        xyz_line = [0 0; 0 0; -2*mag 2*mag; 1 1];

        for j=1:n
            % get coordinate data from the cylinder
            xyz = get(h.joint(j), 'UserData');

            % transform to current joint frame
            xyz = Tn(:,:,j) * xyz;

            % convert to 1-vector
            ncols = numcols(xyz)/2;
            xc = reshape(xyz(1,:), 2, ncols);
            yc = reshape(xyz(2,:), 2, ncols);
            zc = reshape(xyz(3,:), 2, ncols);

            % and update the coordinates of the graphical model
            set(h.joint(j), 'Xdata', xc, 'Ydata', yc, ...
                'Zdata', zc);

            if isfield(h, 'jointaxis')
                % optionally draw the joint axis

                % transform line endpoints to current joint frame
                xyzl = Tn(:,:,j) * xyz_line;

                % and update the graphical line
                set(h.jointaxis(j), 'Xdata', xyzl(1,:), ...
                    'Ydata', xyzl(2,:), ...
                    'Zdata', xyzl(3,:));
                set(h.jointlabel(j), 'Position', xyzl(1:3,1));
            end
        end
    end

    %
    % display the wrist axes and labels
    %
    if isfield(h, 'x')
        %
        % compute the wrist axes, based on final link transformation
        % plus the tool transformation.
        %
        xv = t*[mag;0;0;1];
        yv = t*[0;mag;0;1];
        zv = t*[0;0;mag;1];

        %
        % update the line segments, wrist axis and links
        %
        set(h.x,'xdata',[t(1,4) xv(1)], 'ydata', [t(2,4) xv(2)], ...
            'zdata', [t(3,4) xv(3)]);
        set(h.y,'xdata',[t(1,4) yv(1)], 'ydata', [t(2,4) yv(2)], ...
             'zdata', [t(3,4) yv(3)]);
        set(h.z,'xdata',[t(1,4) zv(1)], 'ydata', [t(2,4) zv(2)], ...
             'zdata', [t(3,4) zv(3)]);
        set(h.xt, 'Position', xv(1:3));
        set(h.yt, 'Position', yv(1:3));
        set(h.zt, 'Position', zv(1:3));
    end

author: Allopart

repo: rbpf-gmapping


text:
%PRVERSION PRtools version number
%
%		[VERSION,STR,DATE] = PRVERSION
%
% OUTPUT
%		VERSION   Version number (double)
%		STR				Version number (string)
%		DATE			Version date (string)
%
% DESCRIPTION
% Returns the numerical version number of PRTools VER (e.g. VER = 3.0205)
% and as a string, e.g. STR = '3.2.5'. In DATE, the version date is returned
% as a string.

% $Id: prversion.m,v 1.2 2006/03/08 22:06:58 duin Exp $

function [version,str,date] = prversion(location)

if nargin < 1 % get it locally
	signature = prtver;
	str = signature{1}.Version;
	date = signature{1}.Date;
  version = str2version(str);
	%version   = str2num(str(1)) + (str2num(str(3))*1000 + str2num(str(5))*10)/10000;
	if nargout == 0
		disp([newline '   PRTools version ' str newline])
		clear version
  end
else % location is url
  if verLessThan('matlab','8.0')
    s = ''; % Cannot read if server is down
  else
    [s,status] = urlread(location,'TimeOut',5);
  end
  %s = urlread(location);
  s = [s ' '];                % make sure there is a space
  s = strrep(s,char(10),' '); % replace newline by space
  n = strfind(lower(s),'version');
  nspace = strfind(s,' ');
  str = s(nspace(1)+1:nspace(2)-1);
  version = str2version(str);
  date = s(nspace(2)+1:nspace(3)-1);
end

function version = str2version(str)
n = strfind(str,'.');
if isempty(n)
  version = str2num(str);
elseif length(n) == 1
  version = str2num(str(1:n-1)) + str2num(str(n+1:end))/100;
else
  version = str2num(str(1:n(1)-1)) + str2num(str(n(1)+1:n(2)-1))/100 + str2num(str(n(2)+1:end))/10000;
end
return;

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
% comp_log_Mel receives a time domain speech signal and produces its
% log Mel filterbank coefficients. It is the Matlab counterpart of the
% feature extraction program provided with AURORA2 database.
% Author: Xiao Xiong
% Created: 18 Jul, 2005
% Last Modified: 28 Jul, 2005
% Inputs:
%       x   1-D time domain signal
% Outputs:
%       hist_log_mel    the log Mel filterbank coefficients
%       hist_mel        the Mel filterbank coefficients
%       hist_abs        the spectral magnitude
%       logE            the log energy item needed in computing MFCC

function [X_mel] = abs2Mel(X_abs, Fs, N_mel)

if nargin < 2
    Fs = 8000;
end
if nargin < 3
    N_mel = 23;             % for AURORA2 database, 23 mel filterbank are used
end

[N_vec,n] = size(X_abs);
FFT_length = 2*n;

% generate the mel window
mel_win = mel_window_FE(N_mel, FFT_length/2, Fs);
% X_mel = zeros(N_vec,N_mel);
% for i=1:N_vec
%     X_mel(i,:) = X_abs(i,:)*mel_win;
% end
% faster implementation
X_mel = X_abs*mel_win;

author: singaxiong

repo: SignalGraph


text:
function [mcD mcN smcN] = depthCuesHelper(pc, pcf, rr, sigmaSpace, qzc, nori, sigmaDisparity)
%function [mcD mcN smcN] = depthCuesHelper(pc, pcf, rr, sigmaSpace, qzc, nori, sigmaDisparity)
% Output:
% 	mcD oriented depth gradient based on the distance between the planes at the point being considered
%	mcN is th eangle betweent he 2 planes
%	smcN is the sign of the N gradient
% Input:
%	pc is the point cloud,
%	pcf is the filled in point cloud,
% 	rr determines the offset at which we want to look up the normals to compute the angle between them
%	sigmaSpace is the sigma for the gaussian to estimate the normals
% 	qzc is used to convert the Z value at a pixel into disparity, and to estimate the error at a particular depth.
%	nori is the number of orientations
%	sigmaDisparity is the value sigma for disparity in plane fitting.

	X = pc(:,:,1); Y = pc(:,:,2); Z = pc(:,:,3);
	Xf = pcf(:,:,1); Yf = pcf(:,:,2); Zf = pcf(:,:,3);

	[h w] = size(Z);
	mcN = zeros([size(Z) nori]);
	mcD = zeros([size(Z) nori]);

	thetaQ = -([0:1:(nori-1)]./nori*pi - pi/2);

	N{1} = NaN([size(Z) 3]);
	N{2} = NaN([size(Z) 3]);
	D{1} = NaN([size(Z) 3]);
	D{2} = NaN([size(Z) 3]);
	cntr{1} = NaN([h w 3]);
	cntr{2} = NaN([h w 3]);
	xyz = cat(3, Xf, Yf, Zf);

	% qZ1 = qzc*ordfilt2(Zf,1,ones(2*rr+1)).^2;
	se = strel(ones(2*rr+1));
	qZ2 = qzc*((-imdilate(-Zf, se)).^2);
	% assert(isequal(qZ1(rr+1:end-rr, rr+1:end-rr), qZ2(rr+1:end-rr, rr+1:end-rr)), 'Oops something is wrong with dilation as a substitue for max filter.');
	qZ = qZ2;

	for k = 1:nori,

		% For computing the signal at different orientations
		[pos(:,:,1) pos(:,:,2)] = ndgrid(1:size(Z,1), 1:size(Z,2));
		[theta, r] = cart2pol(pos(:,:,1), pos(:,:,2));
		theta = theta-thetaQ(k);
		[pos(:,:,1) pos(:,:,2)] = pol2cart(theta, r);

		% Point offsets at which we want to compare the features at
		anchorI = [-(rr+1)/2, (rr+1)/2];
		anchorJ = [0 0];
		[theta, r] = cart2pol(anchorI, anchorJ);
		theta = theta+thetaQ(k);
		[anchorI, anchorJ] = pol2cart(theta, r);
		Rpad = ceil(max(abs([anchorI, anchorJ])))+1;

		pos(:,:,2) = pos(:,:,2)/2; %sqrt(2);

		one_Z = 1./Z;
		X_Z = X./Z;
		Y_Z = Y./Z;
		one = Z; one(~isnan(one)) = 1;
		X_ZZ = X./(Z.*Z);
		Y_ZZ = Y./(Z.*Z);

		AtARaw = cat(3, X_Z.^2, X_Z.*Y_Z, X_Z, Y_Z.^2, Y_Z, one);
		AtbRaw = cat(3, X_ZZ, Y_ZZ, one_Z);
		AtARaw(isnan(AtARaw)) = 0;
		AtbRaw(isnan(AtbRaw)) = 0;

		AtA = filterItChopOffIS(cat(3, AtARaw, AtbRaw), pos, Zf, qzc, sigmaDisparity, sigmaSpace);
		Atb = AtA(:, :, (size(AtARaw,3)+1):end);
		AtA = AtA(:, :, 1:size(AtARaw,3));

		[AtA_1 detAtA] = invertIt(AtA);

		Nboth = mutiplyIt(AtA_1, Atb);
		filteredOne = AtA(:,:,end);

		% to ignore places with very few points.
		badPts = (filteredOne./max(filteredOne(:))) < 0.25;
		Nboth(repmat(badPts, [1 1 3])) = NaN;

		bboth = Nboth(:,:,1);
		bboth(:) = -detAtA;
		bboth = bsxfun(@rdivide, bboth, sqrt(sum(Nboth.^2,3)));
		Nboth = bsxfun(@rdivide, Nboth, sqrt(sum(Nboth.^2,3)));
		%Dboth  = xyz - bsxfun(@times, (sum(Nboth.*xyz, 3) + b), Nboth);

		Nboth = padarray(Nboth, [Rpad Rpad 0], 'replicate', 'both');
		bboth = padarray(bboth, [Rpad Rpad 0], 'replicate', 'both');
		XYZ = padarray(xyz, [Rpad Rpad 0], 'replicate', 'both');

		stI = round(Rpad+1+anchorI(1));
		enI = stI+(h-1);
		stJ = round(Rpad+1+anchorJ(1));
		enJ = stJ+(w-1);
		N{1} = Nboth(stI:enI, stJ:enJ, :);
		b{1} = bboth(stI:enI, stJ:enJ, :);
		D{1} = xyz - bsxfun(@times, (sum(N{1}.*xyz, 3) + b{1}), N{1});
		cntr{1} = XYZ(stI:enI, stJ:enJ, :);

		stI = round(Rpad+1+anchorI(2));
		enI = stI+(h-1);
		stJ = round(Rpad+1+anchorJ(2));
		enJ = stJ+(w-1);
		N{2} = Nboth(stI:enI, stJ:enJ, :);
		b{2} = bboth(stI:enI, stJ:enJ, :);
		D{2} = xyz - bsxfun(@times, (sum(N{2}.*xyz, 3) + b{2}), N{2});
		cntr{2} = XYZ(stI:enI, stJ:enJ, :);

		mcN(:,:,k) = 1-abs(sum(N{1}.*N{2}, 3));
		mcD(:,:,k) = sqrt(sum((D{1} - D{2}).^2, 3));
		fprintf('.');

		% Making the normals consistent.
		for i = 1:2,
			N{i} = N{i}.*repmat(sign(N{i}(:,:,3)),[1 1 3]);
			sn = sign(sum(N{i}.*cat(3, Xf, Yf, Zf),3));
			sn(isnan(sn)) = 1;
			N{i} = repmat(sn,[1 1 3]).*N{i};
		end

		% Computing the orientation for the normal discontinuity
		smcN(:,:,k) = sign(sum((N{1}-N{2}).*(cntr{1}-cntr{2}),3));
		mcD(:,:,k) = mcD(:,:,k).*(mcD(:,:,k) > qZ*1.05);

		fprintf('%d ',k);
	end

	fprintf('\n');
	mcN(isnan(mcN)) = 0;
	mcD(isnan(mcD)) = 0;
end

function fFilt = filterItChopOffIS(f, pos, Zf, qzc, sigmaDisparity, r)
	sp(:,:,1:2) = pos;
	sp(:,:,1:2) = sp(:,:,1:2)/r;
	sp(:,:,4) = (1./(qzc*Zf))/sigmaDisparity;
	fFilt = jointBilateral(sp, f, 1, 1000);
end

function x = mutiplyIt(AtA_1, Atb)
	a = @(k) AtA_1(:,:,k);
	b = @(k) Atb(:,:,k);
	x1 = a(1).*b(1) + a(2).*b(2) + a(3).*b(3);
	x2 = a(2).*b(1) + a(4).*b(2) + a(5).*b(3);
	x3 = a(3).*b(1) + a(5).*b(2) + a(6).*b(3);
	x = cat(3, x1, x2, x3);
end

function [AtA_1 detAtA]= invertIt(AtA)
	a = @(k) AtA(:,:,k);

	A = a(4).*a(6) - a(5).*a(5);
	D = -(a(2).*a(6 )- a(3).*a(5));
	G = a(2).*a(5) - a(3).*a(4);
	E = a(1).*a(6) - a(3).*a(3);
	H = -(a(1).*a(5) - a(2).*a(3));
	K = a(1).*a(4) - a(2).*a(2);

	detAtA = a(1).*A + a(2).*D + a(3).*G;

	AtA_1 = cat(3, A, D, G, E, H, K);
	%AtA_1 = bsxfun(@rdivide, AtA_1, detAtA);
end

author: s-gupta

repo: rgbd

save_path: github-repos/MATLAB/s-gupta-rgbd


text:
function varargout = lvd_EditEnginesGUI(varargin)
% LVD_EDITENGINESGUI MATLAB code for lvd_EditEnginesGUI.fig
%      LVD_EDITENGINESGUI, by itself, creates a new LVD_EDITENGINESGUI or raises the existing
%      singleton*.
%
%      H = LVD_EDITENGINESGUI returns the handle to a new LVD_EDITENGINESGUI or the handle to
%      the existing singleton*.
%
%      LVD_EDITENGINESGUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in LVD_EDITENGINESGUI.M with the given input arguments.
%
%      LVD_EDITENGINESGUI('Property','Value',...) creates a new LVD_EDITENGINESGUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before lvd_EditEnginesGUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to lvd_EditEnginesGUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help lvd_EditEnginesGUI

% Last Modified by GUIDE v2.5 03-Dec-2018 17:11:04

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @lvd_EditEnginesGUI_OpeningFcn, ...
                   'gui_OutputFcn',  @lvd_EditEnginesGUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before lvd_EditEnginesGUI is made visible.
function lvd_EditEnginesGUI_OpeningFcn(hObject, eventdata, handles, varargin)
    % This function has no output args, see OutputFcn.
    % hObject    handle to figure
    % eventdata  reserved - to be defined in a future version of MATLAB
    % handles    structure with handles and user data (see GUIDATA)
    % varargin   command line arguments to lvd_EditEnginesGUI (see VARARGIN)

    % Choose default command line output for lvd_EditEnginesGUI
    handles.output = hObject;

    lvdData = varargin{1};
    setappdata(hObject,'lvdData',lvdData);

    populateGUI(handles, lvdData);

    % Update handles structure
    guidata(hObject, handles);

    % UIWAIT makes lvd_EditEnginesGUI wait for user response (see UIRESUME)
    uiwait(handles.lvd_EditEnginesGUI);

function populateGUI(handles, lvdData)
    set(handles.enginesListBox,'String',lvdData.launchVehicle.getEnginesListBoxStr());

    numEngines = length(get(handles.enginesListBox,'String'));
    if(numEngines <= 1)
        handles.removeEngineButton.Enable = 'off';
    end

% --- Outputs from this function are returned to the command line.
function varargout = lvd_EditEnginesGUI_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
    if(isempty(handles))
        varargout{1} = [];
    else
        close(handles.lvd_EditEnginesGUI);
    end

% --- Executes on button press in closeButton.
function closeButton_Callback(hObject, eventdata, handles)
% hObject    handle to closeButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    uiresume(handles.lvd_EditEnginesGUI);

% --- Executes on selection change in enginesListBox.
function enginesListBox_Callback(hObject, eventdata, handles)
% hObject    handle to enginesListBox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns enginesListBox contents as cell array
%        contents{get(hObject,'Value')} returns selected item from enginesListBox
    if(strcmpi(get(handles.lvd_EditEnginesGUI,'SelectionType'),'open'))
        lvdData = getappdata(handles.lvd_EditEnginesGUI,'lvdData');
        lv = lvdData.launchVehicle;

        selEngine = get(handles.enginesListBox,'Value');
        engine = lv.getEngineForInd(selEngine);
        stageState = getStageStateForEngine(engine, lvdData);
        oldEngineState = stageState.getStateForEngine(engine);
        stageState.removeEngineStateForEngine(engine);

        lvd_EditEngineGUI(engine);

        stageState = getStageStateForEngine(engine, lvdData);
        newEngineState = LaunchVehicleEngineState(stageState);
        newEngineState.engine = engine;

        if(not(isempty(oldEngineState)))
            newEngineState.active = oldEngineState.active;
        end
        stageState.addEngineState(newEngineState);

        set(handles.enginesListBox,'String',lvdData.launchVehicle.getEnginesListBoxStr());
    end

function stageState = getStageStateForEngine(engine, lvdData)
        stage = engine.stage;
        stageStates = lvdData.initStateModel.stageStates;
        stageStateInd = find([stageStates.stage] == stage,1,'first');
        stageState = stageStates(stageStateInd);

% --- Executes during object creation, after setting all properties.
function enginesListBox_CreateFcn(hObject, eventdata, handles)
% hObject    handle to enginesListBox (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: listbox controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes on button press in addEngineButton.
function addEngineButton_Callback(hObject, eventdata, handles)
% hObject    handle to addEngineButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_EditEnginesGUI,'lvdData');

    initStage = lvdData.launchVehicle.getStageForInd(1);

    engine = LaunchVehicleEngine(initStage);
    useEngine = lvd_EditEngineGUI(engine);

    if(useEngine)
%         engine.stage.addEngine(engine);

        stageStates = lvdData.initStateModel.stageStates;
        stageStateInd = find([stageStates.stage] == engine.stage,1,'first');
        stageState = stageStates(stageStateInd);

        newEngineState = LaunchVehicleEngineState(stageState);
        newEngineState.engine = engine;
        stageState.addEngineState(newEngineState);

        set(handles.enginesListBox,'String',lvdData.launchVehicle.getEnginesListBoxStr());

        handles.removeEngineButton.Enable = 'on';
    end

% --- Executes on button press in removeEngineButton.
function removeEngineButton_Callback(hObject, eventdata, handles)
% hObject    handle to removeEngineButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_EditEnginesGUI,'lvdData');
    lv = lvdData.launchVehicle;

    selEngine = get(handles.enginesListBox,'Value');
    engine = lv.getEngineForInd(selEngine);

    tf = engine.isInUse();

    if(tf == false)
        stage = engine.stage;

        lv.removeAllEngineToTanksConnsWithEngine(engine);

        stageStates = lvdData.initStateModel.stageStates;
        stageStateInd = find([stageStates.stage] == stage,1,'first');
        stageState = stageStates(stageStateInd);

        stageState.removeEngineStateForEngine(engine);

        stage.removeEngine(engine);

        listBoxStr = lvdData.launchVehicle.getEnginesListBoxStr();
        set(handles.enginesListBox,'String',listBoxStr);

        numEngines = length(listBoxStr);
        if(selEngine > numEngines)
            set(handles.enginesListBox,'Value',numEngines);
        end

        if(numEngines <= 1)
            handles.removeEngineButton.Enable = 'off';
        end
    else
        warndlg(sprintf('Could not delete the engine "%s" because it is in use as part of an event termination condition, event action, objective function, or constraint.  Remove the engine dependencies before attempting to delete the engine.', engine.name),'Cannot Delete Engine','modal');
    end

% --------------------------------------------------------------------
function copyEngineMenu_Callback(hObject, eventdata, handles)
% hObject    handle to copyEngineMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
    lvdData = getappdata(handles.lvd_EditEnginesGUI,'lvdData');
    lv = lvdData.launchVehicle;

    selEngine = get(handles.enginesListBox,'Value');
    engine = lv.getEngineForInd(selEngine);

    newEngine = engine.copy();
    initStage = newEngine.stage;

	newEngine.stage.addEngine(newEngine);

    stageStates = lvdData.initStateModel.stageStates;
    stageStateInd = find([stageStates.stage] == initStage,1,'first');
    stageState = stageStates(stageStateInd);

    newEngineState = LaunchVehicleEngineState(stageState);
    newEngineState.engine = newEngine;
    stageState.addEngineState(newEngineState);

    set(handles.enginesListBox,'String',lvdData.launchVehicle.getEnginesListBoxStr());

    handles.removeEngineButton.Enable = 'on';

% --------------------------------------------------------------------
function enginesListboxContextMenu_Callback(hObject, eventdata, handles)
% hObject    handle to enginesListboxContextMenu (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% --- Executes on key press with focus on lvd_EditEnginesGUI or any of its controls.
function lvd_EditEnginesGUI_WindowKeyPressFcn(hObject, eventdata, handles)
% hObject    handle to lvd_EditEnginesGUI (see GCBO)
% eventdata  structure with the following fields (see MATLAB.UI.FIGURE)
%	Key: name of the key that was pressed, in lower case
%	Character: character interpretation of the key(s) that was pressed
%	Modifier: name(s) of the modifier key(s) (i.e., control, shift) pressed
% handles    structure with handles and user data (see GUIDATA)
    switch(eventdata.Key)
        case 'return'
            uiresume(handles.lvd_EditEnginesGUI);
        case 'enter'
            uiresume(handles.lvd_EditEnginesGUI);
        case 'escape'
            uiresume(handles.lvd_EditEnginesGUI);
    end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [cc,c0]=lpcpf2cc(pf,np,f)
%LPCPF2CC Convert power spectrum to complex cepstrum CC=(PF,NP)
%
%  Inputs: pf(nf,n)    Power spectrum, uniformly spaced DC to Nyquist
%          np          Number of cepstral coefficients to calculate [n-1]
%          f(1,n)      Frequencies of pf columns [linspace(0,0.5,n)]
%
% Outputs: cc(nf,np)  Complex spectrum from DC to Nyquist
%          c0(nf,1)   The zeroth cepstral coefficient, c(0)
%
% Note since the log spectrum is not normally bandlimited, this conversion is not exact unless n >> np

%      Copyright (C) Mike Brookes 1998-2014
%      Version: $Id: lpcpf2cc.m 5025 2014-08-22 17:07:24Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 2 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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nf,nq]=size(pf);
if nargin<2 || ~numel(np) np=nq-1; end
if nargin<3 || ~numel(f)
    cc=rsfft(log(pf.')).'/(2*nq-2);
    c0=cc(:,1)*0.5;
    cc(:,nq)=cc(:,nq)*0.5;
    if np>nq-1
        cc=[cc(:,2:nq) zeros(nf,np-nq+1)];
    else
        cc=cc(:,2:np+1);
    end
else
    cc=0.5*(log(pf)/cos(2*pi*f(:)*(0:min(np,nq-1))));
    c0=cc(:,1);
    cc=cc(:,2:end);
    if np>nq-1
        cc(1,np)=0;
    end
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function [tmpS] = dicomrt_d2c_buildscanstr(scanInitS)

% Construct the scan structure
tmpS = [];

tmpS.scanArray            = '';
tmpS.scanType             = 'CT';
tmpS.scanInfo             = scanInitS;
tmpS.uniformScanInfo      = '';
tmpS.scanArraySuperior    = '';
tmpS.scanArrayInferior    = '';


text:
function medit_io_test03 ( filename )

%*****************************************************************************80
%
%% MEDIT_IO_TEST03 reads and prints the sizes in a MESH dataset.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 October 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'MEDIT_IO_TEST03\n' );
  fprintf ( 1, '  Read a mesh file and print its sizes.\n' );

  [ dim, vertices, edges, triangles, quadrilaterals, tetrahedrons, ...
    hexahedrons ] = mesh_size_read ( filename );
%
%  Print information.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Header information for "%s":\n', filename );

  mesh_size_print ( dim, vertices, edges, triangles, quadrilaterals, ...
    tetrahedrons, hexahedrons );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function c = set(c, prop_name, val)

   %SET Set properties for correlation object
   %   c = SET(c,prop_name,val)
   %
   %   see help correlation/get.m for valid property names

   % Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   switch upper(prop_name)
      case {'WAVEFORMS' 'WAVEFORM' 'WAVES'}
         c.W = reshape(val,length(val),1); % make 1 column
      case {'TRIG'}
         c.trig = reshape(val,length(val),1);
      case {'CORR'}
         c.corrmatrix = val;
      case {'LAG'}
         c.lags = val;
      case {'STAT'}
         c.stat = val;
      case {'LINK'}
         c.link = val;
      case {'CLUST'}
         c.clust = val;
      otherwise
         warning('can''t understand property name');
         help correlation/set
   end;
end

author: geoscience-community-codes

repo: GISMO


text:
%kkgenloc 'Generate Location Data'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kgenloc.pane file
%
% Parameters:
% InputFile: i 'Input ', optional: 'First Input data object'
% MultiChoice: grid 'Grid Type :', default: 1: 'Grid Type to use when creating Location Data'
%    Choices are:
%   1: 'Uniform'
%   2: 'Rectilinear'
%   3: 'Curvilinear'
% Integer: wsize ' Width ', default: 512: 'Width of location data'
% Integer: hsize ' Height ', default: 512: 'Height of location data'
% Integer: dsize ' Depth ', default: 1: 'Depth of location data'
% OutputFile: o 'Output', required: 'Resulting output data object'
%
% Example: o = kkgenloc(i, {'i','';'grid',1;'wsize',512;'hsize',512;'dsize',1;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kgenloc - Generate Location Data
%
%  DESCRIPTION
% "kgenloc" creates explicit location data in the specified
% destination object (-o). If an input source object is provided (-i)
% then any data segments in that source object will be copied to the
% destination object.
%
% The explicit location data is specified by three sets of begin and end
% values (-xbegin, -xend), (-ybegin, -yend), (-zbegin, -zend) which
% together compose explicit corner markers for volume in (x,y,z) three
% space.  The x-axis corresponds uniformly to the implicit width
% dimension on the polymorphic data model, with the y-axis and z-axis
% corresponding to implicit height, and depth dimensions, respectively.
%
% Data is generated uniformly from those specified corner markers,
% however the explicitness of the generated data can be specified by the
% grid (-grid) option.  A KUNIFORM grid will dictate that the generated
% data consist simply of setting the uniform corner markers on the
% output data set.  A KRECTILINEAR grid will dictate that the generated
% data consist of three arrays of locations, an x-array, a y-array, and
% a z-array, with the data within those arrays changing uniformly over
% the specified (x,y,z) begin and end values.  A KCURVILINEAR grid will
% dictate that the generated data consist of a four dimensional array of
% location data.  Essentially, this data set is a three dimensional
% array of (x,y,z) points.
%
% The size and datatype of the resulting location data segment can be
% specified.  The size of the location segment is controlled with three
% parameters (-wsize, -hsize, -dsize) which correspond to the implicit
% width, height, and depth of the location segment.  The dimension size
% of the destination location segment will always be three corresponding
% to the size of the location points (x,y,z).  The data type of the
% location segment is controlled with a single parameter (-type).
%
% Any size specifications will be overwritten by existing size
% information in the input source object.  If the input source object
% contains any explicit location data, the size and data type of the
% existing location segment will be used to determine the size and data
% type of the destination location segment.  Any explicit location data
% in the source object will be destroyed in favor of the newly generated
% data.  Additionally, if the input source object contains no explicit
% location segment but does contain a value segment or a mask segment,
% the width, height, and depth size of either of those segments will be
% used to set the width, height, and depth size of the destination
% location segment.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% ksetdattr
%
%  RESTRICTIONS
% The data is generated and stored as one block, so generating large
% location data sets are subject to the memory limitations of your
% machine.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkgenloc(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,..] = kkgenloc(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'grid', 1;'wsize', 512;'hsize', 512;'dsize', 1;'o', '__output'};
maxval={1,0,2,2,2,0};
minval={1,0,2,2,2,0};
istoggle=[1,0,1,1,1,0];
was_set=istoggle * 0;
paramtype={'InputFile','MultiChoice','Integer','Integer','Integer','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 'kgenloc"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [k, nErrs] = generateErrorOddsConstant(oddsDesired, nTrials)
if nargin == 0
    oddsDesired = 10;
    nTrials = 1000;
end
if nargin == 1
    nTrials = 1000;
end

% We want to calculate k that is the number to multiply rand() by to
% achieve desired odds. This seems to be a good fomula for
% approximating the odds for all of the various oddsDesired.
if oddsDesired == 0
    c = 0;
elseif oddsDesired < 5
    c = .02;
elseif oddsDesired > 50
    c = .3;
else
    c = .24/oddsDesired;
end
k = oddsDesired/100 + .43 + c;
nErrs = 0;
for ii = 1:nTrials
    r = round(k*rand());
    if r==1
        %fprintf('%d. Error!!\n', ii);
        nErrs = nErrs+1;
    else
        %fprintf('%d\n', ii);
    end
    % pause(.01);
end
fprintf('\nTotal number of errors: %d  ==>  %0.1f%%\n\n', nErrs, 100 * nErrs / nTrials);

author: BUNPC


text:
function reducedPath = tapas_uniqc_genpath_exclude_subdirs(completePath, pathsExcluded)
% works as genpath, but allows to specify excluded subfolders, which are...
% removed incl. their subfolders from the returned path
%
%    reducedPath = tapas_uniqc_genpath_exclude_subdirs(completePath, pathsExcluded)
%
% IN
%   completePath    the input path, as in genpath
%   pathsExcluded   cell(1, nPaths) of excluded sub-folders or paths
%
% OUT
%   reducedPath     string of returned genpath (e.g., use for addpath),
%                   excluding the specified pathsExcluded and their subfolders.
%
% EXAMPLE
%   tapas_uniqc_genpath_exclude_subdirs
%
%   See also

% Author:   Lars Kasper
% Created:  2020-05-22
% Copyright (C) 2020 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public License (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.
%

%create cell of folder names from path string, split at path separator ; (win) or : (unix)
recursivePath = genpath(completePath);
reducedPath = regexp(recursivePath, pathsep, 'split');
for iExcluded = 1:numel(pathsExcluded)
    % note: regexprep might have some issues w/ escape characters, like \., but helps with Windows \ vs Unix /
    idxIncludedFolders = find(cellfun('isempty', regexp(...
        regexprep(reducedPath, '\\', '/'), ...
        regexprep(pathsExcluded{iExcluded}, '\\', '/')  ) ));
    reducedPath = reducedPath(idxIncludedFolders);
end
% add pathseparator to each path, convert back into path string
reducedPath = cellfun(@(x) [x pathsep], reducedPath, 'UniformOutput', false);
reducedPath = [reducedPath{:}];

author: translationalneuromodeling

repo: tapas


text:
function ui = mrViewViewMenu(ui);
%
% ui = mrViewViewMenu(ui);
%
% Attaches a view menu to a mrViewer UI.
%
% ras 07/05.
if ~exist('ui', 'var') | isempty(ui),  ui = get(gcf, 'UserData'); end

ui.menus.view = uimenu(ui.fig, 'Label', 'View');

chkVals = {'off' 'on'}; % will be used to set some menus,  below

% navigation panel toggle callback:
cb = 'ui = mrViewGet; ';
cb = [cb 'mrvPanelToggle(ui.panels.nav, gcbo); clear ui; '];
ui.menus.navToggle = uimenu(ui.menus.view, 'Label', 'Navigation Panel', ...
                            'Checked', 'on', 'Separator', 'off', ...
                            'Accelerator', '1', 'Callback', cb);

% ROI panel toggle callback:
cb = 'ui = mrViewGet; ';
cb = [cb 'mrvPanelToggle(ui.panels.roi, gcbo); clear ui; '];
ui.menus.infoToggle = uimenu(ui.menus.view, 'Label', 'ROI Panel', ...
                            'Checked', 'off', 'Separator', 'off', ...
                            'Accelerator', '2', 'Callback', cb);

% grayscale panel toggle callback:
cb = 'ui = mrViewGet; ';
cb = [cb 'mrvPanelToggle(ui.panels.grayscale, gcbo); clear ui; '];
ui.menus.grayscaleToggle = uimenu(ui.menus.view, 'Label', 'Grayscale Panel', ...
                            'Checked', 'off', 'Separator', 'off', ...
                            'Accelerator', '3', 'Callback', cb);

% info panel toggle callback:
cb = 'ui = mrViewGet; ';
cb = [cb 'mrvPanelToggle(ui.panels.info, gcbo); clear ui; '];
ui.menus.infoToggle = uimenu(ui.menus.view, 'Label', 'Info Panel', ...
                            'Checked', 'off', 'Separator', 'off', ...
                            'Accelerator', '4', 'Callback', cb);

% mesh panel toggle callback:
cb = 'ui = mrViewGet; ';
cb = [cb 'mrvPanelToggle(ui.panels.mesh, gcbo); clear ui; '];
ui.menus.infoToggle = uimenu(ui.menus.view, 'Label', 'Mesh Panel', ...
                            'Checked', 'off', 'Separator', 'off', ...
                            'Accelerator', '5', 'Callback', cb);

%%%%% Dock submenu
ui.menus.dock = uimenu(ui.menus.view, 'Label', 'Dock Panel to Display', ...
                       'Separator', 'on');
% for these submenus, we need to know which panels are currently
% attached to the main figure, and which are not. This subroutine does
% that...
chk = findDockedPanels(ui);

% dock Navigation Panel
cb = 'mrViewDock(gcf, ''nav'', gcbo); ';
ui.menus.navDock = uimenu(ui.menus.dock, 'Label', 'Navigation Panel', ...
                            'Checked', chk{1}, 'Callback', cb);

% dock ROI Panel
cb = 'mrViewDock(gcf, ''roi'', gcbo); ';
ui.menus.roiDock = uimenu(ui.menus.dock, 'Label', 'ROI Panel', ...
                            'Checked', chk{2}, 'Callback', cb);

% dock Grayscale Panel
cb = 'mrViewDock(gcf, ''grayscale'', gcbo); ';
ui.menus.navDock = uimenu(ui.menus.dock, 'Label', 'Grayscale Panel', ...
                            'Checked', chk{3}, 'Callback', cb);

% dock Info Panel
cb = 'mrViewDock(gcf, ''info'', gcbo); ';
ui.menus.infoDock = uimenu(ui.menus.dock, 'Label', 'Info Panel', ...
                            'Checked', chk{4}, 'Callback', cb);

% dock Mesh Panel
cb = 'mrViewDock(gcf, ''mesh'', gcbo); ';
ui.menus.meshDock = uimenu(ui.menus.dock, 'Label', 'Mesh Panel', ...
                            'Checked', chk{5}, 'Callback', cb);

% dock Overlay Panels
cb = 'mrViewDock(gcf, ''overlay'', gcbo); ';
ui.menus.overlayDock = uimenu(ui.menus.dock, 'Label', 'Overlay Panels', ...
                            'Checked', 'on', 'Callback', cb);

% dock Colorbar Panel
cb = 'mrViewDock(gcf, ''colorbar'', gcbo); ';
ui.menus.colorbarDock = uimenu(ui.menus.dock, 'Label', 'Colorbar Panel', ...
                            'Checked', chk{6}, 'Callback', cb);

%%%%% overlays submenu
ui.menus.overlays = uimenu(ui.menus.view, 'Label', 'Overlays', ...
                           'Separator', 'on');

% new overlay
uimenu(ui.menus.overlays, 'Label', 'New Overlay', 'Separator', 'on', ...
       'Accelerator', 'O', 'Callback', 'mrViewAddOverlay;');

% close overlay
uimenu(ui.menus.overlays, 'Label', 'Close Overlay', 'Separator', 'on', ...
       'Accelerator', '9', 'Callback', 'mrViewRemoveOverlay;');

% text for showing current overlays:
uimenu(ui.menus.overlays, 'Label', 'Show Overlay Panels:', ...
                         'ForegroundColor', 'w', 'Separator', 'on');

%%%%% label submenu
ui.menus.label = uimenu(ui.menus.view, 'Label', 'Label', 'Separator', 'on');

% label directions
cb = 'val=umtoggle(gcbo); mrViewSet([], ''labelDirs'', val); mrViewRefresh;';
chk = chkVals{ui.settings.labelDirs+1};
uimenu(ui.menus.label, 'Label', 'Label Directions', ...
       'Checked', chk, 'Separator', 'on', 'Callback', cb);

% label units menu:
cb = 'val=umtoggle(gcbo); mrViewSet([], ''labelAxes'', val); mrViewRefresh;';
chk = chkVals{ui.settings.labelAxes+1};
uimenu(ui.menus.label, 'Label', 'Label Units', ...
       'Checked', chk, 'Separator', 'on', 'Callback', cb);

% preserve aspect menu:
cb = 'val=umtoggle(gcbo); mrViewSet([], ''eqAspect'', val); mrViewRefresh;';
ui.menus.preserveAspect = uimenu(ui.menus.view, 'Label', 'Preserve Aspect', ...
                               'Checked', 'on', 'Separator', 'on', ...
                               'Callback', cb);

% All slices option
cb = ['nslices = mrViewGet(gcf, ''NumSlices''); ' ...
      'nrows = ceil(sqrt(nslices)); ncols = ceil(nslices/nrows); ' ...
      'mrViewSet(gcf, ''MontageRows'', nrows, ''MontageCols'', ncols); ' ...
      'mrViewSet(gcf, ''Slice'', 1, ''DisplayFormat'', 1); ' ...
	  'mrViewRefresh(gcf); ' ...
      'clear nrows ncols nslices; '];
uimenu(ui.menus.view, 'Label', 'All Slices', ...
       'Separator', 'off', 'Callback', cb);

% center of current ROI option
cb = 'mrViewRecenter(get(gcf, ''UserData''), ''roi''); ';
uimenu(ui.menus.view, 'Label', 'Center of Current ROI', ...
       'Separator', 'off', 'Callback', cb);

% toggle for standard MATLAB figure menus
addFigMenuToggle(ui.menus.view);

% refresh options
uimenu(ui.menus.view,  'Label', 'Clean Up Display Panels',  ...
        'Separator',  'on',  'Callback', 'mrViewFixDisplay(gcf);');

uimenu(ui.menus.view,  'Label', 'Refresh',  'Separator',  'off',  ...
       'Accelerator', 'Y',  'Callback', 'mrViewRefresh;');

return
% /--------------------------------------------------------------/ %

% /--------------------------------------------------------------/ %
function chk = findDockedPanels(ui);
% Checkes whether each panel in a mrViewer UI is attached ("docked")
% to the main figure. Returns a cell array of 'on' for attached
% and 'off' for not attached (these are useful for setting the
% 'Checked' property of the dock submenus).
panels = [ui.panels.nav ui.panels.roi ui.panels.grayscale ...
		  ui.panels.info ui.panels.mesh ui.panels.colorbar];
for i = 1:6
	if get(panels(i), 'Parent')==ui.fig
		chk{i} = 'on';
	else
		chck{i} = 'off';
	end
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ X ] = X_Solver_first(D,rho)
  [U,S,V] = svd(D);
  S0 = diag(S);
  r = length(S0);
  P = [ones(r,1), 1-S0, 1/2/rho-S0];
  rt = zeros(r,1);
  for t = 1:r
    p = P(t,:);
    Delta = p(2)^2-4*p(1)*p(3);
    if Delta <= 0
      rt(t) = 0;
    else
      rts = roots(p);
      rts = sort(rts);
      if rts(1)*rts(2)<=0
        rt(t) = rts(2);
      elseif rts(2)<0
        rt(t) = 0;
      else
        funval = log(1+rts(2))+rho*(rts(2)-S0(t)).^2;
        if funval>log(1+0)+rho*(0-S0(t)).^2;
          rt(t) = 0;
        end
      end
    end
  end

  SSS = diag(rt);
  [m,n] = size(D);
  sig = zeros(m,n);
  sig(1:min(m,n),1:min(m,n)) = SSS;

  X = U*sig*V';
end

author: andrewssobral

repo: lrslibrary


text:
% VL_PEGASOS PEGASOS linear SVM solver
%   W = VL_PEGASOS(X, Y, LAMBDA) learns a linear SVM W given training
%   vectors X, their labels Y, and the regularization parameter LAMBDA
%   using the PEGASOS [1] solver. The algorithm finds a minimizer W of
%   the objective function
%
%     LAMBDA/2 |W|^2 + 1/N SUM_i LOSS(W, X(:,i), Y(i))
%
%   where LOSS(W,X,Y) = MAX(0, 1 - Y W'X) is the hinge loss and N is
%   the number of training vectors in X.
%
%   [W B INFO] = VL_SVMPEGASOS(X, Y, LAMBDA) learns a linear SVM W
%   and a bias B given training vectors X, their labels Y, and the
%   regularization parameter LAMBDA using the PEGASOS [1]
%   solver. INFO is a struct containing the input parameters plus
%   diagnostic informations:
%
%   energy::
%     SVM energy value.
%
%   iterations::
%     Number of iterations performed.
%
%   elapseTime::
%     Elapsed time since the start of the SVM learning.
%
%   regulizerTerm::
%     Value of the SVM regulizer term.
%
%   lossPos::
%     Value of loss function only for data points labeled positives.
%
%   lossNeg::
%     Value of loss function onlt for data points labeled negatives.
%
%   hardLossPos::
%     Number of mislabeled positive points.
%
%   hardLossNeg::
%     Number of mislabeled negative points.
%
%   ALGORITHM. PEGASOS is an implementation of stochastic subgradient
%   descent. At each iteration a data point is selected at random, the
%   subgradient of the cost function relative to that data point is
%   computed, and a step is taken in that direction. The step size is
%   inversely proportional to the iteration number. See [1] for
%   details.
%
%   VL_SVMPEGASOS() accepts the following options:
%
%   Epsilon:: [empty]
%     Specify the SVM stopping criterion threshold. If not
%     specified VL_SVMPEGASOS will finish when the maximum number
%     of iterations is reached. The stopping criterion is tested
%     after each ENERGYFREQ iteration.
%
%   MaxIterations:: [10 / LAMBDA]
%     Sets the maximum number of iterations.
%
%   BiasMultiplier:: [0]
%     Appends to the data X the specified scalar value B. This
%     approximates the training of a linear SVM with bias.
%
%   StartingModel:: [null vector]
%     Specify the initial value for the weight vector W.
%
%   StartingIteration:: [1]
%     Specify the iteration number to start from. The only effect
%     is to change the step size, as this is inversely proportional
%     to the iteration number.
%
%   StartingBias:: [0]
%     Specify the inital bias value.
%
%   BiasLearningRate:: [1]
%     Specify the frequency of the bias learning. The default
%     setting updates the bias at each iteration.
%
%   Permutation:: [empty]
%     Specify a permutation PERM to be used to sample the data (this
%     disables random sampling). Specifically, at the T-th iteration
%     the algorithm takes a step w.r.t. the PERM[T']-th data point,
%     where T' is T modulo the number of data samples
%     (i.e. MOD(T'-1,NUMSAMPLES)+1). PERM needs not to be
%     bijective. This allows specifying certain data points more or
%     less frequently, implicitly increasing their relative weight in
%     the error term. A common application is to balance an unbalanced
%     dataset.
%
%   DiagnosticFunction:: [empty]
%     Specify a function handle to be called every ENERGYFREQ iterations.
%
%   DiagnosticCallRef:: [empty]
%     Specify a paramater to be passed to the DIAGNOSTICFUNCTION handle.
%
%   EnergyFreq:: [100]
%     Specify how often the SVM energy is computed.
%
%   HOMKERMAP:: [empty]
%     Specify the use of an Homogeneus Kernel map for the training
%     data (See [2],[3]). The passed value N is such that a 2*N+1
%     dimensional approximated kernel map is computed. Each
%     training data point is expanded online into a vector of
%     dimension 2*N+1.
%
%   KChi2::
%     Compute the map for the Chi2 kernel.
%
%   KINTERS::
%     Compute the map for the intersection kernel.
%
%   KL1::
%     Same as KINTERS, but deprecated as the name is not fully
%     accurate.
%
%   KJS::
%     Compute the map for the JS (Jensen-Shannon) kernel.
%
%   Period:: [automatically tuned]
%     Set the period of the kernel specturm. The approximation is
%     based on periodicizing the kernel specturm. If not specified,
%     the period is automatically set based on the heuristic described
%     in [2].
%
%   Window:: [RECTANGULAR]
%     Set the window used to truncate the spectrum before The window
%     can be either RECTANGULAR or UNIFORM window. See [2] and the API
%     documentation for details.
%
%   Gamma:: [1]
%     Set the homogeneity degree of the kernel. The standard kernels
%     are 1-homogeneous, but sometimes smaller values perform better
%     in applications. See [2] for details.
%
%   Verbose::
%     Be verbose.
%
%   Example::
%     The options StartingModel and StartingIteration can be used
%     to continue training. I.e., the command
%
%       vl_twister('state',0) ;
%       w = vl_pegasos(x,y,lambda,'NumIterations',1000) ;
%
%     produces the same result as the sequence
%
%       vl_twister('state',0) ;
%       w = vl_pegasos(x,y,lambda,'NumIterations',500) ;
%       w = vl_pegasos(x,y,lambda,'NumIterations',500, ...
%                      'StartingIteration', 501, ...
%                      'StartingModel', w) ;
%
%   REFERENCES::
%     [1] S. Shalev-Shwartz, Y. Singer, N. Srebro, and
%     A. Cotter. Pegasos: Primal Estimated sub-GrAdient SOlver for
%     SVM. MBP, 2010.
%
%     [2] A. Vedaldi and A. Zisserman
%     `Efficient Additive Kernels via Explicit Feature Maps',
%     Proc. CVPR, 2010.
%
%     [3] A. Vedaldi and A. Zisserman
%     `Efficient Additive Kernels via Explicit Feature Maps',
%     PAMI, 2011 (submitted).
%
%   See also: VL_HOMKERMAP(), VL_HELP().

% AUTHORIGHTS

author: Cloud-CV

repo: object-proposals


text:
function RGB=showmask(V,M,display_flag);
% showmask(V,M);
%
% M is a nonneg. mask
% Jianbo Shi, 1997

V=V-min(V(:));
V=V/max(V(:));
V=.25+0.75*V; %brighten things up a bit

M=M-min(M(:));
M=M/max(M(:));

H=0.0+zeros(size(V));
S=M;
RGB=hsv2rgb(H,S,V);

%if nargin>2
   image(RGB)
   axis('image')
%end

author: jwyang


text:
% RMSEW - Evaluates weighted root mean squared error.
%
% E = RMSEW(Y,W)
%
% Y'*DIAG(W)*Y / SUM(W)

% Last modified 2011-01-11
% Copyright (c) Jaakko Luttinen (jaakko.luttinen@tkk.fi)

function e = rmsew(Y, W)

if nargin < 2
  W = ones(size(Y));
end

e = sqrt( ((W(:).*Y(:))'*Y(:)) / sum(W(:)) );

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function a_cr = r83_cr_fa ( n, a )

%*****************************************************************************80
%
%% R83_CR_FA decomposes a real tridiagonal matrix using cyclic reduction.
%
%  Discussion:
%
%    The R83 storage format is used for a real tridiagonal matrix.
%    The superdiagonal is stored in entries (1,2:N), the diagonal in
%    entries (2,1:N), and the subdiagonal in (3,1:N-1).  Thus, the
%    original matrix is "collapsed" vertically into the array.
%
%    Once R83_CR_FA has decomposed a matrix A, then R83_CR_SL may be used to solve
%    linear systems A * x = b.
%
%    R83_CR_FA does not employ pivoting.  Hence, the results can be more
%    sensitive to ill-conditioning than standard Gauss elimination.  In
%    particular, R83_CR_FA will fail if any diagonal element of the matrix
%    is zero.  Other matrices may also cause R83_CR_FA to fail.
%
%    R83_CR_FA can be guaranteed to work properly if the matrix is strictly
%    diagonally dominant, that is, if the absolute value of the diagonal
%    element is strictly greater than the sum of the absolute values of
%    the offdiagonal elements, for each equation.
%
%    The algorithm may be illustrated by the following figures:
%
%    The initial matrix is given by:
%
%          D1 U1
%          L1 D2 U2
%             L2 R83 U3
%                L3 D4 U4
%                   L4 D5 U5
%                      L5 D6
%
%    Rows and columns are permuted in an odd/even way to yield:
%
%          D1       U1
%             R83    L2 U3
%                D5    L4 U5
%          L1 U2    D2
%             L3 U4    D4
%                L5       D6
%
%    A block LU decomposition is performed to yield:
%
%          D1      |U1
%             R83   |L2 U3
%                D5|   L4 U5
%          --------+--------
%                  |D2'F3
%                  |F1 D4'F4
%                  |   F2 D6'
%
%    For large systems, this reduction is repeated on the lower right hand
%    tridiagonal subsystem until a completely upper triangular system
%    is obtained.  The system has now been factored into the product of a
%    lower triangular system and an upper triangular one, and the information
%    defining this factorization may be used by R83_CR_SL to solve linear
%    systems.
%
%  Example:
%
%    Here is how a R83 matrix of order 5 would be stored:
%
%       *  A12 A23 A34 A45
%      A11 A22 A33 A44 A55
%      A21 A32 A43 A54  *
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 March 2004
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Roger Hockney,
%    A fast direct solution of Poisson's equation using Fourier Analysis,
%    Journal of the ACM,
%    Volume 12, Number 1, pages 95-113, January 1965.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%    N must be positive.
%
%    Input, real A(3,N), the R83 matrix.
%
%    Output, real A_CR(3,2*N+1), factorization information.
%
  if ( n <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R83_CR_FA - Fatal error!\n' );
    fprintf ( 1, '  Nonpositive N = %d\n', n );
    return
  end

  a_cr = zeros ( 3, 2 * n + 1 );

  if ( n == 1 )
    a_cr(1,1:3) = 0.0;
    a_cr(2,1) = 0.0;
    a_cr(2,2) = 1.0 / a(2,1);
    a_cr(2,3) = 0.0;
    a_cr(3,1:3) = 0.0;
    return
  end
%
%  Zero out the workspace entries.
%
  a_cr(1,1) = 0.0;
  a_cr(1,2:n) = a(1,2:n);
  a_cr(1,n+1:2*n+1) = 0.0;

  a_cr(2,1) = 0.0;
  a_cr(2,2:n+1) = a(2,1:n);
  a_cr(2,n+2:2*n+1) = 0.0;

  a_cr(3,1) = 0.0;
  a_cr(3,2:n) = a(3,1:n-1);
  a_cr(3,n+1:2*n+1) = 0.0;

  il = n;
  ipntp = 0;

  while ( 1 < il )

    ipnt = ipntp;
    ipntp = ipntp + il;
    if ( mod ( il, 2 ) == 1 )
      inc = il + 1;
    else
      inc = il;
    end

    incr = floor ( inc / 2 );
    il = floor ( il / 2 );
    ihaf = ipntp + incr + 1;
    ifulp = ipnt + inc + 2;

    for ilp = incr : -1 : 1
      ifulp = ifulp - 2;
      iful = ifulp - 1;
      ihaf = ihaf - 1;
      a_cr(2,iful+1) = 1.0 / a_cr(2,iful+1);
      a_cr(3,iful+1)  = a_cr(3,iful+1)  * a_cr(2,iful+1);
      a_cr(1,ifulp+1) = a_cr(1,ifulp+1) * a_cr(2,ifulp+2);
      a_cr(2,ihaf+1)  = a_cr(2,ifulp+1) - a_cr(1,iful+1)  * a_cr(3,iful+1) ...
                                  - a_cr(1,ifulp+1) * a_cr(3,ifulp+1);
      a_cr(3,ihaf+1) = -a_cr(3,ifulp+1) * a_cr(3,ifulp+2);
      a_cr(1,ihaf+1) = -a_cr(1,ifulp+1) * a_cr(1,ifulp+2);
    end

  end

  a_cr(2,ipntp+2) = 1.0 / a_cr(2,ipntp+2);

  return
end


text:
function cfg_defaults = cfg_mlbatch_defaults

% function cfg_defaults = cfg_mlbatch_defaults
% This file contains defaults that control the behaviour and appearance
% of matlabbatch.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: cfg_mlbatch_defaults.m 6641 2015-12-11 10:17:11Z volkmar $

rev = '$Rev: 6641 $'; %#ok

try
    % Font definition for cfg_ui user interface
    % cfg_defaults.cfg_ui.Xfont is a param/value list of property values as returned by uisetfont
    % lfont: used in lists, value edit dialogues etc.
    cfg_defaults.cfg_ui.lfont = {'FontAngle','normal',...
        'FontName',get(0,'factoryTextFontName'),...
        'FontSize',12,...
        'FontUnits','points',...
        'FontWeight','normal'};
    % bfont: used for buttons
    cfg_defaults.cfg_ui.bfont = {'FontAngle',get(0, 'factoryUicontrolFontAngle'),...
        'FontName',get(0,'factoryUicontrolFontName'),...
        'FontSize',get(0, 'factoryUicontrolFontSize'),...
        'FontUnits',get(0, 'factoryUicontrolFontUnits'),...
        'FontWeight',get(0, 'factoryUicontrolFontWeight')};
    % Toggle ExpertEdit mode. Value can be 'on' or 'off'
    cfg_defaults.cfg_ui.ExpertEdit = 'off';
catch
    cfg_defaults.cfg_ui = false;
end

% cfg_util
% Parallel execution of independent modules
% Currently, this does not run modules in parallel, but it may reorder
% execution order of modules: all modules without dependencies will be run
% before modules with dependencies will be harvested again. If some modules
% have side effects (e.g. "Change Directory") that are not encoded as
% dependency, this may lead to unwanted results. Disabling parallel
% execution incurs an overhead during job execution because the job
% must be harvested more often.
cfg_defaults.cfg_util.runparallel = false;
% cfg_util('genscript',...) hook to add application specific initialisation
% code to generated scripts. This must be a function handle that takes no
% input arguments and returns two output arguments - the code to be
% inserted as a string array, and a flag indicating whether cfg_util should
% add job execution code (flag == true) or not (flag == false).
% The generated code will be inserted after the for loop which collects the
% input. In the generated code, variables 'jobs' and 'inputs' can be
% referenced. These will hold the jobs and corresponding inputs.
cfg_defaults.cfg_util.genscript_run = [];
% A diary of command window output can be kept for jobs. This is useful for
% debugging various problems related to job execution.
cfg_defaults.cfg_util.run_diary = false;

% Message defaults
cfg_defaults.msgdef.identifier  = 'cfg_defaults:defaultmessage';
cfg_defaults.msgdef.level       = 'info'; % one of 'info', 'warning', 'error'
cfg_defaults.msgdef.destination = 'stdout'; % one of 'none', 'stdout',
                                            % 'stderr', 'syslog'. Errors
                                            % will always be logged to
                                            % the command window, and
                                            % additionally to syslog, if specified
cfg_defaults.msgdef.verbose     = 'off';
cfg_defaults.msgdef.backtrace   = 'off';

cfg_defaults.msgcfg(1)             = cfg_defaults.msgdef;
cfg_defaults.msgcfg(1).identifier  = 'matlabbatch:run:jobfailederr';
cfg_defaults.msgcfg(1).level       = 'error';
cfg_defaults.msgcfg(2)             = cfg_defaults.msgdef;
cfg_defaults.msgcfg(2).identifier  = 'matlabbatch:cfg_util:addapp:done';
cfg_defaults.msgcfg(2).destination = 'none';
cfg_defaults.msgcfg(3)             = cfg_defaults.msgdef;
cfg_defaults.msgcfg(3).identifier  = 'matlabbatch:initialise:invalid';
cfg_defaults.msgcfg(3).level       = 'error';
cfg_defaults.msgcfg(4)             = cfg_defaults.msgdef;
cfg_defaults.msgcfg(4).identifier  = 'matlabbatch:fopen';
cfg_defaults.msgcfg(4).level       = 'error';
cfg_defaults.msgcfg(5)             = cfg_defaults.msgdef;
cfg_defaults.msgcfg(5).identifier  = 'cfg_getfile:notcolumncell';
cfg_defaults.msgcfg(5).level       = 'error';
cfg_defaults.msgcfg(6)             = cfg_defaults.msgdef;
cfg_defaults.msgcfg(6).identifier  = 'matlabbatch:subsref:cfg_dep:multisubs';
cfg_defaults.msgcfg(6).level       = 'warning';
cfg_defaults.msgcfg(6).backtrace   = 'on';

cfg_defaults.msgtpl( 1)             = cfg_defaults.msgdef;
cfg_defaults.msgtpl( 1).identifier  = '^matlabbatch:subsasgn';
cfg_defaults.msgtpl( 1).level       = 'error';
cfg_defaults.msgtpl( 2)             = cfg_defaults.msgdef;
cfg_defaults.msgtpl( 2).identifier  = '^matlabbatch:subsref';
cfg_defaults.msgtpl( 2).level       = 'error';
cfg_defaults.msgtpl( 3)             = cfg_defaults.msgdef;
cfg_defaults.msgtpl( 3).identifier  = '^matlabbatch:constructor';
cfg_defaults.msgtpl( 3).level       = 'error';
cfg_defaults.msgtpl( 4)             = cfg_defaults.msgdef;
cfg_defaults.msgtpl( 4).identifier  = '^matlabbatch:deprecated';
cfg_defaults.msgtpl( 4).destination = 'none';
cfg_defaults.msgtpl( 5)             = cfg_defaults.msgdef;
cfg_defaults.msgtpl( 5).identifier  = '^matlabbatch:usage';
cfg_defaults.msgtpl( 5).level       = 'error';
cfg_defaults.msgtpl( 6)             = cfg_defaults.msgdef;
cfg_defaults.msgtpl( 6).identifier  = '^matlabbatch:setval';
cfg_defaults.msgtpl( 6).destination = 'none';
cfg_defaults.msgtpl( 7)             = cfg_defaults.msgdef;
cfg_defaults.msgtpl( 7).identifier  = '^matlabbatch:run:nomods';
cfg_defaults.msgtpl( 7).level       = 'info';
cfg_defaults.msgtpl( 8)             = cfg_defaults.msgdef;
cfg_defaults.msgtpl( 8).identifier  = '^matlabbatch:cfg_struct2cfg';
cfg_defaults.msgtpl( 8).destination = 'none';
cfg_defaults.msgtpl( 9)             = cfg_defaults.msgdef;
cfg_defaults.msgtpl( 9).identifier  = '^MATLAB:inputdlg';
cfg_defaults.msgtpl( 9).level       = 'error';
cfg_defaults.msgtpl(10)             = cfg_defaults.msgdef;
cfg_defaults.msgtpl(10).identifier  = '^MATLAB:listdlg';
cfg_defaults.msgtpl(10).level       = 'error';
cfg_defaults.msgtpl(11)             = cfg_defaults.msgdef;
cfg_defaults.msgtpl(11).identifier  = '^MATLAB:num2str';
cfg_defaults.msgtpl(11).level       = 'error';
cfg_defaults.msgtpl(12)             = cfg_defaults.msgdef;
cfg_defaults.msgtpl(12).identifier  = '^matlabbatch:ok_subsasgn';
cfg_defaults.msgtpl(12).destination = 'none';
cfg_defaults.msgtpl(13)             = cfg_defaults.msgdef;
cfg_defaults.msgtpl(13).identifier  = 'matlabbatch:checkval:numcheck:transposed';
cfg_defaults.msgtpl(13).destination = 'none';

% value check for cfg_branch/choice/repeat items - set to false after
% configuration has been initialised to speed up job
% initialisation/harvest/run - set this to true if you want to debug some
% configuration or the batch system itself
cfg_defaults.cfg_item.checkval = false;

% verbosity of code generation for cfg_dep objects
% true  - one-line code containing only source information
% false - full code, including source and target specifications
cfg_defaults.cfg_dep.gencode_short = true;

author: spm