text:
%% DEMO_febio_0007_sphere_sliding
% Below is a demonstration for:
%
% * Building geometry for a slab with hexahedral elements, and a
% triangulated sphere.
% * Defining the boundary conditions
% * Coding the febio structure
% * Running the model
% * Importing and visualizing the displacement results

%% Keywords
%
% * febio_spec version 4.0
% * febio, FEBio
% * indentation
% * contact, sliding, sticky, friction
% * rigid body constraints
% * hexahedral elements, hex8
% * triangular elements, tri3
% * slab, block, rectangular
% * sphere
% * static, solid
% * hyperelastic, Ogden
% * displacement logfile
% * stress logfile

%%

clear; close all; clc;

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

%% Control parameters

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

% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=[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

%Specifying dimensions and number of elements for slab
sampleHeight=5; %Height
sampleWidth=sampleHeight*4; %Width
sampleThickness=sampleHeight*1; %Thickness
pointSpacings=[1 1 1]; %Desired point spacing between nodes
numElementsWidth=round(sampleWidth/pointSpacings(1)); %Number of elemens in dir 1
numElementsThickness=round(sampleThickness/pointSpacings(2)); %Number of elemens in dir 2
numElementsHeight=round(sampleHeight/pointSpacings(3)); %Number of elemens in dir 3

%Sphere parameters
numRefineStepsSphere=3;
sphereRadius=sampleHeight/2;
sphereMeshType='tri3'; % tri3 or quad4

%Define applied displacement
sphereIndentationDisplacement=sphereRadius;
sphereSlideDisplacement=sampleWidth-(sphereRadius*2);

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

% FEA control settings
numTimeSteps1=10; %Number of time steps desired
dtmin1=(1/numTimeSteps1)/100; %Minimum time step size
dtmax1=1/numTimeSteps1; %Maximum time step size

numTimeSteps2=10; %Number of time steps desired
dtmin2=(1/numTimeSteps2)/100; %Minimum time step size
dtmax2=1/numTimeSteps2; %Maximum time step size

max_refs=75; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=20; %Optimum number of iterations
max_retries=5; %Maximum number of retires
symmetric_stiffness=0;

runMode='external';% 'internal' or 'external'

%Contact parameters
contactInitialOffset=0.1;
contactPenalty=20;
laugon=0;
minaug=1;
maxaug=10;
fric_coeff=0.3;

%% Creating model geometry and mesh
% A box is created with tri-linear hexahedral (hex8) elements using the
% |hexMeshBox| function. The function offers the boundary faces with
% seperate labels for the top, bottom, left, right, front, and back sides.
% As such these can be used to define boundary conditions on the exterior.

% Create a box with hexahedral elements
beamDimensions=[sampleWidth sampleThickness sampleHeight]; %Dimensions
beamElementNumbers=[numElementsWidth numElementsThickness numElementsHeight]; %Number of elements
outputStructType=2; %A structure compatible with mesh view
[meshStruct]=hexMeshBox(beamDimensions,beamElementNumbers,outputStructType);

%Access elements, nodes, and faces from the structure
E1=meshStruct.elements; %The elements
V1=meshStruct.nodes; %The nodes (vertices)
Fb1=meshStruct.facesBoundary; %The boundary faces
Cb1=meshStruct.boundaryMarker; %The "colors" or labels for the boundary faces
elementMaterialIndices=ones(size(E1,1),1); %Element material indices

%% Creating triangulated sphere surface model

switch sphereMeshType
    case 'tri3'
        [E2,V2,~]=geoSphere(numRefineStepsSphere,sphereRadius);
    case 'quad4'
        [E2,V2]=quadSphere(numRefineStepsSphere,sphereRadius);
end

%Offset indentor
minV2=min(V2,[],1);
minV1=min(V1,[],1);

V2(:,1)=V2(:,1)-minV2(1)+minV1(1);

V2(:,3)=V2(:,3)-minV2(3)+(sampleHeight/2)+contactInitialOffset;

center_of_mass=mean(V2,1);

%%
% Plotting model boundary surfaces and a cut view

hFig=cFigure;

subplot(1,2,1); hold on;
title('Model boundary surfaces and labels','FontSize',fontSize);
gpatch(Fb1,V1,Cb1,'k',faceAlpha1);
gpatch(E2,V2,'kw','k',faceAlpha1);
colormap(gjet(6)); icolorbar;
axisGeom(gca,fontSize);

hs=subplot(1,2,2); hold on;
title('Cut view of solid mesh','FontSize',fontSize);
optionStruct.hFig=[hFig hs];
gpatch(E2,V2,'kw','k',1);
meshView(meshStruct,optionStruct);
axisGeom(gca,fontSize);

drawnow;

%% Joining node sets
V=[V1;V2;]; %Combined node sets
E2=E2+size(V1,1); %Fixed element indices

%%
% Plotting joined geometry
cFigure;
title('Joined node sets','FontSize',fontSize);
xlabel('X','FontSize',fontSize); ylabel('Y','FontSize',fontSize); zlabel('Z','FontSize',fontSize);
hold on;
gpatch(Fb1,V,Cb1,'k',faceAlpha1);
gpatch(E2,V,'kw','k',faceAlpha1);
colormap(gjet(6)); icolorbar;
axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Define contact surfaces

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

% The deformable slave surface of the slab
logicContactSurf1=Cb1==6;
F_contact_primary=Fb1(logicContactSurf1,:);

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

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

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

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

%% Define boundary conditions

%Supported nodes
Fr=Fb1(Cb1==5,:);
bcSupportList=unique(Fr(:));

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

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

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

legend(hl2,{'Rigid body sphere','BC support'});

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

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

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

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

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

%Control section
stepStruct1.Control.analysis='STATIC';
stepStruct1.Control.time_steps=numTimeSteps1;
stepStruct1.Control.step_size=1/numTimeSteps1;
stepStruct1.Control.solver.max_refs=max_refs;
stepStruct1.Control.solver.qn_method.max_ups=max_ups;
stepStruct1.Control.solver.symmetric_stiffness=symmetric_stiffness;
stepStruct1.Control.time_stepper.dtmin=dtmin1;
stepStruct1.Control.time_stepper.dtmax=dtmax1;
stepStruct1.Control.time_stepper.max_retries=max_retries;
stepStruct1.Control.time_stepper.opt_iter=opt_iter;

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

stepStruct2.Control.analysis='STATIC';
stepStruct2.Control.time_steps=numTimeSteps2;
stepStruct2.Control.step_size=1/numTimeSteps2;
stepStruct2.Control.solver.max_refs=max_refs;
stepStruct2.Control.solver.qn_method.max_ups=max_ups;
stepStruct2.Control.solver.symmetric_stiffness=symmetric_stiffness;
stepStruct2.Control.time_stepper.dtmin=dtmin2;
stepStruct2.Control.time_stepper.dtmax=dtmax2;
stepStruct2.Control.time_stepper.max_retries=max_retries;
stepStruct2.Control.time_stepper.opt_iter=opt_iter;

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

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

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

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

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

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

% -> Elements
partName1='Part1';
febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of this part
febio_spec.Mesh.Elements{1}.ATTR.type='hex8'; %Element type
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E1,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E1; %The element matrix

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

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

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

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

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

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

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

%Boundary condition section
% -> Fix boundary conditions
febio_spec.Boundary.bc{1}.ATTR.name='zero_displacement_x';
febio_spec.Boundary.bc{1}.ATTR.type='zero displacement';
febio_spec.Boundary.bc{1}.ATTR.node_set=nodeSetName1;
febio_spec.Boundary.bc{1}.x_dof=1;
febio_spec.Boundary.bc{1}.y_dof=1;
febio_spec.Boundary.bc{1}.z_dof=1;

%Rigid section
% -> Prescribed rigid body boundary conditions
febio_spec.Step.step{1}.Rigid.rigid_bc{1}.ATTR.name='RigidFix_1';
febio_spec.Step.step{1}.Rigid.rigid_bc{1}.ATTR.type='rigid_fixed';
febio_spec.Step.step{1}.Rigid.rigid_bc{1}.rb=2;
febio_spec.Step.step{1}.Rigid.rigid_bc{1}.Rx_dof=1;
febio_spec.Step.step{1}.Rigid.rigid_bc{1}.Ry_dof=1;
febio_spec.Step.step{1}.Rigid.rigid_bc{1}.Ru_dof=1;
febio_spec.Step.step{1}.Rigid.rigid_bc{1}.Rv_dof=1;
febio_spec.Step.step{1}.Rigid.rigid_bc{1}.Rw_dof=1;

febio_spec.Step.step{1}.Rigid.rigid_bc{2}.ATTR.name='RigidPrescribe';
febio_spec.Step.step{1}.Rigid.rigid_bc{2}.ATTR.type='rigid_displacement';
febio_spec.Step.step{1}.Rigid.rigid_bc{2}.rb=2;
febio_spec.Step.step{1}.Rigid.rigid_bc{2}.dof='z';
febio_spec.Step.step{1}.Rigid.rigid_bc{2}.value.ATTR.lc=1;
febio_spec.Step.step{1}.Rigid.rigid_bc{2}.value.VAL=-(sphereIndentationDisplacement+contactInitialOffset);

febio_spec.Step.step{2}.Rigid.rigid_bc{1}.ATTR.name='RigidFix_1';
febio_spec.Step.step{2}.Rigid.rigid_bc{1}.ATTR.type='rigid_fixed';
febio_spec.Step.step{2}.Rigid.rigid_bc{1}.Ry_dof=1;
febio_spec.Step.step{2}.Rigid.rigid_bc{1}.Rz_dof=1;
febio_spec.Step.step{2}.Rigid.rigid_bc{1}.Ru_dof=1;
febio_spec.Step.step{2}.Rigid.rigid_bc{1}.Rv_dof=1;
febio_spec.Step.step{2}.Rigid.rigid_bc{1}.Rw_dof=1;

febio_spec.Step.step{2}.Rigid.rigid_bc{2}.ATTR.name='RigidPrescribe';
febio_spec.Step.step{2}.Rigid.rigid_bc{2}.ATTR.type='rigid_displacement';
febio_spec.Step.step{2}.Rigid.rigid_bc{2}.rb=2;
febio_spec.Step.step{2}.Rigid.rigid_bc{2}.dof='x';
febio_spec.Step.step{2}.Rigid.rigid_bc{2}.value.ATTR.lc=2;
febio_spec.Step.step{2}.Rigid.rigid_bc{2}.value.VAL=sphereSlideDisplacement;
febio_spec.Step.step{2}.Rigid.rigid_bc{2}.relative=1;

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

%LoadData section
% -> load_controller
febio_spec.LoadData.load_controller{1}.ATTR.name='LC_1';
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}.extend='CONSTANT';
febio_spec.LoadData.load_controller{1}.points.pt.VAL=[0 0; 1 1; 2 1];

febio_spec.LoadData.load_controller{2}.ATTR.name='LC_2';
febio_spec.LoadData.load_controller{2}.ATTR.id=2;
febio_spec.LoadData.load_controller{2}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{2}.interpolate='LINEAR';
%febio_spec.LoadData.load_controller{2}.extend='CONSTANT';
febio_spec.LoadData.load_controller{2}.points.pt.VAL=[0 0; 1 0; 2 1];

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

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

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

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

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

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

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

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

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

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

%% Import FEBio results

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

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

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

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

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

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

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

    [CV]=faceToVertexMeasure(E1,V,E_stress_mat(:,:,end));

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    title('$\sigma_{3}$ [MPa]','Interpreter','Latex')
    hp=gpatch(Fb1,V_DEF(:,:,end),CV,'k',1); %Add graphics object to animate
    hp.Marker='.';
    hp.MarkerSize=markerSize2;
    hp.FaceColor='interp';

    hp2=gpatch(E2,V_DEF(:,:,end),'w','none',0.5); %Add graphics object to animate

    axisGeom(gca,fontSize);
    colormap(flipud(gjet(250))); colorbar;
    caxis([min(E_stress_mat(:)) max(E_stress_mat(:))]);
    axis(axisLim(V_DEF)); %Set axis limits statically
    camlight headlight;

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

        [CV]=faceToVertexMeasure(E1,V,E_stress_mat(:,:,qt));

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

