text:
function diagnostic = solvesdp_multiple(varargin);

yalmiptime = clock;
h = varargin{2};
varargin{2} = sum(recover(depends(h)));

if ~is(h,'linear')
    error('Parts of your matrix objective is not linear (multiple solutions can currently only be obtained for linear objectives)');
end

if is(h,'complex')
    error('Parts of your matrix objective is complex-valued (which makes no sense since complex numbers have no natural ordering');
end

if nargin<3
    ops = sdpsettings('saveyalmipmodel',1);
    varargin{3} = ops;
else
    ops = varargin{3};
    varargin{3}.saveyalmipmodel = 1;
end
varargin{3}.pureexport = 1;
yalmiptime = clock;
model = solvesdp(varargin{:});
diagnostic.yalmiptime = etime(clock,yalmiptime);
yalmiptime = clock;

h_variables = getvariables(h);
h_base = getbase(h);
for i = 1:length(h)
    model.f = h_base(i,1);
    model.c = mapObjective(h_variables,model.used_variables,h_base(i,2:end));

    % *************************************************************************
    % TRY TO SOLVE PROBLEM
    % *************************************************************************
    if ops.debug
        eval(['output = ' model.solver.call '(model);']);
    else
        try
            eval(['output = ' model.solver.call '(model);']);
        catch
            output.Primal = zeros(length(model.c),1)+NaN;
            output.Dual  = [];
            output.Slack = [];
            output.solvertime   = nan;
            output.solverinput  = [];
            output.solveroutput = [];
            output.problem = 9;
            output.infostr = yalmiperror(output.problem,lasterr);
        end
    end

    if ops.dimacs
        try
            b = -model.c;
            c = model.F_struc(:,1);
            A = -model.F_struc(:,2:end)';
            x = output.Dual;
            y = output.Primal;
            % FIX this nonlinear crap (return variable type in
            % compilemodel)
            if options.relax == 0 & any(full(sum(model.monomtable,2)~=0))
                if ~isempty(find(sum(model.monomtable | model.monomtable,2)>1))
                    z=real(exp(model.monomtable*log(y+eps)));
                    y = z;
                end
            end

            if isfield(output,'Slack')
                s = output.Slack;
            else
                s = [];
            end

            dimacs = computedimacs(b,c,A,x,y,s,model.K);
        catch
            dimacs = [nan nan nan nan nan nan];
        end
    else
        dimacs = [nan nan nan nan nan nan];
    end

    % ********************************
    % ORIGINAL COORDINATES
    % ********************************
    if isempty(output.Primal)
       output.Primal = zeros(size(model.recoverdata.H,2),1);
    end
    output.Primal(:,i) = model.recoverdata.x_equ+model.recoverdata.H*output.Primal;

    % ********************************
    % OUTPUT
    % ********************************
    diagnostic.yalmiptime = diagnostic.yalmiptime + etime(clock,yalmiptime)-output.solvertime;
    diagnostic.solvertime(:,i) = output.solvertime;
    try
        diagnostic.info = output.infostr;
    catch
        diagnostic.info = yalmiperror(output.problem,model.solver.tag);
    end
    diagnostic.problem(:,i) = output.problem;
    diagnostic.dimacs(:,i) = dimacs;

    % Some more info is saved internally
    solution_internal = diagnostic;
    solution_internal.variables = model.recoverdata.used_variables(:);
    solution_internal.optvar = output.Primal(:,i);

    if ~isempty(model.parametric_variables)
        diagnostic.mpsol = output.solveroutput;
        options.savesolveroutput = actually_save_output;
    end;

    if model.options.savesolveroutput
        diagnostic.solveroutput = output.solveroutput;
    end
    if model.options.savesolverinput
        diagnostic.solverinput = output.solverinput;
    end
    if model.options.saveyalmipmodel
        diagnostic.yalmipmodel = model;
    end

    if ops.warning && warningon && isempty(strfind(output.infostr,'No problems detected'))
        disp(['Warning: ' output.infostr]);
    end

    if ismember(output.problem,ops.beeponproblem)
        try
            beep; % does not exist on all ML versions
        catch
        end
    end

    % And we are done! Save the result
    if ~isempty(output.Primal)
        yalmip('setsolution',solution_internal,i);
    end
    if model.options.saveduals & model.solver.dual
        if isempty(model.Fremoved) | (nnz(model.Q)>0)
            try
                setduals(F,output.Dual,model.K);
            catch
            end
        else
            try
                % Duals related to equality constraints/free variables
                % have to be recovered b-A*x-Ht == 0
                b = -model.oldc;
                A = -model.oldF_struc(1+model.oldK.f:end,2:end)';
                H = -model.oldF_struc(1:model.oldK.f,2:end)';
                x = output.Dual;
                b_equ = b-A*x;
                newdual =  H\b_equ;
                setduals(model.Fremoved + F,[newdual;output.Dual],model.oldK);
            catch
                % this is a new feature...
                disp('Dual recovery failed. Please report this issue.');
            end
        end
    end
end

function newbase = mapObjective(local_vars,global_vars,base)
newbase = spalloc(length(global_vars),1,nnz(base));
for i = 1:length(local_vars)
    j = find(local_vars(i)==global_vars);
    newbase(j) = base(i);
end

function yesno = warningon

s = warning;
yesno = isequal(s,'on');


text:
%% STUDY_01_wrist_scaffold_single_fiber_iFEA_elastic
% Below is a demonstration for:
%
% * Inverse FEA based material parameter optimisation
% * febio_spec version 3.0
% * febio, FEBio
% * hexahedral elements, hex8
% * static, solid
% * isotropic elastic, 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

% FEA control settings
timeEnd=10;
numTimeSteps=10; %Number of time steps desired to reach yield strength
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=60; %Maximum number of retires
dtmin=(timeEnd/numTimeSteps)/100; %Minimum time step size
dtmax=timeEnd/numTimeSteps; %Maximum time step size
tStep=timeEnd/numTimeSteps;%s
min_residual=1e-20;
symmetric_stiffness=1;
lstol=0.9;

%Initial material parameter set
k_factor=100;

%Material model is Ogden 2nd order with initial estimates
c1_ini=57;
m1_ini=2;
c2_ini=57;
m2_ini=-2;
k_ini=(c1_ini+c2_ini)*k_factor;
P=[c1_ini m1_ini c2_ini m2_ini];

% Experimental elastic properties
matname = fullfile(loadPathMatProp,'Elastic_properties.mat');
load(matname);

cFigure;
title('Experimental elastic region','fontSize',fontSize);
hold on;
xlabel('Average displacement, mm','FontSize',fontSize);
ylabel('Average load, N','FontSize',fontSize);
plot(disp_exp,load_exp,'k.-','LineWidth',3);
set(gca,'FontSize',fontSize);
camlight headlight;
displacementMagnitude=max(disp_exp);

%% DEFINE BC's

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

%% Find side face sets
[Nb]=patchNormal(F,V);
zVec=[0 0 1];
d=dot(Nb,zVec(ones(size(Nb,1),1),:),2);
Z=V(:,3);
ZF=mean(Z(F),2);
pointSpacing=mean(patchEdgeLengths(F,V));

logicSides_F1=(d>-1) & ZF>=(max(V(:,3))-2*pointSpacing);
logicSides_F2=(d<1) & ZF<=(min(V(:,3))+2*pointSpacing);

logicSides_FS=logicSides_F1;
logicSides_FL=logicSides_F2;

faceBoundaryMarker(logicSides_FS)=1;
faceBoundaryMarker(logicSides_FL)=2;

%Visualize BC's
cFigure; hold on;
patch('faces',F,'vertices',V,'FaceColor',[0, 0.5, 0],'FaceAlpha',0.1,'EdgeColor','None');

hl(1)=plotV(V(bcRigidList,:),'k.','MarkerSize',markerSize);
hl(2)=plotV(V(bcPrescribeList,:),'.','Color',[0.6350, 0.0780, 0.1840],'MarkerSize',markerSize);

legend(hl,{'BC full support','BC prescribed Z displacement'})
axisGeom;
camlight headlight;
set(gca,'FontSize',fontSize);

axes('position',[.2 .6 .2 .2], 'NextPlot', 'add');
plotV(V(bcPrescribeList,:),'.','Color',[0.6350, 0.0780, 0.1840],'MarkerSize',markerSize*2);
patch('faces',F(faceBoundaryMarker==1,:),'vertices',V,'FaceColor',[0, 0.5, 0],'FaceAlpha',0.5,'EdgeColor','k');
view(3);
axis equal; axis vis3d; axis tight;
set(gca,'XColor','none','YColor','none','ZColor','none');

axes('position',[.2 .3 .2 .2], 'NextPlot', 'add');
plotV(V(bcRigidList,:),'k.','MarkerSize',markerSize*2);
patch('faces',F(faceBoundaryMarker==2,:),'vertices',V,'FaceColor',[0, 0.5, 0],'FaceAlpha',0.5,'EdgeColor','k');
view(3);
axis equal; axis vis3d; axis tight;
set(gca,'XColor','none','YColor','none','ZColor','none');

drawnow;

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

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

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

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

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

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

%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 this part
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;

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

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(:);

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

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

% -> Prescribe boundary conditions
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=1;
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; timeEnd displacementMagnitude];

%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='sz';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.element_data{1}.VAL=1:size(E,1);

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

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

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

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

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

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.runMode='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;

    %%
    % 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 Z force across nodes
    FZ_set=squeeze(FZ_set);

    disp_sim=DZ_set;
    load_sim=FZ_set*NumStrands;

    [~, indUnique] = unique(disp_sim);
    disp_sim=disp_sim(indUnique);
    load_sim=load_sim(indUnique);
    load_sim_exp = interp1(disp_sim,load_sim,disp_exp,'pchip');

    % Visualize curves
    hf=cFigure;

    hold on;
    title('Uniaxial load-displacement curve','FontSize',fontSize);
    xlabel('Displacement','FontSize',fontSize); ylabel('Load','FontSize',fontSize);

    hl(1)=plot(disp_sim,load_sim,'r.-','lineWidth',lineWidth*1.1,'markerSize',markerSize*2);
    hl(2)=plot(disp_exp,load_exp,'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={'c1','m1','c2','m2','k'}; %Parameter names
mat_struct.par_values={c1_ini m1_ini c2_ini m2_ini k_ini}; %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)/10000 -10000 P(3)/10000 -10000]' [P(1)*10000 10000 P(3)*10000 10000]']; %Parameter bounds

febioAnalysis.disp_on=0;
febioAnalysis.disp_log_on=0;

%What should be known to the objective function:
objectiveStruct.h=hl(1);
objectiveStruct.bcPrescribeList=bcPrescribeList;
objectiveStruct.disp_exp=disp_exp;
objectiveStruct.load_exp=load_exp;
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.displacement=fullfile(savePath,febioLogFileName_disp);
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),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),Pn,[],[],OPT_options);
end

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

%% 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('Displacement (mm)','FontSize',fontSize); ylabel('Force (N)','FontSize',fontSize); zlabel('Z','FontSize',fontSize); hold on;

Hn(1)=plot(OPT_stats_out.disp_sim,OPT_stats_out.load_sim,'r.-','lineWidth',lineWidth*1.1,'markerSize',markerSize*2);
Hn(2)=plot(disp_exp,load_exp,'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)

%%

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(1)+P(3))*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:numel(mat_struct.par_names)
    parNameNow=mat_struct.par_names{q};
    parValuesNow=mat_struct.par_values{q};
    febio_spec.Material.material{matId}.(parNameNow)=mat2strIntDouble(parValuesNow);
