text:
function [outputImage, selectedSlices, selectedVolumes] = ...
    select4D(this, varargin)
% Creates new image from selected data range, allows interactive picking of
% volumes/slices via clicking on montage
%
%   Y = MrImage()
%   Y.select4D(inputs)
%
% This is a method of class MrImage.
%
% IN
%
%   varargin    'ParameterName', 'ParameterValue'-pairs for the following
%               properties:
%               'method'            'manual' or 'indexed' (default)
%                                   if 'manual', montage plot is presented
%                                   to select4D slices/volumes
%                                   if 'indexed', extraction options (s.b.)
%                                   are used
%               'fixedWithinFigure' 'slices' or 'volumes' (default)
%                                   Defines which dimension is plotted
%                                   within the same figure when doing
%                                   manual selection
%               'combine'           'and'/'all' or 'or'/'once' (default)
%                                   Specifies whether combination over the
%                                   not-fixed dimension (e.g. slices, if
%                                   fixedWithinFigure = 'volumes') shall be
%                                   performed via an
%                                   and-operation
%                                   i.e. selection only if slice/volume
%                                   selected in all montages
%                                   or an or-operation
%                                   i.e. selection as soon as slice/volume
%                                   selected in at least one montage
%
%
%               Parameters for data extraction:
%
%               'signalPart'        for complex data, defines which signal
%                                   part shall be extracted for plotting
%                                       'all'       - do not change data (default)
%                                       'abs'       - absolute value
%                                       'phase'     - phase of signal
%                                       'real'      - real part of signal
%                                       'imag'      - imaginary part of
%                                                     signal
%               'plotMode',         transformation of data before plotting
%                                   'linear' (default), 'log'
%               'selectedX'         [1, nPixelX] vector of selected
%                                   pixel indices in 1st image dimension
%               'selectedY'         [1, nPixelY] vector of selected
%                                   pixel indices in 2nd image dimension
%               'selectedVolumes'   [1,nVols] vector of selected volumes to
%                                             be displayed
%               'selectedSlices'    [1,nSlices] vector of selected slices to
%                                               be displayed
%                                   choose Inf to display all volumes
%               'sliceDimension'    (default: 3) determines which dimension
%                                   shall be plotted as a slice
%               'exclude'           false (default) or true
%                                   if true, selection will be inverted, i.e.
%                                   selectedX/Y/Slices/Volumes will NOT be
%                                   extracted, but all others in dataset
%               'rotate90'          default: 0; 0,1,2,3; rotates image
%                                   by multiple of 90 degrees AFTER
%                                   flipping slice dimensions
% OUT
%
% EXAMPLE
%   select4D
%
%   See also MrImage MrImage.extract_plot4D_data

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-12-01
% Copyright (C) 2014 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 Licence (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/>.
%

switch nargin
    case 0
        defaults.method     = 'manual';
    otherwise
        defaults.method     = 'indexed';
end

defaults.fixedWithinFigure  = 'slices';
defaults.exclude            = 'false';

[argsUsed, argsExtract] = tapas_uniqc_propval(varargin, defaults);
tapas_uniqc_strip_fields(argsUsed);

if strcmpi(method, 'manual')
    % plot and click montage
    switch fixedWithinFigure

        %% Montage plot and selection per selected slice
        case {'slice', 'slices'}

            %% Take slice selection, if given
            if isfield(argsExtract, 'selectedSlices') && ...
                    ~isinf(argsExtract.selectedSlices)
                selectedSlices = argsExtract.selectedSlices;
            else
                nSlices = this.geometry.nVoxels(3);
                selectedSlices = 1:nSlices;
            end

            argsExtract.fixedWithinFigure = 'slice';

            %% Montage plot per slice
            for indSlice = selectedSlices
                argsExtract.selectedSlices = indSlice;

                this.plot(argsExtract);

                % Now pick volumes from montage via clicking

            end
        case {'volume', 'volumes'}
    end

else
    [dataSelected, displayRange] = this.extract_plot_data(argsExtract);
end

outputImage                     = this.copyobj('exclude', 'data');
outputImage.data                = dataSelected;
nVoxelsOriginal                 = size(dataSelected);
nVoxelsOriginal(end+1:4)        = 1;

% Update nVoxels,FOV; keep resolution
outputImage.update_geometry_dim_info('nVoxels', nVoxelsOriginal);

author: translationalneuromodeling

repo: tapas


text:
% Re-rank the top k candidates by searching a transformation when given camera intrinsics
% pred: the prediction 3D pose
% M_pool: the 2D pose library for matching
% GT_pool: the 3D pose library
% k: the number of extracted nearest candidates
% cameraParams: a camera object compatible with Matlab
function j_NN = kNN_cam_align(M_pool,GT_pool,pred,k,cameraParams)
    temp = zeros(1,28);
    temp_root = 0.5 * (pred(9,:) + pred(12,:));
    for a = 1:14
        temp(1,a*2-1) = pred(a,1) - temp_root(1);
        temp(1,a*2) = pred(a,2) - temp_root(2);
        y_c(a) = temp(1,2*a);
    end
    scale = max(y_c)-min(y_c);
    temp = temp / scale;
    m_idx = knnsearch(M_pool,temp,'k',k);
    e = zeros(1,k);
    for ii=1:k
        j_2d = GT_pool(m_idx(ii),:);
        j_2d = reshape(j_2d,3,14);
        j_2d = double(j_2d');

        %scale = (max(j_2d(:,2))-min(j_2d(:,2)))/(max(pred(:,2))-min(pred(:,2)));
        [r,~] = extrinsics(pred(:,1:2),j_2d,cameraParams);

        p2d = j_2d*r*cameraParams.IntrinsicMatrix';
        [~,Z] = procrustes(p2d,pred(:,1:2));
        su = j_2d(:,1:2)-Z(:,1:2);
        Sum = zeros(size(su,1),1);
        for i = 1:size(su,2)
            Sum = Sum+su(:,i).^2;
        end
        e(ii) = mean(sqrt(Sum));

        %e(ii) = MPJPE_procrus(pred(:,1:2),j_2d(:,1:2));
        %e(ii) = MPJPE_procrus([pred(:,1:2)*scale, j_2d(:,3)] ,j_2d);

    end

    [~,iid] = min(e);
    j_2d = GT_pool(m_idx(iid),:);
    j_2d = reshape(j_2d,3,14);
    j_NN = double(j_2d');
end

author: flyawaychase


text:
function SuperEllipsePlaneRhoTxi
mxMax = 0.99;
myMax = 0.99;
mxMin = 0;
myMin = 0;
c = SuperEllipseParamsRelator.c);
qV = [2 10^6];
color = {'b','r'};
outFile = '/home/alex/Dropbox/PaperStress/AddmissibleSpace/';

% for iq = 1:length(qV)
%     q = qV(iq);
%     rhoMin = 1 - c(q)*mxMax*myMax;
%     rho(:,1) = linspace(rhoMin,1, 100);
%     txiUB(:,1) = atan(mxMax^2*c(q)./(1-rho));
%     txiLB(:,1) = atan((1-rho)/(myMax^2*c(q)));
%     plot(rho,[txiLB txiUB],['-',color{iq}])
%     set(gca,'ytick',[0:pi/8:pi/2]) % where to set the tick marks
%     set(gca,'yticklabels',{'0','\pi/8','\pi/4','3\pi/8','\pi/2'})
%     xlabel('\rho')
% end

myCase = 'A';

switch myCase
    case 'A'

        %figure(1)
        for iq = 1:length(qV)
            q = qV(iq);
            txi = linspace(0,pi/2,1000);
            s.q = q;
            s.txi = txi;
            s.mxMin = mxMin;
            s.myMin = myMin;
            s.mxMax = mxMax;
            s.myMax = myMax;
            rhoBounds = SuperEllipseRhoBoundsComputer(s);
            [rhoMin,rhoMax] = rhoBounds.compute();

            figure()
            h1 = plot(txi,rhoMin,['-',color{iq}]);
            set(h1,'LineWidth',2);
            hold on
            h2 = plot(txi,rhoMax,['-',color{iq}]);
            set(h2,'LineWidth',2);
            set(gca,'xtick',[0:pi/8:pi/2]) % where to set the tick marks
            set(gca,'xticklabels',{'0','\pi/8','\pi/4','3\pi/8','\pi/2'})
            xlabel('\xi')
            ylabel('\rho')
            set(gca,'FontSize',18)
            axis([0 pi/2 0 1])
            print([outFile,'AdmissibleSpaceForQ',num2str(q)],'-dpdf')
        end
    case 'B'
        figure(1)
        hold on
        for iq = 1:length(qV)
            q = qV(iq);
            txi = linspace(0,pi/2,1000);
            rhoMinMx = 1- mxMax^2*c(q)./tan(txi);
            rhoMinMy = 1- myMax^2*c(q)*tan(txi);
            rhoMin = max(rhoMinMx,rhoMinMy);
            h1 = plot(txi,rhoMin,['-',color{iq}]);
            set(h1,'LineWidth',2);
            set(gca,'xtick',[0:pi/8:pi/2]) % where to set the tick marks
            set(gca,'xticklabels',{'0','\pi/8','\pi/4','3\pi/8','\pi/2'})
            xlabel('\xi')
            ylabel('\rho')
            set(gca,'FontSize',18)
            axis([0 pi/2 0 1])
            print([outFile,'AdmissibleSpace',num2str(q)],'-dpdf')
        end
end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function out=strrel(in1,in2,relop)
% function out=strrel(in1,in2,relop)
%
%  Compares strings according to strcmp and issorted according to relop
%  in1 and in2 can be character arrays or cell arrays
%  relop can be <, <=, >, >=, ==, and ~=, or their string equivalents
%    such as 'ne' or 'gt', etc.
%
%  Examples:
%  c = {'How','much','wood','would','a','woodchuck','chuck?'};
%  s = 'wood';
%  r = strrel(s,c,'<')
%  r = strrel(s,c,'ge')

% Author: Ben Barrowes, barrowes@alum.mit.edu, 10/2007
%
% 2007-10-18 bugfix and vectorization as pointed out by Urs (us) Schwarz

out=[];
% empty out if empty in (like []==3)
if isempty(in1) || isempty(in2), return, end

%doesn't handle numerical inputs
if isnumeric(in1) || isnumeric(in2)
 error('arguments to strrel must be char arrays or cell arrays of strings');
end

%promote char arrays to cell arrays
if ischar(in1), in1={in1}; end
if ischar(in2), in2={in2}; end

%repmat smaller cell array to size of larger if needed
if ~all(size(in1)==size(in2))
 if prod(size(in1))==1
  in1=repmat(in1,size(in2));
 elseif prod(size(in2))==1
  in2=repmat(in2,size(in1));
 else
  error('input array size mismatch in strrel');
 end
end

%should both be the same size cell arrays by now
out=false(size(in1));

% Do the logic
isEqual =strcmp(in1,in2);
for i=1:numel(in1)
 %isSorted=issorted({in1{i},in2{i}});
 isSorted=issortedcellchar({in1{i},in2{i}}); %about 3-4 faster
 switch lower(relop)
   case {'<','lt'}
     out(i)=isSorted & ~isEqual(i);
   case {'<=','le'}
     out(i)=isSorted | isEqual(i);
   case {'==','eq'}
     out(i)=isEqual(i);
   case {'~=','ne','!=','/='}
     out(i)=~isEqual(i);
   case {'>','gt'}
     out(i)=~isSorted & ~isEqual(i);
   case {'>=','ge'}
     out(i)=~isSorted | isEqual(i);
   otherwise
     error(['Relational operator, ',relop,', not recognized in strrel.'])
 end
end

%%%% this handle-vectorized version seemed to be somewhat slower
%%%isSorted=cellfun(@(x,y) issorted({x,y}),in1,in2);
%%%isEqual =strcmp(in1,in2);
%%%switch lower(relop)
%%%  case {'<','lt'}
%%%    out=isSorted & ~isEqual;
%%%  case {'<=','le'}
%%%    out=isSorted | isEqual;
%%%  case {'==','eq'}
%%%    out=isEqual;
%%%  case {'~=','ne','!=','/='}
%%%    out=~isEqual;
%%%  case {'>','gt'}
%%%    out=~isSorted & ~isEqual;
%%%  case {'>=','ge'}
%%%    out=~isSorted | isEqual;
%%%  otherwise
%%%    error(['Relational operator, ',relop,', not recognized in strrel.'])
%%%end

author: Sable

repo: mcbench-benchmarks


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

    properties (SetObservable)
        Flags
    end

    properties (SetAccess = private)
        AllowAutomaticModeEntry = false % Determines whether the mode can be entered automatically when the dataset or result is changed
        ExitOnNewPlugin = false % Mode will exit if a new plugin result is called
    end

    properties (Access = private)
        AppDef
        ViewerPanel
        GuiDataset
        Settings
        Reporting
    end

    methods
        function obj = MimMarkerMode(viewer_panel, gui_dataset, app_def, settings, reporting)
            obj.ViewerPanel = viewer_panel;
            obj.GuiDataset = gui_dataset;
            obj.AppDef = app_def;
            obj.Settings = settings;
            obj.Reporting = reporting;
        end

        function exit = ExitOnViewerPanelModeChanged(obj, new_mode)
            exit = ~strcmp(new_mode, MimModes.MarkerMode);
        end

        function EnterMode(obj, current_dataset, plugin_info, current_plugin_name, current_visible_plugin_name, current_context, current_segmentation_name)
            obj.ViewerPanel.SetModes(MimModes.MarkerMode, []);
        end

        function sub_mode = GetSubModeName(obj)
            sub_mode = obj.ViewerPanel.SubMode;
        end

        function ExitMode(obj)
            obj.ViewerPanel.SetControl('W/L');
        end

        function AutoSave(obj)
            obj.GuiDataset.AutoSaveMarkers();
        end

        function OverlayImageChanged(obj, ~, ~)
        end

        function MarkOverlayAsUnchanged(obj)
        end

        function SaveEdit(obj)
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
% SCRIPT TEST FOR THE KINEMATIC PROBLEM FOR SERIAL ROBOTS

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE 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.
%
% ARTE 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 ARTE.  If not, see <http://www.gnu.org/licenses/>.

close all

fprintf('\nTHE DEMO PRESENTS THE DIRECT AND INVERSE KINEMATIC PROBLEM')

%there are eight possible solutions for the inverse kinematic problem for most of these robots
n_solutions = 8;

%Try different configurations beware that, depending on the robot's topology
%not all the eight possible solutions will be feasible for an antropomorphic 6R robot.
q=[0 0 0 0 0 0]

%q=[0.2 0.8 -0.2 0.1 0.1 0.1]

%load robot parameters. You can try different robots
%robot=load_robot('KUKA', 'KR60_3'); n_solutions = 8;
robot=load_robot('MOTOMAN', 'MH6S'); n_solutions = 8;

%adjust 3D view as desired
adjust_view(robot)

%there are just 2 solutions for these robots and 4 DOF
%q = [pi/2 0.2 0.8 pi/4]
%robot=load_robot('kuka', 'KR5_scara_R350_Z200'); n_solutions = 2;
%robot=load_robot('example', 'scara'); n_solutions = 2;
%robot=load_robot('example', '2dofplanar'); n_solutions = 2;

%draw the robot
drawrobot3d(robot, q)

%Now compute direct kinematics for this position q
T = directkinematic(robot, q)

%Set to zero if you want to see the robot transparent
robot.graphical.draw_transparent=0;

%Set to one if you want to see the DH axes
%abb.graphical.draw_axes=1;

%Call the inversekinematic for this robot. All the possible solutions are
%stored at qinv. At least, one of the possible solutions should match q
qinv = inversekinematic(robot, T)

fprintf('\nNOW WE CAN REPRESENT THE DIFFERENT SOLUTIONS TO ACHIEVE THE SAME POSITION AND ORIENTATION\n')
fprintf('\nNote that some solutions may not be feasible since some joints may be out of range.\n')
correct=zeros(1,n_solutions);
%check that all of them are possible solutions!
for i=1:size(qinv,2),

    Ti = directkinematic(robot, qinv(:,i)) %Ti is constant for the different solutions

    % Note that all the solutions may not be feasible. Some of the joints may
    % be out of range. You can test this situation with test_joints
    test_joints(robot, qinv(:,i));

    %now draw the robot to see the solution
    drawrobot3d(robot, qinv(:,i))

    pause(1);

    k=sum(sum((T-Ti).^2));
    if k < 0.01 % a simple threshold to find differences in the solution
        correct(1,i)= 1;
    else
        correct(1,i)= 0; %uncorrect solution
        fprintf('\nERROR: One of the solutions seems to be uncorrect. Sum of errors: %f', i, k);
    end
end

fprintf('\n************** RESULTS **************')

%Display a message if any of the solutions is not correct
if sum(correct)==n_solutions
    fprintf('\nTEST 1--> OK: Every solution in qinv yields the same position/orientation T');
else
    fprintf('\nTEST 1--> ERROR: One or more of the solutions seem to be uncorrect.');
end

%Now, test if any of the solutions in qinv matches q
%find the solution that matches the initial q
%delta is just a squared sum of errors at each of the columns of the matrix
%which store the different solutions of qinv
delta=(repmat(q',[1 n_solutions])-qinv).^2;
i=find(sum(delta,1) < 0.01);
if ~isempty(i)
    fprintf('\nTEST 2--> OK!: Found a matching solution for the initial q.\n');
    solution=qinv(:,i)
else
    error_test2=1
    fprintf('\nTEST 2--> ERROR: Did not find a matching solution for the initial q.');
end

fprintf('\n************** ****** **************\n')

repo: ARTE


text:
function x = p13_root ( i )

%*****************************************************************************80
%
%% P13_ROOT returns a root for problem 13.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the index of the requested root.
%
%    Output, real X, the value of the root.
%
  if ( i == 1 )
    x = 100.0;
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P13_ROOT - Fatal error!\n' );
    fprintf ( 1, '  Illegal root index = %d\n', i );
    error ( 'P13_ROOT - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function printFluxes(model, fluxes, onlyExchange, cutOffFlux, outputFile,outputString,metaboliteList)
% printFluxes
%   Prints reactions and fluxes to the screen or to a file
%
%   model           a model structure
%   fluxes          a vector with fluxes
%   onlyExchange    only print exchange fluxes (opt, default true)
%   cutOffFlux      only print fluxes with absolute values above or equal to this
%                   value (opt, default 10^-8)
%   outputFile      a file to save the print-out to (opt, default is output to
%                   the command window)
%   outputString    a string that specifies the output of each reaction (opt,
%                   default '%rxnID\t(%rxnName):\t%flux\n')
%   metaboliteList  cell array of metabolite names. Only reactions
%                   involving any of these metabolites will be
%                   printed (opt)
%
%   The following codes are available for user-defined output strings:
%
%   %rxnID      reaction ID
%   %rxnName    reaction name
%   %lower      lower bound
%   %upper      upper bound
%   %obj        objective coefficient
%   %eqn        equation
%   %flux       flux
%   %element    equation using the metabolite formulas rather than
%               metabolite names
%   %unbalanced "(*)" if the reaction is unbalanced and "(-)" if it could not
%               be parsed
%   %lumped     equation where the elemental compositions for the left/right
%               hand sides are lumped
%
%   Usage: printFluxes(model, fluxes, onlyExchange, cutOffFlux,
%           outputFile,outputString,metaboliteList)

if nargin<3
    onlyExchange=true;
end
if nargin<4
    cutOffFlux=10^-8;
end
if isempty(cutOffFlux)
    cutOffFlux=10^-8;
end
if nargin<5
    fid=1;
else
    if ~isempty(outputFile)
        outputFile=char(outputFile);
        fid=fopen(outputFile,'w');
    else
        fid=1;
    end
end
if nargin<6 || isempty(outputString)
    outputString='%rxnID\t(%rxnName):\t%flux\n';
else
    outputString=char(outputString);
end
if nargin<7
    metaboliteList={};
else
    metaboliteList=convertCharArray(metaboliteList);
end
if size(fluxes,1)~=numel(model.rxns)
    EM='The number of fluxes and the number of reactions must be the same';
    dispEM(EM);
end

%Only keep reactions involving the defined metabolites
if ~isempty(metaboliteList)
    I=ismember(upper(model.metNames),upper(metaboliteList));
    [~, K]=find(model.S(I,:));

    %Delete all other reactions
    toDelete=true(numel(model.rxns),1);
    toDelete(K)=false;
    model=removeReactions(model,toDelete);
    fluxes(toDelete,:)=[];
end

if onlyExchange==true
    fprintf(fid,'EXCHANGE FLUXES:\n');
else
    fprintf(fid,'FLUXES:\n');
end

%Remove reactions which are below the cut off
toDelete=abs(fluxes)<cutOffFlux;
toDelete=all(toDelete,2);
model=removeReactions(model,toDelete,true,true);
fluxes(toDelete,:)=[];

if any(strfind(outputString,'%eqn'))
    %Construct the equations
    eqn=constructEquations(model);
else
    eqn=cell(numel(model.rxns),1);
    eqn(:)={''};
end
if any(strfind(outputString,'%element'))
    %For printing equations using the composition
    cModel=model;
    cModel.metNames=cModel.metFormulas;
    cModel.metNames(cellfun(@isempty,cModel.metNames))={'?'};
    element=constructEquations(cModel);
else
    element=cell(numel(model.rxns),1);
    element(:)={''};
end

if any(strfind(outputString,'%unbalanced')) || any(strfind(outputString,'%lumped'))
    balanceStructure=getElementalBalance(model);
end

unbalanced=cell(numel(model.rxns),1);
unbalanced(:)={''};
if any(strfind(outputString,'%unbalanced'))
    unbalanced(balanceStructure.balanceStatus==0)={'(*)'};
    unbalanced(balanceStructure.balanceStatus<0)={'(-)'};
end

lumped=cell(numel(model.rxns),1);
lumped(:)={''};
if any(strfind(outputString,'%lumped'))
    for i=1:numel(model.rxns)
        leftGroup='';
        rightGroup='';
        for j=1:numel(balanceStructure.elements.names)
            I=balanceStructure.leftComp(i,j);
            if I~=0
                if I==1
                    leftGroup=[leftGroup balanceStructure.elements.abbrevs{j}];
                else
                    leftGroup=[leftGroup balanceStructure.elements.abbrevs{j} num2str(I)];
                end
            end
            I=balanceStructure.rightComp(i,j);
            if I~=0
                if I==1
                    rightGroup=[rightGroup balanceStructure.elements.abbrevs{j}];
                else
                    rightGroup=[rightGroup balanceStructure.elements.abbrevs{j} num2str(I)];
                end
            end
        end
        if model.rev(i)
            lumped{i}=[leftGroup ' <=> ' rightGroup];
        else
            lumped{i}=[leftGroup ' => ' rightGroup];
        end
    end
end

for i=1:numel(model.rxns)
    %Only print if it's an exchange reaction or if all reactions should be
    %printed. Exchange reactions only have reactants or only products.
    reactants=model.S(:,i)<0;
    products=model.S(:,i)>0;

    %Only print if the absolute value is >= cutOffFlux
    if (onlyExchange==false || (~any(reactants) || ~any(products)))
        printString=outputString;

        %Produce the final string
        printString=strrep(printString,'%rxnID',model.rxns{i});
        printString=strrep(printString,'%eqn',eqn{i});
        printString=strrep(printString,'%rxnName',model.rxnNames{i});
        printString=strrep(printString,'%lower',num2str(model.lb(i)));
        printString=strrep(printString,'%upper',num2str(model.ub(i)));
        printString=strrep(printString,'%obj',num2str(model.c(i)));
        printString=strrep(printString,'%flux',num2str(fluxes(i,:)));
        printString=strrep(printString,'%element',element{i});
        printString=strrep(printString,'%unbalanced',unbalanced{i});
        printString=strrep(printString,'%lumped',lumped{i});
        fprintf(fid,printString);
    end
end

if fid~=1
    fprintf('File successfully saved.\n');
    fclose(fid);
end
end

author: SysBioChalmers

repo: RAVEN


text:
% Van Leer scheme for one-dimensional Euler equations

% Copyright 2001 P. Wesseling
% This program and its subprograms may be freely used, modified and distributed
% under the GNU General Public License: http://www.gnu.org/copyleft/gpl.html

% Theory in Section 10.5 of:

% 	P. Wesseling: Principles of Computational Fluid Dynamics
% 	Springer, Heidelberg, 2000 ISBN 3-5453-0. XII, 642 pp.
% See http://ta.twi.tudelft.nl/nw/users/wesseling/cfdbook.html

% This program makes Figs. 10.15, 10.16 in the book

% Functions called: f, problem_specification, Riemann

global  PRL  CRL MACHLEFT  gamma  pleft  pright  rholeft  rhoright  uleft...
	uright  tend  lambda		% lambda = dt/dx

		% .....................Input............................
gamma = 1.4; 	% Ratio of specific heats
J = 48;		% Number of grid cells
bouncon = 0;	% bouncon chooses outflow boundary conditions
		% = 0: Nothing happens: infinite domain
		% = 1: Solid wall at x = 1 with direct prescription of uwall
		% = 2: Solid wall at x = 1 with reflection b.c.
		% ....................End of input........................

gammab = 1/(gamma - 1); gam1 = gamma-1; gamgam = gamma^gamma;
problem_specification

h = 1/J;  				% Cell size
dt = lambda*h;				% Time step
n = floor(tend/dt);			% Number of time-steps

% 		Definition of grid numbering
%       x=0    					 x=1
% grid   |---o---|---o---|---o---  ...  --|---o---|
%            1   1   2   2   3           J-1  J

xcenter = h*[1:J] - h/2;		% Location of cell centers

press = zeros(size(xcenter));		% Preallocation of pressure,
rhoold = press; uold = press;		%       density and velocity
rhonew = press; mnew = press;		%       momentum
totenew = press;			%	total energy

for j = 1:length(xcenter)		% Initial conditions
  if xcenter(j) < 0.5, press(j) = pleft; rhoold(j) = rholeft;  uold(j) = uleft;
  else,  	     press(j) = pright;  rhoold(j) = rhoright; uold(j) = uright;
  end
end

	% Initialization of cell center variables
totenold = rhoold.*(0.5*uold.*uold + gammab*press./rhoold); % Total energy rho*E
totenleft = totenold(1); totenright = totenold(J);
mold = rhoold.*uold;					    % Momentum m
c = sqrt(gamma*press./rhoold);				    % Sound speed
mach = uold./c;						    % Mach number

% Preallocation of split fluxes
plus1 = c; minus1 = c; plus2 = c; minus2 = c; plus3 = c; minus3 = c;
% Preallocation of van Leer fluxes
flux1 = zeros(J-1,1); flux2 = flux1; flux3 = flux1;

t = 0;
for i = 1:n,  t = t + dt;
  Eflux1 = rhoold.*uold;			% Eflux1,2,3 is Euler flux
  Eflux2 = Eflux1.*uold + (1/gamma)*rhoold.*c.^2;
  Eflux3 = 0.5*Eflux1.*uold.^2 + gammab*Eflux1.*c.^2;

    for j = 1:J
    if mach(j) > 1
      plus1(j) = Eflux1(j);  plus2(j) = Eflux2(j);  plus3(j) = Eflux3(j);
    elseif mach(j) < -1
      plus1(j) = 0;         plus2(j) = 0;         plus3(j) = 0;
    else
      plus1(j) = 0.25*rhoold(j)*c(j)*(1 + mach(j))^2;
      plus2(j) = plus1(j)*c(j)*(2 + gam1*mach(j))/gamma;
      plus3(j) = (plus2(j)^2/plus1(j))*gamma^2*0.5*gammab/(gamma+1);
    end
  end
  minus1 = Eflux1 - plus1; minus2 = Eflux2 - plus2; minus3 = Eflux3 - plus3;

  for j = 1:J-1					% van Leer fluxes
    flux1(j) = plus1(j) + minus1(j+1);
    flux2(j) = plus2(j) + minus2(j+1);
    flux3(j) = plus3(j) + minus3(j+1);
  end

	% Update of state variables
  rhonew(1) = rholeft;		  rhonew(J)  = rhoright;
  mnew(1) = rholeft*uleft; 	  mnew(J)    = rhoright*uright;
  totenew(1) = totenleft; 	  totenew(J) = totenright;
  for j = 2:J-1
    rhonew(j)  = rhoold(j)   - lambda*(flux1(j) - flux1(j-1));
    mnew(j)    = mold(j)     - lambda*(flux2(j) - flux2(j-1));
    totenew(j) = totenold(j) - lambda*(flux3(j) - flux3(j-1));
  end

  if bouncon ~= 0
    if bouncon == 1
      rhowall = rhoold(J); uwall = 0;
      pwall = (gamma-1)*(totenold(J) - 0.5*mold(J)^2/rhoold(J)); % pwall = p(J)
      totenwall = pwall/(gamma-1) + 0.5*rhowall*uwall^2;
    else
      rhowall = rhoold(J); uwall = - mold(J)/rhoold(J);
      pwall = (gamma-1)*(totenold(J) - 0.5*mold(J)^2/rhoold(J)); % pwall = p(J)
      totenwall = pwall/(gamma-1) + 0.5*rhowall*uwall^2;
    end
    wallflux1 = rhowall*uwall;
    pwall = gam1*(totenwall - 0.5*wallflux1*uwall);
    wallflux2 = wallflux1.*uwall + pwall;
    wallflux3 = 0.5*wallflux1*uwall^2 + gammab*gamma*uwall*pwall;
    cwall = sqrt(gamma*pwall/rhowall); machwall = uwall/cwall;
    if machwall > 1
      plus1wall = wallflux1; plus2wall = wallflux2; plus3wall = wallflux3 ;
    elseif machwall < -1
      plus1wall = 0;      plus2wall = 0;      plus3wall = 0;
    else
      plus1wall = 0.25*rhowall*cwall*(1 + machwall)^2;
      plus2wall = plus1wall*cwall*(2 + gam1*machwall)/gamma;
      plus3wall = (plus2wall^2/plus1wall)*gamma^2*0.5*gammab/(gamma+1);
    end
    minus1wall = wallflux1 - plus1wall;
    minus2wall = wallflux2 - plus2wall; minus3wall = wallflux3 - plus3wall;

    flux1wall = plus1(J) + minus1wall;			% van Leer fluxes
    flux2wall = plus2(J) + minus2wall;
    flux3wall = plus3(J) + minus3wall;

    rhonew(J)  = rhoold(J)   - lambda*(flux1wall - flux1(J-1));
    mnew(J)    = mold(J)     - lambda*(flux2wall - flux2(J-1));
    totenew(J) = totenold(J) - lambda*(flux3wall - flux3(J-1));
  end
 	% Update of state variables
  uold = mnew./rhonew; press = gam1*(totenew - 0.5*mnew.*uold);
  rhoold = rhonew; totenold = totenew; mold = mnew;

  c = sqrt(gamma*press./rhoold); mach = uold./c;
end

entropy = log(press./rhoold.^gamma);

figure(1), clf
subplot(2,3,1),hold on,title('DENSITY','fontsize',14),plot(xcenter,rhonew,'o')
subplot(2,3,2),hold on,title('VELOCITY','fontsize',14),plot(xcenter,uold,'o')
subplot(2,3,3),hold on,title('PRESSURE','fontsize',14),plot(xcenter,press,'o')
subplot(2,3,4),hold on,title('MACHNUMBER','fontsize',14),plot(xcenter,mach,'o')
subplot(2,3,5),hold on,title('ENTROPY','fontsize',14),plot(xcenter,entropy,'o')
subplot(2,3,6),axis('off'),hold on,title('Van Leer scheme','fontsize',14)

Riemann		% Plot exact solution

repo: Aero-matlab


text:
function [diffx, diff_hor2]=horizondist(llhref, llhpos)
Ce2n0=llh2dcm_v000(llhref(1:2),[0;1]);
xyz_pos=ecef2geo_v000(llhpos,1);
xyz_ref=ecef2geo_v000(llhref,1);
diffx=Ce2n0*(xyz_pos-xyz_ref);
diff_hor2=sqrt(sum(diffx(1:2).^2));
end

repo: ekfmonoslam


text:
function varargout = drawSquareMesh(nodes, edges, faces, varargin) %#ok<INUSL>
%DRAWSQUAREMESH Draw a 3D square mesh given as a graph
%
%   drawSquareMesh(NODES, EDGES, FACES)
%   Draw the mesh defined by NODES, EDGES and FACES. FACES must be a N-by-4
%   array of vertex indices.
%
%   See Also
%   boundaryGraph, drawGraph
%
%   ---------
%   author: David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 28/06/2004.
%

% input size check up
if size(faces, 2) ~= 4
    error('Requires a face array with 4 columns');
end

% number of faces
Nf = size(faces, 1);

% allocate memory for vertex coordinates
px = zeros(4, Nf);
py = zeros(4, Nf);
pz = zeros(4, Nf);

% initialize vertex coordinates of each face
for f = 1:Nf
    face = faces(f, 1:4);
    px(1:4, f) = nodes(face, 1);
    py(1:4, f) = nodes(face, 2);
    pz(1:4, f) = nodes(face, 3);
end

p = patch(px, py, pz, 'r');

if nargout > 0
    varargout = {p};
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
function [cm] = au2cm(au)
% Convert length from astronomical units to centimeters.
% Chad A. Greene 2012
cm = au*1.49597870691e+13;


text:
function [M] = spm_nwcov (M)
% Get second moments of Normal-Wishart
% FORMAT [M] = spm_nwcov (M)
%
% .mean_prior_cov    Prior covariance of mean
% .sample_prior_cov  Prior covariance of samples
% .mean_post_cov     Posterior covariance of mean
% .sample_pred_cov   Predictive covariance of samples
%
% The latter quantity is also the covariance of the predictive density
% The marginal distributions of the mean and of the samples
% are multivariate-T, not Gaussian.
%
% See J. Bernardo and A. Smith (2000)
% Bayesian Theory, Wiley (page 435)
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny
% $Id: spm_nwcov.m 6548 2015-09-11 12:39:47Z will $

% Get prior covariances
M.mean_prior_cov=M.B0/(M.n0*(M.a0-1));

alpha=2*M.a0-M.P+1;
w_s=(1+1/M.n0)/(0.5*alpha-1);
M.sample_prior_cov=w_s*M.B0;

% Get posterior covariances
M.mean_post_cov=M.BN/(M.nN*(M.aN-1));

alpha=2*M.aN-M.P+1;
w_s=(1+1/M.nN)/(0.5*alpha-1);
M.sample_pred_cov=w_s*M.BN;

author: spm


text:
function updateScatterPolar(obj, plotIndex)

    %-AXIS INDEX-%
    axIndex = obj.getAxisIndex(obj.State.Plot(plotIndex).AssociatedAxis);

    %-PLOT DATA STRUCTURE- %
    plotData = get(obj.State.Plot(plotIndex).Handle);

    %-CHECK FOR MULTIPLE AXES-%
    [xsource, ysource] = findSourceAxis(obj, axIndex);

    %-ASSOCIATE POLAR-AXES LAYOUT-%
    obj.data{plotIndex}.subplot = sprintf('polar%d', xsource+1);

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

    %-parse plot data-%
    rData = plotData.RData;
    thetaData = rad2deg(plotData.ThetaData);

    thetaData(rData<0) = mod(thetaData(rData<0)+180, 360);
    rData = abs(rData);

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

    %-scatterpolar trace setting-%
    obj.data{plotIndex}.type = 'scatterpolar';
    obj.data{plotIndex}.mode = 'markers';
    obj.data{plotIndex}.visible = strcmp(plotData.Visible,'on');
    obj.data{plotIndex}.name = plotData.DisplayName;

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

    %-set scatterpolar data-%
    obj.data{plotIndex}.r = rData;
    obj.data{plotIndex}.theta = thetaData;

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

    %-trace settings-%
    markerStruct = extractScatterMarker(plotData);

    obj.data{plotIndex}.marker = markerStruct;

    if length(markerStruct.size) == 1
        obj.data{plotIndex}.marker.size = markerStruct.size * 0.2;
    end

    if length(markerStruct.line.color) > 1
        obj.data{plotIndex}.marker.line.color = markerStruct.line.color{1};
    end

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

    %-legend setting-%
    leg = get(plotData.Annotation);
    legInfo = get(leg.LegendInformation);

    switch legInfo.IconDisplayStyle
        case 'on'
            obj.data{plotIndex}.showlegend = true;
        case 'off'
            obj.data{plotIndex}.showlegend = false;
    end

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

    %-set polar axes-%
    updatePolaraxes(obj, plotIndex);

    %-------------------------------------------------------------------------%
end

%-------------------------------------------------------------------------%
%
%-SET POLAR AXIS-%
%
%-------------------------------------------------------------------------%

function updatePolaraxes(obj, plotIndex)

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

    %-AXIS INDEX-%
    axIndex = obj.getAxisIndex(obj.State.Plot(plotIndex).AssociatedAxis);

    %-CHECK FOR MULTIPLE AXES-%
    [xsource, ysource] = findSourceAxis(obj, axIndex);

    %-GET DATA STRUCTURES-%
    plotData = get(obj.State.Plot(plotIndex).Handle);
    axisData = get(plotData.Parent);
    thetaAxis = get(axisData.ThetaAxis);
    rAxis = get(axisData.RAxis);

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

    %-set domain plot-%
    xo = axisData.Position(1);
    yo = axisData.Position(2);
    w = axisData.Position(3);
    h = axisData.Position(4);

    polarAxis.domain.x = min([xo xo + w], 1);
    polarAxis.domain.y = min([yo yo + h], 1);

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

    %-setting angular axis-%
    gridColor = sprintf('rgba(%f,%f,%f,%f)', 255*axisData.GridColor, ...
        axisData.GridAlpha);
    gridWidth = axisData.LineWidth;
    thetaLim = thetaAxis.Limits;

    polarAxis.angularaxis.linecolor = gridColor;
    polarAxis.angularaxis.ticklen = mean(thetaAxis.TickLength);

    if isnumeric(thetaLim)
        polarAxis.angularaxis.range = thetaLim;
    else
        polarAxis.angularaxis.autorange = true;
    end

    if strcmp(axisData.ThetaGrid, 'on')
        polarAxis.angularaxis.gridwidth = gridWidth;
        polarAxis.angularaxis.gridcolor = gridColor;
    end

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

    %-set angular axis label-%
    thetaLabel = thetaAxis.Label;

    polarAxis.angularaxis.title.text = thetaLabel.String;
    polarAxis.radialaxis.title.font.family = matlab2plotlyfont(...
        thetaLabel.FontName);
    polarAxis.radialaxis.title.font.size = thetaLabel.FontSize;
    polarAxis.radialaxis.title.font.color = sprintf('rgb(%f,%f,%f)', ...
        255*thetaLabel.Color);

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

    %-setting radial axis-%
    rLim = rAxis.Limits;

    polarAxis.radialaxis.showline = false;
    polarAxis.radialaxis.angle = axisData.RAxisLocation+6;
    polarAxis.radialaxis.tickangle = 90-rAxis.TickLabelRotation;
    polarAxis.radialaxis.ticklen = mean(rAxis.TickLength);

    if isnumeric(rLim)
        polarAxis.radialaxis.range = rLim;
    else
        polarAxis.radialaxis.autorange = true;
    end

    if strcmp(axisData.RGrid, 'on')
        polarAxis.radialaxis.gridwidth = gridWidth;
        polarAxis.radialaxis.gridcolor = gridColor;
    end

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

    %-set radial axis label-%
    rLabel = thetaAxis.Label;

    polarAxis.angularaxis.title.text = 'label';%rLabel.String;
    polarAxis.angularaxis.title.font.family = matlab2plotlyfont(...
        rLabel.FontName);
    polarAxis.angularaxis.title.font.size = rLabel.FontSize;
    polarAxis.angularaxis.title.font.color = sprintf('rgb(%f,%f,%f)', ...
        255*rLabel.Color);

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

    %-angular tick labels settings-%
    tickValues = axisData.ThetaTick;
    tickLabels = axisData.ThetaTickLabel;
    showTickLabels = true;

    try
        if tickValues(1) == 0 && tickValues(end) == 360
            tickValues = tickValues(1:end-1);
        end
    catch
        tickValues = tickValues;
    end

    if isempty(tickValues)
        showTickLabels = false;
        polarAxis.angularaxis.showticklabels = showTickLabels;
        polarAxis.angularaxis.ticks = '';

    elseif isempty(tickLabels)
        polarAxis.angularaxis.tickvals = tickValues;

    else
        polarAxis.angularaxis.tickvals = tickValues;
        polarAxis.angularaxis.ticktext = tickLabels;

    end

    if showTickLabels
        switch thetaAxis.TickDirection
            case 'in'
                polarAxis.angularaxis.ticks = 'inside';
            case 'out'
                polarAxis.angularaxis.ticks = 'outside';
        end

        %-tick font-%
        polarAxis.angularaxis.tickfont.family = matlab2plotlyfont(...
            thetaAxis.FontName);
        polarAxis.angularaxis.tickfont.size = thetaAxis.FontSize;
        polarAxis.angularaxis.tickfont.color = sprintf('rgb(%f,%f,%f)', ...
            255*thetaAxis.Color);
    end

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

    %-radial tick labels settings-%
    tickValues = axisData.RTick;
    tickLabels = axisData.RTickLabel;
    showTickLabels = true;

    if isempty(tickValues)
        showTickLabels = false;
        polarAxis.radialaxis.showticklabels = showTickLabels;
        polarAxis.radialaxis.ticks = '';

    elseif isempty(tickLabels)
        polarAxis.radialaxis.tickvals = tickValues;

    else
        polarAxis.radialaxis.tickvals = tickValues;
        polarAxis.radialaxis.ticktext = tickLabels;
    end

    if showTickLabels
        switch rAxis.TickDirection
            case 'in'
                polarAxis.radialaxis.ticks = 'inside';
            case 'out'
                polarAxis.radialaxis.ticks = 'outside';
        end

        %-tick font-%
        polarAxis.radialaxis.tickfont.family = matlab2plotlyfont(...
            rAxis.FontName);
        polarAxis.radialaxis.tickfont.size = rAxis.FontSize;
        polarAxis.radialaxis.tickfont.color = sprintf('rgb(%f,%f,%f)', ...
            255*rAxis.Color);
    end

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

    %-set polaraxes to layout-%
    obj.layout = setfield(obj.layout, sprintf('polar%d', xsource+1), polarAxis);

    %-------------------------------------------------------------------------%
end

author: plotly


text:
function [inputImages,inputTraces,infoStruct, algorithmStr] = loadNeurodataWithoutBorders(inputFilePath,varargin)
	% Load NWB file containing cell-extraction data.
	% Biafra Ahanonu
	% started: 2020.04.04 [15:02:22]
	% inputs
		% inputFilePath - Str: path to NWB file. If a cell, will only load the first string.
	% outputs
		% inputImages - 3D or 4D matrix containing cells and their spatial information, format: [x y nCells].
		% inputSignals - 2D matrix containing activity traces in [nCells nFrames] format.
		% infoStruct - contains information about the file, e.g. the 'description' property that can contain information about the algorithm.
		% algorithmStr - String of the algorithm name.

	% changelog
		% 2021.02.05 [19:02:44] - Parse the algorithm associated with the NWB signal extraction data.
		% 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.

	%========================
	% DESCRIPTION
	% cellmax, pcaica, cnmf, cnmfe, extract, roi
	options.algorithm = '';
	% String: main name for image mask group
	options.groupImages = '/processing/ophys/ImageSegmentation/PlaneSegmentation';
	options.planeNo = 1;
	% Dataset name for image masks
	options.imagesName = 'image_mask';
	% String: main name for fluorescence series
	options.groupSignalSeries = '/processing/ophys/Fluorescence/RoiResponseSeries';
	% Cell array of strings: alternative names for fluorescence series
	options.groupSignalSeriesAlt = {'/processing/ophys/Fluorescence/Series','/processing/ophys/Fluorescence/RoiResponseSeries'};
	options.signalsName = 'data';
	% Int: user specified series number if needed
	options.signalSeriesNo = [];
	% 1 = load images, 0 = do not load images (e.g. when user only wants the traces)
	options.loadImages = 1;
	% get options
	options = getOptions(options,varargin);
	% display(options)
	% unpack options into current workspace
	% fn=fieldnames(options);
	% for i=1:length(fn)
	% 	eval([fn{i} '=options.' fn{i} ';']);
	% end
	%========================

	try
		% Check that NWB code is downloaded and setup.
		ciapkg.api.setupNwb('checkDependencies',1);

		inputImages = [];
		inputTraces = [];
		infoStruct = struct;
		algorithmStr = '';

		if iscell(inputFilePath)
			inputFilePath = inputFilePath{1};
		end
		% Return if not a string, return as not appropriate input
		if ~ischar(inputFilePath)
			return;
		end
		try
			imagesGroupName = [options.groupImages num2str(options.planeNo) '/' options.imagesName];
			imagesGroupNameAttr = [options.groupImages num2str(options.planeNo)];
			% Check exist else use different default
			h5info(inputFilePath,imagesGroupName);
		catch
			imagesGroupName = [options.groupImages '/' options.imagesName];
			imagesGroupNameAttr = [options.groupImages];
		end

		if options.loadImages==1
			fprintf('Loading: %s file -> %s.\n',inputFilePath,imagesGroupName);
			inputImages = h5read(inputFilePath,imagesGroupName);
		else
			inputImages = [];
		end

		roitestNames = {options.groupSignalSeries,options.groupSignalSeriesAlt{:}};
		tracesGroupNameFinal = '';
		nNums = 2;
		numList = [1:nNums NaN];
		if ~isempty(options.signalSeriesNo)
			numList = [options.signalSeriesNo numList];
		end
		% Iterate until find the correct name for activity series
		for nameNo = 1:length(roitestNames)
			for seriesNo = numList
				try
					if isnan(seriesNo)==1
						tracesGroupName = [roitestNames{nameNo} '/' options.signalsName];
					else
						tracesGroupName = [roitestNames{nameNo} num2str(seriesNo) '/' options.signalsName];
					end
					% Check exist else use different default
					h5info(inputFilePath,tracesGroupName);
					tracesGroupNameFinal = tracesGroupName;
				catch
				end
			end
		end
		fprintf('Loading: %s file -> %s.\n',inputFilePath,tracesGroupNameFinal);
		inputTraces = h5read(inputFilePath,tracesGroupNameFinal);

		disp(['inputImages: ' num2str(size(inputImages))]);
		disp(['inputTraces: ' num2str(size(inputTraces))]);

		% Get description if later need cell-extraction information
		tmp = h5readatt(inputFilePath,imagesGroupNameAttr,'description');
		infoStruct.description = tmp;

		% Check if extraction method in infoStruct then add
		try
			if isfield(infoStruct,'description')
				tmpStr = regexp(infoStruct.description,'Extraction method: \w+','match');
				try
					tmpStr = strsplit(tmpStr{1}{1},': ');
				catch
					tmpStr = strsplit(tmpStr{1},': ');
				end
				algorithmStr = tmpStr{2};
			end
		catch err
			disp(repmat('@',1,7))
			disp(getReport(err,'extended','hyperlinks','on'));
			disp(repmat('@',1,7))
		end
	catch err
		disp(repmat('@',1,7))
		disp(getReport(err,'extended','hyperlinks','on'));
		disp(repmat('@',1,7))
	end
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function roi = roiFromMesh(seg, mapMethod, mr);
%
% roi = roiFromMesh(seg, [mapMethod='all'], [mr]);
%
% Create an ROI based on vertices in the selected mesh contained
% within a given segmentation struct seg, based on the
% sepecified mapping method.
%
% mapMethod can be 'layer1' or 'all' ('all' by default).
%
% If an mr struct is passed as an optional third argument, will
% check the new ROI coords against these coords and return it
% relative to that mr object. (E.g., will xform the ROI from
% segmented volume -> inplanes).
%
% ras, 01/2007.
if nargin < 1, error('Need roi and segmentation inputs.'); end

if notDefined('mapMethod'),     mapMethod = 'all'; end

msh = segGet(seg, 'SelectedMesh');
mrmRoi = mrmGet(msh,'curRoi'); % ROI as defined on mesh
if ~isfield(mrmRoi,'vertices'), error('No ROI defined on mesh!'); end

% initialize empty ROI based on the same anatomy coords as the segmentation
roi = roiCreate('I|P|R');
roi.voxelSize = msh.mmPerVox;

% the main part of this code will be finding the indices I
% from the gray coordinates which are contained within the
% ROI, given the current mapping method.
% We initialize I to include those layer 1 nodes in the ROI:
layer1Nodes = msh.vertexGrayMap(1,mrmRoi.vertices);
curLayer = unique( layer1Nodes(layer1Nodes>0) );
I = curLayer;

if isequal(lower(mapMethod), 'all')
    [nodes edges] = segGet(seg, 'gray');

    % Start with the ROI vertices, which *should* be just layer 1 nodes.
    curLayerNum = 1;
    while ~isempty(curLayer)
        nextLayer = [];
        curLayerNum = curLayerNum+1;
        for ii = 1:length(curLayer)
            offset = nodes(5,curLayer(ii));
            if offset>length(edges), continue; end
            numConnected = nodes(4,curLayer(ii));
            neighbors = edges(offset:offset+numConnected-1);
            nextLayer = [nextLayer, neighbors(nodes(6,neighbors)==curLayerNum)];
        end
        nextLayer = unique(nextLayer);
        I = [I nextLayer];
        curLayer = nextLayer;
    end
end

grayCoords = segGet(seg, 'GrayCoords');
roi.coords = grayCoords(:,I);

% give the ROI a name that describes approx where it is
cen = round( mean(roi.coords, 2) );
roi.name = sprintf('MeshROI ~(%i, %i, %i)', cen(2), cen(1), cen(3));

if exist('mr', 'var') | ~isempty(mr)
    mr = mrParse(mr);
    roi = roiCheckCoords(roi, mr);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function Out   = RiemannLogMap(P,X)

[U Delta] = eig(P);
G = U*sqrt(Delta);
Y = inv(G)*X*inv(G)';
[V Sigma] = eig(Y);
Out = (G*V)*diag(log(diag(Sigma)))*(G*V)';

author: alexandrebarachant

repo: covariancetoolbox


text:
function tutorial_hcp(tutorial_dir)
% TUTORIAL_HCP: Script that reproduces the results of the online tutorial "Human Connectome Project: Resting-state MEG".
%
% CORRESPONDING ONLINE TUTORIALS:
%     https://neuroimage.usc.edu/brainstorm/Tutorials/HCP-MEG
%
% INPUTS:
%     tutorial_dir: Directory where the HCP files have been unzipped

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Author: Francois Tadel, 2017

%% ===== FILES TO IMPORT =====
% You have to specify the folder in which the tutorial dataset is unzipped
if (nargin == 0) || isempty(tutorial_dir) || ~file_exist(tutorial_dir)
    error('The first argument must be the full path to the tutorial dataset folder.');
end
% Subject name
SubjectName = '175237';
% Build the path of the files to import
AnatDir    = fullfile(tutorial_dir, SubjectName, 'MEG', 'anatomy');
Run1File   = fullfile(tutorial_dir, SubjectName, 'unprocessed', 'MEG', '3-Restin', '4D', 'c,rfDC');
NoiseFile  = fullfile(tutorial_dir, SubjectName, 'unprocessed', 'MEG', '1-Rnoise', '4D', 'c,rfDC');
% Check if the folder contains the required files
if ~file_exist(AnatDir) || ~file_exist(Run1File) || ~file_exist(NoiseFile)
    error(['The folder ' tutorial_dir ' does not contain subject #175237 from the HCP-MEG distribution.']);
end

%% ===== CREATE PROTOCOL =====
% The protocol name has to be a valid folder name (no spaces, no weird characters...)
ProtocolName = 'TutorialHcp';
% Start brainstorm without the GUI
if ~brainstorm('status')
    brainstorm nogui
end
% Delete existing protocol
gui_brainstorm('DeleteProtocol', ProtocolName);
% Create new protocol
gui_brainstorm('CreateProtocol', ProtocolName, 0, 0);
% Start a new report
bst_report('Start');

%% ===== IMPORT DATA =====
% Process: Import anatomy folder
bst_process('CallProcess', 'process_import_anatomy', [], [], ...
    'subjectname', SubjectName, ...
    'mrifile',     {AnatDir, 'HCPv3'}, ...
    'nvertices',   15000);

% Process: Create link to raw files
sFilesRun1 = bst_process('CallProcess', 'process_import_data_raw', [], [], ...
    'subjectname',  SubjectName, ...
    'datafile',     {Run1File, '4D'}, ...
    'channelalign', 1);
sFilesNoise = bst_process('CallProcess', 'process_import_data_raw', [], [], ...
    'subjectname',  SubjectName, ...
    'datafile',     {NoiseFile, '4D'}, ...
    'channelalign', 1);
sFilesRaw = [sFilesRun1, sFilesNoise];

%% ===== PRE-PROCESSING =====
% Process: Notch filter: 60Hz 120Hz 180Hz 240Hz 300Hz
sFilesNotch = bst_process('CallProcess', 'process_notch', sFilesRaw, [], ...
    'freqlist',    [60, 120, 180, 240, 300], ...
    'sensortypes', 'MEG, EEG', ...
    'read_all',    1);

% Process: High-pass:0.3Hz
sFilesBand = bst_process('CallProcess', 'process_bandpass', sFilesNotch, [], ...
    'sensortypes', 'MEG, EEG', ...
    'highpass',    0.3, ...
    'lowpass',     0, ...
    'attenuation', 'strict', ...  % 60dB
    'mirror',      0, ...
    'useold',      0, ...
    'read_all',    1);

% Process: Power spectrum density (Welch)
sFilesPsdAfter = bst_process('CallProcess', 'process_psd', sFilesBand, [], ...
    'timewindow',  [0 100], ...
    'win_length',  4, ...
    'win_overlap', 50, ...
    'sensortypes', 'MEG, EEG', ...
    'edit',        struct(...
         'Comment',         'Power', ...
         'TimeBands',       [], ...
         'Freqs',           [], ...
         'ClusterFuncTime', 'none', ...
         'Measure',         'power', ...
         'Output',          'all', ...
         'SaveKernel',      0));

% Mark bad channels
bst_process('CallProcess', 'process_channel_setbad', sFilesBand, [], ...
            'sensortypes', 'A227, A244, A246, A248');

% Process: Snapshot: Frequency spectrum
bst_process('CallProcess', 'process_snapshot', sFilesPsdAfter, [], ...
    'target',         10, ...  % Frequency spectrum
    'modality',       1);      % MEG (All)

% Process: Delete folders
bst_process('CallProcess', 'process_delete', [sFilesRaw, sFilesNotch], [], ...
    'target', 2);  % Delete folders

%% ===== ARTIFACT CLEANING =====
% Process: Select data files in: */*
sFilesBand = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname', 'All');

% Process: Select file names with tag: 3-Restin
sFilesRest = bst_process('CallProcess', 'process_select_tag', sFilesBand, [], ...
    'tag',    '3-Restin', ...
    'search', 1, ...  % Search the file names
    'select', 1);  % Select only the files with the tag

% Process: Detect heartbeats
bst_process('CallProcess', 'process_evt_detect_ecg', sFilesRest, [], ...
    'channelname', 'ECG+, -ECG-', ...
    'timewindow',  [], ...
    'eventname',   'cardiac');

% Process: SSP ECG: cardiac
bst_process('CallProcess', 'process_ssp_ecg', sFilesRest, [], ...
    'eventname',   'cardiac', ...
    'sensortypes', 'MEG', ...
    'usessp',      1, ...
    'select',      1);

% Process: Snapshot: Sensors/MRI registration
bst_process('CallProcess', 'process_snapshot', sFilesRest, [], ...
    'target',         1, ...  % Sensors/MRI registration
    'modality',       1, ...  % MEG (All)
    'orient',         1);  % left

% Process: Snapshot: SSP projectors
bst_process('CallProcess', 'process_snapshot', sFilesRest, [], ...
    'target',         2, ...  % SSP projectors
    'modality',       1);     % MEG (All)

%% ===== SOURCE ESTIMATION =====
% Process: Select file names with tag: task-rest
sFilesNoise = bst_process('CallProcess', 'process_select_tag', sFilesBand, [], ...
    'tag',    '1-Rnoise', ...
    'search', 1, ...  % Search the file names
    'select', 1);  % Select only the files with the tag

% Process: Compute covariance (noise or data)
bst_process('CallProcess', 'process_noisecov', sFilesNoise, [], ...
    'baseline',       [], ...
    'sensortypes',    'MEG', ...
    'target',         1, ...  % Noise covariance     (covariance over baseline time window)
    'dcoffset',       1, ...  % Block by block, to avoid effects of slow shifts in data
    'identity',       0, ...
    'copycond',       1, ...
    'copysubj',       0, ...
    'replacefile',    1);  % Replace

% Process: Compute head model
bst_process('CallProcess', 'process_headmodel', sFilesRest, [], ...
    'sourcespace', 1, ...  % Cortex surface
    'meg',         3);     % Overlapping spheres

% Process: Compute sources [2018]
sSrcRest = bst_process('CallProcess', 'process_inverse_2018', sFilesRest, [], ...
    'output',  2, ...  % Kernel only: one per file
    'inverse', struct(...
         'Comment',        'dSPM: MEG', ...
         'InverseMethod',  'minnorm', ...
         'InverseMeasure', 'dspm2018', ...
         'SourceOrient',   {{'fixed'}}, ...
         'Loose',          0.2, ...
         'UseDepth',       1, ...
         'WeightExp',      0.5, ...
         'WeightLimit',    10, ...
         'NoiseMethod',    'reg', ...
         'NoiseReg',       0.1, ...
         'SnrMethod',      'fixed', ...
         'SnrRms',         1e-06, ...
         'SnrFixed',       3, ...
         'ComputeKernel',  1, ...
         'DataTypes',      {{'MEG'}}));

%% ===== POWER MAPS =====
% Process: Power spectrum density (Welch)
sSrcPsd = bst_process('CallProcess', 'process_psd', sSrcRest, [], ...
    'timewindow',  [0, 100], ...
    'win_length',  4, ...
    'win_overlap', 50, ...
    'clusters',    {}, ...
    'scoutfunc',   1, ...  % Mean
    'edit',        struct(...
         'Comment',         'Power,FreqBands', ...
         'TimeBands',       [], ...
         'Freqs',           {{'delta', '2, 4', 'mean'; 'theta', '5, 7', 'mean'; 'alpha', '8, 12', 'mean'; 'beta', '15, 29', 'mean'; 'gamma1', '30, 59', 'mean'; 'gamma2', '60, 90', 'mean'}}, ...
         'ClusterFuncTime', 'none', ...
         'Measure',         'power', ...
         'Output',          'all', ...
         'SaveKernel',      0));

% Process: Spectrum normalization
sSrcPsdNorm = bst_process('CallProcess', 'process_tf_norm', sSrcPsd, [], ...
    'normalize', 'relative', ...  % Relative power (divide by total power)
    'overwrite', 0);

% Process: Spatial smoothing (3.00)
sSrcPsdNorm = bst_process('CallProcess', 'process_ssmooth_surfstat', sSrcPsdNorm, [], ...
    'fwhm',      3, ...
    'overwrite', 1);

% Screen capture of final result
hFig = view_surface_data([], sSrcPsdNorm.FileName);
set(hFig, 'Position', [200 200 200 200]);
hFigContact = view_contactsheet(hFig, 'freq', 'fig');
bst_report('Snapshot', hFigContact, sSrcPsdNorm.FileName, 'Power');
close([hFig, hFigContact]);

% Save and display report
ReportFile = bst_report('Save', []);
bst_report('Open', ReportFile);

author: brainstorm-tools


text:
function [M, rng] = intRescale(M, type);
%
% [M, rng] = intRescale(M, <type = int16>);
%
% Take a double or other numeric matrix,
% rescale to the range of the specified
% data type (by default 'int16') and return
% M in that type, as well as the range
% [min max] in M before scaling. This
% is intended to be used for dealing
% with large data sets (e.g., time series
% from many scans) which take up tons of
% memory as doubles.
%
% To convert M back to its double form,
% use:
%       M = normalize(M, rng(1), rng(2));
%
% Note that this will introduce slight inaccuracies
% in the reconstructed matrix (see example below). These
% are usually ~1/1000 of the actual values, and if the
% stored numbers are integers, rounding will resolve this.
% But, it can be an issue for applications which require
% great precision (greater than the number of integers in
% the int16 type).
%
% EXAMPLE:
%   [new, rng]  = intRescale(1:10);
%   old = normalize(new, rng(1), rng(2));
%
% ras, 06/05.
if ~exist('type','var') | isempty(type),
    type = 'int16';
end

rng = double([min(M(:)) max(M(:))]);

% get rescaled min/max, based on data type:
% there are functions for this in ML7, but
% in earlier matlab we have to kludge it:
v = version;
if str2num(v(1))>7
    resMin = intmin(type);
    resMax = intmax(type);
else
    resMin = -32768;
    resMax = 32768;
end

if diff(rng)==0
    % if the requested range is 0 (min==max), then we just need to remove the
    % current offset (=rng(1), which also =rng(2)) and apply the new scale
    % and offset.
    M = round((M-rng(1)) * (resMax-resMin) + resMin);
else
    M = round( (M-rng(1)) ./ diff(rng) * (resMax-resMin) + resMin );
end

cmd = sprintf('M = %s(M);',type);
eval(cmd);

return


text:
function U = som_umat(sMap, varargin)

%SOM_UMAT Compute unified distance matrix of self-organizing map.
%
% U = som_umat(sMap, [argID, value, ...])
%
%  U = som_umat(sMap);
%  U = som_umat(M,sTopol,'median','mask',[1 1 0 1]);
%
%  Input and output arguments ([]'s are optional):
%   sMap     (struct) map struct or
%            (matrix) the codebook matrix of the map
%   [argID,  (string) See below. The values which are unambiguous can
%    value]  (varies) be given without the preceeding argID.
%
%   U        (matrix) u-matrix of the self-organizing map
%
% Here are the valid argument IDs and corresponding values. The values which
% are unambiguous (marked with '*') can be given without the preceeding argID.
%   'mask'       (vector) size dim x 1, weighting factors for different
%                         components (same as BMU search mask)
%   'msize'      (vector) map grid size
%   'topol'     *(struct) topology struct
%   'som_topol','sTopol' = 'topol'
%   'lattice'   *(string) map lattice, 'hexa' or 'rect'
%   'mode'      *(string) 'min','mean','median','max', default is 'median'
%
% NOTE! the U-matrix is always calculated for 'sheet'-shaped map and
% the map grid must be at most 2-dimensional.
%
% For more help, try 'type som_umat' or check out online documentation.
% See also SOM_SHOW, SOM_CPLANE.

%%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% som_umat
%
% PURPOSE
%
% Computes the unified distance matrix of a SOM.
%
% SYNTAX
%
%  U = som_umat(sM)
%  U = som_umat(...,'argID',value,...)
%  U = som_umat(...,value,...)
%
% DESCRIPTION
%
% Compute and return the unified distance matrix of a SOM.
% For example a case of 5x1 -sized map:
%            m(1) m(2) m(3) m(4) m(5)
% where m(i) denotes one map unit. The u-matrix is a 9x1 vector:
%    u(1) u(1,2) u(2) u(2,3) u(3) u(3,4) u(4) u(4,5) u(5)
% where u(i,j) is the distance between map units m(i) and m(j)
% and u(k) is the mean (or minimum, maximum or median) of the
% surrounding values, e.g. u(3) = (u(2,3) + u(3,4))/2.
%
% Note that the u-matrix is always calculated for 'sheet'-shaped map and
% the map grid must be at most 2-dimensional.
%
% REFERENCES
%
% Ultsch, A., Siemon, H.P., "Kohonen's Self-Organizing Feature Maps
%   for Exploratory Data Analysis", in Proc. of INNC'90,
%   International Neural Network Conference, Dordrecht,
%   Netherlands, 1990, pp. 305-308.
% Kohonen, T., "Self-Organizing Map", 2nd ed., Springer-Verlag,
%    Berlin, 1995, pp. 117-119.
% Iivarinen, J., Kohonen, T., Kangas, J., Kaski, S., "Visualizing
%   the Clusters on the Self-Organizing Map", in proceedings of
%   Conference on Artificial Intelligence Research in Finland,
%   Helsinki, Finland, 1994, pp. 122-126.
% Kraaijveld, M.A., Mao, J., Jain, A.K., "A Nonlinear Projection
%   Method Based on Kohonen's Topology Preserving Maps", IEEE
%   Transactions on Neural Networks, vol. 6, no. 3, 1995, pp. 548-559.
%
% REQUIRED INPUT ARGUMENTS
%
%  sM (struct) SOM Toolbox struct or the codebook matrix of the map.
%     (matrix) The matrix may be 3-dimensional in which case the first
%              two dimensions are taken for the map grid dimensions (msize).
%
% OPTIONAL INPUT ARGUMENTS
%
%  argID (string) Argument identifier string (see below).
%  value (varies) Value for the argument (see below).
%
%  The optional arguments are given as 'argID',value -pairs. If the
%  value is unambiguous, it can be given without the preceeding argID.
%  If an argument is given value multiple times, the last one is used.
%
%  Below is the list of valid arguments:
%   'mask'      (vector) mask to be used in calculating
%                        the interunit distances, size [dim  1]. Default is
%                        the one in sM (field sM.mask) or a vector of
%                        ones if only a codebook matrix was given.
%   'topol'     (struct) topology of the map. Default is the one
%                        in sM (field sM.topol).
%   'sTopol','som_topol' (struct) = 'topol'
%   'msize'     (vector) map grid dimensions
%   'lattice'   (string) map lattice 'rect' or 'hexa'
%   'mode'      (string) 'min', 'mean', 'median' or 'max'
%                        Map unit value computation method. In fact,
%                        eval-function is used to evaluate this, so
%                        you can give other computation methods as well.
%                        Default is 'median'.
%
% OUTPUT ARGUMENTS
%
%  U   (matrix) the unified distance matrix of the SOM
%               size 2*n1-1 x 2*n2-1, where n1 = msize(1) and n2 = msize(2)
%
% EXAMPLES
%
%  U = som_umat(sM);
%  U = som_umat(sM.codebook,sM.topol,'median','mask',[1 1 0 1]);
%  U = som_umat(rand(10,10,4),'hexa','rect');
%
% SEE ALSO
%
%  som_show    show the selected component planes and the u-matrix
%  som_cplane  draw a 2D unified distance matrix

% Copyright (c) 1997-2000 by the SOM toolbox programming team.
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 1.0beta juuso 260997
% Version 2.0beta juuso 151199, 151299, 200900

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% check arguments

error(nargchk(1, Inf, nargin));  % check no. of input arguments is correct

% sMap
if isstruct(sMap),
  M = sMap.codebook;
  sTopol = sMap.topol;
  mask = sMap.mask;
elseif isnumeric(sMap),
  M = sMap;
  si = size(M);
  dim = si(end);
  if length(si)>2, msize = si(1:end-1);
  else msize = [si(1) 1];
  end
  munits = prod(msize);
  sTopol = som_set('som_topol','msize',msize,'lattice','rect','shape','sheet');
  mask = ones(dim,1);
  M = reshape(M,[munits,dim]);
end
mode = 'median';

% varargin
i=1;
while i<=length(varargin),
  argok = 1;
  if ischar(varargin{i}),
    switch varargin{i},
      % argument IDs
     case 'mask',       i=i+1; mask = varargin{i};
     case 'msize',      i=i+1; sTopol.msize = varargin{i};
     case 'lattice',    i=i+1; sTopol.lattice = varargin{i};
     case {'topol','som_topol','sTopol'}, i=i+1; sTopol = varargin{i};
     case 'mode',       i=i+1; mode = varargin{i};
      % unambiguous values
     case {'hexa','rect'}, sTopol.lattice = varargin{i};
     case {'min','mean','median','max'}, mode = varargin{i};
     otherwise argok=0;
    end
  elseif isstruct(varargin{i}) && isfield(varargin{i},'type'),
    switch varargin{i}(1).type,
     case 'som_topol', sTopol = varargin{i};
     case 'som_map',   sTopol = varargin{i}.topol;
     otherwise argok=0;
    end
  else
    argok = 0;
  end
  if ~argok,
    disp(['(som_umat) Ignoring invalid argument #' num2str(i+1)]);
  end
  i = i+1;
end

% check
[munits dim] = size(M);
if prod(sTopol.msize)~=munits,
  error('Map grid size does not match the number of map units.')
end
if length(sTopol.msize)>2,
  error('Can only handle 1- and 2-dimensional map grids.')
end
if prod(sTopol.msize)==1,
  warning('Only one codebook vector.'); U = []; return;
end
if ~strcmp(sTopol.shape,'sheet'),
  disp(['The ' sTopol.shape ' shape of the map ignored. Using sheet instead.']);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% initialize variables

y = sTopol.msize(1);
x = sTopol.msize(2);
lattice = sTopol.lattice;
shape = sTopol.shape;
M = reshape(M,[y x dim]);

ux = 2 * x - 1;
uy = 2 * y - 1;
U  = zeros(uy, ux);

calc = sprintf('%s(a)',mode);

if size(mask,2)>1, mask = mask'; end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% u-matrix computation

% distances between map units

if strcmp(lattice, 'rect'), % rectangular lattice

  for j=1:y, for i=1:x,
      if i<x,
	dx = (M(j,i,:) - M(j,i+1,:)).^2; % horizontal
	U(2*j-1,2*i) = sqrt(mask'*dx(:));
      end
      if j<y,
	dy = (M(j,i,:) - M(j+1,i,:)).^2; % vertical
	U(2*j,2*i-1) = sqrt(mask'*dy(:));
      end
      if j<y && i<x,
	dz1 = (M(j,i,:) - M(j+1,i+1,:)).^2; % diagonals
	dz2 = (M(j+1,i,:) - M(j,i+1,:)).^2;
	U(2*j,2*i) = (sqrt(mask'*dz1(:))+sqrt(mask'*dz2(:)))/(2 * sqrt(2));
      end
    end
  end

elseif strcmp(lattice, 'hexa') % hexagonal lattice

  for j=1:y,
    for i=1:x,
      if i<x,
	dx = (M(j,i,:) - M(j,i+1,:)).^2; % horizontal
	U(2*j-1,2*i) = sqrt(mask'*dx(:));
      end

      if j<y, % diagonals
	dy = (M(j,i,:) - M(j+1,i,:)).^2;
	U(2*j,2*i-1) = sqrt(mask'*dy(:));

	if rem(j,2)==0 && i<x,
	  dz= (M(j,i,:) - M(j+1,i+1,:)).^2;
	  U(2*j,2*i) = sqrt(mask'*dz(:));
	elseif rem(j,2)==1 && i>1,
	  dz = (M(j,i,:) - M(j+1,i-1,:)).^2;
	  U(2*j,2*i-2) = sqrt(mask'*dz(:));
	end
      end
    end
  end

end

% values on the units

if (uy == 1 || ux == 1),
  % in 1-D case, mean is equal to median

  ma = max([ux uy]);
  for i = 1:2:ma,
    if i>1 && i<ma,
      a = [U(i-1) U(i+1)];
      U(i) = eval(calc);
    elseif i==1, U(i) = U(i+1);
    else U(i) = U(i-1); % i==ma
    end
  end

elseif strcmp(lattice, 'rect')

  for j=1:2:uy,
    for i=1:2:ux,
      if i>1 && j>1 && i<ux && j<uy,    % middle part of the map
	a = [U(j,i-1) U(j,i+1) U(j-1,i) U(j+1,i)];
      elseif j==1 && i>1 && i<ux,        % upper edge
	a = [U(j,i-1) U(j,i+1) U(j+1,i)];
      elseif j==uy && i>1 && i<ux,       % lower edge
	a = [U(j,i-1) U(j,i+1) U(j-1,i)];
      elseif i==1 && j>1 && j<uy,        % left edge
	a = [U(j,i+1) U(j-1,i) U(j+1,i)];
      elseif i==ux && j>1 && j<uy,       % right edge
	a = [U(j,i-1) U(j-1,i) U(j+1,i)];
      elseif i==1 && j==1,              % top left corner
	a = [U(j,i+1) U(j+1,i)];
      elseif i==ux && j==1,             % top right corner
	a = [U(j,i-1) U(j+1,i)];
      elseif i==1 && j==uy,             % bottom left corner
	a = [U(j,i+1) U(j-1,i)];
      elseif i==ux && j==uy,            % bottom right corner
	a = [U(j,i-1) U(j-1,i)];
      else
	a = 0;
      end
      U(j,i) = eval(calc);
    end
  end

elseif strcmp(lattice, 'hexa')

  for j=1:2:uy,
    for i=1:2:ux,
      if i>1 && j>1 && i<ux && j<uy,      % middle part of the map
	a = [U(j,i-1) U(j,i+1)];
	if rem(j-1,4)==0, a = [a, U(j-1,i-1) U(j-1,i) U(j+1,i-1) U(j+1,i)];
	else a = [a, U(j-1,i) U(j-1,i+1) U(j+1,i) U(j+1,i+1)]; end
      elseif j==1 && i>1 && i<ux,        % upper edge
	a = [U(j,i-1) U(j,i+1) U(j+1,i-1) U(j+1,i)];
      elseif j==uy && i>1 && i<ux,       % lower edge
	a = [U(j,i-1) U(j,i+1)];
	if rem(j-1,4)==0, a = [a, U(j-1,i-1) U(j-1,i)];
	else a = [a, U(j-1,i) U(j-1,i+1)]; end
      elseif i==1 && j>1 && j<uy,        % left edge
	a = U(j,i+1);
	if rem(j-1,4)==0, a = [a, U(j-1,i) U(j+1,i)];
	else a = [a, U(j-1,i) U(j-1,i+1) U(j+1,i) U(j+1,i+1)]; end
      elseif i==ux && j>1 && j<uy,       % right edge
	a = U(j,i-1);
	if rem(j-1,4)==0, a=[a, U(j-1,i) U(j-1,i-1) U(j+1,i) U(j+1,i-1)];
	else a = [a, U(j-1,i) U(j+1,i)]; end
      elseif i==1 && j==1,              % top left corner
	a = [U(j,i+1) U(j+1,i)];
      elseif i==ux && j==1,             % top right corner
	a = [U(j,i-1) U(j+1,i-1) U(j+1,i)];
      elseif i==1 && j==uy,             % bottom left corner
	if rem(j-1,4)==0, a = [U(j,i+1) U(j-1,i)];
	else a = [U(j,i+1) U(j-1,i) U(j-1,i+1)]; end
      elseif i==ux && j==uy,            % bottom right corner
	if rem(j-1,4)==0, a = [U(j,i-1) U(j-1,i) U(j-1,i-1)];
	else a = [U(j,i-1) U(j-1,i)]; end
      else
	a=0;
      end
      U(j,i) = eval(calc);
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% normalization between [0,1]

% U = U - min(min(U));
% ma = max(max(U)); if ma > 0, U = U / ma; end

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

author: ilarinieminen

repo: SOM-Toolbox


text:
function model = yalmip2scs(interfacedata);

model.data.A = -interfacedata.F_struc(:,2:end);
model.data.b = full(interfacedata.F_struc(:,1));
model.data.c =  interfacedata.c;
model.cones = interfacedata.K;
model.param = interfacedata.options.scs;

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function [ x, seed ] = sphere_unit_sample2_nd ( dim_num, seed )

%*****************************************************************************80
%
%% SPHERE_UNIT_SAMPLE2_ND picks a random point on the unit sphere in ND.
%
%  Discussion:
%
%    The unit sphere in ND satisfies:
%
%      sum ( 1 <= I <= DIM_NUM ) X(I) * X(I) = 1
%
%    DIM_NUM independent normally distributed random numbers are generated,
%    and then scaled to have unit norm.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the space.
%
%    Input, integer SEED, a seed for the random number generator.
%
%    Output, real X(DIM_NUM), the random point.
%
%    Output, integer SEED, a seed for the random number generator.
%
  [ x, seed ] = r8vec_normal_01 ( dim_num, seed );

  norm = sqrt ( sum ( x(1:dim_num).^2 ) );

  x(1:dim_num) = x(1:dim_num) / norm;

  return
end


text:
function diffNm = maxNorm_ind(V, V0, U, U0)
n = length(V0);
diffV = zeros(n, 1);
diffU = zeros(n, 1);
for i = 1: n
    diffV(i) = norm(V{i} - V0{i}, 'fro');
    diffU(i) = norm(U{i} - U0{i}, 'fro');
end
diffNm = max(mean(diffV), mean(diffU));

end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
%########################################################################
%
%	- PPGI Toolbox -
%   A MATLAB toolbox for Photoplethysmography Imaging (PPGI)
%
% Author   : Christian S. Pilz
% Company  : The Nature of Space of Time
% Date     : 07.05.2019
%
% Contact  : cpi@partofthestars.com
% Web Page : www.partofthestars.com
%
% Version  : beta0.1
%
%########################################################################
%
%	local_group_invariance.m:
%
% Description:
%
%   Implements the Local Group Invariance algorithm.
%   NOTE: the implementation considers the translation invariance solely.
%
% References:
%
%   Christian S. Pilz, S. Zaunseder, J. Krajewski, V. Blazek,
%   Local Group Invariance for Heart Rate Estimation from Face Videos in the Wild,
%   The IEEE Conference on Computer Vision and Pattern Recognition (CVPR) Workshops,
%   pp.1254-1262, Salt Lake City, 2018
%

classdef local_group_invariance

   properties

   end

   methods

       function obj = local_group_invariance()

       end

       function [pulse, obj] = get(obj,skin_pixels)

           pulse=[];
           frames=size(skin_pixels,2);
           C=[];

           for f=1:frames
               x=skin_pixels{f};
               C(f,:)=[mean(x(:,1)); mean(x(:,2)); mean(x(:,3));];
           end

           [U,E,V]=svd(C');

           S=U(:,1)';
           P=eye(3)-S'*S;%rank 1

           F(1:frames,3)=0;
           for f=1:frames
              F(f,:)=(P*C(f,:)')';
           end

           pulse=double(F(:,2));
       end
   end
end

author: partofthestars

repo: PPGI-Toolbox


text:
function [out] = baseflow_7(p1,p2,S,dt)
%baseflow_7

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Non-linear outflow from a reservoir
% Constraints:  f <= S/dt
%               S >= 0
% @(Inputs):    p1   - time coefficient [d-1]
%               p2   - exponential scaling parameter [-]
%               S    - current storage [mm]
%               dt   - time step size [d]

out = min(S/dt,p1.*max(0,S).^p2);

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
function e = dempot(x, mix)
%DEMPOT	Computes the negative log likelihood for a mixture model.
%
%	Description
%	This function computes the negative log of the unconditional data
%	density P(X) for a Gaussian mixture model.  The data structure MIX
%	defines the mixture model, while the matrix X contains the data
%	vectors.
%
%	See also
%	DEMGPOT, DEMHMC1, DEMMET1
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Computes the potential (negative log likelihood)
e = -log(gmmprob(mix, x));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function z = minus( x, y, cheat )

%   Disciplined convex programming information for MINUS:
%      Terms in a difference must have opposite curvature. Real affine
%      expressions are both convex and concave, so they can be involved
%      in a difference with any nonlinear expression. Complex affine (or
%      constant) expressions, however, are neither, so they can only be
%      involved in differences with other affine expressions. So, for
%      example, the following differences are valid:
%         {convex}-{concave}   {concave}-{convex}   {affine}-{affine}
%      The following are not:
%         {convex}-{concave}  {convex}-{complex constant}
%      For vectors, matrices, and arrays, these rules are verified
%      independently for each element.
%
%   Disciplined geometric programming information for MINUS:
%      Non-constant expressions (log-convex or log-concave) may not be
%      involved in a subtraction in disciplined geometric programs.

if nargin < 3, cheat = false; end
z = plus( x, y, true, cheat );

% Copyright 2005-2016 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: cvxr

repo: CVX

save_path: github-repos/MATLAB/cvxr-CVX


text:
%% FUNCTION Logistic_Trace
%   Trace-Norm Regularized Learning with Logistic Loss.
%
%% OBJECTIVE
%   argmin_{W,C} { sum_i^t (- sum(log (1./ (1+ exp(-X{i}*W(:, i) - Y{i} .* C(i)))))/length(Y{i}))
%            + rho1 * \|W\|_*}
%  where \|W\|_* = sum(svd(W, 0)) is the trace norm
%
%% INPUT
%   X: {n * d} * t - input matrix
%   Y: {n * 1} * t - output matrix
%   rho1: trace norm regularization parameter
%
%% OUTPUT
%   W: model: d * t
%   C: model: 1 * t
%   funcVal: function value vector.
%
%% LICENSE
%   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/>.
%
%   Copyright (C) 2011 - 2012 Jiayu Zhou and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 3, 2012.
%
%% RELATED PAPERS
%
%   [1] Ji, S. and Ye, J. An Accelerated Gradient Method for Trace Norm Minimization, ICML 2009
%
%% RELATED FUNCTIONS
%  Least_Trace, init_opts

%% Code starts here
function [W, C, funcVal] = Logistic_Trace(X, Y, rho1, opts)

if nargin <3
    error('\n Inputs: X, Y, and rho1 should be specified!\n');
end
X = multi_transpose(X);

if nargin <4
    opts = [];
end

% initialize options.
opts=init_opts(opts);

if isfield(opts, 'rho_L2')
    rho_L2 = opts.rho_L2;
else
    rho_L2 = 0;
end

task_num  = length (X);
dimension = size(X{1}, 1);
funcVal = [];

%initialize a starting point
C0_prep = zeros(1, task_num);
for t_idx = 1: task_num
    m1 = nnz(Y{t_idx} == 1);
    m2 = nnz(Y{t_idx} == -1);
    if ( m1==0 || m2==0 )
        C0_prep(t_idx) = 0;
    else
        C0_prep(t_idx) = log(m1/m2);
    end
end

if opts.init==2
    W0 = zeros(dimension, task_num);
    C0 = zeros(1, task_num);
elseif opts.init== 0
    W0 = randn(dimension, task_num);
    C0 = C0_prep;
else
    if isfield(opts,'W0')
        W0=opts.W0;
        if (nnz(size(W0)-[dimension, task_num]))
            error('\n Check the input .W0');
        end
    else
        W0 = zeros(dimension, task_num);
    end
    if isfield(opts,'C0')
        C0=opts.C0;
    else
        C0=C0_prep;
    end
end

bFlag=0; % this flag tests whether the gradient step only changes a little

Wz= W0;
Cz= C0;
Wz_old = W0;
Cz_old = C0;

t = 1;
t_old = 0;
iter = 0;
gamma = 1;
gamma_inc = 2;

while iter < opts.maxIter
    alpha = (t_old - 1) /t;

    Ws = (1 + alpha) * Wz - alpha * Wz_old;
    Cs = (1 + alpha) * Cz - alpha * Cz_old;

    % compute function value and gradients of the search point
    [gWs, gCs, Fs ]  = gradVal_eval(Ws, Cs);

    % the Armijo Goldstein line search scheme
    while true
        [Wzp Wzp_tn] = trace_projection(Ws - gWs/gamma, 2 * rho1 / gamma);
        Czp = Cs - gCs/gamma;
        Fzp = funVal_eval  (Wzp, Czp);

        delta_Wzp = Wzp - Ws;
        delta_Czp = Czp - Cs;
        nrm_delta_Wzp = norm(delta_Wzp, 'fro')^2;
        nrm_delta_Czp = norm(delta_Czp, 'fro')^2;
        r_sum = (nrm_delta_Wzp+nrm_delta_Czp)/2;

        %         Fzp_gamma = Fs + trace(delta_Wzp' * gWs)...
        %             + trace(delta_Czp' * gCs)...
        %             + gamma/2 * nrm_delta_Wzp ...
        %             + gamma/2 * nrm_delta_Czp;
        Fzp_gamma = Fs + sum(sum(delta_Wzp .* gWs))...
            + sum(sum(delta_Czp .* gCs))...
            + gamma/2 * nrm_delta_Wzp ...
            + gamma/2 * nrm_delta_Czp;

        if (r_sum <=1e-20)
            bFlag=1; % this shows that, the gradient step makes little improvement
            break;
        end

        if (Fzp <= Fzp_gamma)
            break;
        else
            gamma = gamma * gamma_inc;
        end
    end

    Wz_old = Wz;
    Cz_old = Cz;
    Wz = Wzp;
    Cz = Czp;

    %funcVal = cat(1, funcVal, Fzp + rho1 * sum( svd(Wzp, 0) ));
    funcVal = cat(1, funcVal, Fzp + rho1 * Wzp_tn);

    if (bFlag)
        % fprintf('\n The program terminates as the gradient step changes the solution very small.');
        break;
    end

    % test stop condition.
    switch(opts.tFlag)
        case 0
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <= opts.tol)
                    break;
                end
            end
        case 1
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <=...
                        opts.tol* funcVal(end-1))
                    break;
                end
            end
        case 2
            if ( funcVal(end)<= opts.tol)
                break;
            end
        case 3
            if iter>=opts.maxIter
                break;
            end
    end

    iter = iter + 1;
    t_old = t;
    t = 0.5 * (1 + (1+ 4 * t^2)^0.5);

end

W = Wzp;
C = Czp;

% private functions
    function [grad_W, grad_C, funcVal] = gradVal_eval(W, C)
        grad_W = zeros(dimension, task_num);
        grad_C = zeros(1, task_num);
        lossValVect = zeros (1 , task_num);
        if opts.pFlag
            parfor i = 1:task_num
                [ grad_W(:, i), grad_C(:, i), lossValVect(:, i)] = unit_grad_eval( W(:, i), C(i), X{i}, Y{i});
            end
        else
            for i = 1:task_num
                [ grad_W(:, i), grad_C(:, i), lossValVect(:, i)] = unit_grad_eval( W(:, i), C(i), X{i}, Y{i});
            end
        end
        grad_W = grad_W+ rho_L2 * 2 * W;
        % here when computing function value we do not include
        % l1 norm.
        funcVal = sum(lossValVect) + rho_L2 * norm(W, 'fro')^2;
    end

    function [funcVal] = funVal_eval (W, C )
        funcVal = 0;
        if opts.pFlag
            parfor i = 1: task_num
                funcVal = funcVal + unit_funcVal_eval( W(:, i), C(i), X{i}, Y{i});
            end
        else
            for i = 1: task_num
                funcVal = funcVal + unit_funcVal_eval( W(:, i), C(i), X{i}, Y{i});
            end
        end
        % here when computing function value we do not include
        % l1 norm.
        funcVal = funcVal + rho_L2 * norm(W, 'fro')^2;
    end

end

function [ grad_w, grad_c, funcVal ] = unit_grad_eval( w, c, x, y)
%gradient and logistic evaluation for each task
m = length(y);
weight = ones(m, 1)/m;
weighty = weight.* y;
aa = -y.*(x'*w + c);
bb = max( aa, 0);
funcVal = weight'* ( log( exp(-bb) +  exp(aa-bb) ) + bb );
pp = 1./ (1+exp(aa));
b = -weighty.*(1-pp);
grad_c = sum(b);
grad_w = x * b;
end

function [ funcVal ] = unit_funcVal_eval( w, c, x, y)
%function value evaluation for each task
m = length(y);
weight = ones(m, 1)/m;
aa = -y.*(x'*w + c);
bb = max( aa, 0);
funcVal = weight'* ( log( exp(-bb) +  exp(aa-bb) ) + bb );
end


text:
function  [cpdag] = dag_to_cpdag1(dags)
% 2
% (also works with a cell array of dags, returning a cell array of cpdags)
% DAG_TO_CPDAG produce a N*N matrix which values respect :
%
% 	If the edge is compelled then 1 on the edge.
%	If the edge is reversible then 1 on the edge and 1 in the reverse edge.
%
% Make sure that the entry is a DAG.
%
% See D.M. Chickering: "Learning Equivalence Classes of Bayesian-Network Structures".
%
%
% francois.olivier.c.h@gmail.com, philippe.leray@univ-nantes.fr

if ~iscell(dags)
    dag=cell(1,1);
    dag{1}=dags;
else
    dag=dags;
end

for da=1:length(dag)
    cpdags{da} = abs(label_edges(dag{da}));
end

if ~iscell(dags)
    cpdag=cpdags{1};
else
    cpdag=cpdags;
end

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

function [label] = label_edges(dag)
% LABEL-EDGES produce a N*N matrix which values are
% 	+1 if the edge is compelled or
%	-1 if the edge is reversible.
% Make sure that the entry is a DAG.
%
% francois.olivier.c.h@gmail.com

N=length(dag);
[order xedge yedge] = order_edges(dag);
label = 2*dag; % all edges as unknown

NbEdges = length(xedge) ;
%xedge=x, yedge=y, %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for Edge=1:NbEdges,
    xlow=xedge(Edge);
    ylow=yedge(Edge);
    if label(xlow,ylow)==2
        fin = 0;

        %wcompelled = find(label(:,xlow)~=2);
        wcompelled = find(label(:,xlow)==1);

        parenty = find(label(:,ylow)~=0);
        %sonsy = find(label(ylow,:)~=0);

        for s = 1:length(wcompelled)
           w = wcompelled(s);
           if ~ismember(w,parenty)
                label(xlow,ylow)=1;
                label(ylow,xlow)=0; %
                label(parenty,ylow)=1;
                label(ylow,parenty)=0; %
                %label(ylow,sonsy)=1; %
                %label(sonsy,ylow)=0; %
                fin = 1;

            elseif fin == 0
                label(w,ylow)=1;
                %label(ylow,w)=0; %
            end
        end
        if fin == 0
            parentx = [xlow ; find(label(:,xlow)~=0)];
            if ~isempty(mysetdiff(parenty,parentx))
                %label(xlow,ylow)=1; %
                %label(ylow,xlow)=0; %

                label(find(label(:,ylow)==2),ylow)=1;
                label(ylow,find(label(ylow,:)==2))=1; %
            else
                label(xlow,ylow)=-1;
                label(ylow,xlow)=-1; %
                ttp = find(label(:,ylow)==2);
                label(ttp,ylow)=-1;
                label(ylow,ttp)=-1; %
            end
        end
    end
end

%%%========================================================================================
function [order, x, y] = order_edges(dag)
% ORDER_EDGES produce a total (natural) ordering over the edges in a DAG.
% Make sure that the entry is a DAG.
%
% francois.olivier.c.h@gmail.com
%
% 2 mai 2003

if acyclic(dag)==0
    error('Requires an acyclic graph');
end

N=length(dag);
order = zeros(N,N);

node_order = topological_sort(dag);
[tmp oo] = sort(node_order);

dag2=dag(node_order,node_order);
[x y]=find(flipud(dag2)==1);
nb_edges=length(x);

if nb_edges~=0
  order(sub2ind([N N],N+1-x,y))=1:nb_edges ;
end

order=order(oo,oo);
x=node_order(N+1-x);
y=node_order(y);


text:
classdef ArgumentParser < handle

properties (Constant)
  opt_id = '9kuetrgshwt364759uedhhrytsj$%284' % not set identifier
end

properties

  requireds
  optionals
  keywords
  parameters
end

methods

  function self = ArgumentParser()

    self.requireds = {};
    self.optionals = {};
    self.keywords = {};
    self.parameters = {};
  end

  function addRequired(self, name, checkfh)
    assert(isempty(self.optionals) && isempty(self.keywords) && isempty(self.parameters), ...
      'Required arguments must come first, before positional arguments, keyword arguments, and parameters.');

    entry = struct;
    entry.name = name;
    entry.checkfh = checkfh;

    self.requireds{end+1} = entry;

  end

  function addOptional(self, name, default, checkfh)

    assert(isempty(self.keywords) && isempty(self.parameters), ...
      'Optional arguments must be before keyword arguments and parameters.');

    entry = struct;
    entry.name = name;
    entry.default = default;
    entry.checkfh = checkfh;

    self.optionals{end+1} = entry;
  end

  function addKeyword(self, name, default, checkfh)

    assert(isempty(self.parameters), 'Keyword arguments must be before parameters.');

    entry = struct;
    entry.name = name;
    entry.default = default;
    entry.checkfh = checkfh ;

    self.keywords{end+1} = entry;
  end

  function addParameter(self, name, default, checkfh)
    entry = struct;
    entry.name = name;
    entry.default = default;
    entry.checkfh = checkfh ;

    self.parameters{end+1} = entry;
  end

  function r = parse(self,varargin)

    req_list = self.requireds;
    opt_list = self.optionals;
    key_list = self.keywords;
    param_list = self.parameters;

    ip = inputParser;

    % required arguments
    for k=1:length(req_list)
      entry = req_list{k};
      ip.addRequired(entry.name, entry.checkfh);
    end

    % positional arguments
    for k=1:length(opt_list)
      entry = opt_list{k};
      ip.addOptional(entry.name, entry.default, entry.checkfh);
    end

    % keyword argument as positional argument
    for k=1:length(key_list)
      entry = key_list{k};
      ip.addOptional([entry.name,'Optional'], self.opt_id, @(el) entry.checkfh(el) || strcmp(el, self.opt_id));
    end

    % keyword arguments as parameter
    for k=1:length(key_list)
      entry = key_list{k};
      ip.addParameter(entry.name, entry.default, entry.checkfh);
    end

    % paramter arguments
    for k=1:length(param_list)
      entry = param_list{k};
      ip.addParameter(entry.name, entry.default, entry.checkfh);
    end

    ip.parse(varargin{:});

    % store parsing results in struct
    r = struct;

    % required
    for k=1:length(req_list)
      name = req_list{k}.name;
      r.(name) = ip.Results.(name);
    end

    % optionals
    for k=1:length(opt_list)
      name = opt_list{k}.name;
      r.(name) = ip.Results.(name);
    end

    % keywords
    for k=1:length(key_list)
      name = key_list{k}.name;

      if strcmp(ip.Results.([name,'Optional']), self.opt_id)
        r.(name) = ip.Results.(name);
      else
        r.(name) = ip.Results.([name,'Optional']);
      end
    end

    % parameters
    for k=1:length(param_list)
      name = param_list{k}.name;
      r.(name) = ip.Results.(name);
    end

  end

end

end

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function disp( prob, varargin )
cvx_display( cvx_validate( prob.index_, prob.id_ ), varargin{:} );

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function varargout = createCube()
%CREATECUBE Create a 3D mesh representing the unit cube.
%
%   [V, E, F] = createCube
%   Create a unit cube, as a polyhedra representation.
%   c has the form [V E F], where V is a 8-by-3 array with vertices
%   coordinates, E is a 12-by-2 array containing indices of neighbour
%   vertices, and F is a 6-by-4 array containing vertices array of each
%   face.
%
%   [V, F] = createCube;
%   Returns only the vertices and the face vertex indices.
%
%   MESH = createCube;
%   Returns the data as a mesh structure, with fields 'vertices', 'edges'
%   and 'faces'.
%
%   Example
%   [n, e, f] = createCube;
%   drawMesh(n, f);
%
%   See also
%   meshes3d, drawMesh
%   createOctahedron, createTetrahedron, createDodecahedron
%   createIcosahedron, createCubeOctahedron
%

% ------
% Author: David Legland
% E-mail: david.legland@inra.fr
% Created: 2005-02-10
% Copyright 2005-2022 INRA - TPV URPOI - BIA IMASTE

x0 = 0; dx= 1;
y0 = 0; dy= 1;
z0 = 0; dz= 1;

nodes = [...
    x0 y0 z0; ...
    x0+dx y0 z0; ...
    x0 y0+dy z0; ...
    x0+dx y0+dy z0; ...
    x0 y0 z0+dz; ...
    x0+dx y0 z0+dz; ...
    x0 y0+dy z0+dz; ...
    x0+dx y0+dy z0+dz];

edges = [1 2;1 3;1 5;2 4;2 6;3 4;3 7;4 8;5 6;5 7;6 8;7 8];

% faces are oriented such that normals point outwards
faces = [1 3 4 2;5 6 8 7;2 4 8 6;1 5 7 3;1 2 6 5;3 7 8 4];

% format output
varargout = formatMeshOutput(nargout, nodes, edges, faces);

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [errors, perrors] = ...
    mcmcTestSuperpixelClassifier(spfeatures, imsegs, vclassifier, hclassifier)

spvlabels = {imsegs(:).vert_labels};
sphlabels = {imsegs(:).horz_labels};

% test vertical
[vdata, vlab, imind] = formatData(spfeatures, spvlabels);
vconf = test_boosted_dt_mc(vclassifier, vdata);
vconf = 1 ./ (1+exp(-vconf));
vconf = vconf ./ repmat(sum(vconf, 2), 1, size(vconf, 2));
[tmp, vmax] = max(vconf, [], 2);
verror = mean(vmax~=vlab)

vconferr = 0;
for k = 1:numel(vlab)
    vconferr = vconferr + 1 - vconf(k, vlab(k));
end
vconferr = vconferr / numel(vlab)

% test subclass
[hdata, hlab] = formatData(spfeatures, sphlabels);
hconf = test_boosted_dt_mc(hclassifier, hdata);
hconf = 1 ./ (1+exp(-hconf));
hconf = hconf ./ repmat(sum(hconf, 2), 1, size(hconf, 2));
[tmp, hmax] = max(hconf, [], 2);
herror = mean(hmax~=hlab)

hconferr = 0;
for k = 1:numel(hlab)
    hconferr = hconferr + 1 - hconf(k, hlab(k));
end
hconferr = hconferr / numel(hlab)

% get same_label error
%[serror, sconferr] = sameLabelError(spfeatures, spvlabels, sphlabels, ...
%    {imsegs(:).adjmat}, vclassifier, hclassifier)

errors.verr = verror;
errors.vcerr = vconferr;
errors.herr = herror;
errors.hcerr = hconferr;
%errors.serr = serror;
%errors.scerr = sconferr;

[hdata2, hlab2, imind] = formatData(spfeatures, spvlabels);
hconf2 = test_boosted_dt_mc(hclassifier, hdata2);
hconf2 = 1 ./ (1+exp(-hconf2));
hconf2 = hconf2 ./ repmat(sum(hconf2, 2), 1, size(hconf2, 2));

lab = (vlab==1)*1 + (vlab==2).*(hlab2+1) + (vlab==3)*7;

pEst = [vconf(:, 1) repmat(vconf(:, 2), 1, 5).*hconf2 vconf(:, 3)];
pTrue = (repmat((1:7), numel(lab), 1)==repmat(lab, 1, 7));

[pMean, pVal, estBias] = evaluateProbabilityEstimate(pEst(:), pTrue(:), 0.05);
perrors.pMean = pMean;
perrors.pVal = pVal;
perrors.estBias = estBias;

perrors.pMean2 = 0;
perrors.estBias2 = 0;
snv = 0;
for k = 1:5
    testind = find((imind > (1-k)*50) & (imind < k*50));
    trainind = setdiff((1:numel(imind)),testind)';

    [pMean, pVal, estBias] = evaluateProbabilityEstimate(pEst(trainind), pTrue(trainind), 0.05);
    perrors.pMeanTr{k} = pMean;

    pEstTest = pEst(testind);
    for i = 1:numel(pEstTest)
        [val, maxind] = min(abs(pVal-pEstTest(i)));
        pEstTest(i) = pMean(maxind);
    end
%    disp(num2str(numel(pEstTest)))
    [pMean, pVal, estBias, nv] = evaluateProbabilityEstimate(pEstTest, pTrue(testind), 0.05);
    snv = snv + nv;
    perrors.pMean2 = perrors.pMean2 + pMean.*nv;
    perrors.estBias2 = perrors.estBias2 + estBias/5;
end
perrors.pMean2 = perrors.pMean2 ./ max(snv, 1E-5);

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

function [data, lab, imind] = formatData(features, labels, labels2)
% concatenate data and labels, and give adjacency information

nimages = numel(features);
[tmp, nvars] = size(features{1});
% count superpixels
nsp = 0;
for f = 1:nimages
    nsp = nsp + sum(labels{f}~=0);
end
data = zeros(nsp, nvars);
lab = zeros(nsp, 1);
imind = zeros(nsp, 1);

% concatenate data
vc = 0;
for f = 1:nimages
    ind = find(labels{f}~=0);
    data(vc+1:vc+numel(ind), :) = features{f}(ind, :);
    lab(vc+1:vc+numel(ind)) = labels{f}(ind);
    imind(vc+1:vc+numel(ind)) = f;
    vc = vc + numel(ind);
end

function [serror, sconferr] = sameLabelError(features, vlabels, hlabels, adjmat, vclassifier, hclassifier)

nimages = numel(features);
[tmp, nvars] = size(features{1});
% count superpixels
nsp = 0;
for f = 1:nimages
    nsp = nsp + sum((vlabels{f}~=0) & ((hlabels{f}~=0) | (vlabels{f}~=2)));
end
data = zeros(nsp, nvars);
lab = zeros(nsp, 1);
adj = cell(nsp, 1);
issame = cell(nsp, 1);

% concatenate data
vc = 0;
nadj = 0;
for f = 1:nimages
    ind = find((vlabels{f}~=0) & ((hlabels{f}~=0) | (vlabels{f}~=2)));
    data(vc+1:vc+numel(ind), :) = features{f}(ind, :);
    for k = 1:numel(ind)
        restind = ind(k+1:end);
        adj{vc+k} = restind(find(adjmat{f}(ind(k), restind)));
        issame{vc+k} = zeros(numel(adj{vc+k}), 1);
        for k2 = 1:numel(adj{vc+k})
            issame{vc+k}(k2) = (vlabels{f}(ind(k))==vlabels{f}(adj{vc+k}(k2))) && ...
                (hlabels{f}(ind(k))==hlabels{f}(adj{vc+k}(k2)));
        end
        nadj = nadj  +numel(adj{vc+k});
    end
    vc = vc + numel(ind);
end

vconf = test_boosted_dt_mc(vclassifier, data);
vconf = 1 ./ (1+exp(-vconf));
vconf = vconf ./ repmat(sum(vconf, 2), 1, size(vconf, 2));
hconf = test_boosted_dt_mc(hclassifier, data);
hconf = 1 ./ (1+exp(-hconf));
hconf = hconf ./ repmat(sum(hconf, 2), 1, size(hconf, 2));

sconf = zeros(nadj, 1);
slab = zeros(nadj, 1);
sc = 0;
for k = 1:numel(adj)
    for k2 = 1:numel(adj{k})
        sc = sc + 1;
        slab(sc) = issame{k}(k2);
        sconf(sc) = sum(vconf(k, [1 3]).*vconf(k2, [1 3])) + ...
            sum((vconf(k, 2)*hconf(k, :)) .* (vconf(k2, 2)*hconf(k2, :)));
    end
end

ind1 = find(slab==1);
ind2 = find(slab==0);
disp(num2str([mean(sconf(ind1)) mean(sconf(ind2))]))
[sorted] = sort(sconf);
tval = sorted(floor(numel(sconf)*(1-mean(slab))));
serror = mean(slab~=(sconf>tval));
sconferr = mean(abs(slab-sconf));

author: zouchuhang

repo: LayoutNet