end

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [ result ] = plot_neb_spline( directory  )
%PLOT_NEB_SPLINE Plot the energy along a NEB path.
%   result = PLOT_NEB_SPLINE() plots the energy as a function of
%   hyperdistance using the results of a NEB calculation. Each image
%   directory (01, 02, ... NIMG) must contain an OUTCAR and CONTCAR file.
%   The end points (00 and NIMG+1) must contain an OUTCAR and POSCAR file.
%   The energies are fit with a spline using NEB_SPLINE.
%
%   See also NEB_ENERGIES, NEB_SPLINE, HYPERDISTANCE.

    startDir = pwd();
    if nargin > 0
       cd(directory);
    end
    result = 0;

    [ energy0 x0 ] = neb_energies();
    [ energy x ] = neb_spline();

    figure
    hold on
    plot(x0, energy0, 'ok')
    plot(x, energy, 'k')

    xlabel(['Hyperdistance (' char(197) ')'])
    ylabel('Energy (eV)')

    cd(startDir);
end

author: Sable

repo: mcbench-benchmarks


text:
function [data, xAxis, yAxis, misc] = impload(filename)
% Reads PerkinElmer vis image files files.
% This version is compatible with single image files.
%
% [data, xAxis, yAxis] = visload(filename):
%   data: length(x) * length(y) 2D array
%   xAxis: vector for horizontal axis (e.g. micrometers)
%   yAxis: vector for vertical axis (e.g. micrometers)

% Copyright (C)2007 PerkinElmer Life and Analytical Sciences
% Stephen Westlake, Seer Green
%
% History
% 2007-04-29 SW     Initial version

% Read the bitmap
data = imread(filename, 'bmp');;

fid = fopen(filename,'r');
if fid == -1
    error('Cannot open the file.');
    return
end

% Fixed file header of signature and description
signature = setstr(fread(fid, 2, 'uchar')');
if ~strcmp(signature, 'BM')
    error('This is not a PerkinElmer vis file.');
    return
end
fseek(fid, 0, 'bof');

% Read the trailer
fseek(fid, -4 * 8, 'eof');

% YAxis is reversed to suit image()
xAxis(1) = fread(fid, 1, 'double');     % x1
yAxis(2) = fread(fid, 1, 'double');     % y1
xAxis(2) = fread(fid, 1, 'double');     % x2
yAxis(1) = fread(fid, 1, 'double');     % y2

fclose(fid);


text:
Network torchvision.models.mnasnet {
Layer Conv2d-1 {
Type: CONV
Stride { X: 2, Y: 2 }
Dimensions { K: 32, C: 3, R: 3, S: 3, Y: 224, X: 224 }
}
Layer Conv2d-2 {
Type: DSCONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1, C: 32, R: 3, S: 3, Y: 112, X: 112 }
}
Layer Conv2d-3 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 16, C: 32, R: 1, S: 1, Y: 112, X: 112 }
}
Layer Conv2d-4 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 48, C: 16, R: 1, S: 1, Y: 112, X: 112 }
}
Layer Conv2d-5 {
Type: DSCONV
Stride { X: 2, Y: 2 }
Dimensions { K: 1, C: 48, R: 3, S: 3, Y: 112, X: 112 }
}
Layer Conv2d-6 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 24, C: 48, R: 1, S: 1, Y: 56, X: 56 }
}
Layer Conv2d-7 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 72, C: 24, R: 1, S: 1, Y: 56, X: 56 }
}
Layer Conv2d-8 {
Type: DSCONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1, C: 72, R: 3, S: 3, Y: 56, X: 56 }
}
Layer Conv2d-9 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 24, C: 72, R: 1, S: 1, Y: 56, X: 56 }
}
Layer Conv2d-10 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 72, C: 24, R: 1, S: 1, Y: 56, X: 56 }
}
Layer Conv2d-11 {
Type: DSCONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1, C: 72, R: 3, S: 3, Y: 56, X: 56 }
}
Layer Conv2d-12 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 24, C: 72, R: 1, S: 1, Y: 56, X: 56 }
}
Layer Conv2d-13 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 72, C: 24, R: 1, S: 1, Y: 56, X: 56 }
}
Layer Conv2d-14 {
Type: DSCONV
Stride { X: 2, Y: 2 }
Dimensions { K: 1, C: 72, R: 5, S: 5, Y: 56, X: 56 }
}
Layer Conv2d-15 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 40, C: 72, R: 1, S: 1, Y: 28, X: 28 }
}
Layer Conv2d-16 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 120, C: 40, R: 1, S: 1, Y: 28, X: 28 }
}
Layer Conv2d-17 {
Type: DSCONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1, C: 120, R: 5, S: 5, Y: 28, X: 28 }
}
Layer Conv2d-18 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 40, C: 120, R: 1, S: 1, Y: 28, X: 28 }
}
Layer Conv2d-19 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 120, C: 40, R: 1, S: 1, Y: 28, X: 28 }
}
Layer Conv2d-20 {
Type: DSCONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1, C: 120, R: 5, S: 5, Y: 28, X: 28 }
}
Layer Conv2d-21 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 40, C: 120, R: 1, S: 1, Y: 28, X: 28 }
}
Layer Conv2d-22 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 240, C: 40, R: 1, S: 1, Y: 28, X: 28 }
}
Layer Conv2d-23 {
Type: DSCONV
Stride { X: 2, Y: 2 }
Dimensions { K: 1, C: 240, R: 5, S: 5, Y: 28, X: 28 }
}
Layer Conv2d-24 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 80, C: 240, R: 1, S: 1, Y: 14, X: 14 }
}
Layer Conv2d-25 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 480, C: 80, R: 1, S: 1, Y: 14, X: 14 }
}
Layer Conv2d-26 {
Type: DSCONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1, C: 480, R: 5, S: 5, Y: 14, X: 14 }
}
Layer Conv2d-27 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 80, C: 480, R: 1, S: 1, Y: 14, X: 14 }
}
Layer Conv2d-28 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 480, C: 80, R: 1, S: 1, Y: 14, X: 14 }
}
Layer Conv2d-29 {
Type: DSCONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1, C: 480, R: 5, S: 5, Y: 14, X: 14 }
}
Layer Conv2d-30 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 80, C: 480, R: 1, S: 1, Y: 14, X: 14 }
}
Layer Conv2d-31 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 480, C: 80, R: 1, S: 1, Y: 14, X: 14 }
}
Layer Conv2d-32 {
Type: DSCONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1, C: 480, R: 3, S: 3, Y: 14, X: 14 }
}
Layer Conv2d-33 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 96, C: 480, R: 1, S: 1, Y: 14, X: 14 }
}
Layer Conv2d-34 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 576, C: 96, R: 1, S: 1, Y: 14, X: 14 }
}
Layer Conv2d-35 {
Type: DSCONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1, C: 576, R: 3, S: 3, Y: 14, X: 14 }
}
Layer Conv2d-36 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 96, C: 576, R: 1, S: 1, Y: 14, X: 14 }
}
Layer Conv2d-37 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 576, C: 96, R: 1, S: 1, Y: 14, X: 14 }
}
Layer Conv2d-38 {
Type: DSCONV
Stride { X: 2, Y: 2 }
Dimensions { K: 1, C: 576, R: 5, S: 5, Y: 14, X: 14 }
}
Layer Conv2d-39 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 192, C: 576, R: 1, S: 1, Y: 7, X: 7 }
}
Layer Conv2d-40 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1152, C: 192, R: 1, S: 1, Y: 7, X: 7 }
}
Layer Conv2d-41 {
Type: DSCONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1, C: 1152, R: 5, S: 5, Y: 7, X: 7 }
}
Layer Conv2d-42 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 192, C: 1152, R: 1, S: 1, Y: 7, X: 7 }
}
Layer Conv2d-43 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1152, C: 192, R: 1, S: 1, Y: 7, X: 7 }
}
Layer Conv2d-44 {
Type: DSCONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1, C: 1152, R: 5, S: 5, Y: 7, X: 7 }
}
Layer Conv2d-45 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 192, C: 1152, R: 1, S: 1, Y: 7, X: 7 }
}
Layer Conv2d-46 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1152, C: 192, R: 1, S: 1, Y: 7, X: 7 }
}
Layer Conv2d-47 {
Type: DSCONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1, C: 1152, R: 5, S: 5, Y: 7, X: 7 }
}
Layer Conv2d-48 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 192, C: 1152, R: 1, S: 1, Y: 7, X: 7 }
}
Layer Conv2d-49 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1152, C: 192, R: 1, S: 1, Y: 7, X: 7 }
}
Layer Conv2d-50 {
Type: DSCONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1, C: 1152, R: 3, S: 3, Y: 7, X: 7 }
}
Layer Conv2d-51 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 320, C: 1152, R: 1, S: 1, Y: 7, X: 7 }
}
Layer Conv2d-52 {
Type: CONV
Stride { X: 1, Y: 1 }
Dimensions { K: 1280, C: 320, R: 1, S: 1, Y: 7, X: 7 }
}
Layer Linear-53 {
Type: CONV
Dimensions { K: 1000, C: 1280, R: 1, S: 1, Y: 1, X: 1 }
}
}

author: maestro-project

repo: maestro


text:
function a = polya_fit_m(data, a, weight)
% POLYA_FIT_M   Maximum-likelihood Dirichlet-multinomial (Polya) mean.
%
% POLYA_FIT_M(data,a) returns the MLE (a) for the matrix DATA,
% subject to a constraint on sum(A).
% Each row of DATA is a histogram of counts.
% A is a row vector providing the initial guess for the parameters.
%
% POLYA_FIT_M(data,a,weight) returns the MLE where each histogram is weighted.
% WEIGHT is a column vector of numbers in [0,1] (default all ones).
%
% A is decomposed into S*M, where M is a vector such that sum(M)=1,
% and only M is changed by this function.  In other words, sum(A)
% is unchanged by this function.
%
% The algorithm is a generalized Newton iteration, described in
% "Estimating a Dirichlet distribution" by T. Minka.

% Written by Tom Minka

s = sum(a);
m = a/s;
[N,K] = size(data);

use_weight = (nargin > 2);
row = (rows(a) == 1);

for iter = 1:20
  old_m = m;
  if row
    a = s*m;
    for k = 1:length(m)
      dk = data(:,k);
      vdk = a(k)*di_pochhammer(a(k), dk);
      if use_weight
	vdk = vdk .* weight;
      end
      m(k) = sum(vdk);
    end
  else
    a = repmat(s*m, 1, N);
    vdata = a.*di_pochhammer(a, data);
    if use_weight
      vdata = vdata .* repmat(weight, rows(vdata), 1);
    end
    m = row_sum(vdata);
  end
  m = m ./ sum(m);
  if max(abs(m - old_m)) < 1e-4
    break
  end
end
a = s*m;

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function [ dData, iLC, iEvalInfoMask, drecksMDH, iLCPositions ] = fMeas_main( sFilename, lReadData, measReadPara )
%FMEAS_MAIN main function for extracting Data, Loop counters, evalInfoMask
% and drecksMDH
%
% (c) Christian Wuerslin, Thomas Kuestner
% ---------------------------------------------------------------------

if(nargin < 1)
    [sFile, sPath] = uigetfile({'*.dat', 'ADC meas files (*.dat)'},'Select measurement file', pwd);
    sFilename = [sPath,filesep,sFile];
end

if(nargin < 2)
    lReadData = true;
end

% input read parameter
if(nargin < 3)
    % Set, Seg, Smp
    measReadPara = {[0 inf], [0 64999], [21 inf]}; % for CS_Trufi and CS_Flash, CS_Flash_Retro, CS_EPI
end

[sPath, sName, sExt] = fileparts(sFilename);

if(~exist(fullfile(sPath,[sName,'.mat']),'file'))
    fprintf('Creating look-up table for measdata file\n');
    fMeasCreateLUT_main(sFilename);
end

% load data
if(lReadData)
    [dData, iLC, iEvalInfoMask] = fMeasRead(fullfile(sPath,[sName,sExt]), 'Set', measReadPara{1}, 'Seg', measReadPara{2}, 'Smp', measReadPara{3});
else
    dData = [];
    iLC = [];
    iEvalInfoMask = [];
end
matfile = whos('-file', fullfile(sPath,[sName,'.mat']));
if(ismember('iLCDrecks', {matfile.name}) && ismember('iLCPositions', {matfile.name}))
    % drecksMDH 1.0
    [drecksMDH, iLCPositions] = fMeasReadDrecksMDH(fullfile(sPath,[sName,sExt]));