end
febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode

disp('Done')

%% START FEBio

[runFlag]=runMonitorFEBio(objectiveStruct.febioAnalysis);

pause(0.1);

bcPrescribeList=objectiveStruct.bcPrescribeList;
disp_exp=objectiveStruct.disp_exp;
load_exp=objectiveStruct.load_exp;

if runFlag==1

    %Importing displacement
    [~,N_disp_mat,~]=importFEBio_logfile(objectiveStruct.run_output_names.displacement);
    N_disp_mat=N_disp_mat(:,2:end,:); %Remove index
    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; %With added zeros for t=0

    %Derive applied stretch
    DZ_set=N_disp_mat(bcPrescribeList,3,:); %Final nodal displacements
    DZ_set=mean(DZ_set,1);
    DZ_set=squeeze(DZ_set);

    %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 Z force across nodes
    FZ_set=squeeze(FZ_set);

    disp_sim=DZ_set;
    load_sim=FZ_set*objectiveStruct.NumStrands; %Multiply force by number of strands

    %Remove double entries
    [~, indUnique] = unique(disp_sim);
    disp_sim=disp_sim(indUnique);
    load_sim=load_sim(indUnique);

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

    %Interpolate experiment onto simulated points
    load_sim_exp = interp1(disp_sim,load_sim,disp_exp,'pchip');

    %Derive Fopt
    loadDev=load_exp-load_sim_exp;

    switch objectiveStruct.method
        case 1
            Fopt=sum((loadDev).^2); %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.disp_sim=disp_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 run_writesnap(varargin)
% Send real-time BCI outputs to a SNAP-based stimulus presentation system.
% run_writetcp(Model,SourceStream,OutputHost,OutputPort,OutputForm,MessageFormat,UpdateFrequency,PredictorName,ConnTimeout)
%
% This function runs in the background and processes data from some MATLAB stream (created with some
% other background input plugin, e.g., the BioSemi reader). The processed data is periodically
% forwarded to a program running the SNAP framework.
%
% In:
%   Model : A model data structure (as obtained from bci_train) based on which the predictor shall be
%           created; typically this is a model struct, but for convenience it can be a file name,
%           variable name in the base workspace, or a cell array of {file name, variable name} to
%           refer to a variable inside a .mat file. The model is not modified by this function.
%           (default: 'lastmodel')
%
%   SourceStreamNames : Optional names of stream data structures in the MATLAB base workspace to
%                       consider as possible data sources (previously created with onl_newstream);
%                       if a stream contains all channels that are needed by the predictor, or
%                       alternatively has the right number and type of channels it will be considered
%                       as a potential source stream unless ambiguous. (default: 'laststream')
%
%   OutputHost : destination host name to send results to (computer name, URL or IP address)
%                (default: 'localhost')
%
%   OutputPort : destination port to which the results are sent (default: 7897)
%
%   OutputForm : output data form, see onl_predict (default: 'distribution')
%
%   TargetVariable : name of the variable in the SNAP module which should carry the BCI predictions
%                    (default: 'bci')
%
%   UpdateFrequency : The rate at which new outputs will be computed. (default: 10)
%
%   PredictAt : Predict at markers. If nonempty, this is a cell array of online target markers relative
%               to which predictions shall be made. If empty, predictions are always made on the most recently
%               added sample. (default: {})
%
%   PredictorName : Name of the predictor to be created; a variable of this name will be created in
%                   the MATLAB base workspace to hold the predictor's state. If a variable with this
%                   name already exists it will be overridden. (default: 'lastpredictor')
%
%   Verbose : whether to display verbose outputs (e.g. connection failure) (default: false)
%
% Examples:
%   % open a new BCILAB processing stream, using the previously learned predictive model 'mymodel',
%   % and reading from a previously opened input stream named 'mystream'; send outputs to some IP:port
%   run_writesnap('mymodel','mystream','192.168.1.5')
%
%   % as before, but pass arguments by name
%   run_writetcp('Model','mymodel','SourceStream','mystream','OutputHost','192.168.1.5')
%
%   % as before, but transmit the data in a raw binary encoding
%   run_writetcp('Model','mymodel','SourceStream','mystream','OutputHost','192.168.1.5','TargetVariable','alpha')
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-11-19

% declare the name of this component (shown in the menu)
declare_properties('name','SNAP');

% define arguments
opts = arg_define(varargin, ...
    arg({'pred_model','Model'}, 'lastmodel', [], 'Predictive model. As obtained via bci_train or the Model Calibration dialog.','type','expression'), ...
    arg({'in_stream','SourceStreamNames','SourceStream'}, 'laststream',[],'Input Matlab stream name(s). Optional names of stream data structures in the MATLAB base workspace to consider as possible data sources (previously created with onl_newstream); if a stream contains all channels that are needed by the predictor, or alternatively has the right number and type of channels it will be considered as a potential source stream unless ambiguous.','typecheck',false,'shapecheck',false), ...
    arg({'out_hostname','OutputHost','Host'}, 'localhost',[],'Destination TCP hostname. Can be a computer name, URL, or IP address.'), ...
    arg({'out_port','OutputPort','Port'}, 7897, [],'Destination TCP port. Matches the SNAP default.'), ...
    arg({'out_form','OutputForm','Form'},'expectation',{'expectation','distribution','mode','raw'},'Output form. Can be the expected value (posterior mean) of the target variable, or the distribution over possible target values (probabilities for each outcome, or parametric distribution), or the mode (most likely value) of the target variable.'), ...
    arg({'target_variable','TargetVariable','Target'},'bci',[],'Target Variable. This is the variable name in the SNAP module that should receive the BCI value'), ...
    arg({'update_freq','UpdateFrequency'},10,[0 Inf],'Update frequency. The rate at which new outputs will be computed.'), ...
    arg({'predict_at','PredictAt'}, {},[],'Predict at markers. If nonempty, this is a cell array of online target markers relative to which predictions shall be made. If empty, predictions are always made on the most recently added sample.','type','expression'), ...
    arg({'pred_name','PredictorName'}, 'lastpredictor',[],'Name of new predictor. A variable of this name will be created in the MATLAB base workspace to hold the predictor''s state. If a variable with this name already exists it will be overridden.'), ...
    arg({'verbose_output','Verbose'}, false,[],'Verbose output. Whether to display verbose outputs (e.g., connection failure).'));

% convert format strings to formatting functions
if ~isvarname(opts.target_variable)
    disp('Note: The given target variable is likely not a valid variable name.'); end

% get a fresh connection id
id = sprintf('c%.0f',fresh_id('snapsocket'));

% start background writer job
onl_write_background( ...
    'ResultWriter',@(y)send_message(y,opts.verbose_output,id,opts.out_hostname,opts.out_port,opts.target_variable),...
    'MatlabStream',opts.in_stream, ...
    'Model',opts.pred_model, ...
    'OutputFormat',opts.out_form, ...
    'UpdateFrequency',opts.update_freq, ...
    'PredictorName',opts.pred_name, ...
    'PredictAt',opts.predict_at, ...
    'Verbose',opts.verbose_output, ...
    'StartDelay',0,...
    'EmptyResultValue',[]);

disp('Now writing...');

% background message sending function
function send_message(yy,verbose_output,id,host,port,varname)
persistent conns;
try
    % try to connect if necessary
    if ~isfield(conns,id)
        conns.(id) = connect(host,port); end
    % for each prediction...
    for k=1:size(yy,1)
        y = yy(k,:);
        try
            % send it off
            strm = conns.(id).strm;
            if isscalar(y)
                strm.writeBytes(char([sprintf('setup %s=%.5f',varname,y) 10]));
            elseif ~isempty(y)
                strm.writeBytes(char(['setup ' varname '=(' sprintf('%.5f,',y) ')' 10]));
            end
            strm.flush();
        catch e1
            if strcmp(e1.identifier, 'MATLAB:Java:GenericException')
                % failed to send: try to re-connect...
                conns.(id) = connect(host,port);
            else
                rethrow(e1);
            end
        end
    end
catch e2
    if strcmp(e2.identifier, 'MATLAB:Java:GenericException')
        if verbose_output
            fprintf('Could not connect to %s:%.0f\n',host,port); end
    else
        hlp_handleerror(e2);
    end
end

function newconn = connect(host,port)
import java.io.*
import java.net.*
import java.lang.*
newconn.sock = Socket();
newconn.sock.setTcpNoDelay(1);
newconn.sock.setSoTimeout(5);
newconn.sock.connect(InetSocketAddress(host,port),1);
newconn.strm = DataOutputStream(newconn.sock.getOutputStream());
disp(['Connection to ' host ' established.']);

repo: REST


text:
%  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 should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
clc; zoom on; clf; echo on;

%TFDEMO1 Introduction to the Time-Frequency Toolbox.
%	O. Lemoine - May 1996.
%	Copyright (c) CNRS.

% Welcome to the Time-Frequency Toolbox for MATLAB. This demonstration
% follows the plan of the tutorial, and consider exactly the same examples.
% Therefore, for further information about these illustrations, we advise
% you to refer to the corresponding chapter of the tutorial.
%
% Press any key to continue...

pause; clc;

% First let's create an analytic linear frequency modulated signal, whose
% normalized frequency is changing from 0 to 0.5 :

sig1=fmlin(128,0,0.5);
plot(real(sig1)); axis([1 128 -1 1]);
xlabel('Time'); ylabel('Real part');
title('Linear frequency modulation'); grid

% From this time-domain representation, it is difficult to say what kind
% of modulation is contained in this signal.
%
% Press any key to continue...

pause; clc;

% Now let's consider its energy spectrum :

dsp1=fftshift(abs(fft(sig1)).^2);
plot((-64:63)/128,dsp1);
xlabel('Normalized frequency'); ylabel('Squared modulus');
title('Spectrum'); grid

% We still can not say, from this plot, anything about the evolution in
% time of the frequency content.
%
% In order to have a more informative description of such a signal, it
% would be better to directly represent their frequency content while
% still keeping the time description parameter : this is precisely the
% aim of time-frequency analysis. To illustrate this, let's try the
% Wigner-Ville distribution on this signal :
%
% press any key to continue...

pause;

tfrwv(sig1);

% We can see on this representation that the linear progression of the
% frequency with time, from 0 to 0.5, is clearly shown.
%
% Press any key to continue...

pause; clc;

% If we now add some complex white Gaussian noise on this signal, with
% a 0 dB signal to noise ratio,

sig2=sigmerge(sig1,noisecg(128),0);
Min=min(real(sig2)); Max=max(real(sig2));
clf; plot(real(sig2)); axis([1 128 Min Max]);
xlabel('Time'); ylabel('Real part');
title('Linear frequency modulation plus noise'); grid

% press any key to continue...

pause;

% and consider the spectrum of it :

dsp2=fftshift(abs(fft(sig2)).^2);
plot((-64:63)/128,dsp2);
xlabel('Normalized frequency'); ylabel('Squared modulus');
title('Spectrum'); grid

% It is worse than before to interpret these plots. On the other hand, the
% Wigner-Ville distribution still show quite clearly the linear progression
% of the frequency with time :
%
% press any key to continue...

pause;

tfrwv(sig2);

