text:
function axes_handle = PTKVisualiseIn3D(axes_handle, segmentation, smoothing_size, small_structures, reporting)
    % PTKVisualiseIn3D. Opens a new Matlab figure and visualises the
    % segmentation image in 3D.
    %
    %
    %     Inputs
    %     ------
    %
    %     axes_handle - specify an empty array [] to open a new figure, or the
    %         axes_handle of the existing axes for the visualisation
    %
    %     segmentation - a binary 3D PTKImage containing 1s for the voxels to
    %         visualise
    %
    %     smoothing_size - the amount of smoothing to perform. Higher
    %         smoothing makes the image look better but removes small
    %         structures such as airways and vessels
    %
    %     small_structures - set to true for improved visualisation of
    %         narrow structures such as airways and vessels
    %
    %     reporting - an object implementing CoreReportingInterface
    %         for reporting progress and warnings
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    if nargin < 3
        smoothing_size = 4; % A good value for lobes
    end

    if nargin < 4
        small_structures = false;
    end

    if nargin < 5
        reporting = CoreReportingDefault;
    end

    limit_to_one_component_per_index = false;
    minimum_component_volume_mm3 = 0;
    surface_colormap = CoreSystemUtilities.BackwardsCompatibilityColormap();

    axes_handle = MimVisualiseIn3D(axes_handle, segmentation, smoothing_size, small_structures, limit_to_one_component_per_index, minimum_component_volume_mm3, surface_colormap, reporting);
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [Population,FrontNo,CrowdDis] = EnvironmentalSelection(Population,N,Psi)
% The environmental selection of NSGA-II-conflict

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    Selected = zeros(1,N);
    FrontNo  = zeros(1,N);
    CrowdDis = zeros(1,N);
    PopObj   = Population.objs;
    for i = 1 : length(Psi)
        index = (i-1)*ceil(N/length(Psi))+1 : min(N,i*ceil(N/length(Psi)));
        [Selected(index),FrontNo(index),CrowdDis(index)] = SubSelection(PopObj(:,Psi{i}),length(index));
    end
    Population = Population(Selected);
end

function [Next,FrontNo,CrowdDis] = SubSelection(PopObj,N)
% Environmental selection based on only several objectives

    %% Non-dominated sorting
    [FrontNo,MaxFNo] = NDSort(PopObj,N);
    Next = FrontNo < MaxFNo;

    %% Calculate the crowding distance of each solution
    CrowdDis = CrowdingDistance(PopObj,FrontNo);

    %% Select the solutions in the last front based on their crowding distances
    Last     = find(FrontNo==MaxFNo);
    [~,Rank] = sort(CrowdDis(Last),'descend');
    Next(Last(Rank(1:N-sum(Next)))) = true;

    %% Population for next generation
    Next     = find(Next);
    FrontNo  = FrontNo(Next);
    CrowdDis = CrowdDis(Next);
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function sliderHandle = superSlider(varargin)

% SUPERSLIDER creates a more flexible slider, that allows for
% multiple slides, in place of the uicontrol('Style', 'Slider')
%
%   Developed to imitate a consistant, system independent version
%   of matlab's slider object, with the added capability of creating
%   multiple slides along a single body.
%
%   NOTE: This program was executed using functions rather than the more
%   intuitive use of classes, due to known slow matlab object oriented
%   processing time/ high system overhead. This means that there are few
%   applicable getters and setters, so properties should be handled upon
%   initiation of the function.
%
%   FORMAT: superSlider(handle, 'Property', Value)
%           Calling superSlider alone creates a basic slider with default
%           settings
%
%   PROPERTIES:
%   > handle: figure that the slider is placed in, otherwise it defaults to
%     gcf (use superSlider(figure) to generate a new window + slider)
%   > position: [left bottom width height], uses normalized units
%   > value: an array of length numSliders, containing the starting
%   position of each slide.
%   > tagName: Reference name of the slider object (return val alternative)
%   > min: starting value of the slider
%   > max: ending value of the slider
%   > stepSize: how much space between each arrow/slider click
%   > numSlides: the number of slides to create
%   > discrete: true/false - whether the slide is limitted to certain steps
%   > controlColor: color of the selectable controls [r g b]
%   > backgroundColor: color of the body of the slider [r g b]
%   > callback: a function that gets executed when the slider is used. It
%   is included in single quotations: superSlider('callback',
%   'functionName(Variables)');
%
%   DATA RETRIEVAL:
%   The return value of the function sliderHandle, is the parent to all the
%   slides and contains the values of each one in an array.
%   The value of the slides and the active slide can be found respectively
%   in the rows of the matrix stored in userData, accessed as follows:
%   > infoMatrix = get(sliderHandle, 'UserData')     which returns:
%                       [loc1, loc2, ...]            absolute location
%                       [ 0/1   0/1  ...]            a 1 indicates active
%   To get the handles to the slides (in ascending order):
%   > allSlides = get(sliderHandle, 'Children'); (empty if there was no
%   room for a visible slide in the given slider area)
%
%   Use of this code is encouraged for any interested party, but please
%   keep my name or cite me!
%   Created by: Danielle Ripsman
%               August 27, 2013
%               danielle.ripsman@mail.utoronto.ca

% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

props = []; % Object to contain the slider preferences
parts = []; %Object to contain slider features

% Sort out the handle/parent situation
if nargin > 0
    if ishandle(varargin{1})
        props.parent = varargin{1};
    else
        props.parent = gcf;
    end
end

% Match up inputs to their respective variable names
for ii = 1:nargin - 1
    if ischar(varargin{ii})
        switch lower(varargin{ii}) % Keeping things case insensitive
            case 'position'
                props.position = varargin{ii+1};
            case 'value'
                props.value = varargin{ii+1};
            case 'tagname'
                props.tagName = varargin{ii+1};
            case 'min'
                props.min = varargin{ii+1};
            case 'max'
                props.max = varargin{ii+1};
            case 'stepsize'
                props.stepSize = varargin{ii+1};
            case 'numslides'
                props.numSlides = varargin{ii+1};
            case 'discrete'
                props.discrete = varargin{ii+1};
            case 'controlcolor'
                props.controlColor = varargin{ii+1};
            case 'backgroundcolor'
                props.backgroundColor = varargin{ii+1};
            case 'callback'
                props.callback = varargin{ii+1};
        end
        ii = ii + 1;
    end
end

% Check + complete inputs
props = getDefaults(props);

%Find the size of the figure, in terms of pixels
[parts.pixWidth, parts.pixHeight] = getPixSize(props.parent);

% Determine the location and alignment of the slider components
[props, parts] = sliderLayout(props, parts);

% Body component
propsBod = makeBodySpecs(props, parts);
[sliderHandle, props.largeEnough] = makeSliderComponent(propsBod, false);
clear('propsBod');

% Arrow components
[propsDec, propsInc] = makeArrowSpecs(sliderHandle, props, parts);
decreaseHandle = makeSliderComponent(propsDec, true); % Left arrow button
increaseHandle = makeSliderComponent(propsInc, true);% Right arrow button
clear('propsDec'); clear('propsInc');

% Tagname for the uicontrol, acting as the slider handle
if isfield(props, 'tagName'); set(sliderHandle, 'tag', props.tagName); end

% Make sure the whole slider is deleted when the handle's deleted
set(sliderHandle, 'deleteFcn', {@grabArrows, decreaseHandle, increaseHandle});

% Determining a permissible number of slides
if ~isfield(props, 'numSlides') || ~props.largeEnough || props.numSlides < 1
    props.numSlides = 1;
end
props.numSlides = floor(props.numSlides); % In case a non-integer gets by

% Creating all requested/granted slides dynamically
sliderHandle = makeSlides(sliderHandle, props);

% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

function props = getDefaults(props)
if ~isfield(props, 'parent')
    props.parent = gcf;
end
if ~isfield(props, 'position')
    props.position = [.1 .1 .2 .04];
end
if ~isfield(props, 'controlColor')
    props.controlColor = [.85 .85 .85];
end
if ~isfield(props, 'backgroundColor')
    props.backgroundColor = [1 1 1];
end
if ~isfield(props, 'value');
    props.value = [];
else
    props.value = sort(props.value);
end
if ~isfield(props, 'min')
    props.min = 0;
end
if ~isfield(props, 'max') || props.max <= props.min
    props.max = props.min + 1;
end
if ~isfield(props, 'stepSize') || props.stepSize <= 0
    props.stepSize = props.max - props.min;
end
if ~isfield(props, 'discrete')
    props.discrete = false;
end
if ~isfield(props, 'callback')
    props.callback = '';
end

% Keep the arrows visible, by matching light to dark text/background
if sum(props.controlColor) < 1.5
    props.textColor = [1 1 1];
else
    props.textColor = [0 0 0];
end

% Other auto illustration details
props.borderwidth = 1; %how many pixels wide to make the borders
props.shadow = [.1 .1 .1];
props.highlight = [.65 .65 .65];
props.secondHighlight = [.4 .4 .4];

% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

function [pixWidth, pixHeight] = getPixSize(fig)
% Pixel size, since some built-in features are measured in pixels
preserveUnits = get(fig, 'units');
set(fig, 'units', 'pixels');
figSize = get(fig, 'position');
pixWidth = figSize(3);
pixHeight = figSize(4);
set(fig, 'units', preserveUnits); %return object to initial units

% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

function [props, parts] = sliderLayout(props, parts)
position = props.position;

% Configuration details of the slider, based on the input position
if position(3) >= position(4)
    % Make a HORIZONTAL slider:
    parts.dSymbol = '<';
    parts.iSymbol = '>';
    props.hzn = true;
    constDir = 4; % Height of the component = the slider height
    constPix = parts.pixWidth;
    varDir = 3; % Width is adjustable
    varPix = parts.pixHeight;
else
    % Make a VERTICAL slider:
    parts.dSymbol = 'v';
    parts.iSymbol = '^';
    props.hzn = false;
    constDir = 3; % Width of the component = the slider width
    constPix = parts.pixHeight;
    varDir = 4; % Height is adjustable
    varPix = parts.pixWidth;
end

% Try to create square buttons
if 2*position(constDir)*varPix < position(varDir)*constPix
    size = position(constDir)*varPix/constPix;
else
    size = position(varDir)/2;
end

% Allocate the remaining area properly between the components
posDecreaseArrow = position;
posDecreaseArrow(varDir) = size;
posBody = position;
posBody(varDir-2) = posBody(varDir-2) + size - props.borderwidth/constPix;
posBody(varDir) = posBody(varDir) - 2*size + 2*props.borderwidth/constPix;
posIncreaseArrow = position;
posIncreaseArrow(varDir-2) = posIncreaseArrow(varDir-2)+posIncreaseArrow(varDir)-size;
posIncreaseArrow(varDir) = size;

% Asign to structure for easier passing
parts.posDArrow = posDecreaseArrow;
parts.posBody = posBody;
parts.posIArrow = posIncreaseArrow;

props.constDir = constDir;
props.varDir = varDir;

% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

function [propsBod] = makeBodySpecs(props, parts)
propsBod = props;
propsBod.position = parts.posBody;
propsBod.sym = '';
propsBod.bcolor = props.backgroundColor;
propsBod.color = props.backgroundColor;
propsBod.shadowcolor = props.highlight;
propsBod.callback = '';

% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

function [propsDec, propsInc] = makeArrowSpecs(sliderHandle, props, parts)
% For the decrease arrow:
propsDec = props;
propsDec.position = parts.posDArrow;
propsDec.sym = parts.dSymbol;
propsDec.bcolor = props.secondHighlight;
propsDec.color = props.controlColor;
propsDec.shadowcolor = props.shadow;
propsDec.callback = {@arrowCallback, sliderHandle, props, -1};