elseif(ismember('SDrecksMDH', {matfile.name}))
    % drecksMDH 2.0
    load(fullfile(sPath,[sName,'.mat']), 'SDrecksMDH');
    drecksMDH = SDrecksMDH;
    iLCPositions = [];
else
    % no drecksMDH at all
    drecksMDH = [];
    iLCPositions = [];
end

end

author: thomaskuestner


text:
function [z,a,it,ord,s,fct] = backcor(n,y,ord,s,fct)

% BACKCOR   Background estimation by minimizing a non-quadratic cost function.
%
%   [EST,COEFS,IT] = BACKCOR(N,Y,ORDER,THRESHOLD,FUNCTION) computes and estimation EST
%   of the background (aka. baseline) in a spectroscopic signal Y with wavelength N.
%   The background is estimated by a polynomial with order ORDER using a cost-function
%   FUNCTION with parameter THRESHOLD. FUNCTION can have the four following values:
%       'sh'  - symmetric Huber function :  f(x) = { x^2  if abs(x) < THRESHOLD,
%                                                  { 2*THRESHOLD*abs(x)-THRESHOLD^2  otherwise.
%       'ah'  - asymmetric Huber function :  f(x) = { x^2  if x < THRESHOLD,
%                                                   { 2*THRESHOLD*x-THRESHOLD^2  otherwise.
%       'stq' - symmetric truncated quadratic :  f(x) = { x^2  if abs(x) < THRESHOLD,
%                                                       { THRESHOLD^2  otherwise.
%       'atq' - asymmetric truncated quadratic :  f(x) = { x^2  if x < THRESHOLD,
%                                                        { THRESHOLD^2  otherwise.
%   COEFS returns the ORDER+1 vector of the estimated polynomial coefficients
%   (computed with n sorted and bounded in [-1,1] and y bounded in [0,1]).
%   IT returns the number of iterations.
%
%   [EST,COEFS,IT] = BACKCOR(N,Y) does the same, but run a graphical user interface
%   to help setting ORDER, THRESHOLD and FCT.
%
% For more informations, see:
% - V. Mazet, C. Carteret, D. Brie, J. Idier, B. Humbert. Chemom. Intell. Lab. Syst. 76 (2), 2005.
% - V. Mazet, D. Brie, J. Idier. Proceedings of EUSIPCO, pp. 305-308, 2004.
% - V. Mazet. PhD Thesis, University Henri Poincaré Nancy 1, 2005.
%
% 22-June-2004, Revised 19-June-2006, Revised 30-April-2010,
% Revised 12-November-2012 (thanks E.H.M. Ferreira!)
% Comments and questions to: vincent.mazet@unistra.fr.

% Check arguments
if nargin < 2, error('backcor:NotEnoughInputArguments','Not enough input arguments'); end;
if nargin < 5, [z,a,it,ord,s,fct] = backcorgui(n,y); return; end; % delete this line if you do not need GUI
if ~isequal(fct,'sh') && ~isequal(fct,'ah') && ~isequal(fct,'stq') && ~isequal(fct,'atq'),
    error('backcor:UnknownFunction','Unknown function.');
end;

% Rescaling
N = length(n);
[n,i] = sort(n);
y = y(i);
maxy = max(y);
dely = (maxy-min(y))/2;
n = 2 * (n(:)-n(N)) / (n(N)-n(1)) + 1;
y = (y(:)-maxy)/dely + 1;

% Vandermonde matrix
p = 0:ord;
T = repmat(n,1,ord+1) .^ repmat(p,N,1);
Tinv = pinv(T'*T) * T';

% Initialisation (least-squares estimation)
a = Tinv*y;
z = T*a;

% Other variables
alpha = 0.99 * 1/2;     % Scale parameter alpha
it = 0;                 % Iteration number
zp = ones(N,1);         % Previous estimation

% LEGEND
while sum((z-zp).^2)/sum(zp.^2) > 1e-9,

    it = it + 1;        % Iteration number
    zp = z;             % Previous estimation
    res = y - z;        % Residual

    % Estimate d
    if isequal(fct,'sh'),
        d = (res*(2*alpha-1)) .* (abs(res)<s) + (-alpha*2*s-res) .* (res<=-s) + (alpha*2*s-res) .* (res>=s);
    elseif isequal(fct,'ah'),
        d = (res*(2*alpha-1)) .* (res<s) + (alpha*2*s-res) .* (res>=s);
    elseif isequal(fct,'stq'),
        d = (res*(2*alpha-1)) .* (abs(res)<s) - res .* (abs(res)>=s);
    elseif isequal(fct,'atq'),
        d = (res*(2*alpha-1)) .* (res<s) - res .* (res>=s);
    end;

    % Estimate z
    a = Tinv * (y+d);   % Polynomial coefficients a
    z = T*a;            % Polynomial

end;

% Rescaling
[~,j] = sort(i);
z = (z(j)-1)*dely + maxy;

    a(1) = a(1)-1;
    a = a*dely;% + maxy;

end

% delete lines below if you do not need GUI

function [z,a,it,ord,s,fct] = backcorgui(n,y)

% BACKCORGUI   Graphical User Interface for background estimation.

% Initialization
z = [];
a = [];
it = [];
ord = [];
s = [];
fct = [];

order = 4;
threshold = 0.01;
costfunction = 'atq';

% Main window
hwin = figure('Visible','off','Position',[0 0 750 400],'NumberTitle','off','Name','Background Correction',...
    'MenuBar','none','Toolbar','figure','Resize','on','ResizeFcn',{@WinResizeFcn});
bgclr = get(hwin,'Color');

% Axes
haxes = axes('Units','pixels');

% Buttons OK & Cancel
hok = uicontrol('Style','pushbutton','String','OK','Position',[600,40,80,25],'Callback',{@OKFcn},'BackgroundColor',bgclr);
hcancel = uicontrol('Style','pushbutton','String','Cancel','Position',[510,40,80,25],'Callback',{@CancelFcn},'BackgroundColor',bgclr);

% Cost functions menu
hfctlbl = uicontrol('Style','text','String','Cost function:','HorizontalAlignment','left','BackgroundColor',bgclr);
hfct = uicontrol('Style','popupmenu','Value',4,'BackgroundColor','white','Callback',{@CostFunctionFcn},...
    'String',{'Symmetric Huber function','Asymmetric Huber function','Symmetric truncated quadratic','Asymmetric truncated quadratic'});

% Threshold text
hthresholdlbl = uicontrol('Style','text','String','Threshold:','HorizontalAlignment','left','BackgroundColor',bgclr);
hthreshold = uicontrol('Style','edit','String',num2str(threshold),'BackgroundColor','white','Callback',{@ThresholdFcn});

% Order slider
horderlbl = uicontrol('Style','text','String','Polynomial order:','HorizontalAlignment','left','BackgroundColor',bgclr);
horder = uicontrol('Style','slider','SliderStep',[0.5 0.5],'Min',0,'Max',10,'Value',order,'SliderStep',[0.1 0.1],'Callback',{@OrderFcn});
horderval = uicontrol('Style','text','String',num2str(order),'BackgroundColor',bgclr);

% Move the GUI to the center of the screen
movegui(hwin,'center');

% Plot a first estimation
[ztmp,atmp,ittmp,order,threshold,costfunction] = compute(n,y,order,threshold,costfunction);

% Make the GUI visible
set(hwin,'Visible','on');

% Callback functions

    function CancelFcn(source,eventdata)
        % Just close the window
        uiresume(gcbf);
        close(hwin);
    end

    function OKFcn(source,eventdata)
        % Return the current estimation and close the window
        z = ztmp;
        a = atmp;
        it = ittmp;
        ord = order;
        s = threshold;
        fct = costfunction;
        uiresume(gcbf);
        close(hwin);
    end

    function CostFunctionFcn(source,eventdata)
        % Change cost function
        cf = get(hfct,'Value');
        if cf == 1,
            costfunction = 'sh';
        elseif cf == 2,
            costfunction = 'ah';
        elseif cf == 3,
            costfunction = 'stq';
        elseif cf == 4,
            costfunction = 'atq';
        end
        [ztmp,atmp,ittmp,ord,s,fct] = compute(n,y,order,threshold,costfunction);
    end

    function OrderFcn(source,eventdata)
        % Change order
        order = get(horder,'Value');
        set(horderval,'String',num2str(order));
        [ztmp,atmp,ittmp,ord,s,fct] = compute(n,y,order,threshold,costfunction);
    end

    function ThresholdFcn(source,eventdata)
        % Change threshold
        threshold = get(hthreshold,'String');
        threshold = str2double(threshold);
        [ztmp,atmp,ittmp,ord,s,fct] = compute(n,y,order,threshold,costfunction);
    end

    function [ztmp,atmp,ittmp,order,threshold,costfunction] = compute(n,y,order,threshold,costfunction)
        % Compute and plot an estimation (need to sort the data)
        [ztmp,atmp,ittmp,order,threshold,costfunction] = backcor(n,y,order,threshold,costfunction);
        [~,i] = sort(n);
        plot(n(i),y(i),'b-',n(i),ztmp(i),'r-');
    end

    function WinResizeFcn(source,eventdata)
        % Resize the window
        pos = get(hwin,'Position');
        w = pos(3);
        h = pos(4);
        if w>400 && h>100,
            set(haxes,'Position',[40,40,w-320,h-70]);
        end;
        set(hok,'Position',[w-90,30,80,25]);
        set(hcancel,'Position',[w-180,30,80,25]);
        set(hfctlbl,'Position',[w-240,h-30,220,20]);
        set(hfct,'Position',[w-240,h-50,220,25]);
        set(hthresholdlbl,'Position',[w-240,h-80,220,20]);
        set(hthreshold,'Position',[w-240,h-100,220,20]);
        set(horderlbl,'Position',[w-240,h-130,220,20]);
        set(horder,'Position',[w-210,h-150,190,20]);
        set(horderval,'Position',[w-240,h-150,20,20]);
    end

uiwait(gcf);

end


text:
function test_suite = test_surficial_neighborhood
% tests for cosmo_surficial_neighborhood
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #
    try % assignment of 'localfunctions' is necessary in Matlab >= 2016
        test_functions=localfunctions();
    catch % no problem; early Matlab versions can use initTestSuite fine
    end
    initTestSuite;

function test_surficial_neighborhood_surface_dijkstra
    if cosmo_skip_test_if_no_external('surfing')
        return;
    end
    warning_state=warning();
    warning_resetter=onCleanup(@()warning(warning_state));
    warning('off');

    opt=struct();
    opt.progress=false;

    ds=cosmo_synthetic_dataset('type','surface');

    [vertices,faces]=get_synthetic_surface();

    % dijkstra neighborhood fixed number of voxels
    args={{vertices,faces},'count',4,'metric','dijkstra',opt};
    nh1=cosmo_surficial_neighborhood(ds,args{:});
    assertFalse(isfield(nh1.a,'vol'));
    assert_equal_cell(nh1.neighbors,{ [ 1 2 4 3 ]
                                        [ 2 1 3 5 ]
                                        [ 3 2 6 5 ]
                                        [ 4 1 5 2 ]
                                        [ 5 2 4 6 ]
                                        [ 6 3 5 2 ] });

    assertEqual(nh1.fa.radius,[2 1 sqrt(2) sqrt(2) 1 2]);
    assertEqual(nh1.fa.node_indices,1:6);
    check_area(vertices,faces,nh1);

    args={{vertices,faces},'radius',2.5,'metric','dijkstra',opt};
    nh2=cosmo_surficial_neighborhood(ds,args{:});

    assert_equal_cell(nh2.neighbors,{[1 2 3 4 5];
                                     [1 2 3 4 5 6];
                                     [1 2 3 4 5 6];
                                     [1 2 3 4 5 6];
                                     [1 2 3 4 5 6];
                                     [2 3 4 5 6]; });
    assertEqual(nh2.fa.radius,[2,2,1+sqrt(2),1+sqrt(2),2,2]);
    assertEqual(nh2.fa.node_indices,1:6);
    check_partial_neighborhood(ds,nh2,args);
    check_area(vertices,faces,nh2);

    args{1}{1}([2,6],:)=NaN;

    nh3=cosmo_surficial_neighborhood(ds,args{:});
    assert_equal_cell(nh3.neighbors,{[1 4 5];
                                     [];
                                     [3 4 5];
                                     [1 3 4 5];
                                     [1 3 4 5];
                                     []; });
    assertEqual(nh3.fa.radius,[2,NaN,1+sqrt(2),1+sqrt(2),2,NaN]);
    check_partial_neighborhood(ds,nh3,args);
    check_area(args{1}{:},nh3);

    args{2}='count';
    args{3}=3;

    nh4=cosmo_surficial_neighborhood(ds,args{:});
    assert_equal_cell(nh4.neighbors,{[1 4 5];
                                     [];
                                     [3 5 4];
                                     [4 1 5];
                                     [5 4 3];
                                     []; });
    assertEqual(nh4.fa.radius,[2,NaN,1+sqrt(2),1,sqrt(2),NaN]);
    check_area(args{1}{:},nh4);

    args{1}{1}=vertices;
    args{1}{1}([2,5],:)=NaN; % split in two surfaces
    args{3}=2;
    nh5=cosmo_surficial_neighborhood(ds,args{:});
    assert_equal_cell(nh5.neighbors,{[1 4];
                                     [];
                                     [3 6];
                                     [4 1];
                                     [];
                                     [6 3]; });
    assertEqual(nh5.fa.radius,[1,NaN,1,1,NaN,1]);
    check_area(args{1}{:},nh5);

    % throw error when too many nodes asked for
    aet=@(varargin)assertExceptionThrown(@()...
                    cosmo_surficial_neighborhood(varargin{:}),'');
    args{2}='count';
    args{3}=3;

    aet(ds,args{:})

function test_surficial_neighborhood_surface_direct
    if cosmo_skip_test_if_no_external('surfing')
        return;
    end
    warning_state=warning();
    warning_resetter=onCleanup(@()warning(warning_state));
    warning('off');

    opt=struct();
    opt.progress=false;

    ds=cosmo_synthetic_dataset('type','surface');

    [vertices,faces]=get_synthetic_surface();

    % direct neighborhood
    args={{vertices,faces},'direct',true,opt};
    nh3=cosmo_surficial_neighborhood(ds,args{:});
    assert_equal_cell(nh3.neighbors,{ [ 1 2 4 ]
                                        [ 2 3 1 4 5 ]
                                        [ 3 2 5 6 ]
                                        [ 4 1 2 5 ]
                                        [ 5 3 6 2 4 ]
                                        [ 6 3 5 ] });
    assertElementsAlmostEqual(nh3.fa.radius,sqrt([1 2 2 2 2 1]));
    check_partial_neighborhood(ds,nh3,args);

    args{1}{1}([2 5],:)=NaN;
    nh4=cosmo_surficial_neighborhood(ds,args{:});
    assert_equal_cell(nh4.neighbors,{ [ 1 4 ]
                                        []
                                        [ 3 6 ]
                                        [ 1 4 ]
                                        []
                                        [6 3] });

    check_partial_neighborhood(ds,nh4,args);

function test_surficial_neighborhood_surface_geodesic
    if cosmo_skip_test_if_no_external('fast_marching') || ...
            cosmo_skip_test_if_no_external('surfing')
        return;
    end

    warning_state=warning();
    warning_resetter=onCleanup(@()warning(warning_state));
    warning('off');

    opt=struct();
    opt.progress=false;

    ds=cosmo_synthetic_dataset('type','surface');%,'size','normal');

    [vertices,faces]=get_synthetic_surface();

    args={{vertices,faces},'count',4,opt};
    nh=cosmo_surficial_neighborhood(ds,args{:});
    assert_equal_cell(nh.neighbors,{[ 1 2 4 5 ]
                                    [ 2 1 3 5 ]
                                    [ 3 2 6 5 ]
                                    [ 4 1 5 2 ]
                                    [ 5 2 4 6 ]
                                    [ 6 3 5 2 ] });
    assertEqual(nh.fa.node_indices,1:6);
    assertEqual(nh.fa.radius,[sqrt(.5)+1 1 sqrt(2) sqrt(2) 1 sqrt(.5)+1]);

    vertices2=[NaN NaN NaN; vertices;NaN NaN NaN];
    faces2=[faces+1; 1 1 8];
    args={{vertices2,faces2},'count',4,opt};
    nh2=cosmo_surficial_neighborhood(ds,args{:});

    assertEqual(nh2.neighbors,{ zeros(1,0)
                                 [ 2 3 5 6 ]
                                 [ 3 2 4 6 ]
                                 [ 4 3 6 2 ]
                                 [ 5 2 6 3 ]
                                 [ 6 3 5 4 ] });

function test_surficial_neighborhood_volume_geodesic
    if cosmo_skip_test_if_no_external('fast_marching') || ...
            cosmo_skip_test_if_no_external('surfing')
        return;
    end
    warning_state=warning();
    warning_resetter=onCleanup(@()warning(warning_state));
    warning('off');

    opt=struct();
    opt.progress=false;

    ds=cosmo_synthetic_dataset();
    vertices=[-2 0 2 -2 0 2;...
                -1 -1 -1 1 1 1;...
                -1 -1 -1 -1 -1 -1]';
    faces= [ 3 2 3 2;...
                2 1 5 4;...
                5 4 6 5 ]';

    pial=vertices;
    pial(:,3)=pial(:,3)+1;
    white=vertices;
    white(:,3)=white(:,3)-1;
    nh1=cosmo_surficial_neighborhood(ds,{vertices,[-1 1],faces},...
                                    'count',4,opt);
    nh2=cosmo_surficial_neighborhood(ds,{pial,white,faces},...
                                    'count',4,opt);
    assert_equal_cell(nh1.neighbors,{[ 1 2 4 5 ]
                                        [ 1 2 3 5 ]
                                        [ 2 3 5 6 ]
                                        [ 4 1 5 2 ]
                                        [ 5 2 4 6 ]
                                        [ 6 3 5 2 ] });
    assertEqual(nh1.fa.node_indices,1:6);
    assertEqual(nh1,nh2);

function test_surficial_neighborhood_volume_dijkstra
    if cosmo_skip_test_if_no_external('surfing')
        return;
    end
    warning_state=warning();
    warning_resetter=onCleanup(@()warning(warning_state));
    warning('off');

    opt=struct();
    opt.progress=false;

    ds=cosmo_synthetic_dataset();
    vertices=[-2 0 2 -2 0 2;
                -1 -1 -1 1 1 1
                -1 -1 -1 -1 -1 -1]';
    faces= [ 3 2 3 2
                2 1 5 4
                5 4 6 5 ]';

    pial=vertices;
    pial(:,3)=pial(:,3)+1;
    white=vertices;
    white(:,3)=white(:,3)-1;

    args3={{vertices,[-1 1],faces},'metric','dijkstra','count',4,opt};
    args4={{pial,white,faces},'metric','dijkstra','count',4,opt};
    nh3=cosmo_surficial_neighborhood(ds,args3{:});
    nh4=cosmo_surficial_neighborhood(ds,args4{:});
    assert_equal_cell(nh3.neighbors,{ [ 1 2 4 3 ]
                                        [ 2 1 3 5 ]
                                        [ 3 2 6 5 ]
                                        [ 4 1 5 2 ]
                                        [ 5 2 4 6 ]
                                        [ 6 3 5 2 ] });
    assertEqual(nh3.fa.node_indices,1:6);
    assert_equal_cell(nh4.neighbors,nh3.neighbors);
    assertFalse(isfield(nh3.a,'vol'));

    % TODO
    % check_partial_neighborhood(ds,nh3,args3);
    % check_partial_neighborhood(ds,nh3,args3);

function test_surficial_neighborhood_exceptions
    if cosmo_skip_test_if_no_external('surfing')
        return
    end
    ds=cosmo_synthetic_dataset('type','surface');%,'size','normal');
    [vertices,faces]=get_synthetic_surface();

    aet=@(varargin)assertExceptionThrown(@()...
                    cosmo_surficial_neighborhood(varargin{:},...
                                        'progress',false),'');
    aet(ds,{vertices,faces});

    % need surfaces
    aet(ds,{},'radius',2);

    % center_ids not supported for surface dataset
    aet(ds,{vertices,faces},'radius',2,'center_ids',1);

    % cannot have duplicate feature ids
    ds_double=cosmo_stack({ds,ds},2);
    ds_double.a.fdim.values{1}=[1:6,1:6];
    aet(ds_double,{vertices,faces},'radius',2);

    % outside range
    ds_double.a.fdim.values{1}=[1:5,1:5];
    aet(ds_double,{vertices,faces},'radius',2);

    % cannot have fmri and surface dataset combined
    ds2=cosmo_synthetic_dataset();
    ds2.a.fdim.values{end+1}=ds.a.fdim.values{1};
    ds2.a.fdim.labels{end+1}=ds.a.fdim.labels{1};
    ds2.fa.node_indices=ds.fa.node_indices;
    aet(ds2,{vertices,faces},'radius',2);

    % cannot have MEEG dataset
    ds_meeg=cosmo_synthetic_dataset('type','meeg');
    aet(ds_meeg,{vertices,faces},'radius',2);

    % need positive scalar radius
    aet(ds,{vertices,faces},'radius',-1);
    aet(ds,{vertices,faces},'radius',eye(2));

function check_partial_neighborhood(ds,nh,args)
    % see if when have a partial dataset, the neighborbood reflects
    % that too

    nf=size(ds.samples,2);

    rp=randperm(nf);
    keep_count=round(nf*.7);
    keep_sel=rp(1:keep_count);
    keep_all=[keep_sel keep_sel keep_sel];

    ds_sel=cosmo_slice(ds,keep_all,2);

    fdim=ds_sel.a.fdim.values{1};
    rp_fdim=randperm(numel(fdim));
    ds_sel.a.fdim.values{1}=fdim(rp_fdim);

    nh_sel=cosmo_surficial_neighborhood(ds_sel,args{:});

    assertEqual(numel(nh_sel.neighbors), numel(keep_sel));

    assertEqual(nh_sel.a.fdim.labels,nh.a.fdim.labels);
    assertEqual(nh_sel.a.fdim.values{1},nh.a.fdim.values{1}(rp_fdim));
    assertEqual(numel(nh_sel.a.fdim.values),numel(nh.a.fdim.values));

    assertEqual(ds_sel.a,nh_sel.a);

    opt=cosmo_structjoin(args(2:end));

    if isfield(opt,'radius')
        metric=opt.metric;
        metric_arg=opt.radius;
    elseif isfield(opt,'count')
        metric=opt.metric;
        metric_arg=[10 opt.count];
    elseif isfield(opt,'direct')
        metric='direct';
        if opt.direct
            metric_arg=NaN;
        else
            metric_arg=0;
        end
    else
        assert(false);
    end

    faces=args{1}{2};
    n2f=surfing_invertmapping(faces);

    nodes_ds_sel=ds_sel.a.fdim.values{1}(ds_sel.fa.node_indices);
    nodes_ds=ds.a.fdim.values{1}(ds.fa.node_indices);

    nodes_nh_sel=nh_sel.a.fdim.values{1}(nh_sel.fa.node_indices);

    vertices=args{1}{1};

    nvertices=size(vertices,1);
    nodes_kept=cosmo_match(1:nvertices,nodes_ds_sel);
    vertices(~nodes_kept,:)=NaN;

    node_mask=all(isfinite(vertices),2);

    nodes_removed=setdiff(nodes_ds(:)',nodes_ds_sel(:)');
    assertEqual(setxor(nodes_removed,nodes_ds_sel),1:nf);

    nb_sel=nh_sel.neighbors;
    for k=1:numel(nh_sel.neighbors)
        sel_center_node=nodes_nh_sel(k);
        idx=find(nodes_ds==sel_center_node);
        assert(numel(idx)==1);
        center_node=nodes_ds(idx);

        assertEqual(sel_center_node, center_node);

        switch metric
            case 'direct'
                if node_mask(sel_center_node)
                    direct_neighbors=surfing_surface_nbrs(faces',...
                                                            vertices');
                    around_nodes=direct_neighbors(sel_center_node,:);
                    msk=cosmo_match(around_nodes, ...
                                            find(isfinite(vertices(:,1))));
                    % add node itself
                    around_nodes=[sel_center_node,...
                                    around_nodes(msk & around_nodes>0)];
                else
                    around_nodes=[];
                end
            otherwise
                around_nodes=surfing_circleROI(vertices',faces',...
                            sel_center_node,metric_arg,metric,n2f);
        end

        sel_around_nodes=nodes_ds_sel(nb_sel{k});

        if isempty(sel_around_nodes)
            assertTrue(isempty(around_nodes));
        else
            assertEqual(unique(sel_around_nodes),...
                        setdiff(around_nodes, nodes_removed))
        end
    end

function test_surface_subsampling
    if cosmo_skip_test_if_no_external('surfing')
        return;
    end
    vertices=[0 -1 -2 -1  1  2  1  3  4  3;
          0 -2  0  2  2  0 -2  2  0 -2;
          0  0  0  0  0  0  0  0  0  0]';

    faces=[1 1 1 1 1 1 5 8 6  6;
           2 3 4 5 6 7 8 9 9 10;
           3 4 5 6 7 2 6 6 10 7]';

    % make custom volume
    ds=cosmo_synthetic_dataset('size','small');
    cp=cosmo_cartprod({1:7,1:3})';
    ds.fa.i=cp(1,:);
    ds.fa.j=cp(2,:);
    ds.fa.k=ds.fa.i*0+1;

    ds.a.fdim.values={1:7;1:3;1};
    ds.samples=zeros(numel(ds.sa.targets),numel(ds.fa.i));
    ds.a.vol.dim=cellfun(@numel,ds.a.fdim.values)';
    ds.a.vol.mat(2,4)=-4;
    ds.a.vol.mat(3,4)=-1;
    ds.a.vol.mat(1,1)=1;
    ds.a.vol.mat(2,2)=2;
    ds.a.vol.mat(3,3)=1;

    surfs={vertices,faces,[-4 5]};

    opt=struct();
    opt.progress=false;
    opt.radius=3;
    opt.metric='euclidean';
    nh=cosmo_surficial_neighborhood(ds,surfs,opt);

    assertEqual(nh.neighbors,{ [ 2 4 8 10 12 16 18 ]
                             [ 2 4 8 10 ]
                             [ 2 8 10 16 ]
                             [ 8 10 16 18 ]
                             [ 10 12 16 18 20 ]
                             [ 4 6 10 12 14 18 20 ]
                             [ 2 4 6 10 12 ]
                             [ 12 14 18 20 ]
                             [ 6 12 14 20 ]
                             [ 4 6 12 14 ] });
    assertEqual(nh.origin.fa,ds.fa);
    assertEqual(nh.origin.a,ds.a);

    % test subsampling
    subsample=2;
    surfs={vertices,faces,[-4 5],subsample};
    nh2=cosmo_surficial_neighborhood(ds,surfs,opt);
    assertEqual(nh2.neighbors,{ [ 2 4 8 10 ]
                                 [ 2 8 10 16 ]
                                 [ 8 10 16 18 ]
                                 [ 10 12 16 18 20 ]
                                 [ 2 4 6 10 12 ]
                                 [ 12 14 18 20 ]
                                 [ 6 12 14 20 ]
                                 [ 4 6 12 14 ] });

    assertEqual(nh2.origin.fa,ds.fa);
    assertEqual(nh2.origin.a,ds.a);

    % subsampling with pial surface
    pial=bsxfun(@plus,vertices,[0 0 1]);
    white=bsxfun(@plus,vertices,[0 0 -1]);
    [vo,fo]=surfing_subsample_surface(vertices,faces,2,.2,0);
    surfs={pial,white,faces,vo,fo};
    nh3=cosmo_surficial_neighborhood(ds,surfs,opt);
    assertEqual(nh2,nh3);

    % check center ids options
    slice_ids=[5 3 2];
    nh4=cosmo_surficial_neighborhood(ds,surfs,opt,'center_ids',slice_ids);
    nh4_sl=struct();
    nh4_sl.neighbors=nh3.neighbors(slice_ids);
    nh4_sl.fa=cosmo_slice(nh3.fa,slice_ids,2,'struct');
    nh4_sl.a=nh3.a;

    assertEqual(nh4.a.fdim.values{1}(nh4.fa.node_indices),...
                    nh4_sl.a.fdim.values{1}(nh4_sl.fa.node_indices));
    assertEqual(nh4.neighbors,nh4_sl.neighbors);

    % try with file names
    fn_pial=cosmo_make_temp_filename('pial','.asc');
    fn_white=cosmo_make_temp_filename('white','.asc');
    fn_tiny=cosmo_make_temp_filename('tiny','.asc');

    cleaner1=onCleanup(@()delete(fn_pial));
    cleaner2=onCleanup(@()delete(fn_white));
    cleaner3=onCleanup(@()delete(fn_tiny));

    surfing_write(fn_pial, pial, faces);
    surfing_write(fn_white, white, faces);
    surfing_write(fn_tiny, vo, fo);

    surfs={fn_pial,fn_white,fn_tiny};
    nh5=cosmo_surficial_neighborhood(ds,surfs,opt);
    assertEqual(nh2,nh5);

    % should work with alternative voldef
    ds_bad_vol=ds;
    ds_bad_vol.a.vol.mat(:)=NaN;
    ds_bad_vol.a.vol.dim(:)=NaN;
    nh6=cosmo_surficial_neighborhood(ds_bad_vol,surfs,opt,...
                                            'vol_def',ds.a.vol);
    nh6.origin.a.vol=ds.a.vol;
    assertEqual(nh5,nh6);

    % check exceptions
    aet=@(varargin)assertExceptionThrown(@()...
                    cosmo_surficial_neighborhood(varargin{:}),'');

    surfs={fn_pial,fn_pial,fn_tiny};
    aet(ds,surfs,opt);

    white_bad=white;
    white_bad=white_bad(2:end,:);
    aet(ds,{pial,white_bad,faces},opt);

    % missing faces for output surface
    aet(ds,{fn_pial,fn_white,vo},opt);

    % face mismatch
    faces_bad=faces;
    faces_bad=faces_bad(end:-1:1,:);

    surfing_write(fn_white, white, faces_bad);
    aet(ds,{fn_pial,fn_white},opt);

    % too many surf arguments
    aet(ds,{fn_pial,fn_white,fn_tiny,fn_tiny},opt);
    aet(ds,{pial,white,faces,pial,white,white},opt);
    aet(ds,{pial,white,faces,fn_pial,white},opt);

    % surfs are not a cell
    aet(ds,struct,opt);
    aet(ds,{pial,white,{}});

function check_area(vertices,faces,nh)
    assert(isfield(nh.fa,'area'));
    area=surfing_surfacearea(vertices,faces);

    node_idxs=nh.fa.node_indices;
    n_nodes=numel(node_idxs);
    for k=1:n_nodes
        nbr_idxs=node_idxs(nh.neighbors{k});
        expected_area=sum(area(nbr_idxs));
        if isnan(expected_area)
            assertEqual(expected_area,nh.fa.area(k))
        else
            assertElementsAlmostEqual(expected_area,nh.fa.area(k));
        end
    end

function [vertices,faces]=get_synthetic_surface()
    % return the following surface (face indices in [brackets])
    %
    %  1-----2-----3
    %  |    /|    /|
    %  |[2]/ |[1]/ |
    %  |  /  |  /  |
    %  | /[4]| /[3]|
    %  |/    |/    |
    %  4-----5-----6

    vertices=[0 0 0 1 1 1;
                1 2 3 1 2 3;
                0 0 0 0 0 0]';
    faces= [ 3 2 3 2
                2 1 5 4
                5 4 6 5 ]';

function assert_equal_cell(x,y)
    % small helper
    assertEqual(size(x),size(y))
    for k=1:numel(x)
        xk=x{k};
        yk=y{k};
        if isempty(xk)
            assertTrue(isempty(yk));
        else
            assertEqual(sort(xk),sort(yk));
        end
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [ x, y ] = srot ( n, x, incx, y, incy, c, s )

%*****************************************************************************80
%
%% SROT applies a plane rotation.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 November 2006
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979.
%
%    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
%    Basic Linear Algebra Subprograms for Fortran Usage,
%    Algorithm 539,
%    ACM Transactions on Mathematical Software,
%    Volume 5, Number 3, September 1979, pages 308-323.
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vectors.
%
%    Input, real X(*), one of the vectors to be rotated.
%
%    Input, integer INCX, the increment between successive entries of X.
%
%    Input, real Y(*), one of the vectors to be rotated.
%
%    Input, integer INCY, the increment between successive elements of Y.
%
%    Input, real C, S, parameters (presumably the cosine and
%    sine of some angle) that define a plane rotation.
%
%    Output, real X(*), the rotated vector.
%
%    Output, real Y(*), the rotated vector.
%
  if ( n <= 0 )

    x = [];
    y = [];

  elseif ( incx == 1 & incy == 1 )

    for i = 1 : n
      stemp = c * x(i) + s * y(i);
      y(i) =  c * y(i) - s * x(i);
      x(i) = stemp;
    end

  else

    if ( 0 <= incx )
      ix = 1;
    else
      ix = ( - n + 1 ) * incx + 1;
    end

    if ( 0 <= incy )
      iy = 1;
    else
      iy = ( - n + 1 ) * incy + 1;
    end

    for i = 1 : n
      stemp = c * x(ix) + s * y(iy);
      y(iy) = c * y(iy) - s * x(ix);
      x(ix) = stemp;
      ix = ix + incx;
      iy = iy + incy;
    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


function out = get_scale_sample(im, pos, base_target_sz, scaleFactors, scale_window, scale_model_sz)

% out = get_scale_sample(im, pos, base_target_sz, scaleFactors, scale_window, scale_model_sz)
%
% Extracts a sample for the scale filter at the current
% location and scale.

nScales = length(scaleFactors);

for s = 1:nScales
    patch_sz = floor(base_target_sz * scaleFactors(s));

    xs = floor(pos(2)) + (1:patch_sz(2)) - floor(patch_sz(2)/2);
    ys = floor(pos(1)) + (1:patch_sz(1)) - floor(patch_sz(1)/2);

    % check for out-of-bounds coordinates, and set them to the values at
    % the borders
    xs(xs < 1) = 1;
    ys(ys < 1) = 1;
    xs(xs > size(im,2)) = size(im,2);
    ys(ys > size(im,1)) = size(im,1);

    % extract image
    im_patch = im(ys, xs, :);

    % resize image to model size
    im_patch_resized = imResample(im_patch, scale_model_sz);

    % extract scale features
    temp_hog = fhog(single(im_patch_resized), 4);
    temp = temp_hog(:,:,1:31);

    if s == 1
        out = zeros(numel(temp), nScales, 'single');
    end

    % window
    out(:,s) = temp(:) * scale_window(s);
end


text:
function riverplot(layer1fmri_obj, varargin)
% Make a riverplot of relationships among images, stored in an fmri_data object,
% or two sets of images, stored in two fmri_data objects.
%
% - With one input, uses self
% - With two inputs, uses first as Layer 1 and second as Layer 2 of riverplot
% - Uses obj.image_names for titles of images, so set these before calling this function
%
% :Usage:
% ::
%
%     [list outputs here] = function_name(list inputs here, [optional inputs])
%
% For objects: Type methods(object_name) for a list of special commands
%              Type help object_name.method_name for help on specific
%              methods.
%
% ..
%     Author and copyright information:
%
%     Copyright (C) 2016 Tor Wager
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
% ..
%
% :Inputs:
%
%   **layer1fmri_obj:**
%        An fmri_data object.  Default is to plot associations between
%        images in this object.
%
% :Optional Inputs:
%   **layer2:**
%        A second fmri_data object.  This is layer 2 of the plot.  Plots
%        assocations between images in layer1 object and layer2 object
%
%   **r:**
%        Change similarity metric to correlation.  Default is cosine
%        similarity.
%
%   **dice:**
%        Change similarity metric to dice coefficient.  Default is cosine
%        similarity.
%
%   **pos:**
%        Change similarity matrix thresholding to positive values only.
%        Default is to use original similarity matrix.
%
%   **neg:**
%        Change similarity matrix thresholding to negative values only.
%        Default is to use original similarity matrix.
%
%   **abs:**
%        Change similarity matrix thresholding to absolute values.
%        Default is to use original similarity matrix.
%
%   **reorder:**
%        Automatic reordering of rows of plot to try to improve visibility.
%
%   **recolor:**
%        Automatic color assignment of layer2 colors based on layer1
%        weights
%
%   **layer1order:**
%        Followed by custom order for Layer 1, specified by vector of
%        integers, e.g., [1 3 2 4 5 6]
%
%   **layer2order:**
%        Followed by custom order for Layer 2, specified by vector of
%        integers, e.g., [1 3 2 4 5 6]
%
%   **colors1:**
%        Followed by cell vector of [r g b] colors for each image in Layer 1
%
%   **colors2:**
%        Followed by cell vector of [r g b] colors for each image in Layer 2
%
%   **thin:**
%        Turn off thick/full-rectangle ribbon width coming from Layer 1
%
%
% :Outputs:
%
%   **layer1:**
%        Vector of structures describing rectangles in Layer 1, with
%        handles, positions, etc.
%
%   **layer2:**
%        Vector of structures describing rectangles in Layer 2, with
%        handles, positions, etc.
%
%   **ribbons:**
%        Vector of structures describing ribbons, with
%        handles, positions, etc.
%
%   **layer2colors:**
%        cell array specifying layer2 colors if recolor option is used
%
%
% :Examples:
% ::
%
% layer1fmri_obj = fmri_data(placeboimgs);
% layer2fmri_obj = fmri_data(appraisalimgs);
% riverplot(layer1fmri_obj, 'layer2', layer2fmri_obj);
%
% Add colors and names:
% layer1colors = {[.8 .8 .3] [.3 .3 1]};
% layer2colors = {[0 .2 .8] [1 1 .2] [1 .4 .4] [0 .6 0] [1 .7 0]};  %seaborn_colors(n1);
%
% layer1names = {'Placebo increases', 'Placebo decreases'};  %stats.networknames;
% layer2names = {'Neg Emo' 'Pos Emo' 'Emo Reg' 'Social Cog' 'Self'};
%
% layer1fmri_obj.image_names = char(layer1names{:});
% layer2fmri_obj.image_names = char(layer2names{:});
%
% riverplot(layer1fmri_obj, 'layer2', layer2fmri_obj, 'layer1colors', layer1colors, 'layer2colors', layer2colors);
%
% Example: Plot Buckner lab network maps with themselves:
% [bucknermaps, networknames] = load_image_set('bucknerlab');
% bucknermaps.image_names = char(networknames{:});
%
% layer1colors = bucknerlab_colors(7);
% riverplot(bucknermaps, 'pos', 'reorder', 'layer1colors', layer1colors, 'layer2colors', layer1colors);
%
% riverplot(bucknermaps, 'pos', 'layer1colors', layer1colors, 'layer2colors', layer1colors);
%
%
% :References:
%   Based off of concept in River Plot R package.
%
% :See also:
%   - image_similarity_plot, tor_polar_plot
%

% ..
%    Programmers' notes:
%    List dates and changes here, and author of changes
%    Created July 2016 by Tor Wager
% ..

% BELOW IS A STANDARD TEMPLATE FOR DEFINING VARIABLE (OPTIONAL) INPUT
% ARGUMENTS. MANY FUNCTIONS NEED TO PARSE OPTIONAL ARGS, SO THIS MAY BE
% USEFUL.

% ..
%    DEFAULTS AND INPUTS
% ..

% optional inputs with default values

layer2fmri_obj = layer1fmri_obj;  % Default: Plot object with itself
sim_metric = 'cosine_sim';
threshold = 'none';

nlayer1 = size(layer1fmri_obj.dat, 2);
nlayer2 = size(layer2fmri_obj.dat, 2);

layer1colors = custom_colors([.2 .8 .2], [.8 .8 .2], nlayer1);
layer2colors = custom_colors([.2 .2 .8], [.2 .8 .8], nlayer2);

doreorder = 0;
recolor = 0;

layer1coveragestr = 'layer1fullcoverage';

layer1order = [];

layer2order = [];

for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case {'layer2', 'layer2fmri_obj'},
                layer2fmri_obj = varargin{i+1}; varargin{i+1} = [];
                nlayer2 = size(layer2fmri_obj.dat, 2);
                layer2colors = custom_colors([.2 .2 .8], [.2 .8 .8], nlayer2);

            case 'sim_metric', sim_metric = varargin{i+1}; varargin{i+1} = [];
            case 'threshold', threshold = varargin{i+1}; varargin{i+1} = [];

            case {'colors1', 'layer1colors'}, layer1colors = varargin{i+1}; varargin{i+1} = [];
            case {'colors2', 'layer2colors'}, layer2colors = varargin{i+1}; varargin{i+1} = [];

            case 'reorder', doreorder = 1;
              case 'recolor', recolor = 1;

            case {'r', 'corr', 'correlation'}
                sim_metric = 'r';

            case 'dice'
                sim_metric = 'dice';

            case 'pos'
                threshold = 'pos';

            case 'neg'
                threshold = 'neg';

            case 'abs'
                threshold = 'abs';

            case 'thin'
                layer1coveragestr = ' ';

            case 'layer1order', layer1order = varargin{i+1}; varargin{i+1} = [];
            case 'layer2order', layer2order = varargin{i+1}; varargin{i+1} = [];

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

% -------------------------------------------------------------------------
% Get similarity matrix
% -------------------------------------------------------------------------

switch sim_metric

    case 'cosine_sim' % default

        stats = image_similarity_plot(layer1fmri_obj, 'mapset', layer2fmri_obj, 'noplot', 'cosine_similarity');
        sim_matrix = stats.r;

    case {'r', 'corr', 'correlation'}

        stats = image_similarity_plot(layer1fmri_obj, 'mapset', layer2fmri_obj, 'noplot');
        sim_matrix = stats.r;

    case 'dice'

        dice_coeff = dice_coeff_image(cat(layer1fmri_obj, layer2fmri_obj));

        sim_matrix = dice_coeff(1:nlayer1, nlayer1+1:end)';

    otherwise
        error('Unknown sim_metric. Choose r, cosine_sim, or dice');
end

% -------------------------------------------------------------------------
% Threshold
% -------------------------------------------------------------------------

switch threshold

    case 'none'
        % do nothing - default

    case {'pos', 'positive'}

        % Threshold to look at positive associations
        sim_matrix(sim_matrix < 0) = 0;

    case {'neg', 'negative'}
        sim_matrix = -sim_matrix;
        sim_matrix(sim_matrix < 0) = 0;

    case 'abs'
        sim_matrix = abs(sim_matrix);

    otherwise
        error('Unknown threshold type. Choose pos, neg, or abs');
end

% -------------------------------------------------------------------------
% Names and sizes
% -------------------------------------------------------------------------

layer1names = format_strings_for_legend(layer1fmri_obj.image_names);
layer2names = format_strings_for_legend(layer2fmri_obj.image_names);

[n2, n1] = size(sim_matrix);

if length(layer1names) < n1, error('Layer 1: Not enough names entered. Add names in input layer1 fmri_object.image_names'); end
if length(layer2names) < n2, error('Layer 2: Not enough names entered. Add names in input layer2 fmri_object.image_names'); end

% -------------------------------------------------------------------------
% Reorder rows and columns if asked for
% -------------------------------------------------------------------------

if doreorder

    %  REORDER rows here based on similarity to make plot look cleaner
    [sim_matrix, layer2names, layer2colors] = riverplot_reorder_matrix(sim_matrix, layer2names, layer2colors);

end

if ~isempty(layer1order)

    sim_matrix = sim_matrix(:, layer1order);
    layer1names = layer1names(layer1order);
    layer1colors = layer1colors(layer1order);

end

if ~isempty(layer2order)

    sim_matrix = sim_matrix(layer2order, :);
    layer2names = layer2names(layer2order);
    layer2colors = layer2colors(layer2order);

end

if recolor
    layer2colors=riverplot_recolor_layer2(sim_matrix,layer1colors);
end

% -------------------------------------------------------------------------
% Find start coordinate points for layers
% -------------------------------------------------------------------------

layer2x = max(n1, n2) ./ 2;  % 1/2 largest array size, to keep x proportions balanced with y

% each rect is 1 unit, starting at y = 0 for layer. so if same size, no y_loc offset
% midpoint should be same for both
mid2 = (n2 + 1) ./ 2;
% mid1 should equal mid2  : % (n1 + 1) ./ 2 - y_loc = mid 2
y_loc = -(n1 + 1) ./ 2 + mid2;

% -------------------------------------------------------------------------
% Create the plot
% -------------------------------------------------------------------------

create_figure('riverplot');
set(gca, 'YDir', 'reverse');
layer1 = riverplot_draw_layer(0, n1, 'colors', layer1colors, 'y_loc', y_loc);
layer2 = riverplot_draw_layer(layer2x, n2, 'colors', layer2colors);

ribbons = riverplot_ribbon_matrix(layer1, layer2, sim_matrix, 'colors', layer1colors, layer1coveragestr, 'steepness', 0);

% Turn off lines
riverplot_toggle_lines(ribbons);

% Increase opacity
riverplot_set_ribbon_property(ribbons, 'FaceAlpha', .6);

set(gca, 'XLim', [-3 layer2x+3]);

% Add names
layer1 = riverplot_layer_names(layer1, layer1names);
layer2 = riverplot_layer_names(layer2, layer2names, 'right');

axis off

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [W,G,IM,C] = remesh_planar_patches(V,F,varargin)
  % REMESH_PLANAR_PATCHES  Find nearly planar patches and retriangulate them.
  %   (V,F) should probably not be self-intersecting, at least not near the
  %   planar patches. This will attempt to maintain non-manifold edges
  %   (untested).
  %
  % [W,G] = remesh_planar_patches(V,F)
  % [W,G,IM,C] = remesh_planar_patches(V,F,'ParameterName',ParameterValue,...)
  %
  % Inputs:
  %   V  #V by 3 list of vertices
  %   F  #F by 3 list of triangle indices
  %   Optional:
  %     'Force' followed by whether to remesh even if it will not improve
  %        the number of vertices {false}.
  %     'MinSize' followed by minimum number of facets in group to consider for
  %       remeshing {4}
  %     'MinDeltaAngle' followd by minimum change in angle between neighboring
  %     facets to be considered co-planar {pi-1e-5}
  %     'Except' followed by list of faces not to remesh
  %
  % Outputs:
  %   W  #W by 3 list of output mesh positions
  %   G  #G by 3 list of output triangle indices
  %   IM  indices from final remove_unreferenced
  %   C  #F list of patch indices
  %

  % force remesh even if there are no internal vertices to remove
  force_remesh = false;
  % Minmimum angle for two neighboring facets to be considered coplanar
  % in radians
  min_delta_angle = pi-1e-5;
  %min_delta_angle = pi-1e-3;
  min_size = 4;
  triangle_flags = '';
  except = [];
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Force','MinSize','MinDeltaAngle','TriangleFlags','Except'}, ...
    {'force_remesh','min_size','min_delta_angle','triangle_flags','except'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  % The right way to deal with them would be to make a list of feature edges
  % (aka "exterior edges") on the original mesh and be sure these end up in any
  % planar pataches. They would include the patch boundaries and any other
  % internal non-manifold edges.
  NME = nonmanifold_edges(F);
  C = planar_patches(V,F,'MinDeltaAngle',min_delta_angle,'Except',except);

  [UC,~,ic] = unique(C);
  ucounts = histc(C,UC);
  counts = ucounts(ic);

  %tsurf(F,V,'CData',C);
  %axis equal;
  %drawnow;

  W = V;
  G = F(counts<=min_size,:);
  % loop over connected components
  %for c = reshape(UC(ucounts>min_size),1,[])
  for c = 1:max(C)
    if ucounts(c) > min_size
      Fc = F(C==c,:);
      assert(size(Fc,1)>min_size);
      %S = statistics(V,Fc,'Fast',true);
      %% Not dealing with holes yet
      %if S.num_boundary_loops > 1
      %  warning('Skipping high-genus planar patch...');
      %  G = [G;Fc];
      %  continue;
      %end
      %assert(S.num_boundary_loops == 1);
      uf = unique(Fc);
      % Fit plane to all points
      [~,A] = affine_fit(V(uf,:));
      % Only keep outline points
      O = outline(Fc);

      % Non manifold edges on this patch
      NMEc = NME(all(ismember(NME,uf),2),:);
      E = [O;NMEc];

      uo = unique(E);
      % Nothing to gain by remeshing
      if ~force_remesh && numel(uo) == numel(uf)
        G = [G;Fc];
        continue;
      end
      % only boundary edges projected to plane
      Vuo = V(uo,:)*A;
      J = 1:size(V,1);
      J(uo) = 1:numel(uo);
      % Remap E to Vu
      O = J(O);
      E = J(E);

      % THIS IS INSANELY SLOW
      % DT = delaunayTriangulation(Vuo(:,1),Vuo(:,2));
      % DT.Constraints = E;
      % Gc = DT.ConnectivityList;
      % Triangle is way faster...
      [Wuo,Gc] = triangle(Vuo,E,[],'Quiet','Flags',triangle_flags);
      %tsurf(Gc,Wuo);
      %input('');
      assert(size(Gc,1) >= 1);
      assert(size(Wuo,1) >= size(Vuo,1));
      % easier to remove holes post hoc than pass hole positions to triangle
      AE = adjacency_matrix(E);
      [~,CE] = conncomp(AE);
      if max(CE) > 1
        w = winding_number(Vuo,O,barycenter(Vuo,Gc));
        % should only be 1s and 0s
        Gc = Gc(abs(w)>0.1,:);
      end

      % remap Gc to V
      uo = [uo;(size(W,1)+(1:(size(Wuo,1)-size(Vuo,1))))'];
      Wnew = Wuo((size(Vuo,1)+1):end,:);
      VA = V*A;
      %tsurf(Fc,VA);
      %hold on;
      %scatter(Wnew(:,1),Wnew(:,2));
      %hold off;
      %pause

      I = in_element(VA,Fc,Wnew);
      [~,I] = find(I);
      assert(numel(I) == size(Wnew,1));
      B = barycentric_coordinates( ...
        Wnew,VA((Fc(I,1)),:),VA((Fc(I,2)),:),VA((Fc(I,3)),:));
      Wnew = sum(bsxfun(@times, permute(B,[1 3 2]), ...
        cat(3,V(Fc(I,1),:),V(Fc(I,2),:),V(Fc(I,3),:))),3);
      W = [W;Wnew];
      Gc = uo(Gc);

      %tsurf(Gc,W);
      %pause
      % old mean normal
      oldN = mean(normalizerow(normals(V,Fc)));
      N = mean(normalizerow(normals(W,Gc)));
      flip  = dot(oldN,N) < 0;
      if flip
        Gc = fliplr(Gc);
      end
      %tsurf(Gc,W);
      %input('');
      G = [G;Gc];
      assert(max(G(:))<=size(W,1));
    end
  end

  [W,IM] = remove_unreferenced(W,G);
  G = IM(G);

end

author: alecjacobson

repo: gptoolbox


text:
function s = y_to_s ( y )

%*****************************************************************************80
%
%% Y_TO_S writes a year into a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y, the year.
%
%    Output, string S, a representation of the year.
%
  s = sprintf ( '%d', y );

  return
end


text:
% lsdescent;
% check descent condition
%

cont=max(alist==0); % condition for continue

if cont,
  [fbest,i]=min(flist);
  if alist(i)<0,
    if alist(i)>=4*alist(i+1), cont=0; end;
  elseif alist(i)>0,
    if alist(i)<4*alist(i-1), cont=0; end;
  else
    if i==1, fbest=flist(2);
    elseif i==s, fbest=flist(s-1);
    else fbest=min(flist(i-1),flist(i+1));
    end;
  end;
end;

if cont,
  % force local descent step
  if alist(i)~=0, alp=alist(i)/3;
  elseif i==s, alp=alist(s-1)/3;
  elseif i==1, alp=alist(2)/3;
  else
    % split wider adjacent interval
    if alist(i+1)-alist(i)>alist(i)-alist(i-1), alp=alist(i+1)/3;
    else alp=alist(i-1)/3;
    end;
  end;

  % new function value
  falp=feval(func,data,x+alp*p);
  alist=[alist,alp];flist=[flist,falp];
  lssort;
  if prt>1, disp(['descent check: new point at ',num2str(alp)]); end;
end;

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
  function x = sinc_periodic(t, K)
%|function x = sinc_periodic(t, K)
%| periodic sinc function, obtained by replicates of a sinc:
%| x(t) = \sum_l sinc(t - l K) = ... = sin(pi*t) / tan(pi*t/K) / K for K even.
%| This function is bandlimited and its samples are an impulse train.
%| It is closely related to the Dirichlet function diric() for odd K,
%| but it differs for even K.
%| Copyright 2003-11-2, Jeff Fessler, University of Michigan

if nargin < 1, ir_usage, end
if streq(t, 'test'), sinc_periodic_test, return, end

if ~rem(K,2) % even
	d = tan(pi*t/K);
	j = abs(d) > 1e-12;
	x = ones(size(t));
	t = t(j);
	d = d(j);
	x(j) = sin(pi*t) ./ d / K;
else
	x = nufft_diric(t, K, K, 1);
%	xo = diric(2*pi*t/K,K); % would require matlab's signal toolbox
%	max_percent_diff(xo, x)
end

%
% self test
%
function sinc_periodic_test

Klist = [4 5];
im clf, pl=220;
for kk=1:2
	K = Klist(kk);
	n = [0:(4*K)]';
	t = linspace(0,4*K,401)';
	x = @(t,K) sinc_periodic(t, K);
%	y = @(t,K) diric(2*pi*t/K,K); % would require matlab's signal toolbox
	y = @(t,K) nufft_diric(t,K,K,1);
	if im
		subplot(pl+kk+0)
		plot(t, x(t,K), '-', n, x(n,K), 'o')
		titlef('Sinc-Periodic K=%d', K)
		axis([0 4*K -1 1]), xtick([0:4]*K), grid
		subplot(pl+kk+2)
		plot(t, y(t,K), '-', n, y(n,K), 'o')
		titlef('Dirichlet K=%d', K)
		axis([0 4*K -1 1]), xtick([0:4]*K), grid
		ylabelf('K=%d', K)
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function b = cptsl ( n, d, e, b )

%*****************************************************************************80
%
%% CPTSL solves a Hermitian positive definite tridiagonal linear system.
%
%  Discussion;
%
%    The system does not have to be factored first.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2007
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
%    LINPACK User's Guide,
%    SIAM, 1979,
%    ISBN13: 978-0-898711-72-1,
%    LC: QA214.L56.
%
%  Parameters:
%
%    Input, integer N, the order of the matrix.
%
%    Input, complex D(N), the diagonal of the matrix.
%
%    Input, complex E(N), the superdiagonal of the matrix in E(1:N-1).
%
%    Input, complex B(N), the right hand side.
%
%    Output, complex B(N), the solution.
%

%
%  Check for 1 x 1 case.
%
  if ( n == 1 )
    b(1) = b(1) / d(1);
    return
  end

  nm1 = n - 1;
  nm1d2 = floor ( ( n - 1 ) / 2 );

  if ( n ~= 2 )

    kbm1 = n - 1;
%
%  Zero top half of subdiagonal and bottom half of superdiagonal.
%
    for k = 1 : nm1d2
      t1 = conj ( e(k) ) / d(k);
      d(k+1) = d(k+1) - t1 * e(k);
      b(k+1) = b(k+1) - t1 * b(k);
      t2 = e(kbm1) / d(kbm1+1);
      d(kbm1) = d(kbm1) - t2 * conj ( e(kbm1) );
      b(kbm1) = b(kbm1) - t2 * b(kbm1+1);
      kbm1 = kbm1 - 1;
    end

  end

  kp1 = nm1d2 + 1;
%
%  Clean up for possible 2 x 2 block at center.
%
  if ( mod ( n, 2 ) == 0 )
    t1 = conj ( e(kp1) ) / d(kp1);
    d(kp1+1) = d(kp1+1) - t1 * e(kp1);
    b(kp1+1) = b(kp1+1) - t1 * b(kp1);
    kp1 = kp1 + 1;
  end
%
%  Back solve starting at the center, going towards the top and bottom.
%
  b(kp1) = b(kp1) / d(kp1);

  if ( n ~= 2 )

    k = kp1 - 1;
    ke = kp1 + nm1d2 - 1;

    for kf = kp1 : ke
      b(k) = ( b(k) - e(k) * b(k+1) ) / d(k);
      b(kf+1) = ( b(kf+1) - conj ( e(kf) ) * b(kf) ) / d(kf+1);
      k = k - 1;
    end

  end

  if ( mod ( n, 2 ) == 0 )
    b(1) = ( b(1) - e(1) * b(2) ) / d(1);
  end

  return
end


text:
function y=sigpower(x)

% y=sigpower(x)
%
% Return the average signal power of signal x or mean(abs(x).^2);

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

error(nargchk(1,1,nargin));

y=sum(real(x).*real(x)+imag(x).*imag(x))/length(x);

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function test_issue1184

% MEM 8gb
% WALLTIME 00:30:00
% DEPENDENCY ft_resampledata

global ft_default
ft_default.toolbox.signal = 'matlab'; % do not use the compat version, its resample/firls is too slow
clear ft_defaults % ensure the path setup to be repeated

data = [];
data.time{1} = (1:228864) * 1/1200;
data.trial{1}(1,:) = 1:228864;
data.label = {'time'};
data.fsample = 1200;

cfg = [];
cfg.resamplefs = 100;
resampled = ft_resampledata(cfg, data);

assert(length(resampled.time{1})==19072);

%%

data11 = [];
data11.time{1} = (1:11) * 1/1200;
data11.trial{1}(1,:) = 1:11;
data11.label = {'time'};
data11.fsample = 1200;

cfg = [];
cfg.resamplefs = 100;
resampled11 = ft_resampledata(cfg, data11);

assert(isalmostequal(resampled11.time{1}(1), mean(data11.time{1}), 'abstol', 1e-6));

data12 = [];
data12.time{1} = (1:12) * 1/1200;
data12.trial{1}(1,:) = 1:12;
data12.label = {'time'};
data12.fsample = 1200;

cfg = [];
cfg.resamplefs = 100;
resampled12 = ft_resampledata(cfg, data12);

assert(isalmostequal(resampled12.time{1}(1), mean(data12.time{1}), 'abstol', 1e-6));

%%

% resampled11 and resampled12 have a time axis that is shifted by half a sample (at 1200 Hz)
% since they are constructed from two data structures that are different by one sample
%
% the following code would result in a data structure in which the trials are slightly mis-aligned
% resampledXX = ft_appenddata([], resampled11, resampled12);

% when resampling both trials at the same time, the time axes should remain consistent
dataXX = ft_appenddata([], data11, data12);

for resamplefs=1:1200
  cfg = [];
  cfg.showcallinfo = 'no';
  cfg.resamplefs = resamplefs;
  resampledXX = ft_resampledata(cfg, dataXX);

  % with a resampling rate of 600 or higher, the time axes of the resampled trials will have a different length due to the padding
  % only the first part should be compared
  n = min(length(resampledXX.time{1}), length(resampledXX.time{2}));
  assert(isalmostequal(mean(resampledXX.time{1}(1:n)), mean(resampledXX.time{2}(1:n)), 'abstol', 1e-6));
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [A,B,C,Ypredict] = sparse_bls(X,Y,nhidden,lambda1,lambda2,VERBOSE)

% SPARSE_BLS  Sparse orthogonalized partial least squares

% X: ninput x nsamples input data matrix
% Y: noutput x nsamples output data matrix
% nhidden: number of components [1]
% lambda1: l1 regularization parameter [1]
% lambda2: l2 regularization parameter [0.01]
%
% A: noutput x nhidden weight matrix
% B: ninput x nhidden weight matrix
% C: 1 x nhidden bias vector
% Ypredict: noutput x nsamples predictions

% Parse inputs

if nargin < 3,
    nhidden = 1;
end

if nargin < 4,
	lambda1 = 1;
end

if nargin < 5,
	lambda2 = 0.01;
end

ninput = size(X,1);
noutput = size(Y,1);
nsamples = size(X,2);
options = struct('offset',1,'maxiter',1e4,'tol',1e-2);

if nhidden > 1,
    % Run nhidden times sequentially

    R = Y;
    A = zeros(noutput,nhidden);
    B = zeros(ninput,nhidden);
    C = zeros(1,nhidden);
    for i=1:nhidden,
        [A(:,i),B(:,i),C(i)] = sparse_bls(X,R,1,lambda1,lambda2,VERBOSE);
        if i < nhidden,    % deflate
            Z = B(:,i)'*X + C(i); % hidden activations
            R = R - A(:,i)*Z; % Y activations after deflation
        end
        if VERBOSE,
            fprintf('done %d out of %d; sparsity %g\n',i,nhidden,length(find(B(:,i)))/ninput);
        end
    end
else

    % Run for a single hidden unit

    B = zeros(ninput,1);
    C = 0;
    iter = 0;
    maxiter = 1000;
    tol = nhidden*noutput*(1e-10);

    % Initialize A to first principal component of Y

    optseig.disp = 0;
    if nsamples < noutput,
        [d1,d2] = eigs(Y'*Y,[],1,'LM',optseig);
        A = Y*d1;
        A = A/sqrt(A'*A);
    else
        [A,d2] = eigs(Y*Y',[],1,'LM',optseig);
    end

%   A = randn(size(A));   % random
%   A = A/sqrt(A'*A);

    Aold = A;
    while iter < maxiter,

        if VERBOSE > 1,
            fprintf('   now starting iteration %d\n',iter+1);
        end

        Z = A'*Y;    % reconstruct Z given A from output Y

        % Use elastic net code to fit reconstructed Z given input X

        [B,C] = elastic(X,Z,lambda1,lambda2,options,B,C);

%         opts = glmnetSet;
%         opt.alpha = 1;
%         opts.lambda = lambda1;
%         res = glmnet(X',Z','gaussian',opts);
%         B = res.beta;
%         C = res.a0;
%         m = ft_mv_glmnet('validator',ft_mv_crossvalidator('nfolds',0.66,'metric','coefdet'),'alpha',1,'family','gaussian');
%         m = m.train(X',Z');
%         B = m.weights(1:(end-1));
%         C = m.weights(end);

        Z = B'*X + C;   % reconstruct Z given B and C

        % Find optimal A under constraint A'*A = 1

        Syz = Y*Z'/nsamples;
        denom = sqrt(Syz'*Syz);
        if denom,
            A = Syz/denom;
        else
            A = Aold;
        end
        if any(isnan(A(:))),    % check...
            error('nans!!!\n');
        end

        if sumsqr(A - Aold) < tol,
            if VERBOSE > 1,
                fprintf('   done!!\n',iter+1);
            end
		    iter = maxiter;
        else
	        iter = iter + 1;
            Aold = A;
        end
    end
end

% Parse outputs

if nargout > 3,
    Z = B'*X + C'*ones(1,nsamples);
    Ypredict = A*Z;
end


text:
% draw the frame id (1 based ) and GPS TOW timestamp onto the image
function drawTimeOntoVideo()
close all;
resdir='..\20130808\';
sequencePath =[resdir, 'casio2_3430.MOV'];
imgtimefile=[resdir, 'casio2timestamps.txt'];
initIm =12464; % 415530 GTOW
endIm=12986;
[fimgtime, imgepoch, initIm]=readimgtimeheader(imgtimefile,0, initIm);
hstream=sprintf('%d,%.5f', initIm, imgepoch);
myVid= VideoReader(sequencePath);
figure(1);
tic;
for x = initIm : endIm
    imRGB=read(myVid,x); %78612;
    imwrite(imRGB,strcat('F:\relaylatest\20130808\garage\casio2_',num2str(x),'.jpg'));
    %     if(isempty(imRGB)||imgepoch>416033.2)
    %         break;
    %     else
    %display frmId and timestamp on the image
    %         text(xlbl(:), ylbl(:), lbl(:),'color','w',...
    %     'HorizontalAlignment','center','VerticalAlignment','middle');
    imshow(imRGB);
    text(50, 50, hstream);
    pause(0.02);
    % read in next image epoch
    hstream= fgetl(fimgtime);
    %         if(hstream==-1)
    %             imgepoch=inf;
    %         end
end
toc;
end

repo: ekfmonoslam


text:
function [v,y,w]=nearnonz(x,d)
%NEARNONZ replace each zero element with the nearest non-zero element [V,Y,W]=nearnonz(X,D)
%
%  Inputs:  x         input vector, matrix or larger array
%           d         dimension to apply filter along [default 1st non-singleton]
%
% Outputs:  v         v is the same size as x but with each zero entry replaced by
%                     the nearest non-zero value along dimension d
%                     elements equidistant from two non-zero entries will be taken
%                     from the higher index
%           y         y is the same size as x and gives the index along dimension d
%                     from which the corresponding entry in v was taken
%                     If there are no non-zero entries, then the corresponding
%                     elements of y will be zero.
%           w         w is the same size as x and gives the distance (+ or -) to the
%                     nearest non-zero entry in x

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: nearnonz.m 713 2011-10-16 14:45:43Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
e=size(x);
p=prod(e);
if nargin<2             % if no dimension given, find the first non-singleton
    d=find(e>1,1);
    if ~numel(d)
        d=1;
    end
end
k=e(d);                 % size of active dimension
q=p/k;                  % size of remainder
if d==1
    z=reshape(x,k,q);
else
    z=shiftdim(x,d-1);
    r=size(z);
    z=reshape(z,k,q);
end
xx=z~=0;
cx=cumsum(xx);
[i,j]=find(z);
qq=cx(xx);
pos=full(sparse(qq,j,i,k,q)); % list the positions of non-zero elements in each column
mp=ceil((pos(1:end-1,:)+pos(2:end,:))*0.5); % find the mid point between consecutive non-zero elements
[i2,j2]=find(pos(2:end,:)>0);
zz=1+cumsum(full(sparse(mp(pos(2:end,:)>0),j2,1,k,q)));
y=pos(zz+repmat((0:q-1)*k,k,1));
v=z(max(y,1)+repmat((0:q-1)*k,k,1));
w=y-repmat((1:k)',1,q);
w(y==0)=0;
if d==1
    y=reshape(y,e);
    v=reshape(v,e);
    w=reshape(w,e);
else
    y=shiftdim(reshape(y,r),length(e)+1-d);
    v=shiftdim(reshape(v,r),length(e)+1-d);
    w=shiftdim(reshape(w,r),length(e)+1-d);
end

author: jtkim-kaist

repo: VAD

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


text:
function [newnode,newelem]=highordertet(node,elem,order)
%
% [newnode,newelem]=highordertet(node,elem)
%
% generate high-order straight-edge tetrahedral mesh from
% the 1st order tetrahedral mesh
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%    node: list of nodes
%    elem: list of elements (each row are indices of nodes of each element)
%    order: optional, the order of the generated mesh; if missing, order=2
%
% output:
%    newnode: all new edge-nodes on the output mesh
%    newelem: the indices of the edge nodes for each original tet element
%
%    currently, this function only supports order=2
%    to combine the newnode/newelem with the old mesh, one should use
%
%    elemfull=[elem(:,1:4) newelem+size(node,1)]; % 10-node element
%    nodefull=[node;newnode];
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if(nargin<3)
    order=2;
end
if(order>=3 || order<=1)
    error('currently this function only supports order=2');
end
[edges,idx,newelem]=uniqedges(elem(:,1:4));
newnode=node(edges',1:3);
newnode=reshape(newnode',[3,2,size(edges,1)]);
newnode=squeeze(mean(permute(newnode,[3 2 1]),2));


text:
function d =  testing(a,d)

  eval([a.func ';']);
  d=set_name(d,[get_name(d) ' -> ' get_name(a)]);


text:
function [nstate] = tapas_sampler_mixedlinear_gibbs_node(data, model, ...
    inference, state, node)
%% Samples from a linear multivariate model with fixed and random effects
% and a diagonal covariance matrix.
%
% Input
%
% Output
%
% This level assumes only and unknown mean with known variance for some the
% the elements. The input matrix is required to be simple a weighting of
% prior and mean.
%

% aponteeduardo@gmail.com
% copyright (C) 2016
%

nstate = state;

% Number of subjects and number of chains
[ns, nc] = size(state.graph{node - 1}.y);

% Number of parameters
np = size(state.graph{node}.y{1}.mu, 2);

% Number of regressors
nr = size(state.graph{node}.y{1}.mu, 1);

% First store the values somewhere

for i = 1:nc
    % Get all the parameters at a given temperature
    iy = state.graph{node - 1}.y{i}.mu;

    % Get the grouping of the random variables.
    x = state.graph{node - 1}.u.x;

    % The number of groups is equal to size(x, 2). The ones in a column
    % encode membership to a code. size(x, 1) should be np;

    % Number of groups
    ng = size(x, 2);

    for j = 1:ng
        % Group
        gv = logical(x(:, j));
        ngv = sum(gv);
        % Group mean
        gmu = sum(iy(gv, :), 1) .* state.graph{node - 1}.y{i}.pe;
        gmu = gmu + state.graph{node + 1}.y{i}.pe(j, :) .* ...
            state.graph{node + 1}.y{i}.mu(j, :); % 1 x np

        % Group precision
        gpe = ngv * state.graph{node - 1}.y{i}.pe + ...
            state.graph{node + 1}.y{i}.pe(j, :);

        % Sample using Gibbs
        gmu = gmu./gpe;
        gmu = gmu + sqrt(1./ gpe) .* randn(1, np);

        % Replicate for all the correspondign regressors.
        nstate.graph{node}.y{i}.mu(gv, :) = repmat(gmu, ngv, 1);
    end
end

% Don't need it
%nstate.llh{node} = model.graph{node}.llh(nstate.graph{node}, ...
%    nstate.graph{node + 1}, model.graph{node}.htheta);
%
%nstate.llh{node - 1} = model.graph{node - 1}.llh(...
%    nstate.graph{node - 1}, nstate.graph{node}, model.graph{node - 1}.htheta);

end

author: translationalneuromodeling

repo: tapas


text:
function check = triangle_check ( a, b, c )

%*****************************************************************************80
%
%% TRIANGLE_CHECK checks the parameters of the Triangle CDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, C, the parameters of the PDF.
%    A <= B <= C and A < C.
%
%    Output, logical TRIANGLE_CHECK, is true if the parameters are legal.
%
  if ( b < a )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGLE_CHECK - Fatal error!\n' );
    fprintf ( 1, '  B < A.\n' );
    check = 0;
    return
  end

  if ( c < b )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGLE_CHECK - Fatal error!\n' );
    fprintf ( 1, '  C < B.\n' );
    check = 0;
    return
  end

  if ( a == c )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TRIANGLE_CHECK - Fatal error!\n' );
    fprintf ( 1, '  A == C.\n' );
    check = 0;
    return
  end

  check = 1;

  return
end


text:
function linplus_test153 ( )

%*****************************************************************************80
%
%% TEST153 tests R8BLT_MXV, R8BLT_PRINT, R8BLT_RANDOM, R8BLT_SL, R8BLT_VXM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 March 2009
%
%  Author:
%
%    John Burkardt
%
  ml = 3;
  n = 10;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST153\n' );
  fprintf ( 1, '  For a band matrix in lower triangular storage,\n' );
  fprintf ( 1, '  R8BLT_RANDOM sets a random value;\n' );
  fprintf ( 1, '  R8BLT_SL solves systems;\n' );
  fprintf ( 1, '  R8BLT_MXV computes matrix-vector products;\n' );
  fprintf ( 1, '  R8BLT_VXM computes vector-matrix products;\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
  fprintf ( 1, '  Lower bandwidth ML = %d\n', ml );

  [ a, seed ] = r8blt_random ( n, ml, seed );

  r8blt_print ( n, ml, a, '  The R8BLT matrix:' );

  for job = 0 : 1
%
%  Set the desired solution.
%
    x = r8vec_indicator ( n );
%
%  Compute the corresponding right hand side.
%
    if ( job == 0 )
      b = r8blt_mxv ( n, ml, a, x );
    else
      b = r8blt_vxm ( n, ml, a, x );
    end

    r8vec_print ( n, b, '  The right hand side:' );
%
%  Solve the linear system.
%
    x = r8blt_sl ( n, ml, a, b, job );

    if ( job == 0 )
      r8vec_print ( n, x, '  Solution to the untransposed system:' );
    else
      r8vec_print ( n, x, '  Solution to the transposed system:' );
    end

  end

  return
end