% press any key to continue...

pause; clc;

% The second example we consider is a bat sonar signal, recorded with a
% sampling frequency of 230.4 kHz and an effective bandwidth equal to
% [8 kHz, 80 kHz].
%  First, load the signal from the MAT-file bat.mat :

echo off
DirectoryStr='../data/';
while (exist([DirectoryStr 'bat.mat'])==0),
 fprintf('I can''t find %s\n', [DirectoryStr 'bat.mat']);
 DirectoryStr=input('name of the directory where bat.mat is : ','s');
end;
eval(['load ' DirectoryStr 'bat.mat']);
echo on

Fe=230.4; t0=(0:2047).'/Fe;
clf; plot(t0,bat); xlabel('Time [ms]');
axis([t0(1) t0(2048) -900 800]); grid;

% From this plot, we can not say precisely what is the frequency content
% at each time instant t ; similarly, if we look at its spectrum,
%
% press any key to continue...

pause;
dsp=fftshift(abs(fft(bat)).^2); dsp=dsp/max(dsp);
f0=(-1024:1023)*Fe/2048;
plot(f0,10*log10(dsp)); xlabel('Frequency [kHz]'); ylabel('Squared modulus  (dB)');
axis([0 Fe/2,-30 0])
title('Spectrum  (dB)'); grid; pause

% we can not say at what time the signal is located around 38 kHz, and at
% what time around 40 kHz. Let us now consider a representation called
% the pseudo Wigner-Ville distribution, applied on the most interesting
% part of this signal :
%
% press any key to continue...

tfrpwv(hilbert(bat),8:8:2048,256,tftb_window(127,'hanning'),1);

% We then have a nice description of its spectral content varying with
% time : it is narrow-band signal, whose frequency content is decreasing
% from around 55 kHz to 38kHz, with a non-linear frequency modulation
% (approximately of hyperbolic shape).
%
% press any key to continue...

pause; clc;

% The last introductory example presented here is a transient signal
% embedded in a -5 dB white Gaussian noise. This transient signal is a
% constant frequency modulated by a one-sided exponential amplitude :

trans=amexpo1s(92,1,64).*fmconst(92);
sig=[zeros(100,1) ; trans ; zeros(64,1)];
sign=sigmerge(sig,noisecg(256),-5);
Min=min(real(sign)); Max=max(real(sign));
subplot(211); plot(real(sign)); axis([1 256 Min Max]);
xlabel('Time'); title('Noisy transient signal'); grid
dsp=fftshift(abs(fft(sign)).^2);
subplot(212); plot((-128:127)/256,dsp); grid
xlabel('Normalized frequency'); title('Energy spectrum');

% From these representations, it is difficult to localize precisely the
% signal in the time-domain as well as in the frequency domain.  Now let us
% have a look at the spectrogram of this signal :
%
% press any key to continue...

pause;
tfrsp(sign);

% the transient signal appears distinctly around the normalized frequency
% 0.25, and between time points 125 and 160.
%
% Press any key to return to the main menu.

echo off

author: HeLiangHIT


text:
function indx_max = r8vec_indexed_heap_d_max ( n, a, indx )

%*****************************************************************************80
%
%% R8VEC_INDEXED_HEAP_D_MAX: maximum value in heap descending indexed R8VEC.
%
%  Discussion:
%
%    An R8VEC is a vector of R8's.
%
%    An indexed R8VEC is an R8VEC of data values, and an R8VEC of N indices,
%    each referencing an entry of the data vector.
%
%    This is one of three functions needed to model a priority queue.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 August 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Thomas Cormen, Charles Leiserson, Ronald Rivest,
%    Introduction to Algorithms,
%    MIT Press, 2001,
%    ISBN: 0262032937,
%    LC: QA76.C662.
%
%  Parameters:
%
%    Input, integer N, the number of items in the index vector.
%
%    Input, real A(*), the data vector.
%
%    Input, integer INDX(N), the index vector.
%
%    Output, integer INDX_MAX, the index in A of the maximum value
%    in the heap.
%
  indx_max = indx(1);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function dx = HIVsys_ZURAKOWSKI(t,x,u,p)

try
    run_HIV_params

    dx = zeros(5,1);
    dx(1) = lambda1 - d*x(1) - alpha1*(1-eta*u(1))*x(1)*x(2);
    dx(2) = alpha1*(1-eta*u(1))*x(2)*x(1) - a*x(2) - p1*x(4)*x(2) - p2*x(5)*x(2);
    dx(3) = c2*x(1)*x(2)*x(3) - c2*q*x(2)*x(3) - b2*x(3);
    dx(4) = c1*x(2)*x(4) - b1*x(4);
    dx(5) = c2*q*x(2)*x(3) - h*x(5);
catch
    keyboard
end

author: eurika-kaiser

repo: SINDY-MPC


text:
function obj = misc_history(obj)
%MISC_HISTORY - Records the history of BTB toolbox function calls that
%               affected the object.
%
%Synopsis:
%  MISC_HISTORY(OBJ)
%
%Arguments:
%  OBJ:     struct with neuro data, or montage, or marker
%
%Returns:
%  OBJ:     struct with extra history entry in the substruct .history with
%           the following fields:
%
%  fcn:           function handle
%  fcn_params:    names of named function parameters.
%  varargin:      optional arguments
%  date:          date the function was executed
%
%  All named function arguments also appear as fields except for the object
%  that is modified (e.g. the EEG struct). The latter appears in the
%  parameter list as 'obj'. (note: to avoid confusion proc_ functions
%  should not have another parameter called 'obj')
%
% Note: misc_history saves all arguments passed to functions unless they
% exceed a certain byte size.
%
%See also: misc_applyHistory

% 08-2012 Matthias Treder
global BTB

if ~isempty(BTB.History) && BTB.History==0 ||isempty(obj)
  return
end

misc_checkType(obj,'STRUCT');

if isfield(obj,'history')
  ht = obj.history;
  N = numel(ht)+1;
else
  ht = cell(1,1);
  N = 1;
end
ht{N} = struct();

objname = inputname(1);  % Name of object in caller's workspace

maxSize = 10*1000*1000;   % max byte size of arguments saved in history
% Get function name (of caller)
ST = dbstack('-completenames',1);
ht{N}.fcn = eval(['@' ST(1).name]);

% Find labels of function parameters (arguments)
code = (evalc(['type ' ST(1).name])); % get function code
expr = ['^\s*function s*\[*\s*(\w+|\s*|,)*\]*\s*=\s*' ST(1).name '\s*\((\w+|\s+|,)*\)'];
token = regexpi(code,expr,'tokens');
token = strrep(token{1}{end},' ','');  % remove whitespace
token = regexp(token,',','split');   % split arguments by commas

% Remove params that have not been used (nargin < number of params)
nActualArguments= evalin('caller', 'nargin');
if nActualArguments < numel(token)
  token(nActualArguments+1:end) = [];
end

ht{N}.fcn_params = token(~ismember(token,'varargin','legacy'));