% For the increase arrow:
propsInc = props;
propsInc.position = parts.posIArrow;
propsInc.sym = parts.iSymbol;
propsInc.bcolor = props.secondHighlight;
propsInc.color = props.controlColor;
propsInc.shadowcolor = props.shadow;
propsInc.callback = {@arrowCallback, sliderHandle, props, 1};

% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

function [panel, largeEnough] = makeSliderComponent(props, uiOverlay)
% If the middle of the slide is unused, it shouldn't cause errors.
largeEnough = false;
width = props.borderwidth;

% Panel background for better visual
panel = uipanel('parent', props.parent,...
            'Units','normalized',...
            'backgroundColor', props.bcolor,...
            'bordertype', 'beveledout',...
            'shadowcolor', props.shadowcolor,...
            'highlightcolor', props.highlight,...
            'borderWidth', width,...
            'Position',props.position,...
            'Visible', 'on');

[pixWidth, pixHeight] = getPixSize(panel);

% Set inner region, so the button region doesn't interfere with the frame
innerPos = [0+width/pixWidth, 0+width/pixHeight,...
            1-2*width/pixWidth, 1-2*width/pixHeight];

if innerPos(3) > 0 && innerPos(4)>0
   largeEnough = true;
   if uiOverlay
   % Clickable+text component of the button
   uicontrol('Parent', panel,...
             'Units','normalized',...
             'Style','text',...
             'Enable', 'inactive',...
             'FontUnits', 'Normalized',...
             'FontName', 'Webdings',...
             'FontSize', .6,...
             'HorizontalAlignment', 'center',...
             'Position', innerPos,...
             'ForegroundColor', props.textColor,...
             'background', props.color,...
             'String', props.sym,...
             'FontWeight', 'bold',...
             'Hittest', 'on',...
             'buttonDownFcn', props.callback);
   end
end

% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

function grabArrows(hObject, eventData, ar1, ar2)
% Ensures the entire slider is erased as a single object
if ishandle(ar1)
    delete(ar1);
end
if ishandle(ar2)
    delete(ar2);
end

% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

