text:
%% DEMO_febio_0069_must_points_export
% Below is a demonstration for:
%
% * Building geometry for a cube with hexahedral elements
% * Defining the boundary conditions
% * Coding the febio structure
% * Running the model
% * Importing and visualizing the displacement and stress results

%% Keywords
%
% * febio_spec version 2.5
% * febio, FEBio
% * uniaxial loading
% * compression, tension, compressive, tensile
% * displacement control, displacement boundary condition
% * hexahedral elements, hex8
% * cube, box, rectangular
% * static, solid
% * hyperelastic, Ogden
% * displacement logfile
% * stress logfile

%%

clear; close all; clc;

%% Plot settings
fontSize=20;
faceAlpha1=0.8;
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=fullfile(savePath,[febioFebFileNamePart,'.txt']); %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement
febioLogFileName_stress=[febioFebFileNamePart,'_stress_out.txt']; %Log file name for exporting stress

%Specifying dimensions and number of elements
cubeSize=10;
sampleWidth=cubeSize; %Width
sampleThickness=cubeSize; %Thickness
sampleHeight=cubeSize; %Height
pointSpacings=2*ones(1,3); %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

%Define applied displacement
appliedStrain=0.4; %Linear strain (Only used to compute applied stretch)
loadingOption='tension'; % or 'tension'
switch loadingOption
    case 'compression'
        stretchLoad=1-appliedStrain; %The applied stretch for uniaxial loading
    case 'tension'
        stretchLoad=1+appliedStrain; %The applied stretch for uniaxial loading
end
displacementMagnitude=(stretchLoad*sampleHeight)-sampleHeight; %The displacement magnitude

%Material parameter set
c1=0.7; %Shear-modulus-like parameter
m1=2; %Material parameter setting degree of non-linearity
k_factor=500; %Bulk modulus factor
k=c1*k_factor; %Bulk modulus
formulationType='uncoupled'; %coupled

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

%Define must points
numStepsOutput=11;
t_load_curve=linspace(0,1,numStepsOutput);
dtmax_load_curve=dtmax.*ones(size(t_load_curve));

%% 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
cubeDimensions=[sampleWidth sampleThickness sampleHeight]; %Dimensions
cubeElementNumbers=[numElementsWidth numElementsThickness numElementsHeight]; %Number of elements
outputStructType=2; %A structure compatible with mesh view
[meshStruct]=hexMeshBox(cubeDimensions,cubeElementNumbers,outputStructType);

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

%%
% 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(Fb,V,Cb,'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];
meshView(meshStruct,optionStruct);
axisGeom(gca,fontSize);

drawnow;

%% Defining the boundary conditions
% The visualization of the model boundary shows colors for each side of the
% cube. These labels can be used to define boundary conditions.

%Define supported node sets
logicFace=Cb==1; %Logic for current face set
Fr=Fb(logicFace,:); %The current face set
bcSupportList_X=unique(Fr(:)); %Node set part of selected face

logicFace=Cb==3; %Logic for current face set
Fr=Fb(logicFace,:); %The current face set
bcSupportList_Y=unique(Fr(:)); %Node set part of selected face

logicFace=Cb==5; %Logic for current face set
Fr=Fb(logicFace,:); %The current face set
bcSupportList_Z=unique(Fr(:)); %Node set part of selected face

%Prescribed displacement nodes
logicPrescribe=Cb==6; %Logic for current face set
Fr=Fb(logicPrescribe,:); %The current face set
bcPrescribeList=unique(Fr(:)); %Node set part of selected face

%%
% Visualizing boundary conditions. Markers plotted on the semi-transparent
% model denote the nodes in the various boundary condition lists.

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

gpatch(Fb,V,'kw','k',0.5);

hl(1)=plotV(V(bcSupportList_X,:),'r.','MarkerSize',markerSize);
hl(2)=plotV(V(bcSupportList_Y,:),'g.','MarkerSize',markerSize);
hl(3)=plotV(V(bcSupportList_Z,:),'b.','MarkerSize',markerSize);
hl(4)=plotV(V(bcPrescribeList,:),'k.','MarkerSize',markerSize);

legend(hl,{'BC x support','BC y support','BC z support','BC z prescribe'});

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='2.5';

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

%Control section
febio_spec.Control.time_stepper=rmfield(febio_spec.Control.time_stepper,'dtmax');

febio_spec.Control.analysis.ATTR.type='static';
febio_spec.Control.title='Cube analysis';
febio_spec.Control.time_steps=numTimeStepsAnalysis;
febio_spec.Control.step_size=1/numTimeStepsAnalysis;
febio_spec.Control.time_stepper.dtmin=dtmin;
febio_spec.Control.time_stepper.dtmax.VAL=dtmax;
febio_spec.Control.time_stepper.dtmax.ATTR.lc=2;
febio_spec.Control.time_stepper.max_retries=max_retries;
febio_spec.Control.time_stepper.opt_iter=opt_iter;
febio_spec.Control.max_refs=max_refs;
febio_spec.Control.max_ups=max_ups;

febio_spec.Control.plot_level='PLOT_MUST_POINTS';
febio_spec.Control.output_level='OUTPUT_MUST_POINTS';

%Material section
switch formulationType
    case 'coupled'
        febio_spec.Material.material{1}.ATTR.type='Ogden unconstrained';
        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}.cp=k;
    case 'uncoupled'
        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;
end

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

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

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

febio_spec.Geometry.NodeSet{2}.ATTR.name='bcSupportList_Y';
febio_spec.Geometry.NodeSet{2}.node.ATTR.id=bcSupportList_Y(:);

febio_spec.Geometry.NodeSet{3}.ATTR.name='bcSupportList_Z';
febio_spec.Geometry.NodeSet{3}.node.ATTR.id=bcSupportList_Z(:);

febio_spec.Geometry.NodeSet{4}.ATTR.name='bcPrescribeList';
febio_spec.Geometry.NodeSet{4}.node.ATTR.id=bcPrescribeList(:);

%Boundary condition section
% -> Fix boundary conditions
febio_spec.Boundary.fix{1}.ATTR.bc='x';
febio_spec.Boundary.fix{1}.ATTR.node_set=febio_spec.Geometry.NodeSet{1}.ATTR.name;
febio_spec.Boundary.fix{2}.ATTR.bc='y';
febio_spec.Boundary.fix{2}.ATTR.node_set=febio_spec.Geometry.NodeSet{2}.ATTR.name;
febio_spec.Boundary.fix{3}.ATTR.bc='z';
febio_spec.Boundary.fix{3}.ATTR.node_set=febio_spec.Geometry.NodeSet{3}.ATTR.name;

% -> Prescribe boundary conditions
febio_spec.Boundary.prescribe{1}.ATTR.bc='z';
febio_spec.Boundary.prescribe{1}.ATTR.node_set=febio_spec.Geometry.NodeSet{4}.ATTR.name;
febio_spec.Boundary.prescribe{1}.scale.ATTR.lc=1;
febio_spec.Boundary.prescribe{1}.scale.VAL=1;
febio_spec.Boundary.prescribe{1}.relative=1;
febio_spec.Boundary.prescribe{1}.value=displacementMagnitude;

%LoadData section
febio_spec.LoadData.loadcurve{1}.ATTR.id=1;
febio_spec.LoadData.loadcurve{1}.ATTR.type='linear';
febio_spec.LoadData.loadcurve{1}.point.VAL=[0 0; 1 1];

febio_spec.LoadData.loadcurve{2}.ATTR.id=2;
febio_spec.LoadData.loadcurve{2}.ATTR.type='linear';
febio_spec.LoadData.loadcurve{2}.point.VAL=[t_load_curve(:) dtmax_load_curve(:)];

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

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

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

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

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

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

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

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.disp_log_on=1; %Display convergence information in the command window
febioAnalysis.runMode='external';%'internal';
febioAnalysis.t_check=0.25; %Time for checking log file (dont set too small)
febioAnalysis.maxtpi=1e99; %Max analysis time
febioAnalysis.maxLogCheckTime=10; %Max log file checking time

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