% Get argument values for named arguments
for ii=1:numel(ht{N}.fcn_params)
  if strcmp(ht{N}.fcn_params{ii},objname)
    ht{N}.fcn_params{ii} = 'obj';  % always call the dat/epo/cnt struct obj
  else
    s =evalin('caller',['whos(''' token{ii} ''')']);  % Get size of variable
    if s.bytes > maxSize
       ht{N}.(ht{N}.fcn_params{ii}) = sprintf('variable exceeds maximum size %d kb',maxSize/1000);
    else
       ht{N}.(ht{N}.fcn_params{ii}) = evalin('caller',ht{N}.fcn_params{ii});
    end
  end
end

% Get optional arguments (varargin)
if any(ismember(token,'varargin','legacy'))
  va = evalin('caller','varargin');
  % Check size of variable
  s=evalin('caller','whos(''varargin'')');
  if s.bytes > maxSize
     ht{N}.varargin_note = sprintf('varargin exceeds maximum size %d kb',maxSize/1000);
  else
    for ii=1:numel(va)
        ht{N}.(sprintf('varargin%d',ii)) = va{ii};
    end
  end
end

ht{N}.date= datestr(now,0);
obj.history = ht;

author: bbci


text:
function [x,fval,exitflag,info] = opti_hybrj(fun,grad,x0,opts)
%OPTI_HYBRJ Solve a SNLE using HYBRJ (Powell-Hybrid MINPACK Routine)
%
%   F(x) = 0
%
%   x = opti_hybrj(fun,grad,x0) solves a SNLE where fun is a vector of
%   functions to solve. grad is an optional gradient of the nonlinear
%   equations and x0 is a starting guess.
%
%   x = opti_hybrj(fun,grad,x0,opts) uses opts to pass optiset options to
%   the solver.
%
%   [x,fval,exitflag,info] = opti_hybrj(...) returns the objective value at
%   the solution, together with the solver exitflag, and an information
%   structure.
%
%   THIS IS A WRAPPER FOR HYBRJ + HYBRD
%   See supplied License

%   Copyright (C) 2012 Jonathan Currie (I2C2)

if(nargin < 4), opts = optiset; end
if(nargin < 3), error('HYBRJ requires at least 3 arguments'); end

%Setup display level
opts.display = dispLevel(opts.display);

%Check we have a valid x0
if(isempty(x0) || any(isnan(x0)))
    error('HYBRJ requires an initial guess, x0!');
end

t = tic;
% Run HYBRJ
[x, fval, exitflag, iter] = hybrj(fun,grad,x0,opts);

%Collect Results
info.FuncEvals = iter;
info.Time = toc(t);
info.Algorithm = 'HYBRJ: MINPACK Powell Hybrid';

switch(exitflag)
    case 1
        info.Status = 'Optimal';
    case 0
        info.Status = 'Exceeded Function Evaluations';
    case -1
        info.Status = 'Infeasible / Could not Converge';
    case -2
        info.Status = 'HYBRJ Error';
    case -5
        info.Status = 'User Exited';
    otherwise
        info.Status = 'HYBRJ Error';
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function F = calc_expdecay2(vX, vXdata)
%
% Compute exponential decay combined with x^2 function
%F = vX(1)*exp(-vX(2)*vXdata)+vX(3)*exp(vX(4)*vXdata)
% F = vX(1)+vX(2).*vXdata+vX(3).*(vXdata).^2;
F = vX(1).*exp(vX(2).*vXdata);

author: CelsoReyes


text:
function [out] = coeff2func(alpha,rho,x)
% out = coeff2func(alpha,rho,x)
% Returns the values at x of the function
% \sum_{p} alpha_p e_p(r)
% Where e_p(r) = C_p J_0(rho_p r) (see function Cp.m for a definition of
% C_p).

y = x(:)';
C = Cp(rho);

J0vals = besselj(0,rho(:)*y);
out = J0vals'*(C(:).*alpha(:));

reshape(out,size(x,1),size(x,2));

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [J, grad] = lrCostFunction(theta, X, y, lambda)
%LRCOSTFUNCTION Compute cost and gradient for logistic regression with
%regularization
%   J = LRCOSTFUNCTION(theta, X, y, lambda) computes the cost of using
%   theta as the parameter for regularized logistic regression and the
%   gradient of the cost w.r.t. to the parameters.

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost of a particular choice of theta.
%               You should set J to the cost.
%               Compute the partial derivatives and set grad to the partial
%               derivatives of the cost w.r.t. each parameter in theta
%
% Hint: The computation of the cost function and gradients can be
%       efficiently vectorized. For example, consider the computation
%
%           sigmoid(X * theta)
%
%       Each row of the resulting matrix will contain the value of the
%       prediction for that example. You can make use of this to vectorize
%       the cost function and gradient computations.
%
% Hint: When computing the gradient of the regularized cost function,
%       there're many possible vectorized solutions, but one solution
%       looks like:
%           grad = (unregularized gradient for logistic regression)
%           temp = theta;
%           temp(1) = 0;   % because we don't add anything for j = 0
%           grad = grad + YOUR_CODE_HERE (using the temp variable)
%

addpath('../ex2');
[J, grad] = costFunctionReg(theta, X, y, lambda);

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

end

author: SaveTheRbtz

repo: ml-class


text:
function I = localSum( I, dims, shape, op )
% Fast routine for box, max and min filtering.
%
% Same effect as calling 'C=convn( I, ones(dims), shape)', except more
% efficient. Computes local sums by using running sums. To get sum in
% non-overlapping windows, use shape='block'. Equivalent to doing localSum,
% and then subsampling (except more efficient). If operation op is set to
% 'max' or 'min', computes local maxes or mins instead of sums. Note, that
% when applicable convBox and convMax are significantly faster.
%
% USAGE
%  I = localSum( I, dims, [shape], [op] )
%
% INPUTS
%  I       - matrix to compute sum over
%  dims    - size of volume to compute sum over, can be scalar
%  shape   - ['full'] 'valid', 'full', 'same', or 'block'
%  op      - ['sum'] 'max', or 'min'
%
% OUTPUTS
%  C       - matrix of sums
%
% EXAMPLE - 1
%  A=rand(500,500,1); dim=25; f=ones(dim,1); shape='same'; r=20;
%  tic, for i=1:r, B = localSum(A,dim,shape); end; toc
%  tic, for i=1:r, C = conv2(conv2(A,f,shape),f',shape); end; toc
%  diff=B-C; im(diff), sum(abs(diff(:)))
%
% EXAMPLE - 2
%  load trees; I=ind2gray(X,map); figure(1); im(I);
%  I1=localSum(I,3,'block','sum'); figure(2); im(I1); title('sum')
%  I2=localSum(I,3,'block','max'); figure(3); im(I2); title('max')
%  I3=localSum(I,3,'block','min'); figure(4); im(I3); title('min')
%
% See also convBox, convMax, imShrink
%
% Piotr's Computer Vision Matlab Toolbox      Version 3.22
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<3 || isempty(shape)), shape='full'; end
if( nargin<4 || isempty(op)), op='sum'; end
assert(any(strcmp(shape,{'same','valid','full','block'})));
assert(any(strcmp(op,{'sum','max','min'})));

% standard case (use nlfiltersep)
if( nargin>=3 && strcmp(shape,'block') )
  switch op
    case 'sum', fun = @sum;
    case 'max', fun = @(x,d) max(x,[],d);
    case 'min', fun = @(x,d) min(x,[],d);
  end
else
  switch op
    case 'sum', fun = @nlfiltersep_sum;
    case 'max', fun = @nlfiltersep_max;
    case 'min', error('min not implemented');
  end
end
I = nlfiltersep( I, dims, shape, fun );

end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function update_editfield_value(src,~)
    %update_editfield_value general function to update x.Value when x.String is updated, where x is an editfield
    % use as the callback
    %
    %  uicontrol('Style','edit','string','123','Value',num2str(123),'Callback',@update_editfield_value);
    %
    %  whenever the string (text) is changed, then the value updates, too.
    %
    %  only works for single values, not arrays, or non-numeric text.
    try
        src.Value = str2double(src.String);
    catch
        warning('unable to convert "%s" to a number',src.String);
        src.Value = nan;
    end
end


text:
function obj2 = unitaryopp(f, obj1, varargin);

    % make new memory mapped data file (blank)
    % ----------------------------------------
    newFileName = mmo.getnewfilename;
    data = memmapfile(obj1.dataFile, 'writable', obj1.writable, 'format', { 'single' obj1.dimensions 'x' });
    fid = fopen(newFileName, 'w');
    if obj1.dimensions(end) == 1, obj1.dimensions(end) = []; end;
    ss.type = '()';
    ss.subs(1:length(obj1.dimensions)-1) = { ':' };
    for index = 1:obj1.dimensions(end)
        ss.subs{length(obj1.dimensions)} = index;
        tmpdata = f(subsref(data.Data.x, ss), varargin{:});
        fwrite(fid, tmpdata, 'float');
    end;
    fclose(fid);

    % copy the data
    % -------------
    obj2 = mmo(newFileName, obj1.dimensions, true, obj1.transposed);
    obj2 = updateWorkspace(obj2);

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
classdef MimCopyUID < MimGuiPlugin
    % MimCopyUID. Copy the UID of the current dataset to the clipboard
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Copy UID'
        SelectedText = 'Copy UID'
        ToolTip = 'Copy the UID of the current dataset to the clipboard'
        Category = 'Developer tools'
        Visibility = 'Dataset'
        Mode = 'Segment'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 1

        Location = 103
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            image_info = gui_app.GetImageInfo;
            uid = image_info.ImageUid;
            disp(['Current dataset UID is: ' uid]);
            clipboard('copy', uid);
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.DeveloperMode && gui_app.IsDatasetLoaded;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


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

clear all;
%close all;
clf

num=[1 1];
den=conv([1 0],[0.1 -1]);
w=logspace(-1,2,100);
[m,p]=bode(num,den,w);
subplot(2,1,1),loglog(w,m,'-',w,ones(size(w)),'-');
ylabel('Magnitude');
title('Fig. 6.30 Bode Plot for G=(s+1)/[s(s/10 -1)] (a) magnitude');
bodegrid;
subplot(2,1,2)
if p>0, p=p-360; end      % Matlab quadrant control varies in different versions
semilogx(w,p,'-',w,-180*ones(size(w)),'-');
axis([.1 100 -270 -90])

xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');
title('(b) phase');
bodegrid;

author: Sable

repo: mcbench-benchmarks


text:
vI = cat(3, 255, 102, 51) / 255;
vI = cat(3, 255, 153, 51) / 255;
vI = cat(3, 255, 128, 0) / 255;
vI = cat(3, 204, 127, 102) / 255;

vI = cat(3, 255, 0, 0) / 255;
vI = cat(3, 159, 112, 96) / 255;

% Reds, Yellows, Greens, Cyans, Blues, Magentas
vPhotoshopValues = [25; 75; 35; 65; 45; 55]; %<! Photoshop Neutral

vPhotoshopValues = randi([-200, 300], [6, 1]);
% vPhotoshopValues([1, 2]) = 50;

vPhotoshopValues = [60; 50; 50; 50; 50; 50];

vCoeffValues = (vPhotoshopValues - 50) ./ 50;
mO = ApplyBlackWhiteFilter(vI, vCoeffValues);

round(mO * 255)

% vI = cat(3, 255, 102, 51) / 255;
% round(ConvertRgbToHsl(vI) * 360)
% vI = cat(3, 240, 110, 66) / 255;
% round(ConvertRgbToHsl(vI) * 360)

author: RoyiAvital

repo: StackExchangeCodes


text:
function y = sample_node(CPD, pev)
% SAMPLE_NODE Draw a random sample from P(Xi | x(pi_i), theta_i)  (tabular)
% y = sample_node(CPD, pev)
%
% pev{i} is the value of the i'th parent (if any)

%assert(~any(isemptycell(pev)));

%CPT = CPD_to_CPT(CPD);
%sz = mysize(CPT);
sz = CPD.sizes;
nparents = length(sz)-1;
if nparents > 0
  pvals = cat(1, pev{:});
end
switch nparents
 case 0, T = CPD.CPT;
 case 1, T = CPD.CPT(pvals(1), :);
 case 2, T = CPD.CPT(pvals(1), pvals(2), :);
 case 3, T = CPD.CPT(pvals(1), pvals(2), pvals(3), :);
 case 4, T = CPD.CPT(pvals(1), pvals(2), pvals(3), pvals(4), :);
 otherwise,
  psz = sz(1:end-1);
  ssz = sz(end);
  i = subv2ind(psz, pvals(:)');
  T = reshape(CPD.CPT, [prod(psz) ssz]);
  T = T(i,:);
end

if sz(end)==2
  r = rand(1,1);
  if r > T(1)
    y = 2;
  else
    y = 1;
  end
else
  y = sample_discrete(T);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [b,c,d,x,y,z,qfac] = vox2rasToQform(vox2ras)
% [a,b,c,x,y,z,qfac] = vox2rasToQform(vox2ras)
%
% Converts a vox2ras matrix to NIFTI qform parameters.
% Note: the vox2ras should be 6 DOF. This code mostly just
% follows CH's mriToNiftiQform() in mriio.c.
%
%  hdr.pixdim(1) = qfac;
%  hdr.quatern_b = b;
%  hdr.quatern_c = c;
%  hdr.quatern_d = d;
%  hdr.qoffset_x = x;
%  hdr.qoffset_y = y;
%  hdr.qoffset_z = z;
%  hdr.qform_code = NIFTI_XFORM_SCANNER_ANAT=1;
%

%
% vox2rasToQform.m
%
% Original Author: Doug Greve
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

a = [];
if(nargin ~= 1)
  fprintf('[a,b,c,x,y,z,qfac] = vox2rasToQform(vox2ras)\n');
  return;
end

x = vox2ras(1,4);
y = vox2ras(2,4);
z = vox2ras(3,4);

d = sqrt(sum(vox2ras(:,1:3).^2));
Mdc = vox2ras(1:3,1:3) ./ repmat(d,[3 1]);
if(det(Mdc) == 0.0)
  fprintf('ERROR: vox2ras determinant is 0\n');
  return;
end

r11 = Mdc(1,1);
r21 = Mdc(2,1);
r31 = Mdc(3,1);
r12 = Mdc(1,2);
r22 = Mdc(2,2);
r32 = Mdc(3,2);
r13 = Mdc(1,3);
r23 = Mdc(2,3);
r33 = Mdc(3,3);

if(det(Mdc) > 0.0) qfac = 1.0;
else
  r13 = -r13;
  r23 = -r23;
  r33 = -r33;
  qfac = -1.0;
end

%  /* following mat44_to_quatern() */

a = r11 + r22 + r33 + 1.0;
if(a > 0.5)
  a = 0.5 * sqrt(a);
  b = 0.25 * (r32-r23) / a;
  c = 0.25 * (r13-r31) / a;
  d = 0.25 * (r21-r12) / a;
else
  xd = 1.0 + r11 - (r22+r33);
  yd = 1.0 + r22 - (r11+r33);
  zd = 1.0 + r33 - (r11+r22);
  if(xd > 1.0)
    b = 0.5 * sqrt(xd);
    c = 0.25 * (r12+r21) / b;
    d = 0.25 * (r13+r31) / b;
    a = 0.25 * (r32-r23) / b;
  elseif( yd > 1.0 )
    c = 0.5 * sqrt(yd);
    b = 0.25 * (r12+r21) / c;
    d = 0.25 * (r23+r32) / c;
    a = 0.25 * (r13-r31) / c;
  else
    d = 0.5 * sqrt(zd);
    b = 0.25 * (r13+r31) / d;
    c = 0.25 * (r23+r32) / d;
    a = 0.25 * (r21-r12) / d;
  end
  if(a < 0.0)
    a = -a;
    b = -b;
    c = -c;
    d = -d;
  end
end

return

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


function [x, y] = trans_equi2fisheye(u, v, R, M, D, fe)

[xcam, ycam, zcam] = trans_equi2cam(u, v, fe);

xr = R(1, 1) * xcam + R(1, 2) * ycam + R(1, 3) * zcam ;
yr = R(2, 1) * xcam + R(2, 2) * ycam + R(2, 3) * zcam ;
zr = R(3, 1) * xcam + R(3, 2) * ycam + R(3, 3) * zcam ;

[x, y] = trans_cam2fisheye(xr, yr, zr, M, D);

end


text:
function b = isnumber(str)

if ~ischar(str)
    b = false;
    return;
end
if isempty(str)
    b = false;
    return;
end

b = all( (str>='0' & str<='9') | (str =='e') | (str=='-' | str=='+' | str=='.' | str == ' '));

author: BUNPC


text:
% This script generates a random deployed wireless network, and
% probabilistic broadcast is simulated.
% In probabilistic broadcast, each node re-transmits the packet according
% to a probability.
% A link on probabilistic broadcast:
% http://lsewww.epfl.ch/Documents/acrobat/SCS03.pdf
% A link on DFS algorithm:
% http://en.wikipedia.org/wiki/Depth-first_search

clear;
rand('state', 5);
numOfNodes = 200;   %   number of nodes
envSize=100;        %   envsizeXenvsize environment
txRange = 15;
global floodProb;   %   the pre-defined probability to broadcast;
global savedTransmission;   % The number of saved re-transmissions;

savedTransmission = 0;
global xLocation;   % Array containing the X-coordinations of wireless nodes;
global yLocation;   % Array containing the Y-coordinations of wireless nodes;
xLocation = rand(numOfNodes,1) * envSize;
yLocation = rand(numOfNodes,1) * envSize;   %x,y coords of nodes

% Generate the adjacent matrix to represent the topology graph of the
% randomly deployed wireless networks;
distMatrix = zeros(numOfNodes,numOfNodes);
for i=1:numOfNodes
   for j=1:numOfNodes
      distMatrix(i,j)=sqrt((xLocation(i)-xLocation(j))^2 + (yLocation(i)-yLocation(j))^2); %distance between node pairs
   end
end
% If the Euclidean distance between two nodes is less than the transmission range, there
% exists a link.
global connMatrix;
connMatrix = ( distMatrix < txRange);  %binary connectivity matrix

% The broadcast will start from a sink node;
sinkNode = 1; % sink node;

% Array visited[] stores the boolean value if the broadcast has reach the node.
global visited;
visited = zeros(1, numOfNodes);

% Show the topology;
figure(1);
plot(xLocation, yLocation, '.');
text(xLocation(sinkNode), yLocation(sinkNode), 'sink');
% title(['p = ' num2str(floodProb)]);
% sink node;

visited(sinkNode) = 1;

% recursive Depth first search is adopted to traverse the whole graph;
DFS(sinkNode);
savedTransmission
return;


text:
function aic = gmm_aic(Data,w,Priors,Mu,Sigma,cov_type)
%GMM_AIC Akaike Information criterion
%
%  input ------------------------------------------------------------------
%
%   o Data:     D x N array representing N datapoints of D dimensions.
%
%   o w:        1 x N, set of weights associated with each data point.
%
%   o Priors:   1 x K array representing the prior probabilities of the
%               K GMM components.
%   o Mu:       D x K array representing the centers of the K GMM components.
%
%   o Sigma:    D x D x K array representing the covariance matrices of the
%               K GMM components.
%
%   o cov_type: string, covariance type = 'full','diag' or 'isot'
%
%  output -----------------------------------------------------------------
%
%   o aic   : (1 x 1)
%
%
[D,~]       = size(Data);
K           = length(Priors);
num_param   = K-1 + K * D;

if strcmp(cov_type,'full') == true
    num_param = num_param + K * ( D * ( D - 1)/2 );
elseif strcmp(cov_type,'diag') == true
    num_param = num_param + K * D;
elseif strcmp(cov_type,'isot') == true
    num_param = num_param + K * 1;
else
   error(['no such covariance type: ' cov_type '  only full | diag | isot ']);
end

% compute the loglikelihood of the data given the model
% loglik: (N x 1)
loglik = LogLikelihood_gmm(Data,Priors,Mu,Sigma,w);

aic = - 2 * loglik + 2 * num_param;

end

author: epfl-lasa


text:
function PTKRunTests(tests_or_folder)
    % PTKRunTests. Runs PTK unit tests
    %
    %     tests_to_run - All tests will be run if this argument is omitted
    %
    %
    %
    % Executes all the tests that can be found and reports on the sucess and
    % failures.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    if (nargin == 0) || isempty(tests_or_folder)
        tests_or_folder = PTKDirectories.GetTestSourceDirectory;
    end

    CoreRunTests(tests_or_folder);
end


text:
% COINW  Write a Mathematical File from Matlab Values
%
% THIS IS A LOW LEVEL FUNCTION - USE coinWrite() INSTEAD!
%
% coinW uses the COIN Utilities library to open and write the file.
%
%   coinW(prob,path,filetype)
%
%   Input arguments:
%       prob - problem structure (described below)
%       path - full path to the file on your PC.
%       filetype - one of the following types: [mps,lp]
%
%   Problem Structure Fields (all required):
%       H - quadratic objective matrix (sparse, may be empty)
%       f - linear objective vector
%       A - linear inequality and equality constraint matrix (sparse)
%       rl - linear constraint lower bounds
%       ru - linear constraint upper bounds
%       lb - decision variable lower bounds
%       ub - decision variable upper bounds
%       int - integer variable vector (8 bit integer / char)
%       name - problem name (string)
%       sos_type - SOS types (char array)
%       sos_index - SOS indicies (cell array of double vectors)
%       sos_weight - SOS weights (cell array of double vectors)
%       objbias - Objective bias term (optional, negated internally)
%
%
%   Copyright (C) 2011 Jonathan Currie (I2C2)


text:
%*******************************************************************
%Initialize the BBCI Toolbox:
%Hand over the directory for your (raw) data, which is in
%the Matlab variable 'MyDataDir' (e.g. 'D:\data' under Windows).

%If you have a tmp subfolder in the same directory,
%Start up the Toolbox as follows:
%"startup_bbci_toolbox('DataDir', MyDataDir);"

%If you want to choose a different folder for temporary files (e.g. 'D:\tmp\')
%Start up the Toolbox as follows:
%"startup_bbci_toolbox('DataDir', MyDataDir, 'TmpDir','/tmp/');"
%*******************************************************************

function startup_bbci_toolbox(varargin)

global BTB

% Find directory of the BBCI Toolbox and add it to the path
BBCI_DIR= fileparts(which(mfilename));
addpath(genpath(BBCI_DIR));
rmpath(genpath(fullfile(BBCI_DIR, '.git')));
% The subfolder 'external' is move to the end of the path, see below.

BBCI_PRIVATE_DIR= fullfile(fileparts(BBCI_DIR), 'bbci_private');

BTB= opt_proplistToStruct(varargin{:});

if ~isfield(BTB, 'TypeChecking'),
  BTB.TypeChecking= 1;
end
BTB= opt_setDefaults(BTB, {'DataDir'   '/home/bbci/data/'   'CHAR'});
if ~exist(BTB.DataDir, 'dir'),
  warning('Optional argument ''DataDir'' should specify an existing folder.');
end

% Guess what the location of other directories could be
BBCI_RAW_DIR= fullfile(BTB.DataDir, 'bbciRaw/');
if ~exist(BBCI_RAW_DIR, 'dir'),
  BBCI_RAW_DIR= BTB.DataDir;
end
BBCI_MAT_DIR= fullfile(BTB.DataDir, 'bbciMat/');
if ~exist(BBCI_MAT_DIR, 'dir'),
  BBCI_MAT_DIR= BTB.DataDir;
end
PYFF_DIR= fullfile(fileparts(BBCI_DIR), 'pyff', 'src');
if ~exist(PYFF_DIR, 'dir'),
  PYFF_DIR= '';
end

props= {'Dir'            BBCI_DIR          'CHAR';
        'DataDir'        ''                'CHAR';
        'RawDir'         BBCI_RAW_DIR      'CHAR';
        'MatDir'         BBCI_MAT_DIR      'CHAR';
        'PrivateDir'     BBCI_PRIVATE_DIR  'CHAR';
        'TmpDir'         ''                'CHAR';
        'PyffDir'        PYFF_DIR          'CHAR';
        'FigDir'         ''                'CHAR';
        'Tp'             struct            'STRUCT';
        'Acq'            struct            'STRUCT';
        'History'        1                 '!BOOL';
        'TypeChecking'   1                 '!BOOL'
       };
[BTB, isdefault]= opt_setDefaults(BTB, props);

% import dependencies if not done yet
if ~exist(fullfile(BTB.Dir, 'external', 'auto_import_performed'), 'file')
    inp = input(sprintf(['\nThis is the first startup of the BBCI toolbox - Welcome!\n\n', ...
        'You can now automatically download and import code from other sources.\n', ...
        'This adds additional functionalities such as ICA and SSD to the BBCI toolbox.\n',...
        'Following startups won''t ask for automatic imports. Manual imports are\n',...
        'possible with ""bbci_import_dependencies()"" at any time.\n\n',...
        'Type <yes> to download all dependencies (approx 10MB): ']), 's');
    if strcmp(inp, 'yes') || strcmp(inp, 'y')
        bbci_import_dependencies();
    end
     % write into file
    fileID = fopen(fullfile(BTB.Dir, 'external', 'auto_import_performed'),'w');
    fprintf(fileID, datestr(now()));
    fclose(fileID);
end

% Put external folders to the end of the path, such that they cannot
% overload functions of the BBCI Toolbox:
IMPORT_PATH= genpath(fullfile(BBCI_DIR, 'external'));
rmpath(IMPORT_PATH);
path(path, IMPORT_PATH);

if exist(BTB.PrivateDir, 'dir'),
  private_folders_to_add= {'utils', 'startup'};
  for kk= 1:length(private_folders_to_add),
    folder= fullfile(BTB.PrivateDir, private_folders_to_add{kk});
    if exist(folder, 'dir'),
      addpath(genpath(folder));
    end
  end
end

if isdefault.TmpDir
  BTB.TmpDir= fullfile(BTB.DataDir, 'tmp');
  if ~exist(BTB.TmpDir, 'dir'),
    fprintf('!! Default TEMP dir not existing at\n  %s\n', BTB.TmpDir);
    fprintf('!! Setting TEMP dir to ''''\n');
    BTB.TmpDir= '';
  end
end

% Information about the test person (Tp)
props= {'Dir'    BTB.TmpDir  'CHAR'
        'Code'   ''          'CHAR'
       };
BTB.Tp= opt_setDefaults(BTB.Tp, props);

% Information about data acquistion
props= {'Prefix'          'a'                     'CHAR'
        'StartLetter'     'a'                     'CHAR'
        'Dir'             BTB.TmpDir              'CHAR'
        'Geometry'        []                      'DOUBLE[1 4]'
        'TriggerFcn'      @bbci_trigger_parport   'FCN'
        'TriggerParam'    {}                      'CELL'
       };
BTB.Acq= opt_setDefaults(BTB.Acq, props);

switch computer
 case 'PCWIN'
  BTB.Acq.IoLib= which('inpout32.dll');
 case 'PCWIN64'
  BTB.Acq.IoLib= which('inpoutx64.dll');
end

if isfield(BTB.Acq, 'IoLib') && isfield(BTB.Acq, 'IoAddr'),
  if isempty(BTB.Acq.TriggerParam) && ...
      isequal(BTB.Acq.TriggerFcn, @bbci_trigger_parport),
    BTB.Acq.TriggerParam= {BTB.Acq.IoLib, BTB.Acq.IoAddr};
  end
else
  fprintf('Parport not installed. Triggers will just be printed.\n');
  BTB.Acq.TriggerFcn= @bbci_trigger_print;
end

evalin('base', 'global BTB');


text:
function [imdepthMin, imdepthMax, imdepthCol, contact, x, y, z] = ...
    getDepthRangeForDisplay(bndinfo, glabels, elabels, dt, v0)
% getMinimumDepth(bndinfo, glabels, v0)
%
% Minimum depth is the depth of the foremost occluder.  Maximum depth is
% the depth of the current region, assuming that it is touching ground at
% its lowest point.  If a region is in the foreground, it's minimum depth
% may be a function, and the minimum depth is equal to the maximum depth.
% Otherwise, the minimum and maximum depths are single (different) values.
%
% log(Z) = log(f*y_c) - log(v_0-v_i)

global DO_DISPLAY;

%% Get depth ordering

ordering = occlusion2depthOrdering(bndinfo, elabels, glabels);

%% Set some parameters

MIN_HV_DIST = 0.01;  % minimum distance between horizon and object contact
SKY_DEPTH = 2*(1 ./ MIN_HV_DIST);

%% Get pixel images for ground, vertical, and sky
labim = glabels(bndinfo.wseg);
vim = labim==2;
gim = labim==1;
sim = labim==3;

[imh, imw] = size(vim);

%% Make sure that horizon is above ground and below sky
scol = sum(sim, 2)>0;
gcol = sum(gim, 2)>0;
v0 = (1-v0) * imh;

try
minv = find(scol);  minv = minv(end-1);
maxv = find(gcol);  maxv = maxv(2);

if v0 < minv || v0 > maxv
    if minv > maxv
        v0 = maxv-1;
    else
        v0 = maxv-imh/20;
    end
end
%disp(num2str([v0 maxv minv]))
catch
end

v0 = 1 - v0/imh;
%disp(num2str(v0))

%figure(1), imagesc(cat(3, vim, gim, sim)), axis image

%% Get possible ground-vertical boundary pixels

stats = regionprops(bndinfo.wseg, 'BoundingBox', 'PixelIdxList');
bbox = vertcat(stats.BoundingBox);
idx = {stats.PixelIdxList};

% bbox = [x1 y1 x2 y2]
bbox = [bbox(:, 1) bbox(:, 2) bbox(:, 1)+bbox(:,3) bbox(:,2)+bbox(:,4)];

vfilt = [ones(3, 1) ; zeros(3, 1)];
gfilt = 1-vfilt;
boundaryim = (imfilter(double(vim), vfilt)==sum(vfilt(:))) & ...
    (imfilter(double(gim), gfilt)==sum(gfilt(:)));
boundaryim(end, :) = vim(end, :);
boundaryim = imdilate(boundaryim, ones(3, 1));

%figure(3), imagesc(boundaryim), axis image

bpts = cell(bndinfo.nseg, 1);
edges = bndinfo.edges.indices;
spLR = bndinfo.edges.spLR;
eim = zeros(bndinfo.imsize);
for k = 1:numel(edges)
    sp1 = spLR(k, 1);
    sp2 = spLR(k, 2);
    bndind = edges{k}(boundaryim(edges{k}));
    if ~isempty(bndind)
        if glabels(sp1)==2 % && glabels(sp2)==1
            bpts{sp1} = [bpts{sp1} ; bndind];
        elseif glabels(sp2)==2 % && glabels(sp1)==1
            bpts{sp2} = [bpts{sp2} ; bndind];
        end
    end
    eim(edges{k}) = 1;
end

%% Get a simple estimate of depth (trace columns upwards)

imdepthCol = zeros(bndinfo.imsize(1:2));
maxg = imh*ones(1, imw);
for y = imh:-1:1
    maxg(gim(y, :)) = y;
    imdepthCol(y, :) = maxg;
end
%imy = repmat((1:imh)', [1 imw]);
% imdepthCol = log(1 ./ max(v0 - (imh-imdepthCol)/imh, MIN_HV_DIST));
% imdepthCol(sim==1) = log(SKY_DEPTH);
for k = 1:numel(idx)
    imdepthCol(k) = (imh-median(imdepthCol(idx{k})))/imh;
end
% if DO_DISPLAY
%     figure(3), imagesc(imdepthCol), axis image
% end

%% Get the depth of each object that has a ground boundary

bx = zeros(numel(bpts), 2);
by = zeros(numel(bpts), 2);
bz = zeros(numel(bpts), 2);
bzmin = zeros(numel(bpts), 2);
bzmax = zeros(numel(bpts), 2);
isvalid = false(numel(bpts), 1);
for k = 1:numel(bpts)
    if numel(bpts{k})>5 || bbox(k, 4)+1>=imh

        [py, px] = getPerimeter(bndinfo, k, 50);
        bx(k, :) = bbox(k, [1 3]);

        cu = [];  cv = [];
        if  1 %bx(k, 2)-bx(k, 1) > imw/10  % only get line for larger objects

            keep = true(size(px));
            for tk = 1:numel(px)
                if ~any(boundaryim(py(tk) + imh*(px(tk)-1)))
                    keep(tk) = false;
                end
            end
            %px = px(keep);  py = py(keep);

            %if k==13, keyboard; end

            try
                [cu, cv] = poly2contacts(dt, px, py, bndinfo.imsize, v0, 1.7, [], 0.1);
            catch
            end
            %plot(cu, cv, 'r*')
            cu = round(cu);  cv = round(cv);
            ind = (boundaryim(cv + (cu-1)*imh)); % | (cv==max(cv));
            cu = cu(ind);  cv = cv(ind);
        end

        [cu, ind] = unique(cu);
        cv = cv(ind);

        if ~isempty(cu)
            if numel(cu) == 1
                ypos = max(py);
                imslope = 0;
            elseif numel(cu)==2
                imslope = (cv(2)-cv(1)) / (cu(2)-cu(1));
                ypos = cv(1) - imslope*cu(1);
            else
                rline = robustfit(cu, cv);
                ypos = rline(1);
                imslope = rline(2);
            end
            by(k, :) = imslope*bx(k, :) + ypos;
            bz(k, :) = (1 ./ max(v0 - (imh-by(k, :))/imh, MIN_HV_DIST));
            isvalid(k) = true;
        end

        contact{k} = [cu(:) cv(:)];

    end
end

bzmin(isvalid, :) = bz(isvalid, :);
bzmax(isvalid, :) = bz(isvalid, :);
unknown = (~isvalid) & (glabels==2);

%% Assign a maximum depth to object of unknown depth and sky
missing = unknown;

bzmax(missing, :) = Inf;

for k = find(missing)'
    bx(k, :) = bbox(k, [1 3]);
    by(k, :) = bbox(k, 4);

    tmpmax = 1 ./ max(v0 - (imh-by(k, :))/imh, MIN_HV_DIST);

    if any(tmpmax < bzmax(k, :))
        bzmax(k, :) =tmpmax;
    end
end

% each region must be in back of things that occlude it (as far back as
% possible without being in front of things that it occludes)
missing2 = missing;
changed = 1;
while changed
    lastbz = bzmax;
    %bzmax(missing2, :) = Inf;
    for k = find(missing2)'
        bzmax(k, :) =  1 ./ max(v0 - (imh-by(k, :))/imh, MIN_HV_DIST);
    end
    %changed = 0;
    for k = find(missing)'
        pos1 = find(ordering==k);

        for k2 = 1:numel(ordering)

            pos2 = find(ordering==k2);
            if pos1 < pos2  && glabels(k2)==2 && ~missing2(k2) ...
                    && (any(spLR(spLR(:, 1)==k, 2)==k2) || any(spLR(spLR(:, 2)==k, 1)==k2))
                % k is directly in front of k2
                imslope = (by(k2, 2)-by(k2, 1))./(bx(k2, 2)-bx(k2, 1));
                ypos = by(k2, 1) - imslope*bx(k2, 1);
                by(k, :) = imslope*bx(k, :) + ypos;
                tmpz = 1 ./ max(v0 - (imh-by(k, :))/imh, MIN_HV_DIST);

                if all(tmpz <= bzmax(k, :))
                    bzmax(k, :) = tmpz; %(1 ./ max(v0 - (imh-by(k, :))/imh, 0.01));
                    changed = 1;
                    missing2(k) = false;
                end

            end
        end
    end
    if all(lastbz==bzmax)
        break;
    end
end

bzmax(isinf(bzmax(:, 1)), :) = 1./ MIN_HV_DIST;

%% Assign a minimum depth to objects of unknown depth;

missing = unknown;

% assign depth as maximum depth of object that occludes it
changed = 1;
while changed
    changed = 0;
    isvalid2 = isvalid;
    for k = find(missing)'
        bx(k, :) = bbox(k, [1 3]);

        bzmin(k, :) = 0; %bzmax(k, :);

        ind = find(spLR(:, 1)==k);
        neighbors = spLR(ind, 2);
        validn = isvalid(neighbors) & (glabels(neighbors)==2);
        neighbors = neighbors(validn);
        ind = ind(validn);

        pos1 = find(ordering==k);

        maxov = 0;
        for tk2 = 1:numel(neighbors)
            k2 = neighbors(tk2);
            pos2 = find(ordering==k2);

            ov = min(bbox([k k2], 3)) - max(bbox([k k2], 1));

            if ~missing(k2) && pos2 < pos1  && glabels(k2)==2 & ...
                    (any(spLR(spLR(:, 1)==k, 2)==k2) || any(spLR(spLR(:, 2)==k, 1)==k2))

                 %(ov > maxov) || (ov > 0.25*bbox(k, 3)-bbox(k,1))
                imslope = (by(k2, 2)-by(k2, 1))./(bx(k2, 2)-bx(k2, 1));
                ypos = by(k2, 1) - imslope*bx(k2, 1);
                by(k, :) = imslope*bx(k, :) + ypos;
                tmpz = 1 ./ max(v0 - (imh-by(k, :))/imh, MIN_HV_DIST);

                if any(tmpz > bzmin(k, :))
                    bzmin(k, :) = tmpz; %(1 ./ max(v0 - (imh-by(k, :))/imh, 0.01));
                    missing(k) = false;
                    isvalid2(k) = true;
                    changed = 1;
                end

            end

        end
    end

    isvalid = isvalid2;

end
bzmin(bzmin(:, 1)==0, :) = bzmax(bzmin(:, 1)==0, :);

% check that is not behind anything that it occludes
changed = 1;
while changed
    changed = 0;
    isvalid2 = isvalid;
    for k = find(missing)'
        bx(k, :) = bbox(k, [1 3]);

        ind = find(spLR(:, 1)==k);
        neighbors = spLR(ind, 2);
        validn = isvalid(neighbors) & (glabels(neighbors)==2);
        neighbors = neighbors(validn);
        ind = ind(validn);

        pos1 = find(ordering==k);

        maxov = 0;
        for tk2 = 1:numel(neighbors)
            k2 = neighbors(tk2);
            pos2 = find(ordering==k2);

            if ~missing(k2) && pos1 < pos2  &&  glabels(k2)==2 && ...
                    (any(spLR(spLR(:, 1)==k, 2)==k2) || any(spLR(spLR(:, 2)==k, 1)==k2))

                imslope = (by(k2, 2)-by(k2, 1))./(bx(k2, 2)-bx(k2, 1));
                ypos = by(k2, 1) - imslope*bx(k2, 1);
                by(k, :) = imslope*bx(k, :) + ypos;
                tmpz = 1 ./ max(v0 - (imh-by(k, :))/imh, MIN_HV_DIST);

                if all(tmpz < bzmin(k, :))
                    bzmin(k, :) = tmpz;
                    missing(k) = false;
                    isvalid2(k) = true;
                    changed = 1;
                end

            end

        end
    end

    isvalid = isvalid2;

end

depthim = zeros([imh imw]);
tmpy = zeros(size(depthim));
tmpz = zeros(size(depthim));
for k = 1:bndinfo.nseg
    [y, x] = ind2sub(bndinfo.imsize, idx{k});
    switch glabels(k)
        case 0
            z = bzmin(k, 1) + (x-bx(k, 1)) ./ ...
                (bx(k, 2)-bx(k, 1)) .* (bzmin(k,2)-bzmin(k,1));
%             tmpy(idx{k}) = (y - (x-bx(k, 1)) ./ ...
%                 (bx(k, 2)-bx(k, 1)) .* (by(k,2)-by(k,1))) * 1.7 ./ ((x-bx(k, 1)) ./ ...
%                 (bx(k, 2)-bx(k, 1)) .* (bzmin(k,2)-bzmin(k,1)));
            tmpy(idx{k}) = -(v0 - (imh-(x-bx(k, 1)) ./ ...
                (bx(k, 2)-bx(k, 1)) .* (by(k,2)-by(k,1)))/imh).*z / 1.4 + 1.7;
            z2 = z;
        case 1
            z = 1 ./ max(v0 - (imh-y)/imh, MIN_HV_DIST);
            tmpy(idx{k}) = 0;
            z2 = z*1.4*1.7;
        case 2
            z = bzmin(k, 1) + (x-bx(k, 1)) ./ ...
                (bx(k, 2)-bx(k, 1)) .* (bzmin(k,2)-bzmin(k,1));
            ty = by(k,1) + (x-bx(k, 1)) ./ (bx(k, 2)-bx(k, 1)) .* (by(k,2)-by(k,1));

%             tmpy(idx{k}) = -(v0 - (imh-(x-bx(k, 1)) ./ ...
%                 (bx(k, 2)-bx(k, 1)) .* (by(k,2)-by(k,1)))/imh).*z / 1.4 + 1.7;

            z2 = 1.4*1.7 ./ max(v0 - (imh-ty)/imh, MIN_HV_DIST);
            %(1.7-y(idx{k})) .* 1.4 ./ (v0 - (imh-ty)./imh+1E-5);
        case 3
            z = SKY_DEPTH;
            tmpy(idx{k}) = 500;
            z2 = (1.7-tmpy(idx{k})) .* 1.4 ./ (v0 - (imh-y)./imh+1E-5);
    end
    depthim(idx{k}) = log(z);
    tmpz(idx{k}) = z2;
end

imdepthMin = depthim;

% x = u.*z./f;

if DO_DISPLAY
    figure(1), imagesc(depthim), axis image, colormap jet
end
%print -f1 -djpeg99 ../tmp/gt_mindepth.jpg

%depthim = zeros([imh imw]);
for k = find(glabels==2)'
    [y, x] = ind2sub(bndinfo.imsize, idx{k});
    zk = bzmax(k, 1) + (x-bx(k, 1)) ./ ...
        (bx(k, 2)-bx(k, 1)) .* (bzmax(k,2)-bzmax(k,1));
    depthim(idx{k}) = log(zk);
end

if DO_DISPLAY
    figure(2), imagesc(depthim), axis image, colormap jet
end
%print -f2 -djpeg99 ../tmp/gt_maxdepth.jpg

imdepthMax = depthim;

x = repmat((1:imw)/imw-0.5, [imh 1]).*tmpz/1.4;
y = 1.7 - (v0-repmat(1-(1:imh)'./imh, [1 imw]))/1.4.*tmpz;
z = tmpz;

%% gets perimeter of region
function [py, px] = getPerimeter(bndinfo, r, npts)

imsize = bndinfo.imsize;
spLR = bndinfo.edges.spLR;
spLR = [spLR ; spLR(:, [2 1])];
eind = find(spLR(:, 1)==r);

failed = false;
order = zeros(1, numel(eind));
order(1) = 1;
for k = 2:numel(order)
    laste = eind(order(k-1));
    ismem = ismember(eind, bndinfo.edges.adjacency{laste});
    next = find(ismem);
    if numel(next)~=1
        failed = true;
        break;
    end
    order(k) = next;
end

if ~failed
    eind = eind(order);
end

edges = cell(numel(eind), 1);
for k = 1:numel(eind)
    if eind(k) <= bndinfo.ne
        edges{k} = bndinfo.edges.indices{eind(k)}(2:end-1);
    else
        edges{k} = bndinfo.edges.indices{eind(k)-bndinfo.ne}(end-1:-1:2);
    end
end
[py, px] = ind2sub(imsize, vertcat(edges{:}));

try
    ind = convhull(double(px), double(py));
catch
    ind = (1:5:numel(px))';
end

if npts > numel(ind) && ~failed
    step = ceil(numel(py) / (npts-numel(ind)));
    ind2 = (1:step:numel(py))';
    ind = unique([ind ; ind2]);
end

py = double(py(ind));
px = double(px(ind));

author: Cloud-CV

repo: object-proposals


text:
function evaluate_tDCF_asvspoof19(CM_SCOREFILE, ASV_SCOREFILE)

% Set t-DCF parameters
cost_model.Pspoof       = 0.05;  % Prior probability of a spoofing attack
cost_model.Ptar         = (1 - cost_model.Pspoof) * 0.99; % Prior probability of target speaker
cost_model.Pnon         = (1 - cost_model.Pspoof) * 0.01; % Prior probability of nontarget speaker
cost_model.Cmiss_asv    = 1;     % Cost of ASV system falsely rejecting target speaker
cost_model.Cfa_asv      = 10;    % Cost of ASV system falsely accepting nontarget speaker
cost_model.Cmiss_cm     = 1;     % Cost of CM system falsely rejecting target speaker
cost_model.Cfa_cm       = 10;    % Cost of CM system falsely accepting spoof

% Load organizer's ASV scores
[asv_attacks, asv_key, asv_score] = textread(ASV_SCOREFILE, '%s %s %f');

% Load CM scores (replace these with the scores of your detectors).
[cm_utt_id, cm_attacks, cm_key, cm_score] = textread(CM_SCOREFILE, '%s %s %s %f');

% Extract target, nontarget and spoof scores from the ASV scores
tar_asv     = asv_score(strcmp(asv_key, 'target'))';
non_asv     = asv_score(strcmp(asv_key, 'nontarget'))';
spoof_asv   = asv_score(strcmp(asv_key, 'spoof'))';

% Extract bona fide (real human) and spoof scores from the CM scores
bona_cm     = cm_score(strcmp(cm_key, 'bonafide'));
spoof_cm    = cm_score(strcmp(cm_key, 'spoof'));

% Fix ASV operating point to EER threshold
[eer_asv, asv_threshold] = compute_eer(tar_asv, non_asv);

% Visualize ASV and CM scores
figure(1);
subplot(121);
[h, xi] = ksdensity(tar_asv);   plot(xi, h, 'gr'); hold on;
[h, xi] = ksdensity(non_asv);   plot(xi, h, 'k');
[h, xi] = ksdensity(spoof_asv); plot(xi, h, 're');
h = plot(asv_threshold, 0, 'mo'); set(h, 'MarkerSize', 10);
legend('Target', 'Nontarget', 'Spoof', 'EER threshold');
xlabel('ASV score'); ylabel('PDF');
title('ASV scores', 'FontSize', 12, 'FontWeight', 'bold');

subplot(122);
[h, xi] = ksdensity(bona_cm);   plot(xi, h, 'gr'); hold on;
[h, xi] = ksdensity(spoof_cm);  plot(xi, h, 're');
legend('Bona fide', 'Spoof');
xlabel('CM score'); ylabel('PDF');
title('CM scores', 'FontSize', 12, 'FontWeight', 'bold');

% Obtain the detection error rates of the ASV system
[Pfa_asv, Pmiss_asv, Pmiss_spoof_asv] = obtain_asv_error_rates(tar_asv, non_asv, spoof_asv, asv_threshold);

% Equal error rate of the countermeasure
[eer_cm, ~] = compute_eer(bona_cm, spoof_cm);

% Compute t-DCF
[tDCF_curve, CM_thresholds] = compute_tDCF(bona_cm, spoof_cm, Pfa_asv, Pmiss_asv, Pmiss_spoof_asv, cost_model, true);

% Minimum normalized t-DCF and the corresponding threshold
[min_tDCF, argmin_tDCF] = min(tDCF_curve);
min_tDCF_threshold = CM_thresholds(argmin_tDCF);

fprintf('ASV SYSTEM\n');
fprintf('     EER               \t= %5.5f %%\t\t Equal error rate (target vs. nontarget discrimination)\n',   100 * eer_asv);
fprintf('     Pfa               \t= %5.5f %%\t\t False acceptance rate of nontargets\n',   100 * Pfa_asv);
fprintf('     Pmiss             \t= %5.5f %%\t\t Miss (false rejection) rate of targets\n',   100 * Pmiss_asv);
fprintf('     1-Pmiss,spoof     \t= %5.5f %%\t Spoof false acceptance rate ("NOT miss spoof trial")\n\n', 100 * (1 - Pmiss_spoof_asv));

fprintf('CM SYSTEM\n');
fprintf('     EER               \t= %f %%\t Equal error rate from CM scores pooled across all attacks. \n\n', 100 * eer_cm);

fprintf('TANDEM\n');
fprintf('     min-tDCF          \t= %f\n',   min_tDCF);

% Plot t-DCF as function of the CM threshold.
figure(2);
plot(CM_thresholds, tDCF_curve); hold on;
h = plot(min_tDCF_threshold, min_tDCF, 'mo'); set(h, 'MarkerSize', 10);
xlabel('CM threshold (operating point)');
ylabel('Norm t-DCF');
title('Normalized t-DCF', 'FontSize', 12, 'FontWeight', 'bold');
h = line([min(CM_thresholds) max(CM_thresholds)], [1 1]); set(h, 'LineStyle', '--', 'Color', 'k');
legend('t-DCF', sprintf('min t-DCF (%5.5f)', min_tDCF), 'Arbitrarily bad CM (Norm t-DCF=1)');
axis([min(CM_thresholds) max(CM_thresholds) 0 1.5]);

% Traditional DET plot using Bosaris toolkit
plot_type = Det_Plot.make_plot_window_from_string('old');
plot_obj = Det_Plot(plot_type, 'Detection error trade-off plot of CM');
plot_obj.set_system(bona_cm, spoof_cm, sprintf('CM system (EER=%2.2f %%)', 100* eer_cm));
plot_obj.plot_steppy_det({'b','LineWidth',2},' ');
plot_obj.display_legend();

end

author: nesl


text:
% Reproduces an experiment similar to the one from figure 2 in "A
% sliding-window kernel RLS algorithm and its application to nonlinear
% channel identification".
%
% MSE performance of sliding-window kernel RLS for two different window
% sizes. Execution time: < 1 minute (Intel Pentium Core2 Duo).
%
% S. Van Vaerenbergh, J. Via, and I. Santamaria, "A Sliding-Window Kernel
% RLS Algorithm and Its Application to Nonlinear Channel Identification,"
% 2006 IEEE International Conference on Acoustics, Speech and Signal
% Processing (ICASSP), May 2006,
%
% This file is part of the Kernel Adaptive Filtering Toolbox for Matlab
% https://github.com/steven2358/kafbox/

clear
close all

%% PARAMETERS
opt.chlen = 4; % channel length, and input embedding
opt.N = 1500; % number of training data
opt.N_switch = 500; % iteration after which the channel switch occurs
opt.N_test = 100; % number of test data
opt.SNR = 20; % SNR in dB

num_sim = 10;

setups{1} = struct('M',75,'c',1E-2,'kerneltype','gauss','kernelpar',5);
setups{2} = struct('M',150,'c',1E-2,'kerneltype','gauss','kernelpar',5);

titles = {'SW-KRLS, M=75','SW-KRLS, M=150'};

%% PREPARE DATA
fprintf('Fig. 2 from "A sliding-window kernel RLS algorithm and its \n');
fprintf('application to nonlinear channel identification".\n');

% perform online learning for each algorithm
fprintf('\n')
N = opt.N;
N_switch = opt.N_switch;
num_alg = length(setups);
MSE = zeros(N,num_alg);
MSE_final = zeros(1,num_alg);
for sim_ind = 1:num_sim
    % generate data
    [X,y,y_ref,X_test,y_test] = generate_channel_switch(opt);

    fprintf('SIM %d:\n',sim_ind)
    for algo_ind = 1:num_alg
        t1 = tic;

        kaf = swkrls(setups{algo_ind});

        fprintf('%9s M=%3d: ',upper(class(kaf)),kaf.M);

        mse = zeros(N,1);
        for i=1:N
            if ~mod(i,floor(N/10)), fprintf('.'); end

            y_est = kaf.evaluate(X_test);
            if i<=N_switch
                mse(i) = mean((y_test(:,1)-y_est).^2);
            else
                mse(i) = mean((y_test(:,2)-y_est).^2);
            end

            kaf.train(X(i,:),y(i));
        end
        MSE(:,algo_ind) = MSE(:,algo_ind) + mse/num_sim;
        mse_final = mean(mse(N-500:N));

        fprintf(' %.2fs. Final MSE=%3.2fdB\n',toc(t1),...
            10*log10(mse_final))
    end
    fprintf('\n');
end
fprintf('\n');

%% OUTPUT

figure; hold all
for algo_ind=1:num_alg
    plot(10*log10(MSE(:,algo_ind)),'LineWidth',1)

    axis([0 N 5*floor(min(10*log10(MSE(:)))/5) 0]);
end
xlabel('iteration')
ylabel('MSE (dB)')
legend(titles)

repo: kafbox


text:
function Xr = pcaRandVec( U, mu, vars, k, n, hypershpere, show )
% Generate random vectors in PCA subspace.
%
% Used to generate random vectors from the subspace spanned by the first k
% principal components.  The points generated come from the gaussian
% distribution from within the subspace. Can optionally generate points on
% the subspace that are also on a hypershpere centered on the origin.  This
% may be useful if the original data points were all from a hypershpere --
% for example they were normalized via imNormalize.  Set the optional
% hypershpere flag to 1 to generate points only on the hypersphere.
%
% USAGE
%  Xr = pcaRandVec( U, mu, vars, k, n, [hypershpere], [show] )
%
% INPUTS
%  U           - returned by pca.m
%  mu          - returned by pca.m
%  vars        - returned by pca.m
%  k           - number of principal coordinates to use
%  n           - number of points to generate
%  hypershpere - [0] generate points on hypersphere (see above)
%  show        - [1] figure to use for display (no display if == 0)
%
% OUTPUTS
%  Xr          - resulting randomly generated vectors
%
% EXAMPLE
%
% See also PCA IMNORMALIZE
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.0
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( nargin<6 || isempty(hypershpere) ); hypershpere=0; end
if( nargin<7 || isempty(show) ); show=0; end

siz1 = size(mu);   nd = ndims(mu);
sizX = [siz1, n];  D = prod(siz1);

% generate random vectors inside of subspace.
C = diag( vars(1:k).^-1 );
Yr = C * randn(k,n);
Uk = U(:,1:k);
Xr = Uk * Yr;

if(hypershpere)
  % The final point must lie on hypershpere.  Need to scale each element xr
  % of Xr so that after adding it to mu the resulting vector has length D.
  % So need to find k>0 such that (k*xr + mu) has length D.  To find k
  % simply solve the quadratic equation induced by ||(k*xr + mu)||=D,
  % choosing the root such that k>0.  Note that the mean of resulting
  % vector will not necessarily be 0, but the variance will pretty much be
  % 1! Regardless, renomralize at the end.
  muv = mu(:);  muMag = dot(muv,muv);
  for i=1:n
    xr = Xr(:,i);
    rs = roots( [dot(xr,xr), 2 * dot(xr,muv), muMag-D] );
    Xr(:,i) = muv + max(rs)*xr;
  end
  Xr = reshape( Xr, sizX );
  Xr = fevalArrays( Xr, @imNormalize );
else
  % simply add the mean to reshaped Xr
  Xr = reshape( Xr, sizX );
  muRep = repmat(mu, [ones(1,nd), n ] );
  Xr = Xr + muRep;
end

% optionaly show resulting vectors
if(show && (nd==2 || nd==3)); figure(show); montage2(Xr); end

%%% Little test - see if eigenvectors induced by randomly generated vectors
%%% are the same as the original eigenvectors.  Assumes [U,mu,vars] exist.
%   Xr = pcaRandVec( U, mu, vars, 3, 100 );
%   [ Ur, mur, varsr ] = pca( Xr );
%   ind = 3;
%   Uim = reshape( U(:,ind), [ size(mu,1), size(mu,2) ]  );
%   Uimr = reshape( Ur(:,ind), [ size(mu,1), size(mu,2) ]  );
%   if( sum(abs(Uim-Uimr))>sum(abs(Uim+Uimr))) Uimr=Uimr*-1; end; %sign?
%   clf; subplot(3,1,1); im( Uim ); subplot(3,1,2); im( Uimr );
%   subplot(3,1,3); im( Uim - Uimr);
%%%


text:
function varargout = power(varargin)
%.^	     Pointwise power of a SPHEREFUN.
%
% F.^G returns a SPHEREFUN F to the scalar power G, a scalar F to the
% SPHEREFUN power G, or a SPHEREFUN F to the SPHEREFUN power G.

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

[varargout{1:nargout}] = power@separableApprox(varargin{:});

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%% Creates D-matrix to calculate trends endogenously

function D = TVEregimesDummy(startdate,regimeperiods,T)
% TVEregimesDummy(startdate,regimeperiods,X) creates is a matrix in which
% each row is a selector vector of the active regime

gdates=eval(sprintf(startdate(1:4)))+eval(sprintf(startdate(6)))*0.25:0.25:eval(sprintf(startdate(1:4)))+eval(sprintf(startdate(6)))*0.25-0.25+T/4;
NumberofRegimes=size(regimeperiods,1)+1;

% Create the dummy variable
D=zeros(T,NumberofRegimes);
for it=1:NumberofRegimes-1
    NumberofPeriods=size( regimeperiods(it,:),2)/2;
    for it2=1:NumberofPeriods
        startperiod=eval(sprintf(regimeperiods{it,2*it2-1}(1:4)))+eval(sprintf(regimeperiods{it,2*it2-1}(6)))*0.25;
        endperiod=eval(sprintf(regimeperiods{it,2*it2}(1:4)))+eval(sprintf(regimeperiods{it,2*it2}(6)))*0.25;
        startD=find(gdates==startperiod);
        endD=find(gdates==endperiod);
        D(startD:endD,it+1)=1;
    end
end
if NumberofRegimes>1
    D(:,1)=1-sum(D(:,2:end),2);
else
    D(:,1)=ones(size(D,1),1);
end

end

author: european-central-bank

repo: BEAR-toolbox


text:
function integralHist = integralHistSuperpixels(N)
%comment to be written

N = int16(N);
total_segms = max(max(N));
[height width] = size(N);

integralHist = zeros(height+1,width+1,total_segms);

for sid = 1:total_segms

    superpixelMap = not(N - sid);
    integralHist(:,:,sid) = computeIntegralImage(superpixelMap);
end

end


text:
function [acc_control_learning] = setupAccControlLearning(dataset)
%
% Authors:       Alexander Wischnewski
%
% Description:
%   function used to prepare the variables for the acceleration controller learning algorithms
%
% Input:
%   dataset: Log file with weights used for initialization (if given)

% spread of speed values
acc_control_learning.spread_vx_mps = 10:15:100;
acc_control_learning.spread_ay_req_mps2 = -36:6:36;

% design of basis functions
acc_control_learning.vx_width_mps = 20;
acc_control_learning.ay_width_mps2 = 8;

[VX, AY] = meshgrid(acc_control_learning.spread_vx_mps, acc_control_learning.spread_ay_req_mps2);
acc_control_learning.bf_vx_mps = reshape(VX, numel(VX), 1);
acc_control_learning.bf_ay_req_mps2 = reshape(AY, numel(AY), 1);

% initial weights for lateral acceleration controller learning
% will be updated below if dataset is given
acc_control_learning.w0 = zeros(length(acc_control_learning.bf_vx_mps)+1, 1);

% analyze dataset if available
if(nargin==1)
    % extract weights from log and write to initial weights
    data = load(dataset);
    weights = data.debug_slow.debug_slow_LatAcc_FFweights.Data(end, :);
    acc_control_learning.w0 = weights';
    % calculate learned feedforward part for the end of the log
    visz_vx_mps = 0:5:100;
    visz_ay_req_mps2 = -40:0.5:40;
    [Visz_X, Visz_Y] = meshgrid(visz_vx_mps, visz_ay_req_mps2);
    Visz_Z = zeros(length(visz_ay_req_mps2), length(visz_vx_mps));
    for i = 1:1:length(visz_vx_mps)
        for j = 1:1:length(visz_ay_req_mps2)
            % apply basis functions to each sample point
            dist = ((acc_control_learning.bf_vx_mps - visz_vx_mps(i))./acc_control_learning.vx_width_mps).^2 + ...
                           ((acc_control_learning.bf_ay_req_mps2 - visz_ay_req_mps2(j))./acc_control_learning.ay_width_mps2).^2;
            Visz_Z(j, i) = weights*[exp(-dist); visz_ay_req_mps2(j)];
        end
    end
    figure;
    surface(Visz_X, Visz_Y, Visz_Z);
    xlabel('Velocity in mps');
    ylabel('Lateral acc. target in mps2');
    zlabel('Corrective steering angle in rad');
end

author: TUMFTM