function sliderHandle = makeSlides(sliderHandle, props)
%if there's enough room for slides
if props.largeEnough
    propsSlide = makeSlideSpecs(sliderHandle, props);
    propsSlide.position = zeros(1, 4);
    numStepsAbs = (props.max - props.min)/props.stepSize + 1;
    numSlots = numStepsAbs + props.numSlides - 1;
    size = 1/numSlots;

    varPos = props.varDir - 2; %width/height component of position array (3/4)
    propsSlide.position(props.varDir) = size;
    propsSlide.position(props.constDir) = 1;

    sliderRef = zeros(2, props.numSlides);
    % sliderRef stores slide info - [loc1, loc2, ...]
    %                               [ 0/1   0/1  ...]

    for ii = 1:props.numSlides
        if isempty(props.value)
            % Fit all sliders on in a row, if unspecified
            propsSlide.position(varPos) = (ii - 1)/numSlots;
        else
            dist = (props.value(ii) - props.min)*size/props.stepSize;
            propsSlide.position(varPos) = dist + (ii - 1)*size;
        end
        slideHandle = makeSliderComponent(propsSlide, true);
        loc = get(slideHandle, 'position');
        sliderRef(1, ii) = (loc(varPos) - size*(ii - 1))/size*props.stepSize + props.min;
    end

    % Reference the children in ascending order:
    slideRef = get(sliderHandle, 'Children');
    set(sliderHandle, 'Children', fliplr(slideRef')');
    set(sliderHandle, 'userData', sliderRef);

    slideSelect(get(slideRef(props.numSlides), 'Children'), sliderHandle, props);

    clear('propsSlide');
else
    if isempty(props.value); props.value = props.min; end
    set(sliderHandle, 'userData', [props.value; 1]);
end

% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

function [propsSlide] = makeSlideSpecs(sliderHandle, props)
propsSlide = props;
propsSlide.parent = sliderHandle;
propsSlide.sym = '';
propsSlide.bcolor = props.secondHighlight;
propsSlide.color = props.controlColor;
propsSlide.shadowcolor = props.shadow;
propsSlide.callback = {@slideDrag, sliderHandle, props};

% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

function arrowCallback(hObject, eventData, sliderHandle, props, direction)
set(hObject, 'backgroundcolor', props.textColor);

if props.largeEnough
    slides = get(sliderHandle, 'Children');
    activityMatrix = get(sliderHandle, 'userData');
    activeSlide = slides(activityMatrix(2,:)>0);
    index = find(activityMatrix(2,:));

    location = get(activeSlide, 'Position');

    size = location(props.varDir);
    varPos = props.varDir - 2;

    location(varPos) = location(varPos) + direction*size;

    % Keep the slide in bounds
    location = inBounds(location, varPos, index, size, slides);

    set(activeSlide, 'Position', location);
    sliderRef = get(sliderHandle, 'userData');
    sliderRef(1, index) = (location(varPos) - size*(index - 1))/size*props.stepSize + props.min;
    set(sliderHandle, 'userData', sliderRef);
else
    location = get(sliderHandle, 'userData');
    varPos = 1;
    location(varPos) = location(varPos) + direction*props.stepSize;
    postSlide = props.max;
    preSlide = props.min;
    if location(varPos) < preSlide; location(varPos) = preSlide; end
    if location(varPos) > postSlide; location(varPos) = postSlide; end
    set(sliderHandle, 'userData', location);
end

try % Optional user callback function
   eval(props.callback);
end

pause(.02); %pause for button color effect
set(hObject, 'backgroundcolor', props.controlColor);

% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

function location = inBounds(location, varPos, index, size, slides)
% makes sure the slide doesn't exceed its bounds, or surrounding slides
if index == 1
    preSlide = 0;
else
    presliderPos = get(slides(index - 1), 'Position');
    preSlide = presliderPos(varPos) + size;
end

if index == length(slides)
    postSlide = 1 - size;
else
    postsliderPos = get(slides(index + 1), 'Position');
    postSlide = postsliderPos(varPos) - size;
end

% Keep slider in bounds
if location(varPos) < preSlide; location(varPos) = preSlide; end
if location(varPos) > postSlide; location(varPos) = postSlide; end

% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

function slideDrag(hObject, eventData, sliderHandle, props)
% Tracks where the slide was selected, for consistent dragging
oldUnits = get(gcf, 'units');
set(gcf, 'units', 'pixels');
grabpoint = get(gcf,'CurrentPoint');
set(get(hObject, 'parent'), 'userData', grabpoint);
set(gcf, 'units', oldUnits);

% Handles the slide highlighting/selection update
slideSelect(hObject, sliderHandle, props);

% Uses the figure callback mechanisms to control the interactive slide placement
oldMotionFunction = get(gcf, 'WindowButtonMotionFcn');
set(gcf, 'WindowButtonMotionFcn', {@motionUpdate, sliderHandle, get(hObject, 'parent'), props});
set(gcf, 'WindowButtonUpFcn', {@stopDragging, oldMotionFunction, oldUnits});

% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

function slideSelect(slideHandle, sliderHandle, props)
set(slideHandle, 'backgroundcolor', props.controlColor);
allSlides = get(get(get(slideHandle, 'parent'), 'parent'), 'Children');
sliderRef = zeros(1, length(allSlides));
for ii = 1:length(allSlides) %Goes through all slides, order: largest to smallest
    if get(allSlides(ii), 'Children') ~= slideHandle
        set(get(allSlides(ii), 'Children'), 'backgroundcolor', props.textColor);
    else
        sliderRef(ii) = 1;
    end
    sliderInfo = get(sliderHandle, 'userData');
    sliderInfo(2,:) = sliderRef;
    set(sliderHandle, 'userData', sliderInfo);
end

% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

function motionUpdate(hObject, eventData, sliderHandle, slideHandle, props)
[pixWidth, pixHeight] = getPixSize(sliderHandle);
oldUnits = get(gcf, 'units');
set(gcf, 'units', 'pixels');
newPoint = get(gcf, 'currentPoint');
position = get(slideHandle, 'position');
grabPoint = get(slideHandle, 'userData');

size = position(props.varDir);
varPos = props.varDir - 2;
if props.hzn; convertPix = pixWidth; else convertPix = pixHeight; end

step = (newPoint(varPos) - grabPoint(varPos))/convertPix;
original = position(varPos);

if ~props.discrete %if intermediate points are allowed:
    position(varPos) = position(varPos) + step;
else
    round(step/size);
    position(varPos) = position(varPos) + round(step/size)*size;
end

sliderRef = get(sliderHandle, 'userData');
index = find(sliderRef(2,:));
slides = get(sliderHandle, 'Children');

position = inBounds(position, varPos, index, size, slides);

shift = (position(varPos) - original)*convertPix;
newPoint(varPos) = grabPoint(varPos) + shift;

newPosAbs = (position(varPos) - size*(index - 1))/size*props.stepSize;
sliderRef(1, index) = newPosAbs + props.min;

set(sliderHandle, 'userData', sliderRef);
set(slideHandle, 'position', position);
set(slideHandle, 'userData', newPoint);

set(gcf, 'units', oldUnits);

try % Optional user callback function
    eval(callback);
end

% -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

function stopDragging(hObject, eventData, oldMotionFunction, oldUnits)
set(gcf, 'WindowButtonMotionFcn', oldMotionFunction, 'units', oldUnits);

author: Sable

repo: mcbench-benchmarks


text:
function STMuNet_test(varargin)

if ~isempty(gcp('nocreate'))
    delete(gcp)
end

opts.train.gpus = [ 1 ] ;

opts = cnn_setup_environment(opts);
opts.nSplit =  1 ;

% opts.dataSet = 'hmdb51';
opts.dataSet = 'ucf101';

model = ['ST-MulNet-img50-flow50-final-split=' num2str(opts.nSplit)];
model = ['ST-MulNet-img50-flow152-split=' num2str(opts.nSplit)];

opts.train.memoryMapFile = fullfile(tempdir, 'ramdisk', ['matconvnet' num2str(opts.nSplit ) '.bin']) ;
opts.dataDir = fullfile(opts.dataPath, opts.dataSet) ;
opts.splitDir = 'ucf101_splits'; nClasses = 101;
opts.imdbPath = fullfile(opts.dataDir, [opts.dataSet '_split' num2str(opts.nSplit) 'imdb.mat']);
opts.model = fullfile(opts.modelPath, [opts.dataSet model '.mat']) ;
opts.expDir = fullfile(opts.dataDir, [opts.dataSet '-' model]) ;
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.train.saveAllPredScores = 1;
opts.train.denseEval = 1;
opts.train.plotDiagnostics = 0 ;
opts.train.continue = 1 ;
opts.train.prefetch = 1 ;
opts.train.expDir = opts.expDir ;
opts.train.numAugments = 1;
opts.train.frameSample = 'random';

opts = vl_argparse(opts, varargin) ;

% -------------------------------------------------------------------------
%                                                   Database initialization
% -------------------------------------------------------------------------
if exist(opts.imdbPath)
  imdb = load(opts.imdbPath) ;
  imdb.flowDir = opts.flowDir;
else
  imdb = cnn_setup_data(opts) ;
  save(opts.imdbPath, '-struct', 'imdb', '-v6') ;
end

% -------------------------------------------------------------------------
%                                                    Network initialization
% -------------------------------------------------------------------------
if ~exist(opts.model,'file')
  [~, baseModel] = fileparts(opts.model);
  fprintf('Downloading base model file: %s ...\n', baseModel);
  mkdir(fileparts(opts.model)) ;
  urlwrite(...
  ['http://ftp.tugraz.at/pub/feichtenhofer/st-mul/final/' baseModel '.mat'], ...
    opts.model) ;
end
net = load(opts.model) ;

if isfield(net, 'net'), net=net.net;end
net = dagnn.DagNN.loadobj(net);

opts.train.augmentation = 'f25noCtr';
opts.train.frameSample = 'temporalStrideRandom';
opts.train.nFramesPerVid = 1;
opts.train.temporalStride = 1:15;
opts.train.valmode = 'temporalStrideRandom';
opts.train.numValFrames = 25 ;
opts.train.saveAllPredScores = 1 ;
opts.train.denseEval = 1;
opts.train.temporalFullConvTest = 1;

opts.train.train = NaN;

net.rebuild() ;

net.conserveMemory = 1 ;
fn = getBatchWrapper_rgbflow(net.meta.normalization, opts.numFetchThreads, opts.train) ;
[info] = cnn_train_dag(net, imdb, fn, opts.train) ;

author: feichtenhofer

repo: st-resnet


text:
%% BASIC WAYPOINT EVENT (waypointEvent.m) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This class is designed as a child class to a generic "event", used to
% define the characteristics of a 'waypoint achieved' event.

% Author: James A. Douthwaite

classdef waypointEvent < interactionEvent

%% INITIALISE THE WAYPOINT EVENT SPECIFIC PARAMETERS
    properties
        % Waypoint unique properties
    end
%%  CLASS METHODS
    methods
        % CONSTRUCTION METHOD FOR THE WAYPOINT EVENT
        function obj = waypointEvent(time,objectA,objectB,summaryInfo,enum)
            % INPUTS:
            % time          - Detection time (s)
            % objectA       - Class for the reference object
            % objectIDB     - Class for the second object
            % summaryInfo   - Event description string
            % enum          - Provided alternative enum (nullified?)
            % OUTPUTS:
            % obj           - The waypoint event object

            % INITIALISE THE GENERIC EVENT SUPERCLASS
            obj = obj@interactionEvent(time,'WAYPOINT',objectA,objectB,summaryInfo);

            % ASSIGN AN EVENT NUMERATION
%             if ~exist('enum','var') || ~isa(enum,'uint8')
            if nargin < 5 || ~isa(enum,'uint8')
                obj.type = eventType.waypoint;
            else
                obj.type = enum;
            end
        end
    end
end

repo: OpenMAS


text:
function mtrThresholdResampleSISPathways(fgFilenames,samplerOptsFilename,weightThreshold,dt6Filename,fgResampledFilename,bCalcStats)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% mtrThresholdResampleSISPathways(fgFilenames,samplerOptsFilename,
%%                                numSamplesPerCluster,clusterROI,
%%                                fgResampledFilename)
%%
%% fgFilenames - filenames for metrotrac pathways
%% samplerOptsFilename - parameters file for collecting paths and stats
%% fgResampledFilename - filename for resampled pathways
%% dt6Filename - Filename for dt6.
%% hintNumTotalFibers - Number hopefully greater than the total number of
%% fibers that we will load.  This will speed up processing significantly
%% as it will allow us to preallocate space for validity matrix.
%% bCalcStats - True then calculate the statistics False then they already
%% exist.
%%
%%
%% Author: Anthony Sherbondy
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% MUST MAKE SURE THIS PROGRAM IS COMPILED AND IN RIGHT LOCATION
if(ispc)
    executable = which('updatestats.exe');
else
    error('Not compiled for Linux.');
end

% Get all the needed files
if(ieNotDefined('dt6Filename'))
    [f, p] = uigetfile({'*.mat';'*.*'}, 'Load DT6 file...');
    if(isnumeric(f)), disp('Resampling canceled.'); return; end
    dt6Filename = fullfile(p,f);
end

dt6 = load(dt6Filename);

% Get importance weight file that we need to resample with also load the
% pathfile into matlab
totalFiberCount = 1;
fg_input = dtiNewFiberGroup;
fg_temp = dtiNewFiberGroup;

for ii = 1:length(fgFilenames)
     tic;
    statvec_names{ii} = sprintf('temp%d_statvec_iw.dat',ii);
    if(bCalcStats)
        disp('Calculating Statistics.');
        args = sprintf(' -v %s -i %s -o temp.dat -s 4 -sf %s', samplerOptsFilename, fgFilenames{ii}, statvec_names{ii});
        cmd = [executable args];
        disp(cmd); disp('...')
        [s, ret_info] = system(cmd,'-echo');
        disp('Done.')
    end

    disp('Loading pathway file into matlab...')
    fg_temp = mtrImportFibers(fgFilenames{ii}, dt6.xformToAcPc);
    disp('Done')

    % Keep aggregating all of the fibers together into one group
    if(isempty(fg_input.fibers))
        fg_input.fibers = fg_temp.fibers;
    else
        fg_input.fibers = [fg_input.fibers(:), fg_temp.fibers(:)];
    end
     toc;
end
clear fg_temp;

% Send the valid matrix in with the resample computations
[fg_output, sample_w] = resampleThresholdSISPathways(weightThreshold, statvec_names, fg_input);

% Output the pathways in dtifiberUI space
%dtiWriteFiberGroup(fg_output,fgResampledFilename);

% Output the pathways in metrotrac space
mtrExportFiberGroupToMetrotrac(fgResampledFilename,fg_output,dt6.t1NormParams,dt6.mmPerVox,dt6.xformToAcPc,size(dt6.b0));

% Cleanup temp files
delete('temp.dat');

% Why should I delete these stats?? they are yummy!!
% for ii = 1:length(statvec_names)
%     delete(statvec_names{ii});
% end

return;

function [fg_output, sample_w_out] = resampleThresholdSISPathways(weightThreshold, weightFilenames, fg_input)

numFiles = length(weightFilenames);
sizeFiles = [];

% Setup the output file
fg_output = dtiNewFiberGroup;

% Load multiple importance weight files
iw = [];
for f = 1:numFiles
    iwtemp = mtrLoadStatvec(weightFilenames{f});
    iw = [iw; iwtemp];
    sizeFiles(f) = length(iwtemp);
end
% Now just threshold the paths based on weightThreshold
iwThresholded = find(iw > weightThreshold);

% Get the paths that survive the threshold
fg_output.fibers = fg_input.fibers(iwThresholded);
sample_w_out = ones(size(fg_output.fibers));

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%EDGELIST Return list of edge pixels for region
%
% EG = EDGELIST(IM, SEED) is a list of edge pixels (Nx2) of a region in the
% image IM starting at edge coordinate SEED=[X,Y].  The edgelist has one row per
% edge point coordinate (x,y).
%
% EG = EDGELIST(IM, SEED, DIRECTION) as above, but the direction of edge
% following is specified.  DIRECTION == 0 (default) means clockwise, non
% zero is counter-clockwise.  Note that direction is with respect to y-axis
% upward, in matrix coordinate frame, not image frame.
%
% [EG,D] = EDGELIST(IM, SEED, DIRECTION) as above but also returns a vector
% of edge segment directions which have values 1 to 8 representing W SW S SE E
% NW N NW respectively.
%
% Notes::
% - Coordinates are given assuming the matrix is an image, so the indices are
%   always in the form (x,y) or (column,row).
% - IM is a binary image where 0 is assumed to be background, non-zero
%   is an object.
% - SEED must be a point on the edge of the region.
% - The seed point is always the first element of the returned edgelist.
% - 8-direction chain coding can give incorrect results when used with
%   blobs founds using 4-way connectivty.
%
% Reference::
% - METHODS TO ESTIMATE AREAS AND PERIMETERS OF BLOB-LIKE OBJECTS: A COMPARISON
%   Luren Yang, Fritz Albregtsen, Tor Lgnnestad and Per Grgttum
%   IAPR Workshop on Machine Vision Applications Dec. 13-15, 1994, Kawasaki
%
% See also ILABEL.

% Copyright (C) 1993-2014, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function [e,d] = edgelist(im, P, direction)

    % deal with direction argument
    if nargin == 2
        direction = 0;
    end

    if direction == 0
        neighbours = [1:8]; % neigbours in clockwise direction
    else
        neighbours = [8:-1:1];  % neigbours in counter-clockwise direction
    end

    P = P(:)';
    P0 = P;     % make a note of where we started
    pix0 = im(P(2), P(1));  % color of pixel we start at

    % find an adjacent point outside the blob
    Q = adjacent_point(im, P, pix0);
    if isempty(Q)
        error('no neighbour outside the blob');
    end

    e = P;  % initialize the edge list
    dir = []; % initialize the direction list

    % these are directions of 8-neighbours in a clockwise direction
    dirs = [-1 0; -1 1; 0 1; 1 1; 1 0; 1 -1; 0 -1; -1 -1];

    while 1
        % find which direction is Q
        dQ = Q - P;
        for kq=1:8
            if all(dQ == dirs(kq,:))
                break;
            end
        end

        % now test for directions relative to Q
        for j=neighbours
            % get index of neighbour's direction in range [1,8]
            k = j + kq;
            if k > 8
                k = k - 8;
            end
            dir = [dir; k];

            % compute coordinate of the k'th neighbour
            Nk = P + dirs(k,:);
            try
                if im(Nk(2), Nk(1)) == pix0
                    % if this neighbour is in the blob it is the next edge pixel
                    P = Nk;
                    break;
                end
            end
            Q = Nk;     % the (k-1)th neighbour
        end

        % check if we are back where we started
        if all(P == P0)
            break;
        end

        % keep going, add P to the edgelist
        e = [e; P];
    end

    if nargout > 1
        d = dir;
    end
end

function P = adjacent_point(im, seed, pix0)
    % find an adjacent point not in the region
    dirs = [1 0; 0 1; -1 0; 0 -1];
    for d=dirs'
        P = [seed(1)+d(1), seed(2)+d(2)];
        try
            if im(P(2), P(1)) ~= pix0
                return;
            end
        catch
            % if we get an exception then by definition P is outside the region,
            % since it's off the edge of the image
            return;
        end
    end
    P = [];
end

author: Allopart

repo: rbpf-gmapping


text:
path_to_scatnet = fileparts(mfilename('fullpath'));

addpath(fullfile(path_to_scatnet, 'classification'));
addpath(fullfile(path_to_scatnet, 'convolution'));
addpath(fullfile(path_to_scatnet, 'core'));
addpath(genpath(fullfile(path_to_scatnet, 'demo')));
addpath(fullfile(path_to_scatnet, 'display'));
addpath(fullfile(path_to_scatnet, 'filters'));
addpath(fullfile(path_to_scatnet, 'filters/selesnick'));
addpath(fullfile(path_to_scatnet, 'utils'));
addpath(genpath(fullfile(path_to_scatnet, 'papers')));
addpath(fullfile(path_to_scatnet, 'scatutils'));
addpath(genpath(fullfile(path_to_scatnet, 'unittest')));
addpath(fullfile(path_to_scatnet, 'utils'));
addpath(fullfile(path_to_scatnet, 'reconstruction'));

clear path_to_scatnet;

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
% loss_hsp - hyperbolic secant loss function
%
% Copyright(c) 2009-2011 Ryota Tomioka
%              2009      Stefan Haufe
% This software is distributed under the MIT license. See license.txt
function [floss, gloss]=loss_hsp(zz, bb)

% floss = -sum(log(sech(bb-zz)./pi));
% gloss = tanh(zz-bb);

zz = zz-bb;
mz = abs(zz);

floss = sum(mz + log(exp(zz-mz)+exp(-zz-mz))-log(2/pi));
gloss = (exp(zz-mz)-exp(-zz-mz))./(exp(zz-mz)+exp(-zz-mz));

repo: REST


text:
classdef BodyFixedOrbitVariable < AbstractOrbitModelVariable
    %BodyFixedOrbitVariable Summary of this class goes here
    %   Detailed explanation goes here

    properties
        varObj(1,1) BodyFixedOrbitStateModel = BodyFixedOrbitStateModel.getDefaultOrbitState();

        lb(1,6) double = 0;
        ub(1,6) double = 0;

        varLat(1,1) logical = false;
        varLong(1,1) logical = false;
        varAlt(1,1) logical = false;
        varBFAz(1,1) logical = false;
        varBFEl(1,1) logical = false;
        varBFMag(1,1) logical = false;
    end

    methods
        function obj = BodyFixedOrbitVariable(varObj)
            obj.varObj = varObj;
            obj.varObj.optVar = obj;

            obj.id = rand();
        end

        function x = getXsForVariable(obj)
            x = [];

            if(obj.varLat)
                x(end+1) = obj.varObj.lat;
            end

            if(obj.varLong)
                x(end+1) = obj.varObj.long;
            end

            if(obj.varAlt)
                x(end+1) = obj.varObj.alt;
            end

            if(obj.varBFAz)
                x(end+1) = obj.varObj.vVectNEZ_az;
            end

            if(obj.varBFEl)
                x(end+1) = obj.varObj.vVectNEZ_el;
            end

            if(obj.varBFMag)
                x(end+1) = obj.varObj.vVectNEZ_mag;
            end
        end

        function [lb, ub] = getBndsForVariable(obj)
            useTf = obj.getUseTfForVariable();

            lb = obj.lb(useTf);
            ub = obj.ub(useTf);
        end

        function [lb, ub] = getAllBndsForVariable(obj)
            lb = obj.lb;
            ub = obj.lb;
        end

        function setBndsForVariable(obj, lb, ub)
            if(length(lb) == 6 && length(ub) == 6)
                obj.lb = lb;
                obj.ub = ub;
            else
                useTf = obj.getUseTfForVariable();

                obj.lb(useTf) = lb;
                obj.ub(useTf) = ub;
            end
        end

        function useTf = getUseTfForVariable(obj)
            useTf = [obj.varLat        obj.varLong        obj.varAlt        obj.varBFAz        obj.varBFEl        obj.varBFMag];
        end

        function setUseTfForVariable(obj, useTf)
            obj.varLat  = useTf(1);
            obj.varLong = useTf(2);
            obj.varAlt  = useTf(3);
            obj.varBFAz = useTf(4);
            obj.varBFEl  = useTf(5);
            obj.varBFMag = useTf(6);
        end

        function updateObjWithVarValue(obj, x)
            xInd = 1;

            if(obj.varLat)
                obj.varObj.lat = x(xInd);
                xInd = xInd + 1;
            end

            if(obj.varLong)
                obj.varObj.long = x(xInd);
                xInd = xInd + 1;
            end

            if(obj.varAlt)
                obj.varObj.alt = x(xInd);
                xInd = xInd + 1;
            end

            if(obj.varBFAz)
                obj.varObj.vVectNEZ_az = x(xInd);
                xInd = xInd + 1;
            end

            if(obj.varBFEl)
                obj.varObj.vVectNEZ_el = x(xInd);
                xInd = xInd + 1;
            end

            if(obj.varBFMag)
                obj.varObj.vVectNEZ_mag = x(xInd);
                xInd = xInd + 1; %#ok<NASGU>
            end
        end

        function nameStrs = getStrNamesOfVars(obj, evtNum, varLocType)
            nameStrs = {sprintf('Event %i Latitude', evtNum), ...
                        sprintf('Event %i Longitude', evtNum), ...
                        sprintf('Event %i Altitude', evtNum), ...
                        sprintf('Event %i Body-Fixed Velocity Azimuth', evtNum), ...
                        sprintf('Event %i Body-Fixed Velocity Elevation', evtNum), ...
                        sprintf('Event %i Body-Fixed Velocity Magnitude', evtNum)};

            nameStrs = nameStrs(obj.getUseTfForVariable());
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function tapas_BLT_initCogent(store)
% tapas_BLT_initCogent starts cogent which is needed for the response box
%
% About:
%   Coded by F.Petzschner 19. April 2017
%	Last change: F.Petzschner 19. April 2017

if  ~store.keyboard
config_serial(store.scanner.boxport, 19200, 0, 0, 8);
start_cogent;
end

end

author: translationalneuromodeling

repo: tapas


text:
function result = sphere_unit_07_3d ( func )

%*****************************************************************************80
%
%% SPHERE_UNIT_07_3D approximates surface integrals on the unit sphere in 3D.
%
%  Integration region:
%
%    Points (X,Y,Z) such that:
%
%      X**2 + Y**2 + Z**2 = 1.
%
%  Discussion:
%
%    A 32 point 7-th degree formula is used.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 May 2004
%
%  Author:
%
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    Arthur H Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971.
%
%  Parameters:
%
%    Input, external FUNC, the name of the user supplied
%    function which evaluates F(X,Y,Z), of the form
%      function value = func ( x, y, z )
%
%    Output, real RESULT, the approximate integral of the function.
%
  norder1 = 2;
  norder2 = 4;
  norder3 = 4;
%
%  Set XTAB1 and WATE1.
%
  xtab1(1) = -1.0E+00;
  xtab1(2) =  1.0E+00;
  weight1(1) = 1.0E+00;
  weight1(2) = 1.0E+00;
%
%  Set XTAB2 and WATE2.
%
  for j = 1 : norder2
    angle = pi * ( 2 * j - 1 ) / ( 2 * norder2 );
    xtab2(j) = cos ( angle );
  end

  weight2(1:norder2) = 1.0E+00 / ( 4 * norder2 );
%
%  Set XTAB3 and WATE3.
%
  [ xtab3, weight3 ] = legendre_set ( norder3 );

  quad = 0.0E+00;
  for i = 1 : norder1
    for j = 1 : norder2
      for k = 1 : norder3

        x = xtab1(i) * sqrt ( 1.0E+00 - xtab2(j)^2 ) ...
                     * sqrt ( 1.0E+00 - xtab3(k)^2 );
        y = xtab1(i) * xtab2(j) * sqrt ( 1.0E+00 - xtab3(k)^2 );
        z = xtab1(i) * xtab3(k);

        quad = quad + weight1(i) * weight2(j) * weight3(k) ...
          * feval ( func, x, y, z );

      end
    end
  end

  volume = sphere_unit_area_3d ( );
  result = quad * volume;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [r2,mn,st]= cv_mean(a,r,loss_type,param1)

if nargin<3
  [r2,st]=get_mean2(r);
  [mn,st]=get_mean2(r,[],1);
else
  [r2,st]=get_mean2(r,loss_type);
  [mn,st]=get_mean2(r,loss_type,1);
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [TableChecks, Table_csources, CSourcesTestedRxns, TestSolutionNameOpenSinks,TestSolutionNameClosedSinks] = performSanityChecksonRecon(model,resultsFileName,ExtraCellCompIn,ExtraCellCompOut,runSingleGeneDeletion)
% This function performs various quality control and quality assurance
% tests.
% [TableChecks, Table_csources, CSourcesTestedRxns, TestSolutionNameOpenSinks,TestSolutionNameClosedSinks] = performSanityChecksonRecon(model,resultsFileName,ExtraCellCompIn,ExtraCellCompOut,runSingleGeneDeletion)
%
% INPUT
% model                         model structure
% resultsFileName               File name of the generated output file
% ExtraCellCompIn               [e] compartment by default, if extracellular
%                               uptake compartment is named differently, it
%                               can be specified here
% ExtraCellCompOut              [e] compartment by default, if extracellular
%                               secretion compartment is named differently, it
%                               can be specified here
% runSingleGeneDeletion         if 0 (default): function does not run single gene deletion otw choose 1
%
% OUTPUT
% TableChecks                   Table overview of the performed tests and
%                               their outcomes
% Table_csources                Table the test results for ATP yield from
%                               various carbon scources under aerobic and anaerobic conditions
% CSourcesTestedRxns            List of reactions active when testing the ATP
%                               yield from the various carbon sources

% TestSolutionNameOpenSinks     List of results when testing for 460 metabolic
%                               functions with all sinks open
% TestSolutionNameClosedSinks   List of results when testing for 460 metabolic
%                               functions with all sinks closed
%
% Ines Thiele 2016-2019

global resultsPath
resultsPath = which('MethodSection3.mlx');
resultsPath = strrep(resultsPath,'MethodSection3.mlx',['Results' filesep]);

if ~exist([resultsPath 'OrganChecks'],'dir')
    mkdir([resultsPath 'OrganChecks'])
end

if ~exist('ExtraCellCompIn','var')
    ExtraCellCompIn = '[e]'; % [e] compartment by default
end
if ~exist('ExtraCellCompOut','var')
    ExtraCellCompOut = '[e]'; % [e] compartment by default
end

if ~exist('runSingleGeneDeletion','var')
    runSingleGeneDeletion = 0; % do not run single gene deletion by default
end

global saveDiary
if saveDiary
    %save each diary to PSCM/Results/OrganChecks/
    global resultsPath
    resultsPath = which('MethodSection3.mlx');
    resultsPath = strrep(resultsPath,'MethodSection3.mlx',['Results' filesep]);
    resultsFileName=[resultsFileName '_diary'];
    diary([resultsPath 'OrganChecks' filesep resultsFileName])
end

% if 1
%     changeCobraSolver('tomlab_cplex','lp');
%     changeCobraSolver('tomlab_cplex','qp');
% end
cnt = 1;
tol = 1e-6;

model.rxns(find(ismember(model.rxns,'ATPM')))={'DM_atp_c_'};
model.rxns(find(ismember(model.rxns,'ATPhyd')))={'DM_atp_c_'};
% adds DM_atp to model if not exist

if isempty(strmatch('DM_atp_c_',model.rxns,'exact'))
    [model, rxnIDexists] = addReaction(model,'DM_atp_c_', 'reactionFormula', 'h2o[c] + atp[c]  -> adp[c] + h[c] + pi[c] ');
   % [model, rxnIDexists] = addReaction(model,'DM_atp_c_', 'reactionFormula', '1 atp[c] ->');
end

model.rxns(find(ismember(model.rxns,'EX_biomass_reaction')))={'biomass_reaction'};
model.rxns(find(ismember(model.rxns,'EX_biomass_maintenance')))={'biomass_maintenance'};
model.rxns(find(ismember(model.rxns,'EX_biomass_maintenance_noTrTr')))={'biomass_maintenance_noTrTr'};
model.lb(find(ismember(model.rxns,'biomass_reaction')))=0;
model.lb(find(ismember(model.rxns,'biomass_maintenance_noTrTr')))=0;
model.lb(find(ismember(model.rxns,'biomass_maintenance')))=0;

TestSolutionNameOpenSinks ='';
TestSolutionNameClosedSinks = '';

model.rxns = regexprep(model.rxns,'\(','\[');
model.rxns = regexprep(model.rxns,'\)','\]');
% vanilla leak test

if 1
    modelClosed = model;
    LeakTestRecon;
    TableChecks{cnt,1} = 'fastLeakTest 1';
    if length(LeakMets)>0
        warning('model leaks metabolites!')
        TableChecks{cnt,2} = 'Model leaks metabolites!';
    else
        TableChecks{cnt,2} = 'Leak free!';
    end
    cnt = cnt + 1;
end

if 0
    % test if something leaks when demand are added
    modelClosed = model;
    modelexchanges1 = strmatch('Ex_',modelClosed.rxns);
    modelexchanges4 = strmatch('EX_',modelClosed.rxns);
    modelexchanges2 = strmatch('DM_',modelClosed.rxns);
    modelexchanges3 = strmatch('sink_',modelClosed.rxns);
    selExc = (find( full((sum(abs(modelClosed.S)==1,1) ==1) & (sum(modelClosed.S~=0) == 1))))';

    modelexchanges = unique([modelexchanges1;modelexchanges2;modelexchanges3;modelexchanges4;selExc]);
    modelClosed.lb(find(ismember(modelClosed.rxns,modelClosed.rxns(modelexchanges))))=0;
    [modelClosed,rxnNames] = addDemandReaction(modelClosed,modelClosed.mets);

    [LeakRxns,modelTested] = fastLeakTest(modelClosed,modelClosed.rxns(selExc));

    TableChecks{cnt,1} = 'fastLeakTest 2 - add demand reactions for each metabolite in the model';
    if length(LeakMets)>0
        TableChecks{cnt,2} = 'Model leaks metabolites when demand reactions are added!';
    else
        TableChecks{cnt,2} = 'Leak free when demand reactions are added!';
    end
    cnt = cnt + 1;
end

if 1
    %%model produces energy from water!
    modelClosed = model;

    modelexchanges1 = strmatch('Ex_',modelClosed.rxns);
    modelexchanges4 = strmatch('EX_',modelClosed.rxns);
    modelexchanges2 = strmatch('DM_',modelClosed.rxns);
    modelexchanges3 = strmatch('sink_',modelClosed.rxns);
    selExc = (find( full((sum(abs(modelClosed.S)==1,1) ==1) & (sum(modelClosed.S~=0) == 1))))';

    modelexchanges = unique([modelexchanges1;modelexchanges2;modelexchanges3;modelexchanges4;selExc]);
    modelClosed.lb(find(ismember(modelClosed.rxns,modelClosed.rxns(modelexchanges))))=0;
    modelClosedATP = changeObjective(modelClosed,'DM_atp_c_');
    modelClosedATP = changeRxnBounds(modelClosedATP,'DM_atp_c_',0,'l');
    modelClosedATP = changeRxnBounds(modelClosedATP,'DM_atp_c_',1000,'u');
    modelClosedATP = changeRxnBounds(modelClosedATP,strcat('EX_h2o',ExtraCellCompIn),-1,'l');
    FBA3=optimizeCbModel(modelClosedATP);
    TableChecks{cnt,1} = 'Exchanges, sinks, and demands have  lb = 0, except h2o';
    if abs(FBA3.f) > 1e-6
        TableChecks{cnt,2} = 'model produces energy from water!';
    else
        TableChecks{cnt,2} = 'model DOES NOT produce energy from water!';
    end
    cnt = cnt + 1;
end
%% model produces energy from water and oxygen!
if 1
    modelClosed = model;
    modelexchanges1 = strmatch('Ex_',modelClosed.rxns);
    modelexchanges4 = strmatch('EX_',modelClosed.rxns);
    modelexchanges2 = strmatch('DM_',modelClosed.rxns);
    modelexchanges3 = strmatch('sink_',modelClosed.rxns);
    selExc = (find( full((sum(abs(modelClosed.S)==1,1) ==1) & (sum(modelClosed.S~=0) == 1))))';

    modelexchanges = unique([modelexchanges1;modelexchanges2;modelexchanges3;modelexchanges4;selExc]);
    modelClosed.lb(find(ismember(modelClosed.rxns,modelClosed.rxns(modelexchanges))))=0;
    modelClosedATP = changeObjective(modelClosed,'DM_atp_c_');
    modelClosedATP = changeRxnBounds(modelClosedATP,'DM_atp_c_',0,'l');
    modelClosedATP = changeRxnBounds(modelClosedATP,'DM_atp_c_',1000,'u');
    modelClosedATP = changeRxnBounds(modelClosedATP,strcat('EX_h2o',ExtraCellCompIn),-1,'l');
    modelClosedATP = changeRxnBounds(modelClosedATP,strcat('EX_o2',ExtraCellCompIn),-1,'l');

    FBA6=optimizeCbModel(modelClosedATP);
    TableChecks{cnt,1} = 'Exchanges, sinks, and demands have  lb = 0, except h2o and o2';
    if abs(FBA6.f) > 1e-6
        TableChecks{cnt,2} = 'model produces energy from water and oxygen!';
    else
        TableChecks{cnt,2} = 'model DOES NOT produce energy from water and oxygen!';
    end
    cnt = cnt + 1;
end
%% model produces matter when atp demand is reversed!
if 1
    modelClosed = model;
    modelexchanges1 = strmatch('Ex_',modelClosed.rxns);
    modelexchanges4 = strmatch('EX_',modelClosed.rxns);
    modelexchanges2 = strmatch('DM_',modelClosed.rxns);
    modelexchanges3 = strmatch('sink_',modelClosed.rxns);
    selExc = (find( full((sum(abs(modelClosed.S)==1,1) ==1) & (sum(modelClosed.S~=0) == 1))))';

    modelexchanges = unique([modelexchanges1;modelexchanges2;modelexchanges3;modelexchanges4;selExc]);
    modelClosed.lb(find(ismember(modelClosed.rxns,modelClosed.rxns(modelexchanges))))=0;

    modelClosed = changeObjective(modelClosed,'DM_atp_c_');
    modelClosed.lb(find(ismember(modelClosed.rxns,'DM_atp_c_'))) = -1000;
    modelClosed.ub(selExc)=1000;
    FBA = optimizeCbModel(modelClosed,'min');
    TableChecks{cnt,1} = 'Exchanges, sinks, and demands have  lb = 0, allow DM_atp_c_ to be reversible';
    if abs(FBA.f) > 1e-6
        TableChecks{cnt,2} = 'model produces matter when atp demand is reversed!';
    else
        TableChecks{cnt,2} = 'model DOES NOT produce matter when atp demand is reversed!';
    end
    cnt = cnt + 1;
end
%% model has flux through h[m] demand !
if 1
    modelClosed = model;
    modelexchanges1 = strmatch('Ex_',modelClosed.rxns);
    modelexchanges4 = strmatch('EX_',modelClosed.rxns);
    modelexchanges2 = strmatch('DM_',modelClosed.rxns);
    modelexchanges3 = strmatch('sink_',modelClosed.rxns);
    selExc = (find( full((sum(abs(modelClosed.S)==1,1) ==1) & (sum(modelClosed.S~=0) == 1))))';

    modelexchanges = unique([modelexchanges1;modelexchanges2;modelexchanges3;modelexchanges4;selExc]);
    modelClosed.lb(find(ismember(modelClosed.rxns,modelClosed.rxns(modelexchanges))))=0;
    modelClosed = addDemandReaction(modelClosed,'h[m]');
    modelClosed = changeObjective(modelClosed,'DM_h[m]');
    modelClosed.ub(find(ismember(modelClosed.rxns,'DM_h[m]'))) = 1000;
    modelClosed.ub(selExc)=1000;
    FBA = optimizeCbModel(modelClosed,'max');
    TableChecks{cnt,1} = 'Exchanges, sinks, and demands have  lb = 0, test flux through DM_h[m] (max)';
    if abs(FBA.f) > 1e-6
        TableChecks{cnt,2} = 'model has flux through h[m] demand (max)!';
    else
        TableChecks{cnt,2} = 'model has NO flux through h[m] demand (max)!';
    end
    cnt = cnt + 1;
end
if 0
    modelClosed = model;
    modelexchanges1 = strmatch('Ex_',modelClosed.rxns);
    modelexchanges4 = strmatch('EX_',modelClosed.rxns);
    modelexchanges2 = strmatch('DM_',modelClosed.rxns);
    modelexchanges3 = strmatch('sink_',modelClosed.rxns);
    selExc = (find( full((sum(abs(modelClosed.S)==1,1) ==1) & (sum(modelClosed.S~=0) == 1))))';

    modelexchanges = unique([modelexchanges1;modelexchanges2;modelexchanges3;modelexchanges4;selExc]);
    modelClosed.lb(find(ismember(modelClosed.rxns,modelClosed.rxns(modelexchanges))))=0;
    modelClosed = addDemandReaction(modelClosed,'h[m]');
    modelClosed = changeObjective(modelClosed,'DM_h[m]');
    modelClosed.ub(find(ismember(modelClosed.rxns,'DM_h[m]'))) = 1000;
    modelClosed.lb(find(ismember(modelClosed.rxns,'DM_h[m]'))) = -1000;
    modelClosed.ub(selExc)=1000;
    FBA = optimizeCbModel(modelClosed,'min');
    TableChecks{cnt,1} = 'Exchanges, sinks, and demands have  lb = 0, test flux through DM_h[m] (min)';
    if abs(FBA.f) > 1e-6
        TableChecks{cnt,2} = 'model has flux through h[m] demand (min)!';
    else
        TableChecks{cnt,2} = 'model has NO flux through h[m] demand (min)!';
    end
    cnt = cnt + 1;
end
%% model has flux through h[c] demand !
if 1
    modelClosed = model;
    modelexchanges1 = strmatch('Ex_',modelClosed.rxns);
    modelexchanges4 = strmatch('EX_',modelClosed.rxns);
    modelexchanges2 = strmatch('DM_',modelClosed.rxns);
    modelexchanges3 = strmatch('sink_',modelClosed.rxns);
    selExc = (find( full((sum(abs(modelClosed.S)==1,1) ==1) & (sum(modelClosed.S~=0) == 1))))';

    modelexchanges = unique([modelexchanges1;modelexchanges2;modelexchanges3;modelexchanges4;selExc]);
    modelClosed.lb(find(ismember(modelClosed.rxns,modelClosed.rxns(modelexchanges))))=0;
    modelClosed = addDemandReaction(modelClosed,'h[c]');
    modelClosed = changeObjective(modelClosed,'DM_h[c]');
    modelClosed.ub(find(ismember(modelClosed.rxns,'DM_h[c]'))) = 1000;
    modelClosed.ub(selExc)=1000;
    FBA = optimizeCbModel(modelClosed,'max');
    TableChecks{cnt,1} = 'Exchanges, sinks, and demands have  lb = 0, test flux through DM_h[c] (max)';
    if abs(FBA.f) > 1e-6
        TableChecks{cnt,2} = 'model has flux through h[c] demand (max)!';
    else
        TableChecks{cnt,2} = 'model has NO flux through h[c] demand (max)!';
    end
    cnt = cnt + 1;
end
if 1
    modelClosed = model;
    modelexchanges1 = strmatch('Ex_',modelClosed.rxns);
    modelexchanges4 = strmatch('EX_',modelClosed.rxns);
    modelexchanges2 = strmatch('DM_',modelClosed.rxns);
    modelexchanges3 = strmatch('sink_',modelClosed.rxns);
    selExc = (find( full((sum(abs(modelClosed.S)==1,1) ==1) & (sum(modelClosed.S~=0) == 1))))';

    modelexchanges = unique([modelexchanges1;modelexchanges2;modelexchanges3;modelexchanges4;selExc]);
    modelClosed.lb(find(ismember(modelClosed.rxns,modelClosed.rxns(modelexchanges))))=0;
    modelClosed = addDemandReaction(modelClosed,'h[c]');
    modelClosed = changeObjective(modelClosed,'DM_h[c]');
    modelClosed.ub(find(ismember(modelClosed.rxns,'DM_h[c]'))) = 1000;
    modelClosed.lb(find(ismember(modelClosed.rxns,'DM_h[c]'))) = -1000;
    modelClosed.ub(selExc)=1000;
    modelClosed.ub(find(ismember(modelClosed.rxns,strcat('EX_h',ExtraCellCompOut)))) = 0;
    FBA = optimizeCbModel(modelClosed,'min');
    TableChecks{cnt,1} = 'Exchanges, sinks, and demands have  lb = 0, ub of EX_h[e] = 0, test flux through DM_h[c] (min)';
    if abs(FBA.f) > 1e-6
        TableChecks{cnt,2} = 'model has flux through h[c] demand (min)!';
    else
        TableChecks{cnt,2} = 'model has NO flux through h[c] demand (min)!';
    end
    cnt = cnt + 1;
end
%% model produces too much atp demand from glc -- old remove
if 0
    modelClosed = model;
    modelexchanges1 = strmatch('Ex_',modelClosed.rxns);
    modelexchanges4 = strmatch('EX_',modelClosed.rxns);
    modelexchanges2 = strmatch('DM_',modelClosed.rxns);
    modelexchanges3 = strmatch('sink_',modelClosed.rxns);
    selExc = (find( full((sum(abs(modelClosed.S)==1,1) ==1) & (sum(modelClosed.S~=0) == 1))))';

    modelexchanges = unique([modelexchanges1;modelexchanges2;modelexchanges3;modelexchanges4;selExc]);
    modelClosed.lb(find(ismember(modelClosed.rxns,modelClosed.rxns(modelexchanges))))=0;
    modelClosed.c = zeros(length(modelClosed.rxns),1);
    modelClosed = changeObjective(modelClosed,'DM_atp_c_');
    modelClosed.lb(find(ismember(modelClosed.rxns,strcat('EX_o2',ExtraCellCompIn)))) = -1000;
    modelClosed.lb(find(ismember(modelClosed.rxns,strcat('EX_h2o',ExtraCellCompIn)))) = -1000;
    modelClosed.ub(find(ismember(modelClosed.rxns,strcat('EX_h2o',ExtraCellCompOut)))) = 1000;
    modelClosed.ub(find(ismember(modelClosed.rxns,strcat('EX_co2',ExtraCellCompOut)))) = 1000;
    %modelClosed.lb(find(ismember(modelClosed.rxns,'EX_pi[e]'))) = -1000;
    %modelClosed = addExchangeRxn(modelClosed,{'glc_D[e]'});
    modelClosed.lb(find(ismember(modelClosed.rxns,strcat('EX_glc_D',ExtraCellCompIn)))) = -1;
    modelClosed.ub(find(ismember(modelClosed.rxns,strcat('EX_glc_D',ExtraCellCompIn)))) = -1;

    modelClosed.ub(selExc)=1000;
    %    FBA = optimizeCbModel(modelClosed,'max',1e-6);
    FBA = optimizeCbModel(modelClosed,'max');

    TableChecks{cnt,1} = 'ATP yield ';
    if abs(FBA.f) > 1e-6
        FBA.f
        %  modelClosed = changeObjective(modelClosed,'ATPS4m');
        %  FBA = optimizeCbModel(modelClosed,'max',1e-6);
        FBA.x(find(ismember(modelClosed.rxns,'ATPS4m')))
        % prepare table

        warning('model produces too much atp demand from glc!');
    else
        warning('model DOES NOT produce too much atp demand from glc!');
    end

end

TableChecks{cnt,1} = 'Test metabolic objective functions with open sinks';
if 1 % perform test function
    [TestSolution,TestSolutionNameOpenSinks, TestedRxnsSinks, PercSinks] = Test4HumanFctExtv5(model,'all');
    TableChecks{cnt,2} = strcat('Done. See variable TestSolutionNameOpenSinks for results. The model passes ', num2str(length(find(abs(TestSolution)>tol))),' out of ', num2str(length(TestSolution)), 'tests');
else
    TableChecks{cnt,2} = 'Not performed.';
end
cnt =  cnt + 1;

TableChecks{cnt,1} = 'Test metabolic objective functions with closed sinks (lb)';
if 0 % perform test functions
    [TestSolution,TestSolutionNameClosedSinks, TestedRxnsClosedSinks, PercClosedSinks] = Test4HumanFctExtv5(model,'all',0);
    TableChecks{cnt,2} = strcat('Done. See variable TestSolutionNameClosedSinks for results. The model passes ', num2str(length(find(abs(TestSolution)>tol))),' out of ', num2str(length(TestSolution)), 'tests');
else
    TableChecks{cnt,2} = 'Not performed.';
end
cnt =  cnt + 1;

TableChecks{cnt,1} = 'Compute ATP yield';
if 1 % test ATP yield
    [Table_csources, CSourcesTestedRxns, Perc] = testATPYieldFromCsources(model,[],ExtraCellCompIn,ExtraCellCompOut);
    TableChecks{cnt,2} = 'Done. See variable Table_csources for results.';
else
    TableChecks{cnt,2} = 'Not performed.';
    CSourcesTestedRxns = '';
end
cnt = cnt + 1;

TableChecks{cnt,1} = 'Check duplicated reactions';
if 0
    method='FR';
    removeFlag=0;
    [modelOut,removedRxnInd, keptRxnInd] = checkDuplicateRxn(model,method,removeFlag);
    if isempty(removedRxnInd)
        TableChecks{cnt,2} = 'No duplicated reactions in model.';
    else
        TableChecks{cnt,2} = 'Duplicated reactions in model.';
    end
else
    TableChecks{cnt,2} = 'Not performed.';
end
cnt = cnt + 1;

TableChecks{cnt,1} = 'Check empty columns in rxnGeneMat';
if 1
    E = find(sum(model.rxnGeneMat)==0);
    if isempty(E)
        TableChecks{cnt,2} = 'No empty columns in rxnGeneMat.';
    else
        TableChecks{cnt,2} = 'Empty columns in rxnGeneMat.';
    end
else
    TableChecks{cnt,2} = 'Not performed.';
end
cnt = cnt + 1;

TableChecks{cnt,1} = 'Check that demand reactions have a lb >= 0';
if 1
    DMlb = find(model.lb(strmatch('DM_',model.rxns))<0);
    if isempty(DMlb)
        TableChecks{cnt,2} = 'No demand reaction can have flux in backward direction.';
    else
        TableChecks{cnt,2} = 'Demand reaction can have flux in backward direction.';
    end
else
    TableChecks{cnt,2} = 'Not performed.';
end
cnt = cnt + 1;

TableChecks{cnt,1} = 'Check consistency of model.rev with model.lb';
if 0
    Rev = setdiff(find(model.lb<0), find(model.rev==1));
    if isempty(Rev)
        TableChecks{cnt,2} = 'model.rev and model.lb are consistent.';
    else
        TableChecks{cnt,2} = 'model.rev and model.lb are NOT consistent.';
    end
else
    TableChecks{cnt,2} = 'Not performed.';
end
cnt = cnt + 1;

TableChecks{cnt,1} = 'Check whether singleGeneDeletion runs smoothly';
if runSingleGeneDeletion == 1
    try
        [grRatio,grRateKO,grRateWT,hasEffect,delRxns,fluxSolution] = singleGeneDeletion(model);
        TableChecks{cnt,2} = 'singleGeneDeletion finished without problems.';
    catch
        TableChecks{cnt,2} = 'There are problems with singleGeneDeletion.';
    end
else
    TableChecks{cnt,2} = 'Not performed.';
end
cnt = cnt + 1;

TableChecks{cnt,1} = 'Check for flux consistency';
if 0
    param.epsilon=getCobraSolverParams('LP', 'feasTol')*100;
    param.modeFlag=0;
    %param.method='null_fastcc';
    %param.method='fastcc';
    printLevel = 1;
    [fluxConsistentMetBool,fluxConsistentRxnBool,fluxInConsistentMetBool,fluxInConsistentRxnBool,model] = findFluxConsistentSubset(model,param,printLevel)
    if isempty(find(fluxInConsistentRxnBool))
        TableChecks{cnt,2} = 'Model is flux consistent.';
    else
        TableChecks{cnt,2} = 'Model is NOT flux consistent';
    end
else
    TableChecks{cnt,2} = 'Not performed.';
end
cnt = cnt + 1;

if 0
save([resultsPath 'OrganChecks' filesep resultsFileName,'.mat'],'TableChecks', 'Table_csources', 'CSourcesTestedRxns', 'TestSolutionNameOpenSinks','TestSolutionNameClosedSinks');
end
if saveDiary
    diary off;
end

author: opencobra

repo: cobratoolbox


text:
function dbName = response_demo_database

%RESPONSE_DEMO_DATABASE returns the absolute path to the demo database
%  DBNAME = RESPONSE_DEMO_DATABASE

file = which('response_cookbook');
path = fileparts(file);

if ~exist([path '/demo'])
	error('response_demo_database: demo database not found');
else
	dbName = [path '/demo/plutons'];
end

author: geoscience-community-codes

repo: GISMO


text:
function [A,B,flag] = hmxRSVD(varargin)
%+========================================================================+
%|                                                                        |
%|         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       : hmxRSVD.m                                     |
%|    #    |   VERSION    : 0.52                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 01.01.2019                                    |
%| ( === ) |   SYNOPSIS   : RSVD from Halko et al. 'finding structure with|
%|  `---'  |                randomness'. Adapted from Antoine Liutkus.    |
%+========================================================================+

% Input analysis
if (numel(varargin{2}) == 1)
    % Input
    M   = varargin{1};
    tol = varargin{2};
    rk  = varargin{3};

    % Dimensions
    [m,n] = size(M);

    % Matrix vectot product
    MV = @(V) M * V;
    VM = @(V) V * M;

    % Type
    typ = class(M);

else
    % Input
    A   = varargin{1};
    B   = varargin{2};
    tol = varargin{3};

    % Dimensions
    [m,rk] = size(A);
    n      = size(B,2);

    % Type
    typ = class(A);

    % Matrix vectot product
    MV    = @(V) A * (B * V);
    VM    = @(V) (V * A) * B;
end

% Randomized
p = min(2*rk,n);
X = randn(n,p,typ);
Y = MV(X);
try
    W1 = orth(Y);
catch
    A    = [];
    B    = [];
    flag = 0;
    return
end
B = VM(W1');

% Truncated SVD
try
    [W2,S,V] = svd(B,'econ');
catch
    A    = [];
    B    = [];
    flag = 0;
    return
end

% Product
U = W1*W2;

% Rank with fixed accuracy
if (numel(S) ~= 0)
    I = find(abs(diag(S)/S(1)) >= tol);
else
    I = [];
end

% No values
if isempty(I)
   A    = zeros(m,0);
   B    = zeros(0,n);
   flag = 1;

% Accuracy not reached
elseif (length(I) >= rk)
    A    = [];
    B    = [];
    flag = 0;

% Low-rank representation
else
    A    = U(:,I);
    B    = S(I,I) * V(:,I)';
    flag = 1;
end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text: %% ODFOverview.m

author: mtex-toolbox

repo: mtex

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


text:
%Help file for INTLAB Version 4.1
%
%Program demointlab added.
%
%Few minor changes and corrections.
%


text:
%% (Internal) Mean/Median filtering
%
%   Output = MedianFilt(Input, WinSize, bRobust)
%
% Arguments:
%
%      + Input: the signal
%
%      + WinSize: The size of the window in samples
%
%      + bRobust: use mean (false) or median (true)
%
% Output:
%
%      + Output: filtered output
%
% Example:
%
%     WinSize = round(0.2*SamplingFreq);
%
%     %Baseline estimation.
%     BaselineEstimation = MedianFilt(noisyECG, WinSize );
%
%
% See also BaselineWanderRemovalMedian
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 21/4/2015
% Copyright 2008-2015
%
function Output = MedianFilt(Input, WinSize, bRobust)

if(nargin < 2 || isempty(WinSize) )
    WinSize = 31;
end

if(nargin < 3 || isempty(bRobust) )
    bRobust = true;
end

if(bRobust)
    mean_ptr_func = @nanmedian;
else
    mean_ptr_func = @nanmean;
end

MidPoint = ceil(WinSize/2);

aux_seq = 1:size(Input,1);
laux_seq = length(aux_seq);

each_sample_idx = arrayfun(@(a)(max(1,a):min(laux_seq,a + WinSize-1)), aux_seq - MidPoint+1, ...
                                                'UniformOutput', false);
Output = cellfun(@(a)(mean_ptr_func(Input(a,:),1)), colvec(each_sample_idx), 'UniformOutput', false);

Output = cell2mat(Output);

%old version
% startSample = MidPoint;
% endSample = size(Input,1) - fix(WinSize/2);
% Output = zeros(endSample-startSample+1, size(Input,2));
%
% iCount = 1;
%
% for iSampleIndex = startSample:endSample
%
%     startRange = iSampleIndex-MidPoint+1;
%     iRange = startRange:startRange+WinSize-1;
%
%     Output(iCount,:) = median( Input(iRange,:) );
%
%     iCount = iCount + 1;
%
% end

author: marianux

repo: ecg-kit

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


text:
% This example is from Page.143 of "Probabilistic Networks and Expert Systems",
% Cowell, Dawid, Lauritzen and Spiegelhalter, 1999, Springer.

X = 1; Y = 2; Z = 3;
n = 3;

dag = zeros(n);
dag(X, Y)=1;
dag(Y, Z)=1;

ns = ones(1, n);
dnodes = [];

bnet = mk_bnet(dag, ns, dnodes);
bnet.CPD{X} = gaussian_CPD(bnet, X, 'mean', 0, 'cov', 1);
bnet.CPD{Y} = gaussian_CPD(bnet, Y, 'mean', 0, 'cov', 1, 'weights', 1);
bnet.CPD{Z} = gaussian_CPD(bnet, Z, 'mean', 0, 'cov', 1, 'weights', 1);

engines = {};
engines{end+1} = jtree_inf_engine(bnet);
engines{end+1} = stab_cond_gauss_inf_engine(bnet);
nengines = length(engines);

evidence = cell(1,n);
evidence{Y} = 1.5;

for e=1:nengines
  engines{e} = enter_evidence(engines{e}, evidence);
  margX = marginal_nodes(engines{e}, X);
  assert(approxeq(margX.mu, 0.75))
  assert(approxeq(margX.Sigma, 0.5))

  margZ = marginal_nodes(engines{e}, Z);
  assert(approxeq(margZ.mu, 1.5))
  assert(approxeq(margZ.Sigma, 1))
end

evidence = cell(1,n);
evidence{Z} = 1.5;

for e=1:nengines
  engines{e} = enter_evidence(engines{e}, evidence);
  margX = marginal_nodes(engines{e}, X);
  assert(approxeq(margX.mu, 1/2))
  assert(approxeq(margX.Sigma, 2/3))

  margY = marginal_nodes(engines{e}, Y);
  assert(approxeq(margY.mu, 1))
  assert(approxeq(margY.Sigma, 2/3))
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [nse] = read_neuralynx_nse(filename, begrecord, endrecord)

% READ_NEURALYNX_NSE reads a single electrode waveform file
%
% Use as
%   [nse] = read_neuralynx_nse(filename)
%   [nse] = read_neuralynx_nse(filename, begrecord, endrecord)

% Copyright (C) 2005-2007, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if nargin<2
  begrecord = 1;
end
if nargin<3
  endrecord = inf;
end

% the file starts with a 16*1024 bytes header in ascii, followed by a number of records
hdr = neuralynx_getheader(filename);
fid = fopen_or_error(filename, 'rb', 'ieee-le');

% determine the length of the file
fseek(fid, 0, 'eof');
headersize = 16384;
if isfield(hdr, 'RecordSize')
  recordsize = hdr.RecordSize;
else
  % this is the default, resulting in 32 samples per spike snippet
  recordsize = 112;
end
NSamples = (recordsize-48)/2;
NRecords = floor((ftell(fid) - headersize)/recordsize);

if NRecords>0
  if (ispc), fclose(fid); end
  % read the timestamp from the first and last record
  hdr.FirstTimeStamp = neuralynx_timestamp(filename, 1);
  hdr.LastTimeStamp  = neuralynx_timestamp(filename, inf);
  if (ispc), fid = fopen_or_error(filename, 'rb', 'ieee-le'); end
else
  hdr.FirstTimeStamp = nan;
  hdr.LastTimeStamp  = nan;
end

if begrecord==0 && endrecord==0
  % only read the header
elseif begrecord<1
  ft_error('cannot read before the first record');
elseif begrecord>NRecords
  ft_error('cannot read beyond the last record')
elseif endrecord>NRecords
  endrecord = NRecords;
end

if begrecord>=1 && endrecord>=begrecord
  % rewind to the first record to be read
  fseek(fid, headersize + (begrecord-1)*recordsize, 'bof');

  numrecord    = (endrecord-begrecord+1);
  TimeStamp    = zeros(1,numrecord,'uint64');
  ScNumber     = zeros(1,numrecord);
  CellNumber   = zeros(1,numrecord);
  Param        = zeros(8,numrecord);
  Samp         = zeros(1,NSamples,numrecord);

  k = 1;
  while (begrecord+k-1)<=endrecord
    % read a single electrode record
    TimeStamp(k)   = fread(fid,  1, 'uint64=>uint64');
    ScNumber(k)    = fread(fid,  1, 'int32');
    CellNumber(k)  = fread(fid,  1, 'int32');
    Param(:,k)     = fread(fid,  8, 'int32');
    Samp(1,:,k)    = fread(fid, NSamples, 'int16');
    k = k + 1;
  end

  nse.TimeStamp   = TimeStamp;
  nse.ScNumber    = ScNumber;
  nse.CellNumber  = CellNumber;
  nse.Param       = Param;
  % apply the scaling factor from ADBitVolts and convert to uV
  nse.dat         = Samp * hdr.ADBitVolts * 1e6;
end
fclose(fid);

nse.NRecords = NRecords;
nse.hdr      = hdr;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [dprime, prefCond] = mv_dprime(mv, conditions, verbose);
% [dprime, prefCond] = mv_dprime(mv, [conditions=all], [verbose=1]);
%
%
% Compute a 'd prime' form of selectivity index for a multi-voxel UI, given the selected
% conditions and threshold.
%
% This measure of d-prime performs a series of pairwise comparisons between
% each voxel's 'preferred' condition (the condition which elicits the highest
% response) to each of the other 'nonpreferred' conditions. For each
% comparison, it computes the 'd-prime', separation over spread, between
% those 2 conditions. This is defined (for conditions A and B) as:
%
%    d' = (meanA - meanB) / sqrt((varianceA + varianceB) / 2)
%
% Response amplitudes are computed according the to
% event-related paramter 'ampType': see er_setParams, er_defaultParams.
%
% In addition to coding the degree of response, this also returns a list of
% the preferred condition for each category.
%
% If the verbose flag is set to 1 (default), will report how many voxels prefer each
% condition.
%
% ras 09/06: written.
if notDefined('verbose'), verbose = prefsVerboseCheck; end
if notDefined('conditions'),
    conditions = mv.trials.condNums(mv.trials.condNums>0);
end

mv.params.selConds = conditions;
amps = mv_amps(mv);
amps = amps(:,conditions);
varAll = (mv_stdev(mv)) .^ 2; % variance = sigma^2
nVoxels = size(amps,1);
nConds = size(amps,2);

%%%%% step 1: find preferred cond for each condition
mx = max(amps,[],2); % max values
for i = 1:nVoxels
    if isnan(mx(i)) | sum( amps(i,:)==mx(i) ) > 1
        prefCond(i) = NaN;
        other(i,:) = NaN;
        varA(i,:) = NaN;
        varB(i,:) = NaN;
        continue;
    end

    % preferred condition
    prefCond(i) = find(amps(i,:)==mx(i));

    % amplitudes of other conditions
    other(i,:) = amps(i,setdiff(1:nConds,prefCond(i)));

    % get variance estimates for max, nonmax conds
    varA(i,1) = varAll(i,prefCond(i));
    varB(i,:) = varAll(i,setdiff(1:nConds,prefCond(i)));
end

% remove entries with NaNs
ok = ~isnan(prefCond);
prefCond = prefCond(ok);
other = other(ok,:);
varA = varA(ok,:);
varB = varB(ok,:);
mx = mx(ok,:);

%%%%% step 2: do pairwise comparisons
dprime = [];
for i = 1:size(other, 2)
    dprime(:,i) = [mx - other(:,i)] ./ sqrt([varA + varB(:,i)]./2);
end

% take mean across comparisons for each voxel
dprime = nanmean(dprime, 2);

% report breakdown by condition, if requested
if verbose==1
    for i=1:nConds
        fprintf(1,'cond %i numvoxels %i\n', i, sum(prefCond==i));
    end
end

return


text:
function prob_test0253 ( )

%*****************************************************************************80
%
%% TEST0253 tests BUFFON_PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST0253\n' );
  fprintf ( 1, '  BUFFON_PDF evaluates the Buffon PDF,\n' );
  fprintf ( 1, '  the probability that, on a grid of cells of width A,\n' );
  fprintf ( 1, '  a needle of length L, dropped at random,\n' );
  fprintf ( 1, '  will cross at least one grid line.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      A         L        PDF\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 5
    a = i;

    for k = 0 : 5
      l = k * a / 5.0;
      pdf = buffon_pdf ( a, l );
      fprintf ( 1, '  %8.4f  %8.4f  %14f\n', a, l, pdf );
    end

    fprintf ( 1, '\n' );

  end

  return
end


text:
function eeg_write_scan4_avg(f,filename)

% eeg_write_scan4_avg - Write a Neuroscan 'scan4.x' average file
%
% Useage: eeg_write_scan4_avg(f,filename)
%
% where:    filename is a complete 'path\fileprefix.ext' string;
%           the file will be a neuroscan 4.x avg file.
%
%           f is a structure containing:
%
%               f.header
%               f.electloc
%               f.data
%               f.variance
%               f.tag
%
%           See companion function: eeg_load_scan4_avg.m
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:54 $

% Licence:  GNU GPL, no implied or express warranties
% Created:  08/2000, Sean.Fitzgibbon_at_flinders.edu.au
% Modified: 08/2001, Darren.Weber_at_radiology.ucsf.edu
%                    distribute under GPL, with Sean's permission
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ver = '$Revision: 1.1 $';
fprintf('EEG_WRITE_SCAN4_AVG [v %s]\n',ver(11:15)); tic;

fid = fopen(filename,'w');

fwrite(fid,f.header.rev,'char');
fwrite(fid,f.header.nextfile,'long');
fwrite(fid,f.header.prevfile,'long');
fwrite(fid,f.header.type,'char');
fwrite(fid,f.header.id,'char');
fwrite(fid,f.header.oper,'char');
fwrite(fid,f.header.doctor,'char');
fwrite(fid,f.header.referral,'char');
fwrite(fid,f.header.hospital,'char');
fwrite(fid,f.header.patient,'char');
fwrite(fid,f.header.age,'short');
fwrite(fid,f.header.sex,'char');
fwrite(fid,f.header.hand,'char');
fwrite(fid,f.header.med,'char');
fwrite(fid,f.header.category,'char');
fwrite(fid,f.header.state,'char');
fwrite(fid,f.header.label,'char');
fwrite(fid,f.header.date,'char');
fwrite(fid,f.header.time,'char');
fwrite(fid,f.header.mean_age,'float');
fwrite(fid,f.header.stdev,'float');
fwrite(fid,f.header.n,'short');
fwrite(fid,f.header.compfile,'char');
fwrite(fid,f.header.spectwincomp,'float');
fwrite(fid,f.header.meanaccuracy,'float');
fwrite(fid,f.header.meanlatency,'float');
fwrite(fid,f.header.sortfile,'char');
fwrite(fid,f.header.numevents,'int');
fwrite(fid,f.header.compoper,'char');
fwrite(fid,f.header.avgmode,'char');
fwrite(fid,f.header.review,'char');
fwrite(fid,f.header.nsweeps,'ushort');
fwrite(fid,f.header.compsweeps,'ushort');
fwrite(fid,f.header.acceptcnt,'ushort');
fwrite(fid,f.header.rejectcnt,'ushort');
fwrite(fid,f.header.pnts,'ushort');
fwrite(fid,f.header.nchannels,'ushort');
fwrite(fid,f.header.avgupdate,'ushort');
fwrite(fid,f.header.domain,'char');
fwrite(fid,f.header.variance,'char');
fwrite(fid,f.header.rate,'ushort');
fwrite(fid,f.header.scale,'double');
fwrite(fid,f.header.veogcorrect,'char');
fwrite(fid,f.header.heogcorrect,'char');
fwrite(fid,f.header.aux1correct,'char');
fwrite(fid,f.header.aux2correct,'char');
fwrite(fid,f.header.veogtrig,'float');
fwrite(fid,f.header.heogtrig,'float');
fwrite(fid,f.header.aux1trig,'float');
fwrite(fid,f.header.aux2trig,'float');
fwrite(fid,f.header.heogchnl,'short');
fwrite(fid,f.header.veogchnl,'short');
fwrite(fid,f.header.aux1chnl,'short');
fwrite(fid,f.header.aux2chnl,'short');
fwrite(fid,f.header.veogdir,'char');
fwrite(fid,f.header.heogdir,'char');
fwrite(fid,f.header.aux1dir,'char');
fwrite(fid,f.header.aux2dir,'char');
fwrite(fid,f.header.veog_n,'short');
fwrite(fid,f.header.heog_n,'short');
fwrite(fid,f.header.aux1_n,'short');
fwrite(fid,f.header.aux2_n,'short');
fwrite(fid,f.header.veogmaxcnt,'short');
fwrite(fid,f.header.heogmaxcnt,'short');
fwrite(fid,f.header.aux1maxcnt,'short');
fwrite(fid,f.header.aux2maxcnt,'short');
fwrite(fid,f.header.veogmethod,'char');
fwrite(fid,f.header.heogmethod,'char');
fwrite(fid,f.header.aux1method,'char');
fwrite(fid,f.header.aux2method,'char');
fwrite(fid,f.header.ampsensitivity,'float');
fwrite(fid,f.header.lowpass,'char');
fwrite(fid,f.header.highpass,'char');
fwrite(fid,f.header.notch,'char');
fwrite(fid,f.header.autoclipadd,'char');
fwrite(fid,f.header.baseline,'char');
fwrite(fid,f.header.offstart,'float');
fwrite(fid,f.header.offstop,'float');
fwrite(fid,f.header.reject,'char');
fwrite(fid,f.header.rejstart,'float');
fwrite(fid,f.header.rejstop,'float');
fwrite(fid,f.header.rejmin,'float');
fwrite(fid,f.header.rejmax,'float');
fwrite(fid,f.header.trigtype,'char');
fwrite(fid,f.header.trigval,'float');
fwrite(fid,f.header.trigchnl,'char');
fwrite(fid,f.header.trigmask,'short');
fwrite(fid,f.header.trigisi,'float');
fwrite(fid,f.header.trigmin,'float');
fwrite(fid,f.header.trigmax,'float');
fwrite(fid,f.header.trigdir,'char');
fwrite(fid,f.header.autoscale,'char');
fwrite(fid,f.header.n2,'short');
fwrite(fid,f.header.dir,'char');
fwrite(fid,f.header.dispmin,'float');
fwrite(fid,f.header.dispmax,'float');
fwrite(fid,f.header.xmin,'float');
fwrite(fid,f.header.xmax,'float');
fwrite(fid,f.header.automin,'float');
fwrite(fid,f.header.automax,'float');
fwrite(fid,f.header.zmin,'float');
fwrite(fid,f.header.zmax,'float');
fwrite(fid,f.header.lowcut,'float');
fwrite(fid,f.header.highcut,'float');
fwrite(fid,f.header.common,'char');
fwrite(fid,f.header.savemode,'char');
fwrite(fid,f.header.manmode,'char');
fwrite(fid,f.header.ref,'char');
fwrite(fid,f.header.rectify,'char');
fwrite(fid,f.header.displayxmin,'float');
fwrite(fid,f.header.displayxmax,'float');
fwrite(fid,f.header.phase,'char');
fwrite(fid,f.header.screen,'char');
fwrite(fid,f.header.calmode,'short');
fwrite(fid,f.header.calmethod,'short');
fwrite(fid,f.header.calupdate,'short');
fwrite(fid,f.header.calbaseline,'short');
fwrite(fid,f.header.calsweeps,'short');
fwrite(fid,f.header.calattenuator,'float');
fwrite(fid,f.header.calpulsevolt,'float');
fwrite(fid,f.header.calpulsestart,'float');
fwrite(fid,f.header.calpulsestop,'float');
fwrite(fid,f.header.calfreq,'float');
fwrite(fid,f.header.taskfile,'char');
fwrite(fid,f.header.seqfile,'char');
fwrite(fid,f.header.spectmethod,'char');
fwrite(fid,f.header.spectscaling,'char');
fwrite(fid,f.header.spectwindow,'char');
fwrite(fid,f.header.spectwinlength,'float');
fwrite(fid,f.header.spectorder,'char');
fwrite(fid,f.header.notchfilter,'char');
fwrite(fid,f.header.headgain,'short');
fwrite(fid,f.header.additionalfiles,'int');
fwrite(fid,f.header.unused,'char');
fwrite(fid,f.header.fspstopmethod,'short');
fwrite(fid,f.header.fspstopmode,'short');
fwrite(fid,f.header.fspfvalue,'float');
fwrite(fid,f.header.fsppoint,'short');
fwrite(fid,f.header.fspblocksize,'short');
fwrite(fid,f.header.fspp1,'ushort');
fwrite(fid,f.header.fspp2,'ushort');
fwrite(fid,f.header.fspalpha,'float');
fwrite(fid,f.header.fspnoise,'float');
fwrite(fid,f.header.fspv1,'short');
fwrite(fid,f.header.montage,'char');
fwrite(fid,f.header.eventfile,'char');
fwrite(fid,f.header.fratio,'float');
fwrite(fid,f.header.minor_rev,'char');
fwrite(fid,f.header.eegupdate,'short');
fwrite(fid,f.header.compressed,'char');
fwrite(fid,f.header.xscale,'float');
fwrite(fid,f.header.yscale,'float');
fwrite(fid,f.header.xsize,'float');
fwrite(fid,f.header.ysize,'float');
fwrite(fid,f.header.acmode,'char');
fwrite(fid,f.header.commonchnl,'uchar');
fwrite(fid,f.header.xtics,'char');
fwrite(fid,f.header.xrange,'char');
fwrite(fid,f.header.ytics,'char');
fwrite(fid,f.header.yrange,'char');
fwrite(fid,f.header.xscalevalue,'float');
fwrite(fid,f.header.xscaleinterval,'float');
fwrite(fid,f.header.yscalevalue,'float');
fwrite(fid,f.header.yscaleinterval,'float');
fwrite(fid,f.header.scaletoolx1,'float');
fwrite(fid,f.header.scaletooly1,'float');
fwrite(fid,f.header.scaletoolx2,'float');
fwrite(fid,f.header.scaletooly2,'float');
fwrite(fid,f.header.port,'short');
fwrite(fid,f.header.numsamples,'long');
fwrite(fid,f.header.filterflag,'char');
fwrite(fid,f.header.lowcutoff,'float');
fwrite(fid,f.header.lowpoles,'short');
fwrite(fid,f.header.highcutoff,'float');
fwrite(fid,f.header.highpoles,'short');
fwrite(fid,f.header.filtertype,'char');
fwrite(fid,f.header.filterdomain,'char');
fwrite(fid,f.header.snrflag,'char');
fwrite(fid,f.header.coherenceflag,'char');
fwrite(fid,f.header.continuoustype,'char');
fwrite(fid,f.header.eventtablepos,'long');
fwrite(fid,f.header.continuousseconds,'float');
fwrite(fid,f.header.channeloffset,'long');
fwrite(fid,f.header.autocorrectflag,'char');
fwrite(fid,f.header.dcthreshold,'uchar');

for n = 1:f.header.nchannels
fwrite(fid,f.electloc(n).lab,'char');
fwrite(fid,f.electloc(n).reference,'char');
fwrite(fid,f.electloc(n).skip,'char');
fwrite(fid,f.electloc(n).reject,'char');
fwrite(fid,f.electloc(n).display,'char');
fwrite(fid,f.electloc(n).bad,'char');
fwrite(fid,f.electloc(n).n,'ushort');
fwrite(fid,f.electloc(n).avg_reference,'char');
fwrite(fid,f.electloc(n).clipadd,'char');
fwrite(fid,f.electloc(n).x_coord,'float');
fwrite(fid,f.electloc(n).y_coord,'float');
fwrite(fid,f.electloc(n).veog_wt,'float');
fwrite(fid,f.electloc(n).veog_std,'float');
fwrite(fid,f.electloc(n).snr,'float');
fwrite(fid,f.electloc(n).heog_wt,'float');
fwrite(fid,f.electloc(n).heog_std,'float');
fwrite(fid,f.electloc(n).baseline,'short');
fwrite(fid,f.electloc(n).filtered,'char');
fwrite(fid,f.electloc(n).fsp,'char');
fwrite(fid,f.electloc(n).aux1_wt,'float');
fwrite(fid,f.electloc(n).aux1_std,'float');
fwrite(fid,f.electloc(n).senstivity,'float');
fwrite(fid,f.electloc(n).gain,'char');
fwrite(fid,f.electloc(n).hipass,'char');
fwrite(fid,f.electloc(n).lopass,'char');
fwrite(fid,f.electloc(n).page,'uchar');
fwrite(fid,f.electloc(n).size,'uchar');
fwrite(fid,f.electloc(n).impedance,'uchar');
fwrite(fid,f.electloc(n).physicalchnl,'uchar');
fwrite(fid,f.electloc(n).rectify,'char');
fwrite(fid,f.electloc(n).calib,'float');
end

for i = 1:f.header.nchannels
   fwrite(fid,f.data(i).header,'char');
   fwrite(fid,f.data(i).samples,'float');
end

for j = 1:f.header.nchannels
   fwrite(fid,f.variance(j).samples,'float');
end

fwrite(fid,f.tag,'char');

frewind(fid);
fclose(fid);

t = toc; fprintf('...done (%6.2f sec).\n\n',t);

return

author: PatternRecognition

repo: OpenBMI


text:
%% Script for the Birthday Paradox simulation

% Define paramenters for simulation

  numtrials = 2e5;
  groupsize = 30;

% Preallocate some memory for the matches

  matches = zeros(1, numtrials);

  parfor trial = 1:numtrials

    matches(trial) = birthday(groupsize);

  end

% Probability is the sum of matches divided by number of trials

  prob = sum ( matches ) / numtrials;