%% Import FEBio results

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

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

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

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

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

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

    % Create basic view and store graphics handle to initiate animation
    hf=cFigure; %Open figure
    gtitle([febioFebFileNamePart,': Press play to animate']);
    title('Displacement magnitude [mm]','Interpreter','Latex')
    hp=gpatch(Fb,V_DEF(:,:,end),DN_magnitude,'k',1); %Add graphics object to animate
    hp.Marker='.';
    hp.MarkerSize=markerSize2;
    hp.FaceColor='interp';
    gpatch(Fb,V,0.5*ones(1,3),'k',0.25); %A static graphics object

    axisGeom(gca,fontSize);
    colormap(gjet(250)); colorbar;
    caxis([0 max(DN_magnitude)]); caxis manual;
    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
        DN_magnitude=sqrt(sum(N_disp_mat(:,:,qt).^2,2)); %Current displacement magnitude

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

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

    %Access data
    E_stress_mat=dataStruct.data;

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

    [CV]=faceToVertexMeasure(E,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_{zz}$ [MPa]','Interpreter','Latex')
    hp=gpatch(Fb,V_DEF(:,:,end),CV,'k',1); %Add graphics object to animate
    hp.Marker='.';
    hp.MarkerSize=markerSize2;
    hp.FaceColor='interp';
    gpatch(Fb,V,0.5*ones(1,3),'k',0.25); %A static graphics object

    axisGeom(gca,fontSize);
    colormap(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(E,V,E_stress_mat(:,:,qt));

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

    %%
    % Calculate metrics to visualize stretch-stress curve

    DZ_set=N_disp_mat(bcPrescribeList,end,:); %Z displacements of the prescribed set
    DZ_set=mean(DZ_set,1); %Calculate mean Z displacements across nodes
    stretch_sim=(DZ_set(:)+sampleHeight)./sampleHeight; %Derive stretch
    stress_cauchy_sim=mean(squeeze(E_stress_mat(:,end,:)),1)';

    %%
    % Visualize stress-stretch curve

    cFigure; hold on;
    title('Uniaxial stress-stretch curve','FontSize',fontSize);
    xlabel('$\lambda$ [.]','FontSize',fontSize,'Interpreter','Latex');
    ylabel('$\sigma_{zz}$ [MPa]','FontSize',fontSize,'Interpreter','Latex');

    plot(stretch_sim(:),stress_cauchy_sim(:),'r-','lineWidth',lineWidth);

    view(2); axis tight;  grid on; axis square; box on;
    set(gca,'FontSize',fontSize);
    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


text:
function FatalInitError(str)
%
% function FatalInitError(str)
%
% str: error string
%
% queries user, then blasts the entire thing, deleting mrSESSION and any
% tSeries files
% have been extracted.
%
% djh, 9/26/2001
global HOMEDIR
deleteFlag = questdlg([str,' Start from scratch?'],...
    'Fatal mrInitRet error','Yes','No','Yes');

if strcmp(deleteFlag,'Yes')
    % Delete everything
    sessionFile = fullfile(HOMEDIR,'mrSESSION.mat');
    if exist(sessionFile,'file');
        delete(sessionFile);
    end
    datadir = fullfile(HOMEDIR,'Inplane');
    % Delete anat, if it exists
    delete(fullfile(datadir,'*.mat'));
    % Delete tSeries (if there are any)
    [nscans,scanDirList] = countDirs(fullfile(datadir,'Original','TSeries','Scan*'));
    for s=1:nscans
        delete(fullfile(datadir,'TSeries',scanDirList{s},'*.mat'));
    end
end
clear all
close all

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [doseBinsV, volsHistV] = loadIVHMatrix(IVHNum, planC)
%"loadIVHMatrix"
%   Return the doseBinsV and volsHistV with the binWidth shift included.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.
%
%Usage:
%   function [doseBinsV, volsHistV] = loadIVHMatrix(IVHNum, planC)

indexS = planC{end};

%Calculate width of bins, use to find middle of each bin.
binWidthsV = diff(planC{indexS.IVH}(IVHNum).IVHMatrix(:,1));
lastBinWidth = binWidthsV(end);
binWidthsV(end+1,1) = lastBinWidth;

%Extract dose bin values, adding half to binwidth to get middle.
doseBinsV = planC{indexS.IVH}(IVHNum).IVHMatrix(:,1) + binWidthsV/2;
volsHistV = planC{indexS.IVH}(IVHNum).IVHMatrix(:,2);
return;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [cnmfeAnalysisOutput] = computeCnmfeSignalExtraction(inputMovie,varargin)
    % Wrapper function for CNMF-E, update for most recent versions.
    % Biafra Ahanonu
    % started: 2018.10.20 [16:38:24]
    % Building off of demo_large_data_1p.m in CNMF-E github repo
    % Most recent commit tested on: https://github.com/epnev/ca_source_extraction/commit/187bbdbe66bca466b83b81861b5601891a95b8d1
    % https://github.com/epnev/ca_source_extraction/blob/master/demo_script_class.m
    % inputs
        % inputMovie - a string or a cell array of strings pointing to the movies to be analyzed (recommended).
        % numExpectedComponents - number of expected components
    % outputs
        % cnmfAnalysisOutput - structure containing extractedImages and extractedSignals along with input parameters to the algorithm
    % READ BEFORE RUNNING
        % Get CVX from http://cvxr.com/cvx/doc/install.html
        % Run the below commands in Matlab after unzipping
        % cvx_setup
        % cvx_save_prefs (permanently stores settings)

    % changelog
        % 2016.06.20 - updated to keep in line with recent changes to CNMF functions
        % 2019.04.04 [11:07:21] - Updated to add many of the options to the varargin options structure for function for easier access in parent functions
        % 2021.08.08 [19:30:20] - Updated to handle CIAtah v4.0 switch to all functions inside ciapkg package.
    % TODO
        %

    import ciapkg.api.* % import CIAtah functions in ciapkg package API.

    % ========================
    % OVERALL
    % turn on parallel
    options.nonCNMF.parallel = 1;
    % Binary: 1 = run merging algorithms
    options.runMerge = 1;
    % Binary: 1 = remove false positives using CNMF-E algorithm
    options.runRemoveFalsePositives = 1;
    % ===COMPUTATION
    % Float: GB, memory space you allow to use in MATLAB
    options.memory_size_to_use = 32; %
    % Float: GB, memory space you allow to use in MATLAB
    options.memory_size_per_patch = 1.2; % 0.6
    % Int vector: patch size in pixels
    options.patch_dims = [128, 128]; % [64, 64]
    % ===SPATIAL
    % Int: pixel, gaussian width of a gaussian kernel for filtering the data. 0 means no filtering
    options.gSig = 3;
    % Int: pixel, neuron diameter
    options.gSiz = 11;
    % Int: spatial downsampling factor
    options.ssub = 1;
    % Binary: movie has dendrites?
    options.with_dendrites = true;
    % Int: expand kernel for HALS growing (default: 3) and expansion factor of ellipse (default: 3)
    options.updateA_bSiz = 5;
    % Char: hals, hals_thresh, lars, nnls
    options.spatial_algorithm = 'hals_thresh';
    % ===TEMPORAL
    % Int: temporal downsampling factor
    options.tsub = 1;
    % Int: frame rate
    options.Fs = 10;
    % Float: minimum spike size. When the value is negative, the actual threshold is abs(smin)*noise level
    options.deconv_smin = -5;
    % Int: maximum decay time (unit: frame);
    options.max_tau = 100;
    % Int: detrending the slow fluctuation. usually 1 is fine (no detrending)
    options.nk = 3;
    % ===BACKGROUND
    % Char: model of the background {'ring', 'svd'(default), 'nmf'}
    options.bg_model = 'ring';
    % Int: number of background sources for each patch (only be used in SVD and NMF model)
    options.nb = 1;
    % Int: when the ring model used, it is the radius of the ring used in the background model. otherwise, it's just the width of the overlapping area
    options.ring_radius = 18;
    % Int: downsample background for a faster speed
    options.bg_ssub = 1;
    % ===MERGING
    % Float: 0 to 1, thresholds for merging neurons; [spatial overlap ratio, temporal correlation of calcium traces, spike correlation]
    options.merge_thr = 0.65;
    % Char: method for computing neuron distances {'mean', 'max'}
    options.method_dist = 'max';
    % Int: minimum distances between two neurons. it is used together with merge_thr
    options.dmin = 5;
    % Int: merge neurons if their distances are smaller than dmin_only.
    options.dmin_only = 2;
    % Float vector: merge components with highly correlated spatial shapes (corr=0.8) and small temporal correlations (corr=0.1)
    options.merge_thr_spatial = [0.8, 0.4, -inf];
    % ===INITIALIZATION
    % Int: maximum number of neurons per patch. when K=[], take as many as possible.
    options.K = [];
    % Float: minimum local correlation for a seeding pixel
    options.min_corr = 0.8;
    % minimum peak-to-noise ratio for a seeding pixel
    options.min_pnr = 8;
    options.bd = 0;             % number of rows/columns to be ignored in the boundary (mainly for motion corrected data)
    options.use_parallel = true;    % use parallel computation for parallel computing
    options.show_init = true;   % show initialization results
    options.center_psf = true;  % set the value as true when the background fluctuation is large (usually 1p data)
    % set the value as false when the background fluctuation is small (2p)
    % ===Residual
    % Float: 0 to 1, minimum local correlation for initializing a neuron (default: 0.3)
    options.min_corr_res = 0.7;
    % Float: stands for minimum peak-to-noise ratio to look for a cell
    options.min_pnr_res = 6;
    % Char: method for initializing neurons from the residual. 'auto' or 'manual'
    options.seed_method_res = 'auto';
    % Binary: boolean, update noise level for each pixel
    options.update_sn = true;

    % get options
    options = getOptions(options,varargin);
    % ========================
    options

    % if cvx is not in the path, ask user for file
    runCvxSetup();

    % Make sure consistent
    options.bg_ssub = options.ssub;

    %% clear the workspace and select data
    % clear; clc; close all;

    %% choose data
    inputFilename = inputMovie;
    neuron = Sources2D();

    nam = get_fullname(inputFilename);          % this demo data is very small, here we just use it as an example
    nam = neuron.select_data(nam);  %if nam is [], then select data interactively

    % nams = {inputFilename};          % you can put all file names into a cell array; when it's empty, manually select files
    % nams = neuron.select_multiple_files(nams);  %if nam is [], then select data interactively

    %% parameters
    % -------------------------    COMPUTATION    -------------------------  %
    pars_envs = struct('memory_size_to_use', options.memory_size_to_use, ...   % GB, memory space you allow to use in MATLAB
        'memory_size_per_patch', options.memory_size_per_patch, ...   % GB, space for loading data within one patch
        'patch_dims', options.patch_dims);  %GB, patch size
    % pars_envs = struct('memory_size_to_use', 8, ...   % GB, memory space you allow to use in MATLAB
    %     'memory_size_per_patch', 0.5, ...   % GB, space for loading data within one patch
    %     'patch_dims', [64, 64],...  %GB, patch size
    %     'batch_frames', 1000);           % number of frames per batch

    % -------------------------      SPATIAL      -------------------------  %
    gSig = options.gSig;           % pixel, gaussian width of a gaussian kernel for filtering the data. 0 means no filtering
    gSiz = options.gSiz;          % pixel, neuron diameter
    ssub = options.ssub;           % spatial downsampling factor
    with_dendrites = options.with_dendrites;   % with dendrites or not
    if with_dendrites
        % determine the search locations by dilating the current neuron shapes
        updateA_search_method = 'dilate';  %#ok<UNRCH>
        updateA_bSiz = options.updateA_bSiz;
        updateA_dist = neuron.options.dist;
    else
        % determine the search locations by selecting a round area
        updateA_search_method = 'ellipse'; %#ok<UNRCH>
        updateA_dist = options.updateA_bSiz;
        updateA_bSiz = neuron.options.dist;
    end
    spatial_constraints = struct('connected', true, 'circular', false);  % you can include following constraints: 'circular'
    spatial_algorithm = options.spatial_algorithm;

    % -------------------------      TEMPORAL     -------------------------  %
    Fs = options.Fs;             % frame rate
    tsub = options.tsub;           % temporal downsampling factor
    deconv_options = struct('type', 'ar1', ... % model of the calcium traces. {'ar1', 'ar2'}
        'method', 'foopsi', ... % method for running deconvolution {'foopsi', 'constrained', 'thresholded'}
        'smin', options.deconv_smin, ...         % minimum spike size. When the value is negative, the actual threshold is abs(smin)*noise level
        'optimize_pars', true, ...  % optimize AR coefficients
        'optimize_b', true, ...% optimize the baseline);
        'max_tau', options.max_tau);    % maximum decay time (unit: frame);

    nk = 3;             % detrending the slow fluctuation. usually 1 is fine (no detrending)
    % when changed, try some integers smaller than total_frame/(Fs*30)
    detrend_method = 'spline';  % compute the local minimum as an estimation of trend. method for detrending {'spline', 'local_min'}

    % -------------------------     BACKGROUND    -------------------------  %
    bg_model = options.bg_model;  % model of the background {'ring', 'svd'(default), 'nmf'}
    nb = options.nb;             % number of background sources for each patch (only be used in SVD and NMF model)
    ring_radius = options.ring_radius;  % when the ring model used, it is the radius of the ring used in the background model.
    %otherwise, it's just the width of the overlapping area
    num_neighbors = []; % number of neighbors for each neuron
    bg_ssub = options.bg_ssub;        % downsample background for a faster speed

    % -------------------------      MERGING      -------------------------  %
    show_merge = false;  % if true, manually verify the merging step
    merge_thr = options.merge_thr;     % thresholds for merging neurons; [spatial overlap ratio, temporal correlation of calcium traces, spike correlation]
    method_dist = options.method_dist;   % method for computing neuron distances {'mean', 'max'}
    dmin = options.dmin;       % minimum distances between two neurons. it is used together with merge_thr
    dmin_only = options.dmin_only;  % merge neurons if their distances are smaller than dmin_only.
    merge_thr_spatial = options.merge_thr_spatial;  % merge components with highly correlated spatial shapes (corr=0.8) and small temporal correlations (corr=0.1)

    % -------------------------  INITIALIZATION   -------------------------  %
    K = options.K;             % maximum number of neurons per patch. when K=[], take as many as possible.
    min_corr = options.min_corr;     % minimum local correlation for a seeding pixel
    min_pnr = options.min_pnr;       % minimum peak-to-noise ratio for a seeding pixel
    min_pixel = gSig^2;      % minimum number of nonzero pixels for each neuron
    bd = options.bd;             % number of rows/columns to be ignored in the boundary (mainly for motion corrected data)
    frame_range = [];   % when [], uses all frames
    save_initialization = false;    % save the initialization procedure as a video.
    use_parallel = options.use_parallel;    % use parallel computation for parallel computing
    show_init = options.show_init;   % show initialization results
    choose_params = false; % manually choose parameters
    center_psf = options.center_psf;  % set the value as true when the background fluctuation is large (usually 1p data)
    % set the value as false when the background fluctuation is small (2p)

    % -------------------------  Residual   -------------------------  %
    min_corr_res = options.min_corr_res;
    min_pnr_res = options.min_pnr_res;
    seed_method_res = options.seed_method_res;  % method for initializing neurons from the residual
    update_sn = options.update_sn;

    % ----------------------  WITH MANUAL INTERVENTION  --------------------  %
    with_manual_intervention = false;

    % -------------------------  FINAL RESULTS   -------------------------  %
    save_demixed = true;    % save the demixed file or not
    kt = 3;                 % frame intervals

    % -------------------------    UPDATE ALL    -------------------------  %
    neuron.updateParams(...
        'save_intermediate',false,...
        'gSig', gSig, ...       % -------- spatial --------
        'gSiz', gSiz, ...
        'ring_radius', ring_radius, ...
        'ssub', ssub, ...
        'search_method', updateA_search_method, ...
        'bSiz', updateA_bSiz, ...
        'dist', updateA_bSiz, ...
        'spatial_constraints', spatial_constraints, ...
        'spatial_algorithm', spatial_algorithm, ...
        'tsub', tsub, ...                       % -------- temporal --------
        'deconv_options', deconv_options, ...
        'nk', nk, ...
        'detrend_method', detrend_method, ...
        'background_model', bg_model, ...       % -------- background --------
        'nb', nb, ...
        'ring_radius', ring_radius, ...
        'num_neighbors', num_neighbors, ...
        'bg_ssub', bg_ssub, ...
        'merge_thr', merge_thr, ...             % -------- merging ---------
        'dmin', dmin, ...
        'method_dist', method_dist, ...
        'min_corr', min_corr, ...               % ----- initialization -----
        'min_pnr', min_pnr, ...
        'min_pixel', min_pixel, ...
        'bd', bd, ...
        'center_psf', center_psf);
    neuron.Fs = Fs;

    %% distribute data and be ready to run source extraction
    neuron.getReady(pars_envs);
    % neuron.getReady_batch(pars_envs);

    %% initialize neurons from the video data within a selected temporal range
    if choose_params
        % change parameters for optimized initialization
        [gSig, gSiz, ring_radius, min_corr, min_pnr] = neuron.set_parameters();
    end

    [center, Cn, PNR] = neuron.initComponents_parallel(K, frame_range, save_initialization, use_parallel);
    neuron.compactSpatial();
    if show_init
        figure();
        ax_init= axes();
        imagesc(Cn, [0, 1]); colormap gray;
        hold on;
        plot(center(:, 2), center(:, 1), '.r', 'markersize', 10);
    end
    subfxnDisplayState();

    %% estimate the background components
    neuron.update_background_parallel(use_parallel);
    neuron_init = neuron.copy();
    subfxnDisplayState();

    if options.runMerge==1
        %%  merge neurons and update spatial/temporal components
        neuron.merge_neurons_dist_corr(show_merge);
        neuron.merge_high_corr(show_merge, merge_thr_spatial);
        subfxnDisplayState();
    end

    %% update spatial components

    %% pick neurons from the residual
    [center_res, Cn_res, PNR_res] = neuron.initComponents_residual_parallel([], save_initialization, use_parallel, min_corr_res, min_pnr_res, seed_method_res);
    if show_init
        axes(ax_init);
        plot(center_res(:, 2), center_res(:, 1), '.g', 'markersize', 10);
    end
    neuron_init_res = neuron.copy();
    subfxnDisplayState();

    %% udpate spatial&temporal components, delete false positives and merge neurons
    % update spatial
    if update_sn
        neuron.update_spatial_parallel(use_parallel, true);
        udpate_sn = false;
    else
        neuron.update_spatial_parallel(use_parallel);
    end
    subfxnDisplayState();
    if options.runMerge==1
        % merge neurons based on correlations
        neuron.merge_high_corr(show_merge, merge_thr_spatial);
        subfxnDisplayState();
    end

    for m=1:2
        % update temporal
        neuron.update_temporal_parallel(use_parallel);

        if options.runRemoveFalsePositives==1
            % delete bad neurons
            neuron.remove_false_positives();
        end
        subfxnDisplayState();
        if options.runMerge==1
            % merge neurons based on temporal correlation + distances
            neuron.merge_neurons_dist_corr(show_merge);
            subfxnDisplayState();
        end
    end

    %% add a manual intervention and run the whole procedure for a second time
    neuron.options.spatial_algorithm = 'nnls';
    if with_manual_intervention
        show_merge = true;
        neuron.orderROIs('snr');   % order neurons in different ways {'snr', 'decay_time', 'mean', 'circularity'}
        neuron.viewNeurons([], neuron.C_raw);
        if options.runMerge==1
            % merge closeby neurons
            neuron.merge_close_neighbors(true, dmin_only);
        end

        % delete neurons
        tags = neuron.tag_neurons_parallel();  % find neurons with fewer nonzero pixels than min_pixel and silent calcium transients
        ids = find(tags>0);
        if ~isempty(ids)
            neuron.viewNeurons(ids, neuron.C_raw);
        end
        subfxnDisplayState();
    end
    %% run more iterations
    neuron.update_background_parallel(use_parallel);
    neuron.update_spatial_parallel(use_parallel);
    neuron.update_temporal_parallel(use_parallel);
    subfxnDisplayState();

    K = size(neuron.A,2);
    % find neurons with fewer nonzero pixels than min_pixel and silent calcium transients
    tags = neuron.tag_neurons_parallel();

    if options.runRemoveFalsePositives==1
        neuron.remove_false_positives();
        subfxnDisplayState();
    end
    if options.runMerge==1
        neuron.merge_neurons_dist_corr(show_merge);
        neuron.merge_high_corr(show_merge, merge_thr_spatial);
        subfxnDisplayState();
    end

    if K~=size(neuron.A,2)
        neuron.update_spatial_parallel(use_parallel);
        neuron.update_temporal_parallel(use_parallel);
        if options.runRemoveFalsePositives==1
            neuron.remove_false_positives();
        end
        subfxnDisplayState();
    end

    %%
    % Get the folder path string
    [PATHSTR,NAME,EXT] = fileparts(inputFilename);
    [~,folderName,~] = fileparts(PATHSTR);

    results = neuron.obj2struct();
    cnmfeAnalysisOutput.success = 1;
    cnmfeAnalysisOutput.params = results.options;
    cnmfeAnalysisOutput.inputOptions = options;
    cnmfeAnalysisOutput.movieList = inputFilename;
    cnmfeAnalysisOutput.extractedImages = reshape(full(results.A),[neuron.options.d1 neuron.options.d2 size(results.C,1)]);
    % cnmfeAnalysisOutput.extractedImages = reshape(full(results.A),[size(results.P.sn) size(results.C,1)]);
    cnmfeAnalysisOutput.extractedSignals = results.C;
    cnmfeAnalysisOutput.extractedSignalsEst = results.C_raw;
    cnmfeAnalysisOutput.extractedPeaks = results.S;
    cnmfeAnalysisOutput.Cn = results.Cn;
    cnmfeAnalysisOutput.P = results.P;
    % save([PATHSTR filesep folderName '_cnmfeAnalysis.mat'],'cnmfeAnalysisOutput');

    % cnmfAnalysisOutput = cnmfeAnalysisOutput;
    % save([PATHSTR filesep folderName '_cnmfAnalysis.mat'],'cnmfAnalysisOutput');

    %% save the workspace for future analysis
    neuron.orderROIs('snr');
    try
        cnmfe_path = neuron.save_workspace();
    catch err
        display(repmat('@',1,7))
        disp(getReport(err,'extended','hyperlinks','on'));
        display(repmat('@',1,7))
    end
    %% show neuron contours
    % Coor = neuron.show_contours(0.6);

    %% create a video for displaying the
    % amp_ac = 140;
    % range_ac = 5+[0, amp_ac];
    % multi_factor = 10;
    % range_Y = 1300+[0, amp_ac*multi_factor];

    % avi_filename = neuron.show_demixed_video(save_demixed, kt, [], amp_ac, range_ac, range_Y, multi_factor);

    %% save neurons shapes
    % neuron.save_neurons();
    function subfxnDisplayState()
        disp(['A (space):' num2str(size(neuron.A)) ' | C (time): ' num2str(size(neuron.C)) ' | b (background): ' num2str(size(neuron.b))])
    end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function [ output_args ] = create_html_file_gm( filename, rx_lat, rx_long, hyperbola_lat, hyperbola_long, heatmap_cell, heatmap_threshold )
%CREATE_HTML_FILE   Generates the html code for google maps showing RX
%                   positions and the hyperbolas with a heatmap
%
%   rx_lat:         list of latitudes for RX positions (for displaying receiver positions)
%   rx_long:        list of longitudes for RX positions
%   hyperbola_lat:  cell array of latitudes for RX positions (one cell array elements correspond to one hyperbola)
%   hyperbola_lat:  cell array of longitudes for RX positions (one cell array elements correspond to one hyperbola)
%   heatmap_cell:   cell array, that contains {1} = long vector, {2} = lat vector and {3} = magnitude of heatmap points vector

    disp('writing html for google maps... ');

    % consistency checks
    if ~iscell(heatmap_cell)
        error('Parameter heatmap_cell needs to be a cell array with long, lat and magnitude');
    end

    if size(heatmap_cell) ~= 3
        error('Parameter heatmap_cell needs to be a cell array with long, lat and magnitude');
    end

    if ~iscell(hyperbola_lat) || ~iscell(hyperbola_long)
        error('Parameter hyperbola_lat, hyperbola_long need to be cell arrays');
    end

    if length(hyperbola_lat) ~= length(hyperbola_long)
        error('Length of hyperbola latitude and longitude values do not match');
    end

    for ii=1:length(hyperbola_lat)
        if length(cell2mat(hyperbola_lat(ii))) ~= length(cell2mat(hyperbola_long(ii)))
            error(['Length of hyperbola latitude and longitude in cell array doesn not for hyperbola ' num2str(ii)]);
        end
    end

    if size(rx_lat) ~= size(rx_long)
        error('Dimensions of rx latitude and longitude values do not match');
    end

    num_rx_positions = length(rx_lat);

    num_hyperbolas = length(hyperbola_lat);  % length of cell array

    num_hyperb_points = zeros(1, num_hyperbolas);
    for ii = 1:num_hyperbolas
        num_hyperb_points(ii) = length(hyperbola_lat{ii}); % length of each vector in the cell array
    end

    %% generate html code
    disp(['write html to file: ' filename]);

    fid = fopen(filename,'w');

    fprintf(fid, [...
        '<!DOCTYPE html>\n'...
        '\n'...
        '\n'...
        '\n'...
        '\n'...
        '

\n'...
        '\n'...
        '\n'...
        '\n'...
        '\n'...
        '\n'...
        '\n'...
        '\n'...
        '\n'...
        ]);

    fclose(fid);

    disp('writing html, done!');
end

repo: tdoa-evaluation-rtlsdr


text:
% Create folder if it doesn't exist
%
function mkfolder(folder)
    if exist(folder, 'dir') == 0
        mkdir(folder);
        if exist(folder, 'dir') == 0
            error('Failed to create output directory. Check the argument or create the directory manually');
        end
    end
end

author: kavli-ntnu


text:
function [ a, b, c, d, e, f ] = rs_to_xy_map ( t )

%*****************************************************************************80
%
%% RS_TO_XY_MAP returns the linear map from reference to physical triangle.
%
%  Location:
%
%    http://people.sc.fsu.edu/~jburkardt/m_src/triangle_poly_integral/rs_to_xy_map.m
%
%  Discussion:
%
%    This function returns the coefficients of the linear map that sends
%    the vertices of the reference triangle, (0,0), (1,0) and (0,1), to
%    the vertices of a physical triangle T, of the form:
%
%      X = A + B * R + C * S;
%      Y = D + E * R + F * S.
%
%  Reference Element:
%
%    |
%    1  3
%    |  |\
%    |  | \
%    S  |  \
%    |  |   \
%    |  |    \
%    0  1-----2
%    |
%    +--0--R--1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real T(2,3), the coordinates of the vertices.  The vertices are
%    assumed to be the images of (0,0), (1,0) and (0,1) respectively.
%
%    Output, real A, B, C, D, E, F, the mapping coefficients.
%
  a = t(1,1);
  b = t(1,2) - t(1,1);
  c = t(1,3) - t(1,1);

  d = t(2,1);
  e = t(2,2) - t(2,1);
  f = t(2,3) - t(2,1);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function res = vl_ffdnet_matlab(net, input)

%% If you did not install the matconvnet package, you can use this for testing.

global sigmas;
n = numel(net.layers);
res = struct('x', cell(1,n+1));
res(1).x = input;

for i = 1 : n
    l = net.layers{i};
    switch l.type

        case 'conv'
            disp(['Processing ... ',int2str(i),'/',int2str(n)]);
            for noutmaps = 1 : size(l.weights{1},4)
                z = zeros(size(res(i).x,1),size(res(i).x,2),'single');
                for ninmaps = 1 : size(res(i).x,3)
                    z = z + convn(res(i).x(:,:,ninmaps), rot90(l.weights{1}(:,:,ninmaps,noutmaps),2),'same'); % 180 degree rotation for kernel
                end
                res(i+1).x(:,:,noutmaps) = z + l.weights{2}(noutmaps);
            end

        case 'relu'
            res(i+1).x = max(res(i).x,0);

        case 'concat'
            if size(sigmas,1)~=size(res(i).x,1)
                sigmaMap   = bsxfun(@times,ones(size(res(i).x,1),size(res(i).x,2),1,size(res(i).x,4),'single'),permute(sigmas,[3 4 1 2]));
                res(i+1).x = cat(3,res(i).x,sigmaMap);
            else
                res(i+1).x = cat(3,res(i).x,sigmaMap);
            end

        case 'SubP'
            res(i+1).x = vl_nnSubP(res(i).x, [],'scale',l.scale);

    end
    res(i).x = [];
end

end

author: cszn

repo: DnCNN

save_path: github-repos/MATLAB/cszn-DnCNN


text:
function [exp_apx]=exp_approx3(dboun, exp_par,show_res)

%convert i
nind=find(exp_par(:,2)<=0);
exp_par(nind,2)=-1*exp_par(nind,2);

dind1=dboun(1,1):dboun(1,2);
nd1=dboun(1,2)-dboun(1,1)+1;
dind2=dboun(2,1):dboun(2,2);
nd2=dboun(2,2)-dboun(2,1)+1;
dind3=dboun(3,1):dboun(3,2);
nd3=dboun(3,2)-dboun(3,1)+1;
dind4=dboun(4,1):dboun(4,2);
nd4=dboun(4,2)-dboun(4,1)+1;
nd=nd1+nd2+nd3+nd4;

H=zeros(nd,5);
Y=zeros(nd,1);

H(1:nd1,1)=dind1;
H(1:nd1,3)=1;
H((nd1+1):(nd1+nd2),2)=dind2;
H((nd1+1):(nd1+nd2),4)=1;
H((nd1+nd2+1):(nd1+nd2+nd3),1)=dind3;
H((nd1+nd2+1):(nd1+nd2+nd3),5)=1;
H((nd1+nd2+nd3+1):nd,2)=dind4;
H((nd1+nd2+nd3+1):nd,5)=1;

Rinv1=1./(linspace(dboun(1,1),dboun(1,2),length(dind1)));
Rinv2=1./(linspace(dboun(2,1),dboun(2,2),length(dind2)));
Rinv3=1./(linspace(dboun(3,1),dboun(3,2)^2,length(dind3)));
Rinv4=1./(linspace(dboun(4,1),dboun(4,2)^2,length(dind4)));
HR(1:nd1,1)=dind1.*Rinv1;
HR(1:nd1,3)=Rinv1;
HR((nd1+1):(nd1+nd2),2)=dind2.*Rinv2;
HR((nd1+1):(nd1+nd2),4)=Rinv2;
HR((nd1+nd2+1):(nd1+nd2+nd3),1)=dind3.*Rinv3;
HR((nd1+nd2+1):(nd1+nd2+nd3),5)=Rinv3;
HR((nd1+nd2+nd3+1):nd,2)=dind4.*Rinv4;
HR((nd1+nd2+nd3+1):nd,5)=Rinv4;

Y(1:nd1)=log(exp_par(1,2)*exp_par(1,1).^dind1);
Y((nd1+1):(nd1+nd2))=log(exp_par(2,2)*exp_par(2,1).^dind2);
Y((nd1+nd2+1):(nd1+nd2+nd3))=log(exp_par(3,2)*exp_par(3,1).^dind3);
Y((nd1+nd2+nd3+1):nd)=log(exp_par(4,2)*exp_par(4,1).^dind4);

res=inv(HR'*H)*HR'*Y;
exp_apx=[exp(res(1)) exp(res(3));exp(res(2)) exp(res(4));exp(res(1)) exp(res(5));exp(res(2)) exp(res(5))];
exp_apx(nind,2)=-1*exp_apx(nind,2);

if (show_res==1)
    exp_par(nind,2)=-1*exp_par(nind,2);
    figure;
    for i=1:4
        ind=0:dboun(i,2);
        plot(ind,exp_par(i,1).^ind*exp_par(i,2))
        hold on
        plot(ind,exp_apx(i,1).^ind*exp_apx(i,2),'r')
    end
    grid;
end

repo: ekfmonoslam


text:
function A = vandermonde_interp_1d_matrix ( n, x )

%*****************************************************************************80
%
%% VANDERMONDE_INTERP_1D_MATRIX computes a Vandermonde 1D interpolation matrix.
%
%  Discussion:
%
%    We assume the interpolant has the form
%
%      p(x) = c1 + c2 * x + c3 * x^2 + ... + cn * x^(n-1).
%
%    We have n data values (x(i),y(i)) which must be interpolated:
%
%      p(x(i)) = c1 + c2 * x(i) + c3 * x(i)^2 + ... + cn * x(i)^(n-1) = y(i)
%
%    This can be cast as an NxN linear system for the polynomial
%    coefficients:
%
%      [ 1 x1 x1^2 ... x1^(n-1) ] [  c1 ] = [  y1 ]
%      [ 1 x2 x2^2 ... x2^(n-1) ] [  c2 ] = [  y2 ]
%      [ ...................... ] [ ... ] = [ ... ]
%      [ 1 xn xn^2 ... xn^(n-1) ] [  cn ] = [  yn ]
%
%    and if the x values are distinct, the matrix A is theoretically
%    invertible (though in fact, generally badly conditioned).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 July 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of data points.
%
%    Input, real X(N,1), the data values.
%
%    Output, real A(N,N), the Vandermonde matrix for X.
%
  A = ( x(1:n,1) * ones(1,n) ) .^ ( ones(n,1) * (0:n-1) );

  return
end


text:
function mvroi_mdsfig_plot2(CLUSTER,SPEC,sigmatavg,sigmatdif,titlestr)
% ::
%
%    mvroi_mdsfig_plot2(CLUSTER,SPEC,DATA.CORRELS.AVGSTATS.sigmat_uncorrected,DATA.CORRELS.DIFSTATS.sigmat_uncorrected,titlestr)
%
% :Example:
% ::
%
%    mvroi_mdsfig_plot2(DATA.CLUSTER,DATA.SPEC,DATA.CORRELS.AVGSTATS.sigmat_uncorrected,DATA.CORRELS.DIFSTATS.sigmat_uncorrected,'Uncorrected');
%
scnsize = get(0,'ScreenSize');
figure('position',[50 50 scnsize(3)-100 scnsize(4)/2],'color','white')

subplot(1,2,1);
set(gca,'FontSize',14)
mdsfig(CLUSTER.Gs(:,1:2),CLUSTER.names,CLUSTER.classes,sigmatavg);
title(['Average across conditions, ', titlestr])

subplot(1,2,2);
set(gca,'FontSize',14)
mdsfig(CLUSTER.Gs(:,1:2),CLUSTER.names,CLUSTER.classes,sigmatdif);
title(sprintf(['Contrast [' repmat('%3.0f ',1,size(SPEC.comps,2)) '] across conditions, %s'],SPEC.comps(1,:),SPEC.comptitle))

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function Bh = hmxMldivide(Mh,Bh)
%+========================================================================+
%|                                                                        |
%|         OPENHMX - LIBRARY FOR H-MATRIX COMPRESSION AND ALGEBRA         |
%|           openHmx is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : hmxMldivide.m                                 |
%|    #    |   VERSION    : 0.51                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.01.2019                                    |
%| ( === ) |   SYNOPSIS   : Solve Mh x = Bh with the rule                 |
%|  `---'  |                Compr > Full > H-Matrix                       |
%+========================================================================+

% Check dimensions
if (size(Mh,1) ~= size(Bh,1))
    error('hmxMldivide.m : matrix dimensions must agree.')
end

% Check input data
if ~isa(Mh,'hmx')
    error('hmxMldivide.m : unavailable case.')
end

%%% Lower H-Matrix \ H-Matrix -> H-Matrix
if islower(Mh) && isa(Bh,'hmx')
    % H-Matrix \ H-Matrix -> H-Matrix (recursion)
    if (Mh.typ == 0) && (Bh.typ == 0)
        % X11 -> L11 \ B11
        Bh.chd{1} = hmxMldivide(Mh.chd{1},Bh.chd{1});

        % X12 -> L11 \ B12
        Bh.chd{2} = hmxMldivide(Mh.chd{1},Bh.chd{2});

        % X21 -> L22 \ (B21 - L21*X11)
        Bh.chd{3} = Bh.chd{3} - Mh.chd{3} * Bh.chd{1};
        Bh.chd{3} = hmxMldivide(Mh.chd{4},Bh.chd{3});

        % X22 -> L22 \ (B22 - L21 * X12)
        Bh.chd{4} = Bh.chd{4} - Mh.chd{3} * Bh.chd{2};
        Bh.chd{4} = hmxMldivide(Mh.chd{4},Bh.chd{4});

        % Fusion
        Bh = hmxFusion(Bh);

    % H-Matrix \ Compr -> Compr
    elseif (Mh.typ == 0) && (Bh.typ == 1)
        Bh.dat = { hmxMldivide(Mh,Bh.dat{1}) , Bh.dat{2} };

    % H-Matrix \ Full -> Full
    elseif (Mh.typ == 0) && (Bh.typ == 2)
        Bh.typ = 2;
        Bh.dat = hmxMldivide(Mh,Bh.dat);

    % Compr \ --- -> ---
    elseif (Mh.typ == 1)
        error('hmxMldivide : unavailable case')

    % Full \ H-Matrix -> Unavailable
    elseif (Mh.typ == 2) && (Bh.typ == 0)
        error('hmxMldivide : unavailable case')

    % Full \ Compr -> Compr
    elseif (Mh.typ == 2) && (Bh.typ == 1)
        Bh.dat = { Mh.dat\Bh.dat{1} , Bh.dat{2} };

    % Full \ Full -> Full
    elseif (Bh.typ == 2)
        Bh.dat = Mh.dat \ Bh.dat;

    else
        error('hmxMldivide : unavailable case')
    end

%%% Upper H-Matrix \ H-Matrix -> H-Matrix
elseif isupper(Mh) && isa(Bh,'hmx')
    % H-Matrix \ H-Matrix -> H-Matrix (recursion)
    if (Mh.typ == 0) && (Bh.typ == 0)
        % X22 -> U22 \ B22
        Bh.chd{4} = hmxMldivide(Mh.chd{4},Bh.chd{4});

        % X21 -> U22 \ B21
        Bh.chd{3} = hmxMldivide(Mh.chd{4},Bh.chd{3});

        % X12 -> U11 \ (B12 - U12*X22)
        Bh.chd{2} = Bh.chd{2} - Mh.chd{2} * Bh.chd{4};
        Bh.chd{2} = hmxMldivide(Mh.chd{1},Bh.chd{2});

        % X11 -> U11 \ (B11 - U12 * X21)
        Bh.chd{1} = Bh.chd{1} - Mh.chd{2} * Bh.chd{3};
        Bh.chd{1} = hmxMldivide(Mh.chd{1},Bh.chd{1});

        % Fusion
        Bh = hmxFusion(Bh);

    % H-Matrix \ Compr -> Compr
    elseif (Mh.typ == 0) && (Bh.typ == 1)
        Bh.dat = { hmxMldivide(Mh,Bh.dat{1}) , Bh.dat{2} };

    % H-Matrix \ Full -> Full
    elseif (Mh.typ == 0) && (Bh.typ == 2)
        Bh.typ = 2;
        Bh.dat = hmxMldivide(Mh,Bh.dat);

    % Compr \ --- -> ---
    elseif (Mh.typ == 1)
        error('hmxMldivide : unavailable case')

    % Full \ H-Matrix -> Unavailable
    elseif (Mh.typ == 2) && (Bh.typ == 0)
        error('hmxMldivide : unavailable case')

    % Full \ Compr -> Compr
    elseif (Mh.typ == 2) && (Bh.typ == 1)
        Bh.dat = { Mh.dat\Bh.dat{1} , Bh.dat{2} };

    % Full \ Full -> Full
    elseif (Bh.typ == 2)
        Bh.dat = Mh.dat \ Bh.dat;

    else
        error('hmxMldivide : unavailable case')
    end

%%% Lower H-Matrix \ Matrix -> Matrix
elseif islower(Mh)
    if (size(Bh,2) > 0)
        % H-Matrix (recursion)
        if (Mh.typ == 0)
            % X1 -> L11 \ B1
            X1 = hmxMldivide(Mh.chd{1},Bh(Mh.row{1},:));

            % X2 -> L22 \ (B2 - L21*X1)
            X2 = Bh(Mh.row{3},:) - Mh.chd{3}*X1;
            X2 = hmxMldivide(Mh.chd{4},X2);

            % Bh = [X1 X2]
            if issparse(X1)
                Bh = sparse(size(Mh,1),size(Bh,2));
            else
                Bh = zeros(size(Mh,1),size(Bh,2),class(X1));
            end
            Bh(Mh.col{1},:) = X1;
            Bh(Mh.col{2},:) = X2;

        % Compressed leaf
        elseif (Mh.typ == 1)
            error('hmxMldivide : unavailable case')

        % Full leaf
        elseif (Mh.typ == 2)
            Bh = Mh.dat \ Bh;

        % Unknown type
        else
            error('hmxMldivide : unavailable case')
        end
    else
        Bh = zeros(size(Mh,1),0);
    end

%%% Upper H-Matrix \ Matrix -> Matrix
elseif isupper(Mh)
    if (size(Bh,2) > 0)
        % H-Matrix (recursion)
        if (Mh.typ == 0)
            % X2 -> U22 \ B2
            X2 = hmxMldivide(Mh.chd{4},Bh(Mh.row{4},:));

            % X1 -> U11 \ (B1 - U12*X2)
            X1 = Bh(Mh.row{1},:) - Mh.chd{2}*X2;
            X1 = hmxMldivide(Mh.chd{1},X1);

            % Bh = [X1 ; X2]
            if issparse(X1)
                Bh = sparse(size(Mh,1),size(Bh,2));
            else
                Bh = zeros(size(Mh,1),size(Bh,2),class(X1));
            end
            Bh(Mh.col{1},:) = X1;
            Bh(Mh.col{2},:) = X2;

        % Compressed leaf
        elseif (Mh.typ == 1)
            error('hmxMldivide : unavailable case')

        % Full leaf
        elseif (Mh.typ == 2)
            Bh = Mh.dat \ Bh;

        % Unknown type
        else
            error('hmxMldivide : unavailable case')
        end
    else
        Bh = zeros(size(Mh,1),0);
    end

%%% H-Matrix \ (H-)Matrix -> (H-)Matrix
else
    [Lh,Uh] = lu(Mh);
    Bh      = Uh\(Lh\Bh);
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function timing_matlab_commands
% TIMING_MATLAB_COMMANDS  Testing for speed different MATLAB commands.
%
% Main conclusion: RESHAPE and * (i.e. MTIMES) are very quick!

% Paolo de Leva
% University of Rome, Foro Italico, Rome, Italy
% 2008 Dec 24

clear all

% Checking whether needed software exists
if ~exist('bsxfun', 'builtin')
    message = sysrequirements_for_testing('bsxmex', 'timeit');
else
    message = sysrequirements_for_testing('timeit');
end
if message
    disp ' ', error('timing_matlab_commands:Missing_subfuncs', message)
end

disp ' '
disp '---------------------------------- Experiment 1 ----------------------------------'
N = 10000; P = 3; Q = 3; R = 1;
timing(N,P,Q,R);

disp '---------------------------------- Experiment 2 ----------------------------------'
N = 1000; P = 3;  Q = 30; R = 1;
timing(N,P,Q,R);

disp '---------------------------------- Experiment 3 ----------------------------------'
N = 1000; P = 9; Q = 10;  R = 3;
timing(N,P,Q,R);

disp '---------------------------------- Experiment 4 ----------------------------------'
N = 100; P = 9; Q = 100;  R = 3;
timing(N,P,Q,R);

disp '---------------------------------- Experiment 5 ----------------------------------'
disp ' '
timing2(4, 10000);
timing2(200, 200);
timing2(10000, 4);

disp '---------------------------- Experiment 6 ----------------------------'
disp ' '
a = rand(4096, 4096);
fprintf ('Size of A:  %0.0f x %0.0f\n', size(a))
disp ' '
disp '   SUM(A,1)  SUM(A,2)'
f1 = @() sum(a, 1);
f2 = @() sum(a, 2);
disp ([timeit(f1), timeit(f2)])

clear all
b = rand(256, 256, 256);
fprintf ('Size of B:  %0.0f x %0.0f x %0.0f\n', size(b))
disp ' '
disp '   SUM(B,1)  SUM(B,2)  SUM(B,3)'
f1 = @() sum(b, 1);
f2 = @() sum(b, 2);
f3 = @() sum(b, 3);
disp ([timeit(f1), timeit(f2), timeit(f3)])

disp '---------------------------- Experiment 7 ----------------------------'
disp ' '
a = rand(101,102,103);
fprintf ('Size of A:  %0.0f x %0.0f x %0.0f\n', size(a))
disp ' '
disp 'Moving last dimension to first dimension:'
disp 'PERMUTE(A,[3 2 1])  PERMUTE(A,[3 1 2])  SHIFTDIM(A,2)'
disp '(SWAPPING)          (SHIFTING)          (SHIFTING)'
f1 = @() permute(a, [3 2 1]);
f2 = @() permute(a, [3 1 2]);
f3 = @() shiftdim(a, 2);
fprintf(1, '%8.2g            ', [timeit(f1), timeit(f2), timeit(f3)])
disp ' ', disp ' '
a2 = f1(); s = size(a2);
a2 = f2(); s(2,:) = size(a2);
a2 = f3(); s(3,:) = size(a2);
disp (s)

disp 'Moving first dimension to last dimension:'
disp 'PERMUTE(A,[3 2 1])  PERMUTE(A,[2 3 1])  SHIFTDIM(A,1)'
disp '(SWAPPING)          (SHIFTING)          (SHIFTING)'
f1 = @() permute(a, [3 2 1]);
f2 = @() permute(a, [2 3 1]);
f3 = @() shiftdim(a, 1);
fprintf(1, '%8.2g            ', [timeit(f1), timeit(f2), timeit(f3)])
disp ' ', disp ' '
a2 = f1(); s = size(a2);
a2 = f2(); s(2,:) = size(a2);
a2 = f3(); s(3,:) = size(a2);
disp (s)

disp ' '
a = rand(21,22,23,24,25);
fprintf ('Size of A:  %0.0f x %0.0f x %0.0f x %0.0f x %0.0f\n', size(a))
disp ' '
disp 'Moving 4th dimension to 1st dimension:'
disp 'PERMUTE(A,[4 2 3 1 5])  PERMUTE(A,[4 1 2 3 5])  PERMUTE(A,[4 5 1 2 3])'
disp '(SWAPPING)              (PARTIAL SHIFTING)      (SHIFTING)'
f1 = @() permute(a, [4 2 3 1 5]);
f2 = @() permute(a, [4 1 2 3 5]);
f3 = @() permute(a, [4 5 1 2 3]);
fprintf(1, '%8.2g                ', [timeit(f1), timeit(f2), timeit(f3)])
disp ' ', disp ' '
a2 = f1(); s = size(a2);
a2 = f2(); s(2,:) = size(a2);
a2 = f3(); s(3,:) = size(a2);
disp (s)

disp 'Moving 2nd dimension to 5th dimension:'
disp 'PERMUTE(A,[1 5 3 4 2])  PERMUTE(A,[1 3 4 5 2])  PERMUTE(A,[3 4 5 1 2])'
disp '(SWAPPING)              (PARTIAL SHIFTING)      (SHIFTING)'
f1 = @() permute(a, [1 5 3 4 2]);
f2 = @() permute(a, [1 3 4 5 2]);
f3 = @() permute(a, [3 4 5 1 2]);
fprintf(1, '%8.2g                ', [timeit(f1), timeit(f2), timeit(f3)])
disp ' ', disp ' '
a2 = f1(); s = size(a2);
a2 = f2(); s(2,:) = size(a2);
a2 = f3(); s(3,:) = size(a2);
disp (s)

disp '---------------------------- Experiment 8 ----------------------------'
disp ' '
a =rand(101,102,103);
order = [1 2 3];
shape = [101,102,103];
f1 = @() perm(a,order);
f2 = @() ifpermute(a,order);
f3 = @() ifpermute2(a,order);
f4 = @() resh(a,shape);
f5 = @() ifreshape(a,shape);
f6 = @() ifreshape2(a,shape);
disp 'COMPARING STATEMENTS THAT DO NOTHING!'
disp ' '
fprintf ('Size of A:  %0.0f x %0.0f x %0.0f\n', size(a))
disp ' '
disp 'ORDER = [1 2 3]       % (keeping same order)'
disp 'SHAPE = [101,102,103] % (keeping same shape)'
disp ' '
fprintf (1,'PERMUTE(A,ORDER) ..........................................  %0.4g\n', timeit(f1))
fprintf (1,'IF ~ISEQUAL(ORDER,1:LENGTH(ORDER)), A=PERMUTE(A,ORDER); END  %0.4g\n', timeit(f2))
fprintf (1,'IF ~ISEQUAL(ORDER,1:3),             A=PERMUTE(A,ORDER); END  %0.4g\n', timeit(f3))
disp ' '
fprintf (1,'RESHAPE(A,SHAPE) ..........................................  %0.4g\n', timeit(f4))
fprintf (1,'IF ~ISEQUAL(SHAPE,SIZE(A)), A=RESHAPE(A,SHAPE); END .......  %0.4g\n', timeit(f5))
fprintf (1,'IF ~ISEQUAL(SHAPE,SHAPE),   A=RESHAPE(A,SHAPE); END .......  %0.4g\n', timeit(f5))
disp ' '

function a=perm(a, order)
a=permute(a, order);
function a=resh(a,shape)
a=reshape(a,shape);
function a=ifpermute(a, order)
if ~isequal(order, 1:length(order)), a=permute(a,order); end
function a=ifreshape(a, shape)
if ~isequal(shape, size(a)), a=reshape(a,shape); end
function a=ifpermute2(a, order)
if ~isequal(order, 1:3), a=permute(a,order); end
function a=ifreshape2(a, shape)
if ~isequal(shape, shape), a=reshape(a,shape); end

function timing(N,P,Q,R)

a0 = rand(1, P, Q);
b0 = rand(1, Q, R);
a = a0(ones(1,N),:,:); % Cloning along first dimension
b = b0(ones(1,N),:,:); % Cloning along first dimension
[n1 p q1] = size(a); % reads third dim even if it is 1.
[n2 q2 r] = size(b); % reads third dim even if it is 1.
disp ' '
disp        'Array  Size      Size               Number of elements'
fprintf (1, 'A      Nx(PxQ)   %0.0f x (%0.0f x %0.0f) %8.0f\n', [n1 p q1  numel(a)])
fprintf (1, 'B      Nx(QxR)   %0.0f x (%0.0f x %0.0f) %8.0f\n', [n2 q2 r  numel(b)])
f1 = @() permute(a, [2 3 1]);
f2 = @() permute(a, [1 3 2]);
f3 = @() permute(a, [2 1 3]);
f4 = @() permute(a, [1 2 3]);
f5 = @() permute(b, [2 3 1]);
f6 = @() permute(b, [1 3 2]);
f7 = @() permute(b, [2 1 3]);
f8 = @() permute(b, [1 2 3]);
disp ' '
disp '   PERMUTE(A,[2 3 1])  PERMUTE(A,[1 3 2])  PERMUTE(A,[2 1 3])  PERMUTE(A,[1 2 3])'
fprintf(1, '%20.5f', [timeit(f1), timeit(f2), timeit(f3), timeit(f4)])
disp ' '
disp '   PERMUTE(B,[2 3 1])  PERMUTE(B,[1 3 2])  PERMUTE(B,[2 1 3])  PERMUTE(B,[1 2 3])'
fprintf(1, '%20.5f', [timeit(f5), timeit(f6), timeit(f7), timeit(f8)])
disp ' '
disp ' '
disp '   RESHAPE(A,[N*P Q])  RESHAPE(B,[N R Q])  RESHAPE(B,[N 1 R Q])'
f1 = @() reshape(a, [N*P Q]);
f2 = @() reshape(b, [N R Q]);
f3 = @() reshape(b, [N 1 R Q]);
fprintf(1, '%20.5f', [timeit(f1), timeit(f2), timeit(f3)])
disp ' '
f1 = @() a .* a;
f2 = @() bsxfun(@times, a, a);
f3 = @() b .* b;
f4 = @() bsxfun(@times, b, b);
disp ' '
disp '              A .* A   BSXFUN(@TIMES,A,A)'
fprintf(1, '%20.5f%20.5f\n', [timeit(f1), timeit(f2)])
disp '              B .* B   BSXFUN(@TIMES,B,B)'
fprintf(1, '%20.5f%20.5f\n', [timeit(f3), timeit(f4)])
if R==1
    disp ' '
    disp '   NOTE: If R=1 then RESHAPE(B,[N R Q]) is equivalent to'
    disp '                     PERMUTE(B,[1 3 2]) but much faster!'
    disp '                     (at least on my system)'
end
disp ' '

function timing2(P,Q)

a = rand(P, Q);
b = rand(Q, 1);
fprintf ('Size of A:  %0.0f x %0.0f\n', size(a))
fprintf ('Size of B:  %0.0f x %0.0f\n', size(b))
disp ' '
disp '        A * B   TONY''S TRICK     BSXFUN'
f1 = @() a * b;
f2 = @() clone_multiply_sum(a, b', P);
f3 = @() sum(bsxfun(@times, a, b'), 2);
fprintf(1, '%13.5f', [timeit(f1), timeit(f2), timeit(f3)])
disp ' '
disp ' '
c = f1() - f2();
d = max(c(:));
if  d > eps*20
    disp 'There is an unexpected output difference:';
    disp (d);
end

function c = clone_multiply_sum(a,b,P)
c = sum(a .* b(ones(1,P),:), 2);

repo: Aero-matlab


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

%%
clear; close all; clc;

%% Syntax
% |[varargout]=importFaceFile_tetGen(fileName);|

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

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function DataSet = prtDataGenProstate(prostateFile)
%   prtDataGenIris Read in data from the UCI PROSTATE data file
%
%   DataSet = prtDataGenIris reads in the data from the classic Prostate.UCI
%   file used in many pattern recognition texts and papers. For more
%   information about this data, see the following URL:
%
%   http://lib.stat.cmu.edu/S/Harrell/data/descriptions/prostate.html
%
%   Example:
%
%   ds = prtDataGenProstate;
%   ds.plot;
%
%   See also: prtDataSetClass, prtDataGenBiModal, prtDataGenIris,
%   prtDataGenMary, prtDataGenNoisySinc, prtDataGenOldFaithful,
%   prtDataGenSpiral, prtDataGenUnimodal, prtDataGenUnimodal, prtDataGenXor

%   prtDataGenProstate Read in data from a CSV formatted PROSTATE-like data file
%
% DataSet = prtDataGenProstate; Read in the data from the Prostate.UCI file
% available with the prt and on the web.
%
% DataSet = prtDataGenProstate(ProstateFile);  Read in the CSV formatted,
% file matching the specifications of ProstateFile
% (see: ProstateFile in <<prtRoot>>\data\)
%
% http://lib.stat.cmu.edu/S/Harrell/data/descriptions/prostate.html

if nargin == 0
    prostateFile = fullfile(prtRoot,'dataGen','dataStorage','prostate','prostate.csv');
end

fid = fopen(prostateFile,'r');
%patno,stage,rx,dtime,status,age,wt,pf,hx,sbp,dbp,ekg,hg,sz,sg,ap,bm,sdate
topLineSpec = repmat('%s ',1,18);
%1,3,0.2 mg estrogen,72,alive,75,76,normal activity,0,15,9,heart strain,13.79882813,2,8,0.299987793,0,2778
dataSpec = '%d %d %s %d %s %d %d %s %d %d %d %s %f %d %d %f %d %f';

varNames = textscan(fid,topLineSpec,1,'whitespace',' \b\t""','Delimiter',',');
cellData = textscan(fid,dataSpec,'whitespace',' \b\t"','Delimiter',',');
fclose(fid);

%handle the data
cellData{3} = prostateFixEstrogen(cellData{3});
cellData{3} = cat(1,cellData{3}{:});
[cellData{5},col5Labels] = prtUtilStringsToClassNumbers(cellData{5});
[cellData{8},col8Labels] = prtUtilStringsToClassNumbers(cellData{8});
[cellData{12},col12Labels] = prtUtilStringsToClassNumbers(cellData{12});

targetInd = 5;
featureInd = setdiff(2:length(cellData),targetInd);
Y = cellData{targetInd};
X = double(cat(2,cellData{featureInd}));  %drop the counting var, and the Y
colNames = varNames(featureInd);

varNames = cellfun(@(c)c{1},varNames,'uniformoutput',false);
colNames = varNames(featureInd);

userData = struct('column6IntegerExplanations',{col8Labels},...
    'column10IntegerExplanations',{col12Labels});

DataSet = prtDataSetClass(X,Y,'name','UCI Prostate Data');
DataSet = DataSet.setFeatureNames(colNames);
DataSet.userData = userData;
DataSet = DataSet.setClassNames(col5Labels);

function cellData = prostateFixEstrogen(cellData)
for i = 1:length(cellData)
    switch lower(cellData{i})
        case 'placebo'
            cellData{i} = 0;
        otherwise
            cellData{i} = str2double(cellData{i}(1:3));
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [ know, x ] = p07_sol ( )

%*****************************************************************************80
%
%% P07_SOL returns the solution for problem 7.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 February 2002
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer KNOW.
%    If KNOW is 0, then the solution is not known.
%    If KNOW is positive, then the solution is known, and is returned in X.
%
%    Output, real X, the solution, if known.
%
  know = 1;

  x = -0.6795786599525;

  return
end


text:
function P=compute_constraint_distance_3param_6eq_6unk(m1,m2,m3)

% Copyright (C) <2007>  <Francesc Moreno-Noguer, Vincent Lepetit, Pascal Fua>
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the version 3 of the GNU General Public License
% as published by the Free Software Foundation.
%
% 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/>.
%
% Francesc Moreno-Noguer, CVLab-EPFL, September 2007.
% fmorenoguer@gmail.com, http://cvlab.epfl.ch/~fmoreno/

%redefine variables name, for compatibility with maple
m1_1=m1(1);
m1_2=m1(2);
m1_3=m1(3);
m1_4=m1(4);
m1_5=m1(5);
m1_6=m1(6);
m1_7=m1(7);
m1_8=m1(8);
m1_9=m1(9);
m1_10=m1(10);
m1_11=m1(11);
m1_12=m1(12);

m2_1=m2(1);
m2_2=m2(2);
m2_3=m2(3);
m2_4=m2(4);
m2_5=m2(5);
m2_6=m2(6);
m2_7=m2(7);
m2_8=m2(8);
m2_9=m2(9);
m2_10=m2(10);
m2_11=m2(11);
m2_12=m2(12);

m3_1=m3(1);
m3_2=m3(2);
m3_3=m3(3);
m3_4=m3(4);
m3_5=m3(5);
m3_6=m3(6);
m3_7=m3(7);
m3_8=m3(8);
m3_9=m3(9);
m3_10=m3(10);
m3_11=m3(11);
m3_12=m3(12);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t1 = (m1_2 ^ 2);
t4 = (m1_6 ^ 2);
t5 = (m1_3 ^ 2);
t6 = (m1_5 ^ 2);
t11 = (m1_4 ^ 2);
t12 = (m1_1 ^ 2);
t20 = m1_4 * m2_4;
t21 = m1_3 * m2_3;
t22 = m1_5 * m2_5;
t23 = m1_2 * m2_2;
t24 = m1_6 * m2_6;
t25 = m1_1 * m2_1;
t26 = (-m2_4 * m1_1 - m2_5 * m1_2 - m2_6 * m1_3 - m1_6 * m2_3 - m1_4 * m2_1 - m1_5 * m2_2 + t20 + t21 + t22 + t23 + t24 + t25);
t27 = m1_6 * m3_6;
t29 = m1_5 * m3_5;
t30 = m1_4 * m3_4;
t33 = m1_3 * m3_3;
t35 = m1_1 * m3_1;
t38 = m1_2 * m3_2;
t39 = (t27 - m1_6 * m3_3 + t29 + t30 - m1_4 * m3_1 - m3_6 * m1_3 + t33 - m3_5 * m1_2 + t35 - m3_4 * m1_1 - m1_5 * m3_2 + t38);
t40 = (m2_4 ^ 2);
t41 = (m2_2 ^ 2);
t42 = (m2_5 ^ 2);
t43 = (m2_1 ^ 2);
t44 = (m2_6 ^ 2);
t45 = (m2_3 ^ 2);
t53 = m2_4 * m3_4;
t56 = m2_5 * m3_5;
t57 = m2_2 * m3_2;
t60 = m2_1 * m3_1;
t62 = m2_6 * m3_6;
t63 = m2_3 * m3_3;
t65 = (t53 - m2_4 * m3_1 - m3_4 * m2_1 + t56 + t57 - m2_5 * m3_2 - m2_6 * m3_3 + t60 - m3_6 * m2_3 + t62 + t63 - m3_5 * m2_2);
t66 = (m3_5 ^ 2);
t69 = (m3_4 ^ 2);
t70 = (m3_3 ^ 2);
t71 = (m3_2 ^ 2);
t72 = (m3_6 ^ 2);
t75 = (m3_1 ^ 2);
t81 = (m1_9 ^ 2);
t82 = (m1_8 ^ 2);
t83 = (m1_7 ^ 2);
t90 = m1_8 * m2_8;
t91 = m1_7 * m2_7;
t93 = m1_9 * m2_9;
t98 = (-m1_9 * m2_3 + t21 + t90 + t23 + t91 + t25 - m2_7 * m1_1 + t93 - m2_9 * m1_3 - m1_8 * m2_2 - m1_7 * m2_1 - m2_8 * m1_2);
t101 = m1_9 * m3_9;
t102 = m1_7 * m3_7;
t106 = m1_8 * m3_8;
t108 = (-m3_8 * m1_2 + t38 + t33 - m3_9 * m1_3 + t101 + t102 - m1_7 * m3_1 - m1_9 * m3_3 + t35 - m3_7 * m1_1 + t106 - m1_8 * m3_2);
t113 = (m2_8 ^ 2);
t116 = (m2_9 ^ 2);
t117 = (m2_7 ^ 2);
t119 = m2_8 * m3_8;
t122 = m2_9 * m3_9;
t125 = m2_7 * m3_7;
t128 = (t119 - m2_7 * m3_1 - m3_8 * m2_2 + t57 + t63 + t60 + t122 - m2_8 * m3_2 - m3_7 * m2_1 + t125 - m3_9 * m2_3 - m2_9 * m3_3);
t129 = (m3_7 ^ 2);
t130 = (m3_8 ^ 2);
t133 = (m3_9 ^ 2);
t141 = (m1_10 ^ 2);
t142 = (m1_11 ^ 2);
t143 = (m1_12 ^ 2);
t151 = m1_10 * m2_10;
t152 = m1_12 * m2_12;
t154 = m1_11 * m2_11;
t158 = (-m2_10 * m1_1 - m2_12 * m1_3 + t151 + t23 + t25 + t152 + t21 - m1_10 * m2_1 + t154 - m1_11 * m2_2 - m2_11 * m1_2 - m1_12 * m2_3);
t160 = m1_12 * m3_12;
t164 = m1_10 * m3_10;
t165 = m1_11 * m3_11;
t168 = (-m3_10 * m1_1 + t160 - m3_11 * m1_2 + t38 + t33 - m1_12 * m3_3 - m3_12 * m1_3 + t164 + t35 + t165 - m1_11 * m3_2 - m1_10 * m3_1);
t169 = (m2_12 ^ 2);
t170 = (m2_10 ^ 2);
t171 = (m2_11 ^ 2);
t179 = m2_10 * m3_10;
t181 = m2_12 * m3_12;
t185 = m2_11 * m3_11;
t188 = (t57 + t60 + t179 - m2_10 * m3_1 + t181 - m2_12 * m3_3 - m3_12 * m2_3 - m3_10 * m2_1 + t63 + t185 - m2_11 * m3_2 - m3_11 * m2_2);
t191 = (m3_12 ^ 2);
t192 = (m3_10 ^ 2);
t193 = (m3_11 ^ 2);
t212 = (t22 + t91 + t93 - m1_7 * m2_4 + t20 - m2_7 * m1_4 + t90 - m2_8 * m1_5 - m1_8 * m2_5 - m1_9 * m2_6 + t24 - m2_9 * m1_6);
t219 = (t101 + t30 - m3_9 * m1_6 - m3_8 * m1_5 - m1_8 * m3_5 - m3_7 * m1_4 + t27 - m1_9 * m3_6 - m1_7 * m3_4 + t102 + t106 + t29);
t233 = (t53 - m2_8 * m3_5 - m2_9 * m3_6 - m2_7 * m3_4 + t125 + t122 - m3_8 * m2_5 + t62 + t119 + t56 - m3_7 * m2_4 - m3_9 * m2_6);
t254 = (-m2_11 * m1_5 + t154 + t151 - m1_12 * m2_6 + t20 - m1_10 * m2_4 - m2_12 * m1_6 - m2_10 * m1_4 - m1_11 * m2_5 + t152 + t24 + t22);
t261 = (t30 - m3_12 * m1_6 - m1_10 * m3_4 - m1_11 * m3_5 + t160 + t27 + t164 + t165 - m3_11 * m1_5 - m3_10 * m1_4 - m1_12 * m3_6 + t29);
t275 = (-m3_10 * m2_4 + t56 - m2_10 * m3_4 + t62 - m3_12 * m2_6 - m2_11 * m3_5 + t53 - m3_11 * m2_5 - m2_12 * m3_6 + t179 + t181 + t185);
t296 = (-m2_12 * m1_9 + t152 + t91 + t90 - m1_11 * m2_8 + t151 - m2_11 * m1_8 + t93 - m1_10 * m2_7 - m1_12 * m2_9 - m2_10 * m1_7 + t154);
t303 = (t164 + t102 - m3_10 * m1_7 + t160 - m3_12 * m1_9 - m1_10 * m3_7 - m3_11 * m1_8 + t101 - m1_12 * m3_9 + t106 + t165 - m1_11 * m3_8);
t317 = (t125 + t119 - m3_12 * m2_9 - m3_10 * m2_7 - m2_11 * m3_8 - m3_11 * m2_8 + t122 + t179 - m2_12 * m3_9 + t181 - m2_10 * m3_7 + t185);
P(1,1) = t1 - 2 * m1_4 * m1_1 + t4 + t5 + t6 - 2 * m1_5 * m1_2 - 2 * m1_6 * m1_3 + t11 + t12;
P(1,2) = 2 * t26;
P(1,3) = 2 * t39;
P(1,4) = t40 + t41 + t42 + t43 + t44 + t45 - 2 * m2_5 * m2_2 - 2 * m2_6 * m2_3 - 2 * m2_4 * m2_1;
P(1,5) = 2 * t65;
P(1,6) = t66 - 2 * m3_4 * m3_1 + t69 + t70 + t71 + t72 - 2 * m3_6 * m3_3 + t75 - 2 * m3_5 * m3_2;
P(2,1) = -2 * m1_8 * m1_2 + t12 + t81 + t5 + t82 + t83 + t1 - 2 * m1_9 * m1_3 - 2 * m1_7 * m1_1;
P(2,2) = 2 * t98;
P(2,3) = 2 * t108;
P(2,4) = -2 * m2_8 * m2_2 - 2 * m2_9 * m2_3 + t113 - 2 * m2_7 * m2_1 + t116 + t117 + t41 + t43 + t45;
P(2,5) = 2 * t128;
P(2,6) = t75 + t70 + t129 + t71 + t130 - 2 * m3_9 * m3_3 + t133 - 2 * m3_8 * m3_2 - 2 * m3_7 * m3_1;
P(3,1) = -2 * m1_11 * m1_2 + t141 + t142 + t12 + t1 + t5 + t143 - 2 * m1_10 * m1_1 - 2 * m1_12 * m1_3;
P(3,2) = 2 * t158;
P(3,3) = 2 * t168;
P(3,4) = t169 + t41 + t43 + t45 + t170 + t171 - 2 * m2_12 * m2_3 - 2 * m2_10 * m2_1 - 2 * m2_11 * m2_2;
P(3,5) = 2 * t188;
P(3,6) = t71 - 2 * m3_12 * m3_3 + t75 + t191 + t70 + t192 + t193 - 2 * m3_10 * m3_1 - 2 * m3_11 * m3_2;
P(4,1) = -2 * m1_9 * m1_6 + t11 + t4 + t81 + t82 + t6 + t83 - 2 * m1_7 * m1_4 - 2 * m1_8 * m1_5;
P(4,2) = 2 * t212;
P(4,3) = 2 * t219;
P(4,4) = t117 + t113 - 2 * m2_8 * m2_5 - 2 * m2_7 * m2_4 - 2 * m2_9 * m2_6 + t116 + t40 + t42 + t44;
P(4,5) = 2 * t233;
P(4,6) = t129 - 2 * m3_9 * m3_6 + t69 + t66 + t72 + t133 - 2 * m3_8 * m3_5 - 2 * m3_7 * m3_4 + t130;
P(5,1) = t4 + t143 + t11 - 2 * m1_12 * m1_6 - 2 * m1_11 * m1_5 - 2 * m1_10 * m1_4 + t6 + t141 + t142;
P(5,2) = 2 * t254;
P(5,3) = 2 * t261;
P(5,4) = t170 + t171 + t169 - 2 * m2_10 * m2_4 - 2 * m2_11 * m2_5 - 2 * m2_12 * m2_6 + t40 + t42 + t44;
P(5,5) = 2 * t275;
P(5,6) = t69 + t66 + t72 - 2 * m3_12 * m3_6 - 2 * m3_10 * m3_4 + t193 - 2 * m3_11 * m3_5 + t192 + t191;
P(6,1) = t142 - 2 * m1_10 * m1_7 + t141 + t83 + t143 + t82 - 2 * m1_12 * m1_9 + t81 - 2 * m1_11 * m1_8;
P(6,2) = 2 * t296;
P(6,3) = 2 * t303;
P(6,4) = t171 + t170 - 2 * m2_12 * m2_9 - 2 * m2_10 * m2_7 - 2 * m2_11 * m2_8 + t113 + t169 + t116 + t117;
P(6,5) = 2 * t317;
P(6,6) = -2 * m3_11 * m3_8 + t193 + t133 + t129 + t192 + t130 - 2 * m3_12 * m3_9 + t191 - 2 * m3_10 * m3_7;

author: cvlab-epfl

repo: EPnP

save_path: github-repos/MATLAB/cvlab-epfl-EPnP


text:
function [selectionImage, selectionIndexArray, unusedVarargin] = select(this, varargin)
% Prototype for select-method of n-dimensional image data using dimInfo
%
%   Y = MrImage()
%   selectionImage = Y.select('type', 'index'/'label', 'invert', true/false, ...
%              'removeDims', true/false
%               'dimLabel1', arrayIndicesDim1/samplingPointsDim1, ...,
%               ...
%              'dimLabelK', arrayIndicesDimK/samplingPointsDimK, ...)
%
%
% This is a method of class MrImage.
%
% NOTE: Y.select([]) returns a COPY of the current object Y.
%
%   PropertyName/Value pairs
%   General parameters
%   'type'      'index' (default) or 'label'
%               defines how selection elements should be
%               interpreted as array indices or sampling points
%   'invert'    true or false (default)
%               if true, selection is inverted and given arrayIndices are
%               excluded from returned selection, based on
%               dimInfo-information about whole array
%   'removeDims' true or false (default)
%               if true, singleton dimensions (1 or less dimensions) will be removed
%
% OUT
%   selectionIndexArray     indexArray of selected samples in original
%                           dimInfo
%   selectionImage          dimInfo of specified selection, derived as
%                           subset from dimInfo
%   unusedVarargin          if specified, returns all input arguments that
%                           did not match a dimLabel/Value pair of dimInfo
%                           Note: If this output is omitted, select returns
%                           an error for every unknown dimension requested
%
% EXAMPLE
%   select
%
%   See also MrImage MrDimInfo.select

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

returnUnusedVarargin = nargout >=3;
% extract removeDims info
defaults.removeDims = 0;
[args, argsSelect] = tapas_uniqc_propval(varargin, defaults);
tapas_uniqc_strip_fields(args);

if returnUnusedVarargin
    [selectionDimInfo, selectionIndexArray, unusedVarargin] = ...
        this.dimInfo.select(argsSelect{:});
else
    % be strict about wrong varargin, return error for unused select dims
    [selectionDimInfo, selectionIndexArray] = ...
        this.dimInfo.select(argsSelect{:});
end

selectionImage = this.copyobj;
selectionImage.data = selectionImage.data(selectionIndexArray{:});
selectionImage.dimInfo = selectionDimInfo.copyobj();

if removeDims
    selectionImage = selectionImage.remove_dims();
end

author: translationalneuromodeling

repo: tapas


text:
function [groupTable,out_no_group] = arrangerules(system,intTable)
%   Arranges the rules into groups function
%
%
%
%   Neelamugilan Gobalakrishnan
%   2006
%   $Revision:1.0 $ $Date: 15/03/2006 $

if size(system.output, 2) == 1
    out_no_group = 1;
else
    out_no_group = input('Enter the output number which is used to group the integer table : ');
end

row = 1;
for i=1:size(system.output(out_no_group).mf, 2)
    for j=1:size(intTable, 1)
        if intTable(j,size(system.input, 2)+out_no_group) == i
            grouped_intTable(row,1:size(system.input, 2))=intTable(j,1:size(system.input, 2));
            grouped_intTable(row,size(system.input, 2)+1)=intTable(j,size(system.input, 2)+out_no_group);
            row = row + 1;
        end
    end
end

%===========================================================
%   Display grouped integer tables with choosen output number.
%===========================================================

fprintf('\n');
disp('-------------------------------------------------------------');
fprintf(' ');
disp('                    Grouped integer table                    ');
fprintf(' ');
disp('-------------------------------------------------------------');

fprintf('\n');
disp('-------------------------------------------------------------');
fprintf(' ');

for i=1:size(system.input, 2)
    fprintf('Input%d ', i);
end

fprintf('Output%d ', out_no_group);

fprintf('\n');
disp('-------------------------------------------------------------');
fprintf('\n');
disp(grouped_intTable);
disp('-------------------------------------------------------------');

fprintf('\n\n');

groupTable = grouped_intTable;

author: Sable

repo: mcbench-benchmarks


text:
function y = sigmoidabTransform(x, transform, varargin)

% SIGMOIDABTRANSFORM Constrains a parameter to be between A and B
% by a scaled logistic sigmoid function.
%
% FORMAT
%
% DESC contains commands to constrain parameters to be between A
% and B via the sigmoid function, y=A+(B-A)/(1+exp(-x)).
%
% ARG x : input argument.
%
% ARG transform : type of transform, 'atox' maps a value into
% the transformed space (i.e. makes it between A and B). 'xtoa'
% maps the parameter back from transformed space to the original
% space. 'gradfact' gives the factor needed to correct gradients
% with respect to the transformed parameter, that is, it gives
% the gradient dx/da where x is the transformed parameter and a
% is the untransformed parameter.
%
% ARG transformsettings (first element of varargin): vector [A B]
% giving the minimum and maximum values A and B for the transformed
% parameter. If not given, assume A=0 and B=1 as in the function
% sigmoidTransform.
%
% OUTPUT y : return argument.
%
% SEEALSO : negLogLogitTransform, expTransform
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005, 2006, 2007
%
% COPYRIGHT : Jaakko Peltonen, 2011
%
% COPYRIGHT : Antti Honkela, 2012

% OPTIMI

% Get the transformation settings from the varargin structure
if ~isempty(varargin),
  transformsettings=varargin{1};
else
  transformsettings=[];
end;

identitymode=0;
if identitymode==1,
  % DEBUG ONLY, turns transformation to identity!
  fprintf(1,'Using identity mode for transform %s of parameter value %f\n',transform, x)
  switch transform
   case 'atox'
    y=x;
   case 'xtoa'
    y=x;
   case 'gradfact'
    y=0*x+1;
  end
else
  % normal mode, logistic sigmoid between A and B

  if isempty(transformsettings),
    warning(sprintf('sigmoidabTransform(%f,%s) called without transformation settings\n',x,transform));
    A=0; B=1;
  else
    A=transformsettings(1);
    B=transformsettings(2);
  end;

  limVal = 36;
  minval_sigmoid=A+(B-A)*eps;
  maxval_sigmoid=A+(B-A)*(1-eps);

  y = zeros(size(x));
  switch transform
   case 'atox'
    %x
    %A
    %B
    %fprintf(1,'transforming x %f, A %f, B %f\n',x,A,B);

    I1 = x < -limVal;
    y(I1) = A+(B-A)*eps;

    I2 = x > limVal;
    y(I2) = A+(B-A)*(1-eps);

    I3 = ~I1 & ~I2;
    y(I3) = A+(B-A)*sigmoid(x(I3));

   case 'xtoa'
    %  [x A B]
    %fprintf(1,'inverse-transforming x %f, A %f, B %f\n',x,A,B);
    I1 = x<=minval_sigmoid;
    y(I1) = -limVal;

    I2 = x>=maxval_sigmoid;
    y(I2) = limVal;

    I3 = ~I1 & ~I2;
    y(I3) = invSigmoid((x(I3)-A)/(B-A));

   case 'gradfact'
    %fprintf(1,'gradfact-transforming x %f, A %f, B %f\n',x,A,B);
    %y = (B-A)*x.*(1-x);
    y = (x-A).*(B-x)/(B-A);
    %y=0*x+1; % DEBUG ONLY, turns factors off!
  end
end;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function noiseBand = GetNoiseBand(view, scan)
%
% noiseBand = GetNoiseBand(view, scan)
%
% Gets the noise band, the frequency-domain indices of that
% portion of the power spectrum to be used to estimate the
% noise power. Returns value of zero to indicate default
% behavior: use ALL power-spectral samples to estimate noise.
%
% Ress, 8/30/02
global dataTYPES;
if ~exist('scan', 'var')
    if isfield(view, 'ui')
        scan = getCurScan(view);
    else
        scan = 1;
    end
end

noiseBand = 0;
if isfield(dataTYPES(view.curDataType).blockedAnalysisParams(scan), 'noiseBand')
  noiseBand = dataTYPES(view.curDataType).blockedAnalysisParams(scan).noiseBand;
end


text:
function [p,CI] = compareAUC(X1,X2,Y)
% -------------------------------------------------------------------------
% function [p,CI] = compareAUC(X1,X2,Y)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function calculates the significance (p-value) for the difference
% between two CORRELATED AUCs (i.e. computed from the same samples)
% computed for two different prediction situations. The implementation of
% this function is based on ref. [1].
% -------------------------------------------------------------------------
% REFERENCE:
% [1] DeLong, E. R., DeLong, D. M. & Clarke-Pearson, D. L. Comparing the
%     areas under two or more correlated receiver operating characteristic
%     curves: a nonparametric approach. Biometrics 44, 837–845 (1988).
% -------------------------------------------------------------------------
% INPUTS:
% 1. X1: [N X 1] vector of probabilities of outcome Y for situation 1.
% 2. X2: [N X 1] vector of probabilities of outcome Y for situation 2.
% 3. Y:  [N X 1] vector defining if each nth intance has the outcome or not
%        (1's or 0's), where N is the total number of patients.
% -------------------------------------------------------------------------
% OUTPUTS:
% 1. p: p-value for the difference in AUC.
% 2. CI: 95 % confidence interval for the difference in AUCs (diff), to be
%    subsequently applied as diff - CI and diff + CI.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2017
% -------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015-2017  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

% SANITY CHECKS
if numel(X1) ~= numel(Y) || numel(X2) ~= numel(Y)
    error('The number of instances in vectors X1, X2, and Y must be the same')
end
Y = logical(Y);

% PARAMETERS OR INTEREST
sigma = 1.96; % For 95 % confidence interval
L = [1,-1]; % Contrast vector. We are always comparing 2 variables at a time, not more. This is hardcoded for now and should stay as is.

% INITIALIZATION
X = [X1,X2]; % Combining the variables into a single matrix. We are always comparing 2 variables at a time, not more.
nVar = size(X,2); % Number of variables we compare (referred to "r" in ref. [1], so here, "v" --> "r" in ref. [1])
nInst = numel(Y); nPos = sum(Y); nNeg = nInst - nPos; % Total number of instances, number of positive instances ("nPos" --> "m" in ref. [1]) and number of negative instances ("nNeg" --> "n" in ref. [1]), respectively
Vpos = zeros(nPos,nVar); Vneg = zeros(nNeg,nVar); % Matrices "Vpos" --> "V10" and "Vneg" --> "V01" in ref. [1], respectively.
auc = zeros(1,nVar); % AUC estimates ("auc" --> "theta" in ref. [1])

% COMPUTATION OF AUC ESTIMATES ("auc" --> "theta"), VPOS ("Vpos" --> "V10") and VNEG ("Vneg" --> "V01")
Xpos = X(Y,:); Xneg = X(~Y,:); % Separating positive and negative instances ("Xpos" --> "X" and Xneg --> "Y" in ref. [1]. Here, our "Y" refers to the outcomes, or targets, to defined if an instance as the condition or not)
for v = 1:nVar
    for i = 1:nPos
        val = Xpos(i,v);
        Vpos(i,v) = sum(Xneg(:,v) < val) + 0.5*sum(Xneg(:,v) == val);
    end
    for i = 1:nNeg
        val = Xneg(i,v);
        Vneg(i,v) = sum(Xpos(:,v) > val) + 0.5*sum(Xpos(:,v) == val);
    end
end
auc = sum(Vpos)/(nNeg*nPos);
Vpos = Vpos/nNeg; Vneg = Vneg/nPos;

% COMPUTATION OF COVARIANCE MATRICES ("Spos" --> "S10" and "Sneg" --> "S01" in ref. [1])
Spos = (1/nPos) * ((Vpos'*Vpos) - nPos*(auc*auc'));
Sneg = (1/nNeg) * ((Vneg'*Vneg) - nNeg*(auc*auc'));
S = (1/nPos)*Spos + (1/nNeg)*Sneg;

% COMPUTATION OF CHI-SQUARE STATISTICS
LSL = L*S*L';
val_chi2 = (auc*L')*(inv(LSL))*(L*auc')*sqrt(2);
df_chi2 = rank(LSL);
p = 1-chi2cdf(val_chi2,df_chi2);
CI = sigma * sqrt(LSL);

end

author: mvallieres

repo: radiomics


text:
function val = isOctave()
%
%
%        val = isOctave()
%
%
%        Description: it checks if the running environment is MATLAB
%                     or Octave.
%
%        Input:
%
%        Output:
%           -val: a boolean value: 1 if the environment is Octave, otherwise 0.
%
%     Copyright (C) 2015  Francesco Banterle
%
%     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/>.
%

tmp = ver();

val = 0;

if(isfield(tmp, 'Name'))
    try
        val = (strcmp(tmp(1).Name, 'Octave') == 1);
    catch expr
        disp(expr);
    end
end

end

author: banterle


text:
classdef mtexFigure < handle
%
% a class to handles figures equaly sized axes with the following features
%
% * nicely resize axes
% * global or individual colorbars
% *
%
% Syntax
%   newMTEXFigure
%
%
% Class Properties
%  parent            - handle of the parent figure
%  children          - handles to all axes
%  cBarAxis          - handles to all colorbar axes
%  innerPlotSpacing  -
%  keepAspectRatio   -
%  nrows = 1         - number of rows
%  ncols = 1         - number of columns
%  axisWidth         - width of an individual axis
%  axisHeight        - height of an individual axis
%  cbx = 0           - colorbar width
%  cby = 0           - colorbar height
%  tightInset        - is added to axisSize
%  figTightInset     - is added to figSize
%  layoutMode        - 'auto' or 'user'
%
% Dependent Class Properties
%   currentAxes      - handle of the current axis
%   currentId        - id of te current axis
%   axesWidth        -
%   axesHeight       -
%   outerPlotSpacing -
%   dataCursorMenu   - handle of the data cursor context menu
%
% Description
%
% The calculation of the layout is initiated by the command
% <mtexFigure.drawNow.html drawNow>. This involves calls to the following
% functions:
%
%  drawNow
%    |  |
%    V  |
%  calcTightInset -> compute width of boundary around each axis
%       |
%       V
%  updateLayout
%       |
%       V
%  calcPartition -> compute partition (nrows, ncols)
%       |
%       V
%  calcAxesSize  -> compute axes size
%
% A mtexFigure may have the following children
%
% * <mapPlot.html mapPlot> -> micronbar
% * <sphericalPlot.html sphericalPlot> (stored in appdata of axes handle)
% * pfPlot [CS,SS,h]
% * ipdfPlot [CS]
% * MillerPlot [r,SS]
%
% See also
%

  properties
    parent            % the parent figure
    children          % the axes
    cBarAxis          % the colorbar axes
    innerPlotSpacing  %
    keepAspectRatio   %
    nrows = 1         % number of rows
    ncols = 1         % number of columns
    axisWidth         % width of an individual axis
    axisHeight        % height of an individual axis
    cbx = 0           % colorbar width
    cby = 0           % colorbar height
    tightInset = [0,0,0,0] % is added to axisSize
    figTightInset = [10,10,10,10] % is added to figSize
    layoutMode = 'auto' % set to user to fix it
    figSizeFactor = 0 % relative to the full screen
  end

  properties (Dependent = true)
    currentAxes       % current axis
    currentId         % current axis id
    axesWidth         %
    axesHeight        %
    outerPlotSpacing  %
    dataCursorMenu    % handle of the data cursor context menu
  end

  properties (Access = protected)
    cBarShift
  end

  methods

    function mtexFig = mtexFigure(varargin)

      % clear figure and set it up

      clf('reset');
      rmallappdata(gcf);
      set(gcf,'Visible','on');

      % set figure name
      if ~isempty(get_option(varargin,'FigureTitle'))
        set(gcf,'Name',get_option(varargin,'FigureTitle'));
      end

      % set custom resize function
      set(gcf,'ResizeFcn',@(src,evt) updateLayout(mtexFig));

      mtexFig.parent = gcf;

      mtexFig.outerPlotSpacing = get_option(varargin,'outerPlotSpacing',...
        getMTEXpref('outerPlotSpacing'));
      mtexFig.innerPlotSpacing = get_option(varargin,'innerPlotSpacing',...
        getMTEXpref('innerPlotSpacing'));
      mtexFig.keepAspectRatio = get_option(varargin,'keepAspectRatio',true);

      colrow = get_option(varargin,'layout',[1 1]);
      mtexFig.nrows = get_option(varargin,'nrows',colrow(1));
      mtexFig.ncols = get_option(varargin,'ncols',colrow(2));

      if check_option(varargin,{'nrows','ncols','layout'})
        mtexFig.layoutMode = 'user';
      end

      switch get_option(varargin,'figSize','','char')
        case 'huge'
          mtexFig.figSizeFactor = 1;
        case 'large'
          mtexFig.figSizeFactor = 0.8;
        case {'normal','medium'}
          mtexFig.figSizeFactor = 0.5;
        case 'small'
          mtexFig.figSizeFactor = 0.35;
        case 'tiny'
          mtexFig.figSizeFactor =  0.25;
        otherwise
          mtexFig.figSizeFactor = get_option(varargin,'figSize',0,'double');
      end

      set(mtexFig.parent,'color',[1 1 1],'nextPlot','replaceChildren');
      setappdata(mtexFig.parent,'mtexFig',mtexFig);

      varargin = delete_option(varargin,'color',1);
      varargin = delete_option(varargin,'position',1);
      optiondraw(mtexFig.parent,varargin{:});

      % set data cursor
      if check_option(varargin,'datacursormode')
        dcm_obj = datacursormode(mtexFig.parent);
        set(dcm_obj,'SnapToDataVertex','off')
        set(dcm_obj,'UpdateFcn',ensurecell(get_option(varargin,'datacursormode')));
        if ~check_option(varargin,'3d')
          datacursormode on;
        end
      end

      set(mtexFig.parent,'DefaultAxesCreateFcn', @updateChildren);
      set(mtexFig.parent,'DefaultAxesDeleteFcn', @deleteChildren);

      MTEXFigureMenu(mtexFig,varargin{:});

      h = findall(mtexFig.parent,'ToolTipString','Insert Colorbar');
      set(h,'ClickedCallback',@(a,b) mtexColorbar(mtexFig,a,b));

      function updateChildren(a,b)

        % prevent that this is called by colorbar
        x = dbstack;
        if any(strcmpi(x(2).name,{'colorbar','legend'})), return;end

        mtexFig.children = ...
        flipud(findobj(mtexFig.parent,'type','axes',...
        '-not','tag','Colorbar','-and','-not','tag','legend'));
      end

      function deleteChildren(a,b)
        mtexFig.children(mtexFig.children==a) = [];
      end

    end

    function ax = gca(mtexFig)
      % return current axis if exist otherwise create a new one

      if mtexFig.currentId == 0, mtexFig.currentId = 1;end

      ax = mtexFig.currentAxes;

    end

    % ---------------------------------------------------------

    function ax = get.currentAxes(mtexFig)
      ax = get(mtexFig.parent,'CurrentAxes');
    end

    function set.currentAxes(mtexFig,ax)
      try
        set(mtexFig.parent,'CurrentAxes',ax);
      end
    end

    function id = get.currentId(mtexFig)
      if isempty(mtexFig.currentAxes)
        id = 0;
      else
        id = find(mtexFig.currentAxes == mtexFig.children);
      end
    end

    function set.currentId(mtexFig,id)

      if id > numel(mtexFig.children) % create a nex axes
        for k = numel(mtexFig.children)+1:id
          axes('visible','off','parent',mtexFig.parent);
        end
        set(mtexFig.parent,'nextplot','add');

        if numel(mtexFig.children) > mtexFig.ncols * mtexFig.nrows
          mtexFig.ncols = ceil(numel(mtexFig.children) / mtexFig.nrows);
        end

      else
        set(mtexFig.parent,'CurrentAxes',mtexFig.children(id));
      end
    end

    function w = get.outerPlotSpacing(mtexFig)
      w = min(mtexFig.figTightInset);
    end

    function set.outerPlotSpacing(mtexFig,w)
      mtexFig.figTightInset = w * [1,1,1,1];
    end

    function aw = get.axesWidth(mtexFig)
      % the width of all axes is the number of columns times the width of
      % each single axis + inner and outer spacing

      aw = mtexFig.ncols * (mtexFig.axisWidth + sum(mtexFig.tightInset([1,3])))...
        + sum(mtexFig.figTightInset([1,3])) + ...
        (mtexFig.ncols-1) * mtexFig.innerPlotSpacing;
    end

    function ah = get.axesHeight(mtexFig)
      % the height of all axes is the number of rows times the height of
      % each single axis + inner and outer spacing

      ah = mtexFig.nrows * (mtexFig.axisHeight + sum(mtexFig.tightInset([2,4])))...
        + sum(mtexFig.figTightInset([2,4])) + ...
        (mtexFig.nrows-1) * mtexFig.innerPlotSpacing;
    end

    function setCamera(mtexFig,varargin)

      for a = 1:numel(mtexFig.children)
        setCamera(mtexFig.children(a),varargin{:});
      end
      mtexFig.drawNow;
    end

    function dcm = get.dataCursorMenu(mtexFig)

     dcm_obj = datacursormode(mtexFig.parent);
     dcm = dcm_obj.UIContextMenu;

    end

  end
end

author: mtex-toolbox

repo: mtex

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