text:
classdef PTKFissureness < PTKPlugin
    % PTKFissureness. Plugin to detect fissures
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %     This is an intermediate stage towards lobar segmentation.
    %
    %     PTKFissureness computes the fissureness by combining two components
    %     generated from the two plugins PTKFissurenessHessianFactor and
    %     PTKFissurenessVesselsFactor.
    %
    %     For more information, see
    %     [Doel et al., Pulmonary lobe segmentation from CT images using
    %     fissureness, airways, vessels and multilevel B-splines, 2012]
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Fissureness'
        ToolTip = 'Fissureness filter for detecting plane-like points with suppression of points close to vessels'
        Category = 'Fissures'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = false
        PTKVersion = '1'
        ButtonWidth = 6
        ButtonHeight = 2
        GeneratePreview = true
        Visibility = 'Developer'
    end

    methods (Static)

        function results = RunPlugin(dataset, reporting)
            fissureness_from_hessian = single(dataset.GetResult('PTKFissurenessHessianFactor').RawImage)/100;
            fissureness_from_vessels = single(dataset.GetResult('PTKFissurenessVesselsFactor').RawImage)/100;
            results = dataset.GetTemplateImage(PTKContext.LungROI);
            results.ChangeRawImage(100*fissureness_from_vessels.*fissureness_from_hessian);
            results.ImageType = PTKImageType.Scaled;
        end

    end
end


text:
function[b]=choose(n,k)
%CHOOSE  Binomial coefficient: CHOOSE(N,K) = N!K!/(N-K)!
%
%   CHOOSE(N,K) returns the binomial coefficient "N choose K", written
%
%         ( N )         N!
%         (   )   =  --------
%         ( K )      K!(N-K)!
%
%   with a little imagination for the left-hand side notation.
%
%   N and K can each either be an array or a scalar, and must of course
%   be the same size if they are both arrays.
%
%   Usage: b=choose(n,k);
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2007--2015 J.M. Lilly --- type 'help jlab_license' for details

if length(k)==1&&length(n)~=1
    k=k+0*n;
elseif length(k)~=1&&length(n)==1
    n=n+0*k;
end
bool=k>n;
k(k>n)=0;

% b=zeros(size(n));
% n,k
% for i=1:length(n)
%     b(i)=(factorial(n(i)))./(factorial(k(i)).*factorial(n(i)-k(i)));
% end

b=(factorial(n))./(factorial(k).*factorial(n-k));
b(bool)=nan;

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
%DEMOIMBOXFILTER  One-line description here, please.
%
%   output = demoImBoxFilter(input)
%
%   Example
%   demoImBoxFilter
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inra.fr
% Created: 2018-11-22,    using Matlab 8.6.0.267246 (R2015b)
% Copyright 2018 INRA - Cepia Software Platform.

% reference image of fixed size
img = zeros([150 150 150], 'double');
img(150:100, 150:100, 150:100) = 255;

% kernel sizes
kernelSizes = 3:4:25;
nSizes = length(kernelSizes);

resRef = zeros(1, nSizes);
resBox = zeros(1, nSizes);

for i = 1:nSizes
    disp(sprintf('iter %d/%d', i, nSizes));

    dim = kernelSizes(i);
    kernel = ones([dim dim dim]) / (dim*dim*dim);
    tic; resf = imfilter(img, kernel); resRef(i) = toc;

    tic; resf2 = imBoxFilter(img, [dim dim dim]); resBox(i) = toc;
end

figure; hold on;
plot(kernelSizes, resRef, 'k');
plot(kernelSizes, resBox, 'b');
legend({'imfilter', 'imBoxFilter'});

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
%POSTPLOT Postprocessing and visualization function.
%
%   [ H ] = POSTPLOT( PROB, VARARGIN ) Postprocessing function to plot and
%   visualize FEATool problems. PROB is a valid finite element problem
%   struct. Accepts the following property/value pairs and optionally
%   returns an array H of plot object handles.
%
%       Property    Value/{Default}           1D 2D 3D Description
%       -----------------------------------------------------------------------------------
%       surfexpr    string expression         x  x  x  Surface plot expression
%       surfhexpr   string expression            x     Expression for surface height plot
%       isoexpr     string expression            x  x  Contour/iso-surface plot expression
%       isolev      scalar{10}|vector            x  x  Number/vector of contour/iso-levels
%       isores      {[20 20 20]}                    x  Isosurface interpolation resolution
%       isomap      {jet}                        x  x  Colormap for iso-lines/surfaces
%       arrowexpr   n_sdim string expressions    x  x  Cell with arrow plot expressions
%       arrowcolor  {r}                          x  x  Arrow color
%       arrowspacing vector {[10 10 (10)]}       x  x  Number of arrows in each direction
%       arrowscale  scalar {1}                   x  x  Arrow size scaling
%       streamexpr  n_sdim string expressions    x  x  Streamline plot expressions
%       streamres   {[20 20 (20)]}               x  x  Streamline interpolation resolution
%       streamstart n_sdim x n_p array {20}      x  x  Streamline start coordinates
%       streamcolor {r}                          x  x  Streamline colors
%       sliceexpr   string expression               x  Slice plot expression
%       slicex      [x-dir. middle of domain]       x  Slice x-position coordinates
%       slicey      [y-dir. middle of domain]       x  Slice y-position coordinates
%       slicez      [z-dir. middle of domain]       x  Slice z-position coordinates
%       deformexpr  n_sdim string expressions    x  x  Cell with deformation expressions
%       deformscale scalar {auto}                x  x  Deformation scaling (geom fraction)
%       solnum      scalar {n_sols}           x  x  x  Solution number/time to plot
%       evalstyle   exact|{average}           x  x     Shape function evaluation type
%       gradrec     scalar {1}                x  x  x  Use gradient recovery for deriv.
%                                                      up to i:th order shape functions
%       axis        off|{on}                  x  x  x  Show axes
%       grid        on |{off}                 x  x  x  Show grid
%       view        {[-35 20]}                      x  3D view setting
%       axequal     off|{on}                  x  x  x  Axis equal setting
%       edges       on |{off}                       x  Show internal edges for surface plot
%       boundary    off|{on}                     x  x  Show geometry boundary
%       bbox        {0.05}                    x  x  x  Size of bounding box (0=off)
%       alpha       {1}                             x  Transparency level
%       facecolor   {interp}                     x  x  Cell face style or fixed color
%       linestyle   {1/3D:-}|{2D:none}        x  x  x  Line style
%       linewidth   {1}                       x  x     Line width
%       edgecolor   {k}                          x     Cell edge line color
%       color       {k}                       x        Line color
%       marker      {none}                    x        Marker style
%       markersize  {5}                       x        Marker size
%       markercolor {k}                       x        Marker color
%       cellabels   on|{off}                  x  x  x  Show cell numbers
%       nodelabels  on|{off}                  x  x  x  Show vertex/node numbers
%       fontsize    {axes default}            x  x  x  Font size used in text
%       colormap    {jet}                        x  x  Colormap for surface/slice plots
%       caxis       {auto}                    x  x  x  Surface/contour/iso/slice plot range
%       colorbar    {1D:off}|{2/3D:on}        x  x  x  Colorbar toggle
%       location    {eastoutside}             x  x  x  Colorbar location
%       selexpr     string expression         x  x  x  Expression to select cells to plot
%       selcells    {all}                     x  x  x  Index vector to cells to plot
%       selsubd     {all}                     x  x  x  Index vector to subdomains to plot
%       title       string                    x  x  x  Plotly plot title
%       renderer    {empty string}|plotly     x  x  x  Render in MATLAB or Plotly
%       plotlyfile  {default}                 x  x  x  Filename to use for plotly export
%       parent      {gca}                     x  x  x  Plot axes handle
%       setaxes     {on}|off                  x  x  x  Use axes modification settings

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function rsom_show(sMap, D, varargin)

%RSOM_SHOW visualizes information about a given RSOM, currently only the
%U-matrix.
%
% rsom_show(sMap, D);
%
% Input and output arguments:
%   sMap     (struct) RSOM map struct
%   D        (matrix) Dissimilarity matrix used for training of the RSOM
%
% For more help, try 'type rsom_lininit' or check out online documentation.
% See also RSOM_UMAT.

%%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% rsom_show
%
% PURPOSE
%
% Visualizes information concerning the given RSOM, currently only the
% U-matrix.
%
% SYNTAX
%
%  sM = rsom_show(sMap, D);
%
% DESCRIPTION
%
% Visualizes the U-matrix defined by the SOM.
%
% REFERENCES
%
%   Barbara Hammer, Alexander Hasenfuss: Topographic Mapping of Large
%   Dissimilarity Data Sets. Neural Computation 22(9): 2229-2284 (2010)
%
% REQUIRED INPUT ARGUMENTS
%
%   sMap     (struct) RSOM map struct
%   D        (matrix) Dissimilarity matrix used for training of the RSOM
%
% EXAMPLES
%
%   sM = rsom_show(sMap, D);
%
% SEE ALSO
%
%   rsom_umat    Calculates the U-matrix

% Contributed to SOM Toolbox vs2, December 7th, 2012 by Alexander Schulz
% Copyright (c) Alexander Schulz
% http://www.cis.hut.fi/projects/somtoolbox/

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

h_axes = subplot(1,1,1);

lattice = sMap.topol.lattice;
msize   = sMap.topol.msize;

u = rsom_umat(sMap, D);
u = u(:);

tmp_h=som_cplane([lattice 'U'],msize,u);
set(tmp_h,'EdgeColor','none');
set(h_axes,'Tag','Uplane');
h_label=xlabel('U-matrix');
set(h_label,'Visible','on','verticalalignment','top');
set(gca,'plotboxaspectratio',[msize(2) msize(1) msize(1)]);

set(h_label,'interpreter','none');
colormap(flipud(gray(64).^2));
colorbar;

author: ilarinieminen

repo: SOM-Toolbox


text:
function varargout = drawEdgeLabels(p, e, value)
%DRAWEDGELABELS Draw values associated to graph edges.
%
%   usage:
%   drawEdgeLabels(NODES, EDGES, VALUES);
%   NODES: array of double, containing x and y values of nodes
%   EDGES: array of int, containing indices of in and out nodes
%   VALUES is an array the same length of EDGES, containing values
%   associated to each edges of the graph.
%
%   The function computes the center of each edge, and puts the text with
%   associated value.
%
%   H = drawEdgeLabels(...) return array of handles to each text structure,
%   making possible to change font, color, size
%

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

if length(p) > 1 && length(e) > 1
    h = zeros(length(e), 1);
    hold on;
    for l=1:length(e)
        % indices of source and target nodes
        n1 = e(l, 1);
        n2 = e(l, 2);

        % node coordinates
        x1 = p(n1, 1);
        y1 = p(n1, 2);
        x2 = p(n2, 1);
        y2 = p(n2, 2);

        % display the edge
        line([x1 x2], [y1 y2]);

        % coordinates of edge label
        xm = (x1 + x2)/2;
        ym = (y1 + y2)/2;

        % display label
        h(l) = text(xm, ym, sprintf('%3d', floor(value(l))));
    end
end

if nargout == 1
    varargout = {h};
end

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
%CODEGENERATOR.GENFRICTION Generate code for joint friction
%
% F = cGen.genfriction() is the symbolic vector (1xN) of joint friction
% forces.
%
% Notes::
% - Side effects of execution depends on the cGen flags:
%   - saveresult: the symbolic expressions are saved to
%     disk in the directory specified by cGen.sympath
%   - genmfun: ready to use m-functions are generated and
%     provided via a subclass of SerialLink stored in cGen.robjpath
%   - genslblock: a Simulink block is generated and stored in a
%     robot specific block library cGen.slib in the directory
%     cGen.basepath
%   - genccode: generates C-functions and -headers in the directory
%     specified by the ccodepath property of the CodeGenerator object.
%   - mex: generates robot specific MEX-functions as replacement for the
%     m-functions mentioned above. Access is provided by the SerialLink
%     subclass. The MEX files rely on the C code generated before.
%
% Author::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
% See also CodeGenerator.CodeGenerator, CodeGenerator.geninvdyn, CodeGenerator.genfdyn.

% Copyright (C) 2012-2014, by Joern Malzahn
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
%
% The code generation module emerged during the work on a project funded by
% the German Research Foundation (DFG, BE1569/7-1). The authors gratefully
% acknowledge the financial support.

function [ F ] = genfriction( CGen )

%% Derivation of symbolic expressions
CGen.logmsg([datestr(now),'\tDeriving joint friction model']);

[~,qd] = CGen.rob.gencoords;
F = CGen.rob.friction(qd);

CGen.logmsg('\t%s\n',' done!');

%% Save symbolic expressions
if CGen.saveresult
    CGen.logmsg([datestr(now),'\tSaving symbolic friction expression']);

    CGen.savesym(F,'friction','friction.mat')

    CGen.logmsg('\t%s\n',' done!');
end

%% M-Functions
if CGen.genmfun
    CGen.genmfunfriction;
end

%% Embedded Matlab Function Simulink blocks
if CGen.genslblock
    CGen.genslblockfriction;
end

%% C-Code
if CGen.genccode
    CGen.genccodefriction;
end

%% MEX
if CGen.genmex
    CGen.genmexfriction;
end

end

author: petercorke

repo: robotics-toolbox-matlab


text:
function [bndinfo, bndinfo_all] = im2boundaries(im, pbim, gconf, dtBnd, dtBnd_fast, dtCont, ...
    gdata, gclassifiers, thresh)
%
% bndinfo = im2boundaries(im, pbim, gconf, dtBnd, dtBnd_fast, dtCont, gdata, gclassifiers)
%
% Finds occluding contours for im.
%
% Input:
%   im: hxwx3 RGB image
%   pbim:  hxwx4 low level probability of boundary (for four orientations)
%   gconf: hxwx7 geometric context likelihoods
%   dtBnd: 1x3 boundary classifiers for each stage
%   dtBnd_fast: 1x3 fast boundary classifiers for each stage
%   dtCont: 1x3 continuity classifiers for each stage (first is empty)
%   gdata.imsegs: the superpixel structure for the geometric context classifers
%   gclassifiers: the geometric context classifiers
%   thresh: 1x3 thresholds for region merging (default = [0.105 0.25 0.6])
%
% Output:
%   bndinfo.wseg: the superpixels
%   bndinfo.edges.indices: the indices for each edge
%   bndinfo.result.edgeProb: probability of edge being on
%   bndinfo.result.geomProb: probabilty of geometric label (gnd, vert, sky)
%   bndinfo.result.boundaries: most likely ege label (0=off, 1=on)
%
% Notes:
%   - length of result.boundaries is twice length of edges.indices.  If
%   label in first half is "on", left side occludes; if label in second
%   half is "on", right side occludes.
%

global DO_DISPLAY;

DO_DISPLAY = 0;

if ~exist('thresh', 'var') || isempty(thresh)
    thresh = [0.105 0.25 0.6];
end

%% Create bndinfo structure
wseg = pb2wseg(pbim, 100000);

[edges, juncts, neighbors, wseg] = seg2fragments(double(wseg), im, 25);
bndinfo = processBoundaryInfo(wseg, edges, neighbors);

gdata.imsegs = bndinfo2imsegs(bndinfo, gdata.imsegs);

%% Do iter 1 (min merge from local boundary classification)
bndinfo_all{1} = bndinfo;
s = 1;
X = getFeatures(bndinfo, im, pbim, gconf);
[result, pbnd] = mergeStageMin(X, bndinfo, dtBnd(s), dtBnd_fast(s), 2, thresh(s));
bndinfo = updateBoundaryInfo(bndinfo, result, im);
bndinfo_all{1}.pbnd = pbnd{1};
bndinfo_all{2} = bndinfo;

%% Do iter 2 (min merge after CRF inference)

s = 2;
X = getFeatures(bndinfo, im, pbim, gconf);
[result, pbnd] = mergeStageBp(X, bndinfo, dtBnd(s), dtBnd_fast(s), dtCont(s), thresh(s));
bndinfo = updateBoundaryInfo(bndinfo, result, im);
bndinfo_all{2}.pbnd = pbnd{1};
bndinfo_all{3} = bndinfo;

%% Do iter 3..N (min merge after CRF inference, including geometry)

s = 3;

bndinfotmp = bndinfo;
bndinfotmp.labels = (1:bndinfo.nseg);
bndinfotmp = transferSuperpixelLabels(bndinfotmp, wseg);
segmaps = bndinfotmp.labels;

[bndinfo, lab, plab_e, plab_g, pbnd] = mergeStageFinalBpGeometry(...
    bndinfo, dtBnd(s), dtBnd_fast(s), dtCont(s), thresh(s), ...
    im, pbim, gconf, gdata, gclassifiers, segmaps, wseg);

bndinfo.edges.boundaryType = lab;
bndinfo.result.edgeProb = plab_e;
bndinfo.result.geomProb = plab_g;
bndinfo.result.boundaries = lab;
bndinfo.result.thresh = thresh;

bndinfo_all{3}.pbnd = pbnd;
bndinfo_all{4} = bndinfo;

author: Cloud-CV

repo: object-proposals


text:
function  [h, compUp] =  lfmComputeH4JA(gamma1_p, gamma1_m, sigma2, t1, ...
    preFactor, preExp, mode)

% LFMCOMPUTEH4JA Helper function for computing part of the LFMJA kernel.
% FORMAT
% DESC computes a portion of the LFMJA kernel.
% ARG gamma1 : Gamma value for first system.
% ARG gamma2 : Gamma value for second system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG mode: indicates in which way the vectors t1 and t2 must be transposed
% RETURN h : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : Mauricio Alvarez, 2010
%
% SEEALSO : lfmComputeH4.m, lfmComputeH4AA.m

% KERN

if mode==0
    if nargout > 1
        compUp{1} = lfmjpComputeUpsilonVector(gamma1_p,sigma2, t1);
        compUp{2} = lfmjpComputeUpsilonVector(gamma1_m,sigma2, t1);
        h =  compUp{1}*(preExp(:,1)/preFactor(1) - preExp(:,2)/preFactor(2)).' ...
            + compUp{2}*(preExp(:,2)/preFactor(3) - preExp(:,1)/preFactor(4)).';
    else
        h =  lfmjpComputeUpsilonVector(gamma1_p,sigma2, t1)*(preExp(:,1)/preFactor(1) - preExp(:,2)/preFactor(2)).' ...
            + lfmjpComputeUpsilonVector(gamma1_m,sigma2, t1)*(preExp(:,2)/preFactor(3) - preExp(:,1)/preFactor(4)).';
    end
else
    if nargout > 1
        compUp{1} = lfmapComputeUpsilonVector(gamma1_p,sigma2, t1);
        compUp{2} = lfmapComputeUpsilonVector(gamma1_m,sigma2, t1);
        h =  compUp{1}*(preExp(:,1)/preFactor(1) - preExp(:,2)/preFactor(2)).' ...
            + compUp{2}*(preExp(:,2)/preFactor(3) - preExp(:,1)/preFactor(4)).';
    else
        h =  lfmapComputeUpsilonVector(gamma1_p,sigma2, t1)*(preExp(:,2)/preFactor(2) - preExp(:,1)/preFactor(1)).' ...
            + lfmapComputeUpsilonVector(gamma1_m,sigma2, t1)*(preExp(:,1)/preFactor(4) - preExp(:,2)/preFactor(3)).';
    end
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function kern = velotransKernExpandParam(kern, params)

% VELOTRANSKERNEXPANDPARAM Create kernel structure from VELOTRANS kernel's parameters.
% FORMAT
% DESC returns a velocity translate kernel structure filled with the
% parameters in the given vector. This is used as a helper function to
% enable parameters to be optimised in, for example, the NETLAB
% optimisation functions.
% ARG kern : the kernel structure in which the parameters are to be
% placed.
% ARG param : vector of parameters which are to be placed in the
% kernel structure.
% RETURN kern : kernel structure with the given parameters in the
% relevant locations.
%
% SEEALSO : velotransKernParamInit, velotransKernExtractParam, kernExpandParam
%
% COPYRIGHT : Neil D. Lawrence, 2011

% KERN

endVal = length(params)-kern.inputDimension+1;
kern = cmpndKernExpandParam(kern, params(1, 1:endVal));
startVal = endVal + 1;
kern.velocity = params(1, startVal:end);


text:
function linplus_test47 ( )

%*****************************************************************************80
%
%% TEST47 tests R8PBU_FA, R8PBU_RANDOM, R8PBU_SL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  n = 50;
  mu = 1;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST47\n' );
  fprintf ( 1, '  For a banded positive definite symmetric matrix,\n' );
  fprintf ( 1, '  R8PBU_FA computes the LU factors.\n' );
  fprintf ( 1, '  R8PBU_RANDOM sets a random value.\n' );
  fprintf ( 1, '  R8PBU_SL solves a linear system.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Matrix order N = %d\n', n );
  fprintf ( 1, '  Upper bandwidth MU = %d\n', mu );
  fprintf ( 1, '\n' );
%
%  Set the matrix values.
%
  [ a, seed ] = r8pbu_random ( n, mu, seed );
%
%  Set the desired solution.
%
  x = r8vec_indicator ( n );
%
%  Compute the right hand side.
%
  b = r8pbu_mxv ( n, mu, a, x );
%
%  Factor the matrix.
%
  [ a_lu, info ] = r8pbu_fa ( n, mu, a );
%
%  Solve the linear system.
%
  x = r8pbu_sl ( n, mu, a_lu, b );

  r8vec_print_some ( n, x, 1, 10, '  Solution:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function h=ciplot(lower,upper,x,colour,rotate,varargin);

% ciplot(lower,upper)
% ciplot(lower,upper,x)
% ciplot(lower,upper,x,colour)
%
% Plots a shaded region on a graph between specified lower and upper confidence intervals (L and U).
% l and u must be vectors of the same length.
% Uses the 'fill' function, not 'area'. Therefore multiple shaded plots
% can be overlayed without a problem. Make them transparent for total visibility.
% x data can be specified, otherwise plots against index values.
% colour can be specified (eg 'k'). Defaults to blue.

% Raymond Reynolds 24/11/06

if length(lower)~=length(upper)
    error('lower and upper vectors must be same length')
end

if nargin<4
    colour='b';
end

if nargin<3
    x=1:length(lower);
end

if nargin<4
    rotate=0;
end

% convert to row vectors so fliplr can work
if find(size(x)==(max(size(x))))<2
x=x'; end
if find(size(lower)==(max(size(lower))))<2
lower=lower'; end
if find(size(upper)==(max(size(upper))))<2
upper=upper'; end

if rotate
    h=fill([upper fliplr(lower)],[x fliplr(x)],colour,varargin{:});
else
    fill([x fliplr(x)],[upper fliplr(lower)],colour)
end

repo: REST


text:
function pdf = exponential_pdf ( x, a, b )

%*****************************************************************************80
%
%% EXPONENTIAL_PDF evaluates the Exponential PDF.
%
%  Discussion:
%
%    PDF(X)(A,B) = ( 1 / B ) * EXP ( ( A - X ) / B )
%
%    The time interval between two Poisson events is a random
%    variable with the Exponential PDF.  The parameter B is the
%    average interval between events.
%
%    In another context, the Exponential PDF is related to
%    the Boltzmann distribution, which describes the relative
%    probability of finding a system, which is in thermal equilibrium
%    at absolute temperature T, in a given state having energy E.
%    The relative probability is
%
%      Boltzmann_Relative_Probability(E,T) = exp ( - E / ( k * T ) ),
%
%    where k is the Boltzmann constant,
%
%      k = 1.38 * 10**(-23) joules / degree Kelvin
%
%    and normalization requires a determination of the possible
%    energy states of the system.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%    A <= X
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < B.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x < a )
    pdf = 0.0;
  else
    pdf = ( 1.0 / b ) * exp ( ( a - x ) / b );
  end

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of the
%   example 3 DOF strange, non-planar. robot.
%
%   Author: Arturo Gil. Universidad Miguel Hern�ndez de Elche.
%   email: arturo.gil@umh.es date:   05/10/2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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/>.
function robot = parameters()

robot.name='Example 3DOF strange arm';
L0=1;
L1=1;
L2=1;
L3=0.2;
L4=1;

%kinematic data DH parameters
robot.DH.theta='[q(1) q(2)-pi/2 q(3)+pi/2]';
robot.DH.d='[1 1 1]';
robot.DH.a='[.2  0  0]';
robot.DH.alpha='[pi/2  pi/2  0]';

%number of degrees of freedom
robot.DOF = 3;

%Jacobian matrix
robot.J=[];
robot.kind=['R' 'R' 'R'];

%Function name to compute inverse kinematic
robot.inversekinematic_fn = 'inversekinematic_3dofstrange(robot, T)';
robot.directkinematic_fn = 'directkinematic(robot, q)';

%minimum and maximum rotation angle in rad
robot.maxangle =[deg2rad(-180) deg2rad(180); %Axis 1, minimum, maximum
                deg2rad(-180) deg2rad(180);
                deg2rad(-180) deg2rad(180)]; %Axis 2, minimum, maximum

%maximum absolute speed of each joint rad/s or m/s
robot.velmax = []; %empty, not available

robot.accelmax=robot.velmax/0.1; % 0.1 is here an acceleration time
% end effectors maximum velocity
robot.linear_velmax = 0; %m/s, example, not available

%base reference system
robot.T0 = [0 0 1 0; 0 -1 0 0; 1 0 0 1; 0 0 0 1];

%INITIALIZATION OF VARIABLES REQUIRED FOR THE SIMULATION
%position, velocity and acceleration
robot=init_sim_variables(robot);
robot.path = pwd;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GRAPHICS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%read graphics files
robot.graphical.has_graphics=0;
robot.graphical.color = [255 20 40]./255;
%for transparency
robot.graphical.draw_transparent=0;
%draw DH systems
robot.graphical.draw_axes=1;
%DH system length and Font size, standard is 1/10. Select 2/20, 3/30 for
%bigger robots
robot.graphical.axes_scale=2;
%adjust for a default view of the robot
robot.axis=[-3.5 3.5 -3.5 3.5 0 1];
robot = read_graphics(robot);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DYNAMIC PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
robot.has_dynamics=1;

%consider friction in the computations
robot.dynamics.friction=0;

%link masses (kg)
robot.dynamics.masses=[1 1 1];

%COM of each link with respect to own reference system
robot.dynamics.r_com=[-0.5      0         0; %(rx, ry, rz) link 1
                      -0.5      0         0;
                      -0.5      0         0];%(rx, ry, rz) link 2

%link masses
m1 = robot.dynamics.masses(1);
m2 = robot.dynamics.masses(2);
m3 = robot.dynamics.masses(3);

a=eval(robot.DH.a);
L1 = a(1);
L2 = a(2);
L3 = a(3);

%Momentos de inercia de cada eslabon.
% Ixx	Iyy	Izz	Ixy	Iyz	Ixz, por cada fila
robot.dynamics.Inertia=[0   m1*L1^2/3   m1*L1^2/3    0	0	0;
                        0   m2*L2^2/3   m2*L2^2/3    0	0	0;
                        0   m3*L3^2/3   m3*L3^2/3    0	0	0];

%Actuator rotor inertia
robot.motors.Inertia=[0 0 0];
%Reduction ratio motor/joint speed
robot.motors.G=[1 1 1];
%consider friction
robot.friction=0;
%Viscous friction factor, motor referred
%robot.B = [1e-3  1e-3 1e-3];
robot.motors.Viscous = [10  10 10];
%Coulomb friction, motor referred
robot.motors.Coulomb = [0 0;
            0 0;
            0 0];

repo: ARTE


text:
% Psychoacoustic scaling function from [1]
%
% Description
% Based of the experiment in [1]_ wherein participants adjusted a
% second tone until it was half the pitch of the first. The functional
% approximation to the scale is implemented with the formula from
% [2]_ (being honest, from Wikipedia):
%
% ..math:: s = 1127 \ln \left(1 + \frac{f}{700} \right)
%
% Where `s` is the scale and `f` is the frequency in Hertz.
%
% References
% [1] S. S. Stevens, J. Volkmann & E. B. Newman (1937). A Scale for
%    the Measurement of the Psychological Magnitude Pitch. The Journal
%    of the Acoustical Society of America, 8, 185-190.
% [2] O'Shaughnessy, D. (1987). Speech communication: human and
%    machine. Addison-Wesley Pub. Co.
%
%
% Copyright (c) 2018 Department of Computer Science,
%                    University of Toronto, Canada,
%                    Vector Institute, Canada
%
% License
% This file is under the LGPL license,  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. This file 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.
%
% This function is part of the Covarep project: http://covarep.github.io/covarep
%
% Author
%  Yingxue Wang <yingxue@cs.toronto.edu>
%  Sean Robertson <sdrobert@cs.toronto.edu>
%

classdef MelScaling < ScalingFunction
    methods (Static)
        function hertz = scale_to_hertz(scale)
            hertz = 700 * (exp(scale / 1127) - 1);
        end

        function scale = hertz_to_scale(hertz)
            scale = 1127 * log(1 + hertz / 700);
        end
    end
end

author: covarep

repo: covarep

save_path: github-repos/MATLAB/covarep-covarep


text:
function [ac] = acc(label_pre,label_gt)
%Measure percentage of Accuracy and the Rand index of clustering results
% The number of class must equal to the number cluster

%Output
% Acc = Accuracy of clustering results
% rand_index = Rand's Index,  measure an agreement of the clustering results
% match = 2xk mxtrix which are the best match of the Target and clustering results

%Input
% T = 1xn target index
% idx =1xn matrix of the clustering results

% EX:
% X=[randn(200,2);randn(200,2)+6,;[randn(200,1)+12,randn(200,1)]]; T=[ones(200,1);ones(200,1).*2;ones(200,1).*3];
% idx=kmeans(X,3,'emptyaction','singleton','Replicates',5);
%  [Acc,rand_index,match]=Acc_measure(T,idx)

res = bestMap(label_gt(:), label_pre(:));
ac = length(find(label_gt(:) == res(:)))/length(label_gt(:));

return;

k=max(T);
n=length(T);
for i=1:k
    temp=find(T==i);
    a{i}=temp; %#ok<AGROW>
end

b1=[];
t1=zeros(1,k);
for i=1:k
    tt1=find(idx==i);
    for j=1:k
       t1(j)=sum(ismember(tt1,a{j}));
    end
    b1=[b1;t1]; %#ok<AGROW>
end
    Members=zeros(1,k);

P = perms((1:k));
    Acc1=0;
for pi=1:size(P,1)
    for ki=1:k
        Members(ki)=b1(P(pi,ki),ki);
    end
    if sum(Members)>Acc1
        match=P(pi,:);
        Acc1=sum(Members);
    end
end

rand_ss1=0;
rand_dd1=0;
for xi=1:n-1
    for xj=xi+1:n
        rand_ss1=rand_ss1+((idx(xi)==idx(xj))&&(T(xi)==T(xj)));
        rand_dd1=rand_dd1+((idx(xi)~=idx(xj))&&(T(xi)~=T(xj)));
    end
end
rand_index=200*(rand_ss1+rand_dd1)/(n*(n-1));
Acc=Acc1/n*100;
match=[1:k;match];

author: jwyang


text:
% Optimising k-Wave Performance Example
%
% This example demonstrates how to increase the computational performance
% of k-Wave using optional input parameters and data casting. A separate
% standardised benchmarking script benchmark is also included within the
% k-Wave toolbox to allow computational times to be compared across
% different computers and GPUs.
%
% author: Bradley Treeby
% date: 15th July 2009
% last update: 25th August 2014
%
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2014 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave 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.
%
% k-Wave 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 Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

clear all;

example_number = 1;
% 1: default input options
% 2: nearest neighbour Cartesian interpolation and plotting switched off
% 3: as above with 'DataCast' set to 'single'
% 4: as aobve with 'DataCast' set to 'gpuArray-single'

% change scale to 2 to increase the computational time
scale = 1;

% =========================================================================
% SIMULATION
% =========================================================================

% assign the grid size and create the computational grid
Nx = 256*scale;     % number of grid points in the x (row) direction
Nz = 256*scale;     % number of grid points in the y (column) direction
x = 10e-3;          % grid size in the x direction [m]
z = 10e-3;          % grid size in the y direction [m]
dx = x/Nx;          % grid point spacing in the x direction [m]
dz = z/Nz;          % grid point spacing in the y direction [m]
kgrid = makeGrid(Nx, dx, Nz, dz);

% load the initial pressure distribution from an image and scale
p0_magnitude = 2;   % [au]
source.p0 = p0_magnitude*loadImage('EXAMPLE_source_two.bmp');
source.p0 = resize(source.p0, [Nz, Nx]);

% define the properties of the propagation medium
medium.sound_speed = 1500;	% [m/s]

% define a centered Cartesian circular sensor
sensor_radius = 4.5e-3;     % [m]
num_sensor_points = 100;
sensor.mask = makeCartCircle(sensor_radius, num_sensor_points);

% create the time array
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);

% set the input options
switch example_number
    case 1

        % set input arguments
        input_args = {};

    case 2

        % set input arguments
        input_args = {'PlotSim', false};

        % convert Cartesian sensor mask to binary mask
        sensor.mask = cart2grid(kgrid, sensor.mask);

    case 3

        % set input arguments
        input_args = {'PlotSim', false, 'DataCast', 'single'};

        % convert Cartesian sensor mask to binary mask
        sensor.mask = cart2grid(kgrid, sensor.mask);

    case 4

        % set input arguments
        input_args = {'PlotSim', false, 'DataCast', 'gpuArray-single'};

        % convert Cartesian sensor mask to binary mask
        sensor.mask = cart2grid(kgrid, sensor.mask);

        % run GPU warmup
        tic; fft2(gpuArray(rand(128, 128, 'single'))); toc;
        tic; fft2(gpuArray(rand(128, 128, 'single'))); toc;
        tic; fft2(gpuArray(rand(128, 128, 'single'))); toc;

end

% run the simulation
kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});

repo: Aero-matlab


text:
function sts = subsasgn_check_funhandle(val)

% function sts = subsasgn_check_funhandle(val)
% Return true if val is either empty, or a function or function handle.
% One could also check for nargin == 1 and nargout == 1.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: subsasgn_check_funhandle.m 1764 2008-05-30 13:09:40Z volkmar $

rev = '$Rev: 1764 $'; %#ok

sts = isempty(val) || isa(val, 'function_handle') || ...
    (ischar(val) && (any(exist(val) == 2:6) || ~isempty(which(val))));
if sts && isa(val, 'function_handle')
    try
        f = functions(val);
    catch
        % fail silently if "functions" evaluation fails, keep sts == true
        return;
    end;
    try
        if ~strcmp(f.type,'anonymous')
            % file name should not be empty for file functions
            sts = ~isempty(f.file) || ~isempty(which(f.function));
        end;
    catch
        % fail silently if something is wrong with f, keep sts == true
        return;
    end;
end;

author: spm


text:
function [RMS_position, RMS_orientation, NEES_pose, NEES_orientation] = REKF_plot_rms_nees( estimation_results, data, do_vis )
% plot rms and nees for R-EKF
N  = size(estimation_results, 2);

T = 2:N;
RMS_position=[];
RMS_orientation=[];
for i = T
    position = estimation_results{i}.position;
    ap = data.poses.position(:,i);

    RMS_position = [RMS_position norm(position-ap)];
    RMS_orientation=  [RMS_orientation norm(so3_log(estimation_results{i}.orientation*(data.poses.orientation(3*i-2:3*i,1:3))'))];

end

fprintf( 'mean(RMS:Position) = %f\n', mean(RMS_position) );
fprintf( 'mean(RMS:Orientation) = %f\n', mean(RMS_orientation) );

if do_vis == 1
    figure;

    subplot(2,2,1);plot(T, RMS_position);
    title('RMS:position(meter)');xlim([0,N]);
    subplot(2,2,2); plot(T,RMS_orientation);
    title('RMS:orientation(radius)');xlim([0,N]);
end

NEES_pose=[];
NEES_orientation=[];
for i = T
    position = estimation_results{i}.position;
    ap = data.poses.position(:,i);

    dw=so3_log(estimation_results{i}.orientation*(data.poses.orientation(3*i-2:3*i,1:3))');
    dv=(jaco_r(-dw))\( (estimation_results{i}.position-estimation_results{i}.orientation*data.poses.orientation(3*i-2:3*i,1:3)'*data.poses.position(:,i)));
    %dv=( (estimation_results{i}.position-estimation_results{i}.orientation*data.poses.orientation(3*i-2:3*i,1:3)'*data.poses.position(:,i)));

     cov_o=estimation_results{i}.cov(1:3,1:3);
     cov_pose=estimation_results{i}.cov(1:6,1:6);

     invcov_o=eye(3)/cov_o;
     invcov_pose=eye(6)/cov_pose;
    dP=[dw;dv];

    NEES_orientation = [NEES_orientation dw'*invcov_o*dw/3];
    NEES_pose=  [NEES_pose dP'*invcov_pose*dP/6];

end

fprintf( 'mean(NEES:Position) = %f\n', mean(NEES_pose(2:end)) );
fprintf( 'mean(NEES:Orientation) = %f\n', mean(NEES_orientation(2:end)) );

if do_vis == 1
    subplot(2,2,3);
    plot(T,NEES_orientation);
    title('NEES:orientation');xlim([0,N]);
    subplot(2,2,4);
    plot(T,NEES_pose);
    title('NEES:pose');xlim([0,N]);
end

author: RomaTeng

repo: EKF-SLAM-on-Manifold


text:
function [dfofMatrix, inputMovieF0, inputMovieStd] = dfofMovie(inputMovie, varargin)
	% Does deltaF/F and other relative fluorescence changes calculations for a movie using bsxfun for faster processing.
	% Biafra Ahanonu
	% started 2013.11.09 [09:12:36]
	% inputs
		% inputMovie - either a [x y t] matrix or a char string specifying a HDF5 movie.
	% outputs
		%

	[dfofMatrix, inputMovieF0, inputMovieStd] = ciapkg.movie_processing.dfofMovie(inputMovie,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function image3(mat)
%image3
% usage:
% clf
% mat = 0.1*rand(5,6,7);
% mat(2,3,4) = 1;
% image3(mat)
[nx,ny,nz] = size(mat);
c = [1 0 0];
[yrng_org, xrng_org, zrng_org] = meshgrid(1:ny,1:nx,1:nz);
xrng(1,:,:,:) = xrng_org-0.5;
xrng(2,:,:,:) = xrng_org+0.5;
yrng(1,:,:,:) = yrng_org-0.5;
yrng(2,:,:,:) = yrng_org+0.5;
zrng(1,:,:,:) = zrng_org-0.5;
zrng(2,:,:,:) = zrng_org+0.5;
tbl = [1 1 1 2;1 1 2 2; 1 2 2 2; 1 2 1 2];
xtbl = tbl(:,[3 3 4 1 3 3]);
ytbl = tbl(:,[2 2 3 3 4 1]);
ztbl = tbl(:,[4 1 2 2 2 2]);
x_tot = zeros(4,6,nx,ny,nz);
y_tot = x_tot;
z_tot = x_tot;
for vertex = 1:4
    for face = 1:6
        x_tot(vertex,face,:,:,:) = xrng(xtbl(vertex,face),:,:,:);
        y_tot(vertex,face,:,:,:) = yrng(ytbl(vertex,face),:,:,:);
        z_tot(vertex,face,:,:,:) = zrng(ztbl(vertex,face),:,:,:);
    end
end
x_tot = reshape(x_tot, 4, nx*ny*nz*6);
y_tot = reshape(y_tot, 4, nx*ny*nz*6);
z_tot = reshape(z_tot, 4, nx*ny*nz*6);

h = patch(x_tot,y_tot,z_tot,c);
set(h, 'FaceAlpha', 'flat');
alphamat = (mat(:)*ones(1,6))';
set(h,'FaceVertexAlphaData',alphamat(:))
set(h, 'EdgeAlpha', 0);
xlabel('x');
ylabel('y');
zlabel('z');
axis([0.5 nx+0.5 0.5 ny+0.5 0.5 nz+0.5])
view([41,28]);

author: Sable

repo: mcbench-benchmarks


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% vert2.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [x1,x2,f1,f2] = vert2(j,x,z,f,x1,x2,f1,f2)
% subprogram used by vertex.m

function [x1,x2,f1,f2] = vert2(j,x,z,f,x1,x2,f1,f2)
if j == 1
  j1 = 2;
else
  j1 = 1;
end
if x1 == Inf
  x1 = z(j);
  f1 = f1 + f(j);
  if x ~= z(j1)
    x2 = z(j1);
    f2 = f2 + f(j1);
  end
elseif x2 == Inf & x1 ~= z(j)
  x2 = z(j);
  f2 = f2 + f(j);
elseif x2 == Inf
  x2 = z(j1);
  f2 = f2 + f(j1);
end

author: lacerbi

repo: optimviz

save_path: github-repos/MATLAB/lacerbi-optimviz


text:
function [cond1_uptake, cond2_uptake, cond1_secretion, cond2_secretion, slope_Ratio] = defineUptakeSecretionProfiles(input_A,input_B, data_RXNS, tol,essAA_excl,exclude_upt,exclude_secr,add_secr, add_upt)
% This function calculated the slope ratios and gives out uptake&secretion profiles for condition 1 and condition 2
%
% USAGE:
%
%    [cond1_uptake, cond2_uptake, cond1_secretion, cond2_secretion, slope_Ratio] = defineUptakeSecretionProfiles(input_A, input_B, data_RXNS, tol, essAA_excl, exclude_upt, exclude_secr, add_secr, add_upt)
%
% INPUTS:
%    input_A:             Matrix, 4 columns
%
%                           1. Control TP1,
%                           2. Control TP2,
%                           3. Condition 1 TP1,
%                           4. Condition 1 TP 2
%    input_B:             Matrix, 4 columns
%
%                           1. Control2 TP1,
%                           2. Control2 TP2,
%                           3. Condition 2 TP1,
%                           4. Condition 2 TP 2
%    data_RXNS:           Exchange reactions same order as `input_A` and `input_B`
%    tol:                 minimal threshold to call uptake or secretion (Default 0.05)=5%
%
% OPTIONAL INPUTS:
%    essAA_excl:          If essential amino acids should be excluded from the secretion profile 1(yes), or 0(no) (Default = 0);
%                         `essAAs` = {`EX_his_L(e)`; `EX_ile_L(e)`; `EX_leu_L(e)`; `EX_lys_L(e)`; `EX_met_L(e)`; `EX_phe_L(e)`; `EX_thr_L(e)`; `EX_trp_L(e)`; `EX_val_L(e)`}
%    exclude_upt:         Exclude uncertain metabolites from uptake (e.g., metabolites from GlutaMax, e.g., `EX_gln_L(e)`, `EX_cys_L(e)`, `EX_ala_L(e)`)
%    exclude_secr:        Exclude uncertain metabolites from secretion (e.g., metabolites from GlutaMax, e.g., `EX_gln_L(e)`, `EX_cys_L(e)`, `EX_ala_L(e)`)
%    add_secr:            Due to mising data points automatic analysis might do wrong assignment of metabolites to secretion
%    add_upt:             Due to mising data points automatic analysis might do wrong assignment of metabolites to uptake
%
% OUTPUTS:
%    cond1_uptake:        List of exchanges that specify consumed metabolites in condition 1
%    cond2_uptake:        List of exchanges that specify consumed metabolites in condition 2
%    cond1_secretion:     List of exchanges that specify released metabolites in condition 1
%    cond2_secretion:     List of exchanges that specify released metabolites in condition 2
%    slope_Ratio:         For length of inputA/B, relative difference
%
% .. Author: - Maike K. Aurich 27/05/15

if ~exist('tol','var') || isempty(tol)
    tol = 0.05;
end

if ~exist('essAA_excl','var') || isempty(essAA_excl)
    essAA_excl = 0;
end

%% calculate slope for control (TP2/TP1) and condition (TP2/TP1) of both conditions (cell types, or cell lines), and slope ration between conditions

for i=1:length(input_A)

    control1(i,1) = input_A(i,2)/input_A(i,1);
    cond1(i,1) = input_A(i,4)/input_A(i,3);
    slope_cond1(i,1) = cond1(i,1)/ control1(i,1);

    control2(i,1) = input_B(i,2)/input_B(i,1);
    cond2(i,1) = input_B(i,4)/input_B(i,3);
    slope_cond2(i,1) = cond2(i,1)/ control2(i,1);

    slope_Ratio(i,1) = slope_cond2(i,1)/slope_cond1(i,1);

end

%% Make uptake and secretion profiles
%%define uptake and secretion for each condition based on the slope

cond1_uptake = {};
cond2_uptake = {};
cond1_secretion = {};
cond2_secretion = {};

cond1_uptake = data_RXNS(find(slope_cond1<1-tol));
cond2_uptake = data_RXNS(find(slope_cond2<1-tol));

cond1_secretion = data_RXNS(find(slope_cond1>1+tol));
cond2_secretion = data_RXNS(find(slope_cond2>1+tol));

%% exclude secretion of essAAs and modify as specified

if exist('essAA_excl','var')

    cond1_secretion(find(ismember(cond1_secretion, essAA_excl)))='';
    cond2_secretion(find(ismember(cond2_secretion, essAA_excl)))='';

end

if exist('exclude_secr','var')
    cond1_secretion(find(ismember(cond1_secretion, exclude_secr)))='';
    cond2_secretion(find(ismember(cond2_secretion, exclude_secr)))='';
end

if exist('exclude_upt','var')
    cond1_uptake(find(ismember(cond1_uptake, exclude_upt)))='';
    cond2_uptake(find(ismember(cond2_uptake, exclude_upt)))='';

end

if exist('add_secr','var')
    cond1_secretion = [cond1_secretion;add_secr];
    cond2_secretion = [cond2_secretion;add_secr];

end

if exist('add_upt','var')
    cond1_uptake = [cond1_uptake;add_upt];
    cond2_uptake = [cond2_uptake;add_upt];

end

%% make sure vectors are unique
cond1_secretion = unique(cond1_secretion);
cond2_secretion = unique(cond2_secretion);
cond1_uptake = unique(cond1_uptake);
cond2_uptake = unique(cond2_uptake);

 end

author: opencobra

repo: cobratoolbox


text:
function Y = minmaxfilt3(X,filter,varargin)
%  MAXFILT3    Three-dimensional max filter
%
%     Y = MINMAXFILT3(X,'max',[N1 N2 N3]) performs three-dimensional maximum
%     Y = MINMAXFILT3(X,'min',[N1 N2 N3]) performs three-dimensional minimum
%     filtering on the image X using an N1-by-N2-by-N3 window. The result
%     Y contains the maximun value in the N1-by-N2-by-N3 neighborhood around
%     each pixel in the original image.
%     This function uses the van Herk algorithm for max filters.
%
%     Y = MINMAXFILT3(X,M) is the same as Y = MINMAXFILT3(X,[M M M])
%
%     Y = MINMAXFILT3(X) uses a 3-by-3-by-3 neighborhood.
%
%     Y = MINMAXFILT3(..., 'shape') returns a subsection of the 3D
%     filtering specified by 'shape' :
%        'full'  - Returns the full filtering result,
%        'same'  - (default) Returns the central filter area that is the
%                   same size as X,
%        'valid' - Returns only the area where no filter elements are outside
%                  the image.
%
%     See also : MINFILT2, VANHERK
%

% Initialization
[S, shape] = parse_inputs(varargin{:});

% filtering
dim = size(X);
for k=1:dim(3);
	im=squeeze(X(:,:,k));
	im = vanherk(im,S(1),filter,shape);
	im = vanherk(im,S(2),filter,'col',shape);
	X2(:,:,k) = im;
end

if dim(3) < S(3)
    f = str2func(filter);
    im = f(X2,[],3);
    for k=1:dim(3)
        Y(:,:,k) = im;
    end
else
    dim = size(X2);
    for k=1:dim(2);
        im = squeeze(X2(:,k,:));
        im = vanherk(im,S(3),filter,shape);
        Y(:,:,k) = im;
    end
    Y = permute(Y,[1 3 2]);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [S, shape] = parse_inputs(varargin)
shape = 'same';
flag = [0 0]; % size shape

for i = 1 : nargin
   t = varargin{i};
   if strcmp(t,'full') & flag(2) == 0
      shape = 'full';
      flag(2) = 1;
   elseif strcmp(t,'same') & flag(2) == 0
      shape = 'same';
      flag(2) = 1;
   elseif strcmp(t,'valid') & flag(2) == 0
      shape = 'valid';
      flag(2) = 1;
   elseif flag(1) == 0
      S = t;
      flag(1) = 1;
   else
      error(['Too many / Unkown parameter : ' t ])
   end
end

if flag(1) == 0
   S = [3 3 3];
end
if length(S) == 1;
   S(2) = S(1); S(3) = S(1);
end
if length(S) ~= 3
   error('Wrong window size parameter.')
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function lp = create_lp_solve_native_model(A,b,f,xint,LB,UB,e,options);

[m,n] = size(A);
lp = mxlpsolve('make_lp', m, n);

mxlpsolve('set_mat', lp, A);
mxlpsolve('set_rh_vec', lp, b);
mxlpsolve('set_obj_fn', lp, f);
mxlpsolve('set_maxim', lp); % default is solving minimum lp.
for i = 1:length(e)
    if e(i) < 0
        con_type = 1;
    elseif e(i) == 0
        con_type = 3;
    else
        con_type = 2;
    end
    mxlpsolve('set_constr_type', lp, i, con_type);
end
for i = 1:length(LB)
    %if ~isinf(LB(i))
        mxlpsolve('set_lowbo', lp, i, LB(i));
    %end
end
for i = 1:length(UB)
   if ~isinf(UB(i))
        mxlpsolve('set_upbo', lp, i, UB(i));
   end
end
for i = 1:length(xint)
    mxlpsolve('set_int', lp, xint(i), 1);
end

if options.lpsolve.scalemode~=0
    mxlpsolve('set_scaling', lp, scalemode);
end

% for i = 1:length(sos)
%     mxlpsolve('add_SOS', lp, ['dummy' num2str(i)], 1, i, sos{i}, 1:length(sos{i}));
% end

switch options.verbose
    case 0
        mxlpsolve('set_verbose', lp, 0);%options.verbose)
    case 1
        mxlpsolve('set_verbose', lp, 4);%options.verbose)
    case 2
        mxlpsolve('set_verbose', lp, 5);%options.verbose)
    otherwise
        mxlpsolve('set_verbose', lp, 6);%options.verbose)
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function [ imdb ] = get_mnist( opts )

% --------------------------------------------------------------------
% Preapre the imdb structure, returns image data with mean image subtracted
files = {'train-images-idx3-ubyte', ...
         'train-labels-idx1-ubyte', ...
         't10k-images-idx3-ubyte', ...
         't10k-labels-idx1-ubyte'} ;

if ~exist(opts.dataDir, 'dir')
  mkdir(opts.dataDir) ;
end

for i=1:4
  if ~exist(fullfile(opts.dataDir, files{i}), 'file')
    url = sprintf('http://yann.lecun.com/exdb/mnist/%s.gz',files{i}) ;
    fprintf('downloading %s\n', url) ;
    gunzip(url, opts.dataDir) ;
  end
end

f=fopen(fullfile(opts.dataDir, 'train-images-idx3-ubyte'),'r') ;
x1=fread(f,inf,'uint8');
fclose(f) ;
x1=permute(reshape(x1(17:end),28,28,60e3),[2 1 3]) ;

f=fopen(fullfile(opts.dataDir, 't10k-images-idx3-ubyte'),'r') ;
x2=fread(f,inf,'uint8');
fclose(f) ;
x2=permute(reshape(x2(17:end),28,28,10e3),[2 1 3]) ;

f=fopen(fullfile(opts.dataDir, 'train-labels-idx1-ubyte'),'r') ;
y1=fread(f,inf,'uint8');
fclose(f) ;
y1=double(y1(9:end)')+1 ;

f=fopen(fullfile(opts.dataDir, 't10k-labels-idx1-ubyte'),'r') ;
y2=fread(f,inf,'uint8');
fclose(f) ;
y2=double(y2(9:end)')+1 ;

set = [ones(1,numel(y1)) 3*ones(1,numel(y2))];
data = single(reshape(cat(3, x1, x2),28,28,1,[]));
dataMean = mean(data(:,:,:,set == 1), 4);
data = bsxfun(@minus, data, dataMean) ;

imdb.images.data = data ;
imdb.images.data_mean = dataMean;
imdb.images.labels = cat(2, y1, y2) ;
imdb.images.set = set ;
imdb.meta.sets = {'train', 'val', 'test'} ;
imdb.meta.classes = arrayfun(@(x)sprintf('%d',x),0:9,'uniformoutput',false) ;

end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet


text:
function [U, S] = pca(X)
%PCA Run principal component analysis on the dataset X
%   [U, S, X] = pca(X) computes eigenvectors of the covariance matrix of X
%   Returns the eigenvectors U, the eigenvalues (on diagonal) in S
%

% Useful values
[m, n] = size(X);

% You need to return the following variables correctly.
U = zeros(n);
S = zeros(n);

% ====================== YOUR CODE HERE ======================
% Instructions: You should first compute the covariance matrix. Then, you
%               should use the "svd" function to compute the eigenvectors
%               and eigenvalues of the covariance matrix.
%
% Note: When computing the covariance matrix, remember to divide by m (the
%       number of examples).
%

Sigma = (1 / m) * X' * X;
[U, S] = svd(Sigma);

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

end

repo: machine-learning


text:
%MODELS Summarise and search available robot models
%
% MODELS() lists keywords associated with each of the models in Robotics Toolbox.
%
% MODELS(QUERY) lists those models that match the keyword QUERY.  Case is
% ignored in the comparison.
%
% M = MODELS(QUERY) as above but returns a cell array (Nx1) of the names of the
% M-files that define the models.
%
% Examples::
%         models
%         models('modified_DH')  % all models using modified DH notation
%         models('kinova')       % all Kinova robot models
%         models('6dof')         % all 6dof robot models
%         models('redundant')    % all redundant robot models, >6 DOF
%         models('prismatic')    % all robots with a prismatic joint
%
% Notes::
% - A model is a file mdl_*.m in the models folder of the RTB directory.
% - The keywords are indicated by a line '% MODEL: ' after the main comment
%   block.
%

% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
function name_ = models(query)

    % first find the path to the models
    pth = mfilename('fullpath');
    pth = fileparts(pth);
    pth = fullfile(pth, 'models');
    assert(exist(pth, 'dir') > 0, 'RTB:models:nomodel', 'models folder not found');

    models = dir(fullfile(pth, 'mdl_*.m'));
    info = {};
    name = {};
    for model=models'
        fid = fopen( fullfile(pth, model.name), 'r');
        while true
            line = fgetl(fid);
            if line == -1
                break;
            end
            if (length(line) > 11) && (strcmp(line(1:11), '% Copyright:') == 1)
                % an empty line, done with the header block
                break;
            end
            [p,n,e] = fileparts(model.name);

            if (length(line) > 8) && (strcmp(line(1:8), '% MODEL:') == 1)
                % we have a model description line

                line = line(10:end);  % trim off the tag

                if nargin > 0
                    % we have a  query
                    if strfind(lower(line), lower(query))
                        info = [info; [line ' (' n ')' ]];
                        name = [name; n];
                    end
                else
                    % no query, report all
                    info = [info; [line ' (' n ')' ]];
                    name = [name n];
                end
                break
            end
        end
        fclose(fid);
    end

    % optionally return a list of model names
    if nargout == 1
        name_ = name';
    else
        % now sort and print the matching models
        for i = sort(info)'
            fprintf('%s\n', i{1});
        end
    end

end

repo: robotics-toolbox-matlab


text:
function adj = triangulation_order6_adj_set ( node_num, tri_num, ...
  triangle_node, triangle_neighbor, adj_num, adj_col )

%*****************************************************************************80
%
%% TRIANGULATION_ORDER6_ADJ_SET sets adjacencies in a triangulation.
%
%  Discussion:
%
%    This routine is called to count the adjacencies, so that the
%    appropriate amount of memory can be set aside for storage when
%    the adjacency structure is created.
%
%    The triangulation is assumed to involve 6-node triangles.
%
%    Two nodes are "adjacent" if they are both nodes in some triangle.
%    Also, a node is considered to be adjacent to itself.
%
%    This routine can be used to create the compressed column storage
%    for a quadratic triangle finite element discretization of
%    Poisson's equation in two dimensions.
%
%  Diagram:
%
%       3
%    s  |\
%    i  | \
%    d  |  \
%    e  6   5  side 2
%       |    \
%    3  |     \
%       |      \
%       1---4---2
%
%         side 1
%
%    The local node numbering
%
%
%   21-22-23-24-25
%    |\    |\    |
%    | \   | \   |
%   16 17 18 19 20
%    |   \ |   \ |
%    |    \|    \|
%   11-12-13-14-15
%    |\    |\    |
%    | \   | \   |
%    6  7  8  9 10
%    |   \ |   \ |
%    |    \|    \|
%    1--2--3--4--5
%
%    A sample grid.
%
%
%    Below, we have a chart that lists the nodes adjacent to each node, with
%    an asterisk to indicate the adjacency of the node to itself
%    (in some cases, you want to count this self adjacency and in some
%    you don't).
%
%    N   Adjacencies
%
%    1:  *  2  3  6  7 11
%    2:  1  *  3  6  7 11
%    3:  1  2  *  4  5  6  7  8  9 11 12 13
%    4:  3  *  5  8  9 13
%    5:  3  4  *  8  9 10 13 14 15
%    6:  1  2  3  *  7 11
%    7:  1  2  3  6  *  8 11 12 13
%    8:  3  4  5  7  *  9 11 12 13
%    9:  3  4  5  8  * 10 13 14 15
%   10:  5  9  * 13 14 15
%   11:  1  2  3  6  7  8  * 12 13 16 17 21
%   12:  3  7  8 11  * 13 16 17 21
%   13:  3  4  5  7  8  9 10 11 12  * 14 15 16 17 18 19 21 22 23
%   14:  5  9 10 13  * 15 18 19 23
%   15:  5  9 10 13 14  * 18 19 20 23 24 25
%   16: 11 12 13  * 17 21
%   17: 11 12 13 16  * 18 21 22 23
%   18: 13 14 15 17  * 19 21 22 23
%   19: 13 14 15 18  * 20 23 24 25
%   20: 15 19  * 23 24 25
%   21: 11 12 13 16 17 18  * 22 23
%   22: 13 17 18 21  * 23
%   23: 13 14 15 17 18 19 20 21 22  * 24 25
%   24: 15 19 20 23  * 25
%   25: 15 19 20 23 24  *
%
%    Below, we list the number of adjancencies to lower-indexed nodes, to
%    the node itself, to higher-indexed nodes, the total number of
%    adjacencies for this node, and the location of the first and last
%    entries required to list this set of adjacencies in a single list
%    of all the adjacencies.
%
%    N   Below  Self   Above   Total First  Last
%
%   --      --    --      --      --   ---     0
%    1:      0     1       5       6     1     6
%    2:      1     1       4       6     7    12
%    3:      2     1       9      12    13    24
%    4:      1     1       4       6    25    30
%    5:      2     1       6       9    31    39
%    6:      3     1       2       6    40    45
%    7:      4     1       4       9    46    54
%    8:      4     1       4       9    55    63
%    9:      4     1       4       9    62    72
%   10:      2     1       3       6    73    78
%   11:      6     1       5      12    79    90
%   12:      4     1       4       9    91    99
%   13:      9     1       9      19   100   118
%   14:      4     1       4       9   119   127
%   15:      5     1       6      12   128   139
%   16:      3     1       2       6   140   145
%   17:      4     1       4       9   146   154
%   18:      4     1       4       9   155   163
%   19:      4     1       4       9   164   172
%   20:      2     1       3       6   173   178
%   21:      6     1       2       9   179   187
%   22:      4     1       1       6   188   193
%   23:      9     1       2      12   194   205
%   24:      4     1       1       6   206   211
%   25:      5     1       0       6   212   217
%   --      --    --      --      --   218   ---
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 August 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer TRI_NUM, the number of triangles.
%
%    Input, integer TRIANGLE_NODE(6,TRI_NUM), lists the nodes that
%    make up each triangle.  The first three nodes are the vertices,
%    in counterclockwise order.  The fourth value is the midside
%    node between nodes 1 and 2; the fifth and sixth values are
%    the other midside nodes in the logical order.
%
%    Input, integer TRIANGLE_NEIGHBOR(3,TRI_NUM), for each side of
%    a triangle, lists the neighboring triangle, or -1 if there is
%    no neighbor.
%
%    Input, integer ADJ_NUM, the number of adjacencies.
%
%    Input, integer ADJ_COL(NODE_NUM+1).  Information about column J is stored
%    in entries ADJ_COL(J) through ADJ_COL(J+1)-1 of ADJ.
%
%    Output, integer ADJ(ADJ_NUM), the adjacency information.
%
  triangle_order = 6;
  adj(1:adj_num) = -1;
  adj_copy(1:node_num) = adj_col(1:node_num);
%
%  Set every node to be adjacent to itself.
%
  for node = 1 : node_num
    adj(adj_copy(node)) = node;
    adj_copy(node) = adj_copy(node) + 1;
  end
%
%  Examine each triangle.
%
  for triangle = 1 : tri_num

    n1 = triangle_node(1,triangle);
    n2 = triangle_node(2,triangle);
    n3 = triangle_node(3,triangle);
    n4 = triangle_node(4,triangle);
    n5 = triangle_node(5,triangle);
    n6 = triangle_node(6,triangle);
%
%  For sure, we add the adjacencies:
%    43 / (34)
%    51 / (15)
%    54 / (45)
%    62 / (26)
%    64 / (46)
%    65 / (56)
%
    adj(adj_copy(n3)) = n4;
    adj_copy(n3) = adj_copy(n3) + 1;
    adj(adj_copy(n4)) = n3;
    adj_copy(n4) = adj_copy(n4) + 1;

    adj(adj_copy(n1)) = n5;
    adj_copy(n1) = adj_copy(n1) + 1;
    adj(adj_copy(n5)) = n1;
    adj_copy(n5) = adj_copy(n5) + 1;

    adj(adj_copy(n4)) = n5;
    adj_copy(n4) = adj_copy(n4) + 1;
    adj(adj_copy(n5)) = n4;
    adj_copy(n5) = adj_copy(n5) + 1;

    adj(adj_copy(n2)) = n6;
    adj_copy(n2) = adj_copy(n2) + 1;
    adj(adj_copy(n6)) = n2;
    adj_copy(n6) = adj_copy(n6) + 1;

    adj(adj_copy(n4)) = n6;
    adj_copy(n4) = adj_copy(n4) + 1;
    adj(adj_copy(n6)) = n4;
    adj_copy(n6) = adj_copy(n6) + 1;

    adj(adj_copy(n5)) = n6;
    adj_copy(n5) = adj_copy(n5) + 1;
    adj(adj_copy(n6)) = n5;
    adj_copy(n6) = adj_copy(n6) + 1;
%
%  Add edges (1,2), (1,4), (2,4) if this is the first occurrence,
%  that is, if the edge (1,4,2) is on a boundary (TRIANGLE2 <= 0)
%  or if this triangle is the first of the pair in which the edge
%  occurs (TRIANGLE < TRIANGLE2).
%
%  Maybe add
%    21 / 12
%    41 / 14
%    42 / 24
%
    triangle2 = triangle_neighbor(1,triangle);

    if ( triangle2 < 0 | triangle < triangle2 )
      adj(adj_copy(n1)) = n2;
      adj_copy(n1) = adj_copy(n1) + 1;
      adj(adj_copy(n2)) = n1;
      adj_copy(n2) = adj_copy(n2) + 1;
      adj(adj_copy(n1)) = n4;
      adj_copy(n1) = adj_copy(n1) + 1;
      adj(adj_copy(n4)) = n1;
      adj_copy(n4) = adj_copy(n4) + 1;
      adj(adj_copy(n2)) = n4;
      adj_copy(n2) = adj_copy(n2) + 1;
      adj(adj_copy(n4)) = n2;
      adj_copy(n4) = adj_copy(n4) + 1;
    end
%
%  Maybe add
%    32 / 23
%    52 / 25
%    53 / 35
%
    triangle2 = triangle_neighbor(2,triangle);

    if ( triangle2 < 0 | triangle < triangle2 )
      adj(adj_copy(n2)) = n3;
      adj_copy(n2) = adj_copy(n2) + 1;
      adj(adj_copy(n3)) = n2;
      adj_copy(n3) = adj_copy(n3) + 1;
      adj(adj_copy(n2)) = n5;
      adj_copy(n2) = adj_copy(n2) + 1;
      adj(adj_copy(n5)) = n2;
      adj_copy(n5) = adj_copy(n5) + 1;
      adj(adj_copy(n3)) = n5;
      adj_copy(n3) = adj_copy(n3) + 1;
      adj(adj_copy(n5)) = n3;
      adj_copy(n5) = adj_copy(n5) + 1;
    end
%
%  Maybe add
%    31 / 13
%    61 / 16
%    63 / 36
%
    triangle2 = triangle_neighbor(3,triangle);

    if ( triangle2 < 0 | triangle < triangle2 )
      adj(adj_copy(n1)) = n3;
      adj_copy(n1) = adj_copy(n1) + 1;
      adj(adj_copy(n3)) = n1;
      adj_copy(n3) = adj_copy(n3) + 1;
      adj(adj_copy(n1)) = n6;
      adj_copy(n1) = adj_copy(n1) + 1;
      adj(adj_copy(n6)) = n1;
      adj_copy(n6) = adj_copy(n6) + 1;
      adj(adj_copy(n3)) = n6;
      adj_copy(n3) = adj_copy(n3) + 1;
      adj(adj_copy(n6)) = n3;
      adj_copy(n6) = adj_copy(n6) + 1;
    end

  end
%
%  Ascending sort the entries for each node.
%
  for node = 1 : node_num
    k1 = adj_col(node);
    k2 = adj_col(node+1)-1;
    adj(k1:k2) = i4vec_sort_heap_a ( k2+1-k1, adj(k1:k2) );
  end

  return
end


text:
function y = vl_nnbilinear(x, varargin)
% VL_NNBILINEAR  extracts bilinear features at each location
% This function differs from vl_nnbilinerpool whichs pools features across
% locations

% input:
% forward pass:
% x: input featre of size [hight, width, channels, batches]
% backward pass:
% x: input featre of size [hight, width, channels, batches]
% dzdy: the gradient with respect to output y.

% output:
% forward pass:
% y: self outer product of x at each pixel location. The output size if [hight, width, channels*channels, batches]
% backward pass:
% y: graident with respect to x. y will have the same size as input x.

backMode = numel(varargin) > 0 && ~isstr(varargin{1}) ;
if backMode
  dzdy = varargin{1} ;
end

[h, w, ch, bs] = size(x);
gpuMode = isa(x, 'gpuArray');

if backMode

    if gpuMode
        y = gpuArray(zeros(size(x), 'single'));
    else
        y = zeros(size(x), 'single');
    end
    for b=1:bs
        for yy=1:h
            for xx=1:w
                dzdy_b = reshape(dzdy(yy,xx,:,b), [ch, ch]);
                a = squeeze(x(yy,xx,:,b));
                y(yy, xx, :, b) = reshape(a'*dzdy_b, [h, w, ch]);
            end
        end
    end
else

    if gpuMode
        y = gpuArray(zeros([h, w, ch*ch, bs], 'single'));
    else
        y = zeros([h, w, ch*ch, bs], 'single');
    end
    for b = 1:bs,
        for yy = 1:h,
            for xx = 1:w,
                a = squeeze(x(yy,xx,:,b));
                y(yy,xx,:, b) = reshape(a*a', [1 ch*ch]);
            end
        end
    end
end

repo: DBCNN


text:
function [LSM LSM_R] =lsmat(H,R, mtd)
if (mtd==0)
    R_inv=inv(R);
    LSM_R=inv(H'*R_inv*H);
    LSM=LSM_R*H'*R_inv;
elseif (mtd==1)
    [u,s,v]=svd(H);
    T2=u(:,1:size(H,2))';
    T1=u(:,size(H,2)+1:size(H,1))';

    mx_a=inv(T1*R*T1');
    mx_b=inv(T2*H);
    LSM=mx_b*(T2-T2*R*T1'*mx_a*T1);
    LSM_R=LSM*R*LSM';
elseif(mtd==2)
    [u,s,v]=svd(R);
    Ht=u'*H;
    Rt=u'*R*u;
    in=find(diag(Rt)<1e-7, 1, 'first');
    if (isempty(in))
        in=size(H,1);
    else
        in=in-1;
    end
    Htr=Ht(1:in,:);
    Rtr_inv=inv(Rt(1:in,1:in));
    LSM_R=inv(Htr'*Rtr_inv*Htr);
    LSM=LSM_R*Htr'*Rtr_inv*u(:,1:in)';
end

repo: ekfmonoslam


text:
function srf = nrbcoons(u1, u2, v1, v2)
%
% Function Name:
%
%   nrbcoons - Construction of a Coons patch.
%
% Calling Sequence:
%
%   srf = nrbcoons(ucrv1, ucrv2, vcrv1, vcrv2)
%
% Parameters:
%
%  ucrv1	: NURBS curve defining the bottom U direction boundary of
% 		the constructed NURBS surface.
%
%   ucrv2	: NURBS curve defining the top U direction boundary of
% 		the constructed NURBS surface.
%
%   vcrv1	: NURBS curve defining the bottom V direction boundary of
% 		the constructed NURBS surface.
%
%   vcrv2	: NURBS curve defining the top V direction boundary of
% 		the constructed NURBS surface.
%
%   srf		: Coons NURBS surface patch.
%
% Description:
%
%   Construction of a bilinearly blended Coons surface patch from four NURBS
%   curves that define the boundary.
%
%   The orientation of the four NURBS boundary curves.
%
%          ^ V direction
%          |
%          |     ucrv2
%          ------->--------
%          |              |
%          |              |
%    vcrv1 ^   Surface    ^ vcrv2
%          |              |
%          |              |
%          ------->-----------> U direction
%                ucrv1
%
%
% Examples:
%
%   // Define four NURBS curves and construct a Coons surface patch.
%   pnts = [ 0.0  3.0  4.5  6.5 8.0 10.0;
%            0.0  0.0  0.0  0.0 0.0  0.0;
%            2.0  2.0  7.0  4.0 7.0  9.0];
%   crv1 = nrbmak(pnts, [0 0 0 1/3 0.5 2/3 1 1 1]);
%
%   pnts= [ 0.0  3.0  5.0  8.0 10.0;
%           10.0 10.0 10.0 10.0 10.0;
%           3.0  5.0  8.0  6.0 10.0];
%   crv2 = nrbmak(pnts, [0 0 0 1/3 2/3 1 1 1]);
%
%   pnts= [ 0.0 0.0 0.0 0.0;
%           0.0 3.0 8.0 10.0;
%           2.0 0.0 5.0 3.0];
%   crv3 = nrbmak(pnts, [0 0 0 0.5 1 1 1]);
%
%   pnts= [ 10.0 10.0 10.0 10.0 10.0;
%           0.0   3.0  5.0  8.0 10.0;
%           9.0   7.0  7.0 10.0 10.0];
%   crv4 = nrbmak(pnts, [0 0 0 0.25 0.75 1 1 1]);
%
%   srf = nrbcoons(crv1, crv2, crv3, crv4);
%   nrbplot(srf,[20 20],220,45);

%  D.M. Spink
%  Copyright (c) 2000.

if nargin ~= 4
  error('Incorrect number of input arguments');
end

r1 = nrbruled(u1, u2);
r2 = nrbtransp(nrbruled(v1, v2));
t  = nrb4surf(u1.coefs(:,1), u1.coefs(:,end), u2.coefs(:,1), u2.coefs(:,end));

% raise all surfaces to a common degree
du = max([r1.order(1), r2.order(1), t.order(1)]);
dv = max([r1.order(2), r2.order(2), t.order(2)]);
r1 = nrbdegelev(r1, [du - r1.order(1), dv - r1.order(2)]);
r2 = nrbdegelev(r2, [du - r2.order(1), dv - r2.order(2)]);
t  = nrbdegelev(t,  [du - t.order(1),  dv - t.order(2)]);

% merge the knot vectors, to obtain a common knot vector

% U knots
k1 = r1.knots{1};
k2 = r2.knots{1};
k3 = t.knots{1};
k = unique([k1 k2 k3]);
n = length(k);
kua = [];
kub = [];
kuc = [];
for i = 1:n
  i1 = length(find(k1 == k(i)));
  i2 = length(find(k2 == k(i)));
  i3 = length(find(k3 == k(i)));
  m = max([i1, i2, i3]);
  kua = [kua k(i)*ones(1,m-i1)];
  kub = [kub k(i)*ones(1,m-i2)];
  kuc = [kuc k(i)*ones(1,m-i3)];
end

% V knots
k1 = r1.knots{2};
k2 = r2.knots{2};
k3 = t.knots{2};
k = unique([k1 k2 k3]);
n = length(k);
kva = [];
kvb = [];
kvc = [];
for i = 1:n
  i1 = length(find(k1 == k(i)));
  i2 = length(find(k2 == k(i)));
  i3 = length(find(k3 == k(i)));
  m = max([i1, i2, i3]);
  kva = [kva k(i)*ones(1,m-i1)];
  kvb = [kvb k(i)*ones(1,m-i2)];
  kvc = [kvc k(i)*ones(1,m-i3)];
end

r1 = nrbkntins(r1, {kua, kva});
r2 = nrbkntins(r2, {kub, kvb});
t  = nrbkntins(t,  {kuc, kvc});

% combine coefficient to construct Coons surface
coefs(1,:,:) = r1.coefs(1,:,:) + r2.coefs(1,:,:) - t.coefs(1,:,:);
coefs(2,:,:) = r1.coefs(2,:,:) + r2.coefs(2,:,:) - t.coefs(2,:,:);
coefs(3,:,:) = r1.coefs(3,:,:) + r2.coefs(3,:,:) - t.coefs(3,:,:);
coefs(4,:,:) = r1.coefs(4,:,:) + r2.coefs(4,:,:) - t.coefs(4,:,:);
srf = nrbmak(coefs, r1.knots);


text:
function [out]=fcnSensitivityRun_3D(Objfun, BCcoeff, setts, param)
% performes sensitivity corelation check by pairs
% for any function
% Objfun- objective function with one scalar return
% BCcoeff- base case parameter values
% setts -settings e.g. range,step
% param - additional parameters for the Objfun function
% returns X,Y,Z cells for each parameter pair.

% ver. 1.0
% by Levente L. SIMON - ETH Zuerich
% email: l.simon@chem.ethz.ch
% =============================================================

Rangesize = setts.Range;
Step = setts.step;

% ---------------------------------------------------------
Srange = (-1) * Rangesize :Step: Rangesize%  set step
for i =1: size(BCcoeff,2)
    Coeffrange(i,:)  = BCcoeff(i) + (Srange./100 .* BCcoeff(i)) ; % calculate the interval for each coefficient
end

% ---------------------------------------------------------
CoeffComb2 = combnk(1:size(BCcoeff,2), 2); % calculation of combinations based on the parameter number
% sorts matrix from back to front , v 1.1
CoeffComb = CoeffComb2(end:-1:1,:);

% Plot settings
MaxSubPlot = ceil(size(CoeffComb,1) / 3);

for n=1:size(CoeffComb,1)

    % get interval for each data pair
    X = Coeffrange (CoeffComb(n,1), :);
    Y = Coeffrange (CoeffComb(n,2), :);
   clear Z
    for i=1:size(Y,2) % have allways same size, simetric

        for j=1:size(X,2)
            % update vector of BC coefficients
            Scoeff = BCcoeff; Scoeff(CoeffComb(n,1)) = X(j); Scoeff(CoeffComb(n,2)) = Y(i);

            Z(j,i)=feval(Objfun, Scoeff, param);
          end
   end

    % -------- normalize the Z  0..100  --
    Zz = fcnNormalize0_1(Z);

      Xcell{n} = X; Ycell{n} = Y; Zcell{n} = Zz;
      xlabel (strcat ('Parameter No. ', num2str(CoeffComb(n,1))));
      ylabel (strcat ('Parameter No. ', num2str(CoeffComb(n,2))));
      XlabelCell{n} = strcat ('Parameter No. ', num2str(CoeffComb(n,1)));
      YlabelCell{n} = strcat ('Parameter No. ', num2str(CoeffComb(n,2)));

      disp(strcat('Z matrix calculated for pair:',num2str(n)));

end

% Function output structure
out.Xcell = Xcell;
out.Ycell =Ycell;
out.Zcell=Zcell;
out.XlabelCell=XlabelCell;
out.YlabelCell=YlabelCell;
out.MaxSubPlot=MaxSubPlot;
out.setts=setts;


text:
function Y = slapplypca(S, X)
%SLAPPLYPCA Applies PCA model to samples
%
% $ Syntax $
%   Y = slapplypca(S, X)
%
% $ Arguments $
%   - S:        the struct representing the PCA model
%   - X:        the sample matrix
%   - Y:        the principal component vectors of the samples
%
% $ Description $
%   - Y = slapplypca(S, X) applies the PCA model S to reduce the
%     vector dimensions of samples X. It outputs the PCA features by Y.
%     The formula for transform for each sample x in X is:
%     y = S.P' * (x - S.vmean)
%
% $ History $
%   - Created by Dahua Lin on May 1st, 2006
%   - Modified by Dahua Lin, on Sep 10, 2006
%       - change sladd to sladdvec
%

%% parse and verify input arguments
if nargin < 2
    raise_lackinput('slapplypca', 2);
end

if ndims(X) ~= 2
    error('sltoolbox:invaliddims', ...
        'The sample matrix X should be a 2D matrix');
end

if size(X, 1) ~= S.sampledim
    error('sltoolbox:sizmismatch', ...
        'The sample dimension does not match that of PCA model');
end

%% compute
Y = S.P' * sladdvec(X, -S.vmean, 1);

author: lmthang


text:
%EVALEXPRP Evaluates an expression in the grid points.
%
%   [ VEVAL ] = EVALEXPRP( S_EXPR, PROB, SOLNUM, IND ) Evaluates the
%   expression S_EXPR in the points specified in PROB.GRID.P (in
%   global x/y/z-coordinates). PROB is a valid finite element problem
%   struct. Optional parameters are SOLNUM which indicates the chosen
%   solution, and IND for indexing which grid points are to be
%   evaluated for.
%
%       Input       Value/[Size]           Description
%       -----------------------------------------------------------------------------------
%       s_expr      string                 String expression to evaluate
%       prob        struct                 Problem definition struct
%       solnum      scalar {n_sols}        Solution number/time to evaluate
%       ind         vector {1:n_p}         Index to grid points for which to evaluate
%                                                                                .
%       Output      Value/(Size)           Description
%       -----------------------------------------------------------------------------------
%       vEval       [n_ind,1]              Output vector of evaluated values

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
function y = acquire(video_file)

if ischar(video_file),
    display(['Loading file ' video_file]);
    v = VideoReader(video_file);
else
    v = video_file;
end

numFrames = v.NumberOfFrames;

display(['Total frames: ' num2str(numFrames)]);

y = zeros(1, numFrames);
for i=1:numFrames,
    display(['Processing ' num2str(i) '/' num2str(numFrames)]);
    frame = read(v, i);
    redPlane = frame(:, :, 1);
    y(i) = sum(sum(redPlane)) / (size(frame, 1) * size(frame, 2));
end

display('Signal acquired.');
display(' ');
display(['Sampling rate is ' num2str(v.FrameRate) '. You can now run process(your_signal_variable, ' num2str(v.FrameRate) ')']);

end

author: uavster


text:
function a = c8vec_indicator ( n )

%*****************************************************************************80
%
%% C8VEC_INDICATOR sets a C8VEC to the indicator vector.
%
%  Discussion:
%
%    X(1:N) = ( 1-1i, 2-2i, 3-3i, 4-4i, ... )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of elements of A.
%
%    Output, complex A(N), the array to be initialized.
%
  a = zeros ( n, 1 );

  for j = 1 : n
    a(j) = j - j * i;
  end

  return
end

repo: jburkardt-m


text:
function pass = test_isempty()
% Test spherefun isempty() command.

f = spherefun( );
pass(1) = isempty( f );

f = spherefun( @(x,y,z) cos(x) );
pass(2) = ~isempty( f );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% eml_psca_test.m
% compare aspire and matlab versions of E-ML-PSCA
% Copyright Apr 2000, Jeff Fessler, The University of Michigan

%
% generate data
%
if ~isvar('yi'),	disp 'setup eml_psca_test'
	if has_aspire
		f.dir	= test_dir;
		f.wtf	= [f.dir 't,g.wtf'];
		f.yi	= [f.dir 'yi.fld'];
		f.ci	= [f.dir 'ci.fld'];
		f.ri	= [f.dir 'ri.fld'];
		f.mask	= [f.dir 'mask.fld'];
	end
	ig = image_geom('nx', 16, 'ny', 14, 'fov', 500);
	em_test_setup
prompt
end

%
% matlab iterations
%
if ~isvar('xmat'),	disp 'matlab E-ML-PSCA'
	f.niter = 9;

	xinit = ig.ones;
	x = xinit(ig.mask);

	xmat = zeros(sum(ig.mask(:)), f.niter);
	xmat(:,1) = x;
	for ii=2:f.niter
		x = eml_psca(x, G, yi, ci, ri);
		xmat(:,ii) = x;
	end
	xmat = ig.embed(xmat);

	im clf, im(xmat)
prompt
end

if ~has_aspire, return, end

%
% aspire iterations
%
if ~isvar('xasp'),	disp 'aspire E-ML-PSCA'

	f.init	= [f.dir 't,init.fld'];
	fld_write(f.init, xinit, 'check', 0)
	f.out = [f.dir 't,out.fld'];
	if (exist(f.out) == 2), delete(f.out), end

	f.scaleinit = 0;

	f.saver = 'stack,1';
	f.alg = 'psca,1,raster1';
	f.penal = sprintf('%g,quad,0,-', -100);
	f.method = sprintf('@%d@%s@%s', f.niter-1, f.alg, f.penal);
%	f.fitype = ['2z@' f.wtr '@-'];
%	f.com = sprintf(['i -chat 0 empl3 %s %s  %s %s 1 %s 1 %s -' ...
%			' %s %s 1 1 1e30 0 -'], ...
%		f.out, f.init, f.yi, f.ci, f.ri, f.fitype, ...
%		f.method, f.saver);

	f.com = sprintf(['i -chat 0 empl2 %s %s  %s %s %s 1 %s -' ...
			' %s %s 0 1e30 0 -'], ...
		f.out, f.init, f.yi, f.ci, f.ri, f.wtf, ...
		f.method, f.saver);

		os_run(f.com)
		xasp = double(fld_read(f.out));
end

	t = vcorrcoef(xasp, xmat);
	printf('corr. %g,%g', t, t-1)

	t1 = eql_obj(xmat, G, yi(:), ci(:), ri(:), [], ig.mask);
	t2 = eql_obj(xasp, G, yi(:), ci(:), ri(:), [], ig.mask);

if 1 && im
	im clf, im(221, xmat, 'xhat matlab'), cbar
	im(222, xasp, 'xhat aspire'), cbar
	im(223, (xasp-xmat)/max(xmat(:)), 'aspire-matlab'), cbar

	subplot(224)
	plot(0:f.niter-1, t1-t1(1), '-o', 0:f.niter-1, t2-t1(1), '-x')
	legend('mat', 'asp', 4)
prompt
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function varargout = process_delete( varargin )
% PROCESS_DELETE: Delete files, subject, or condition.
%
% USAGE:     sProcess = process_delete('GetDescription')
%                       process_delete('Run', sProcess, sInputs)

% @=============================================================================
% 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.
% =============================================================================@
%
% Authors: Francois Tadel, 2012-2019

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Delete files';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'File';
    sProcess.Index       = 1023;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/SelectFiles?highlight=%28Delete+files%29';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'raw', 'data', 'results', 'timefreq', 'matrix', 'pdata', 'presults', 'ptimefreq', 'pmatrix'};
    sProcess.OutputTypes = {'raw', 'data', 'results', 'timefreq', 'matrix', 'pdata', 'presults', 'ptimefreq', 'pmatrix'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Definition of the options
    % === TARGET
    sProcess.options.target.Comment = {'Delete selected files', 'Delete folders', 'Delete subjects'};  % , 'Delete raw file'
    sProcess.options.target.Type    = 'radio';
    sProcess.options.target.Value   = 1;
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = [sProcess.options.target.Comment{sProcess.options.target.Value}];
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
    % Get options
    DeleteTarget = sProcess.options.target.Value;
    % Returned files: NONE
    OutputFiles = {};

    % Group files in different ways: by subject, by condition, or all together
    switch (DeleteTarget)
        % === DATA FILES ===
        case 1
            % Delete all the files
            file_delete(file_fullpath(unique({sInputs.FileName})), 1);
            % Reload studies
            db_reload_studies(unique([sInputs.iStudy]));

        % === CONDITIONS ===
        case 2
            % Remove all the special studies (intra, inter...)
            specialNames = {bst_get('DirDefaultStudy'), bst_get('DirAnalysisIntra'), bst_get('DirAnalysisInter')};
            iSpecial = find(cellfun(@(c)ismember(c,specialNames), {sInputs.Condition}));
            if ~isempty(iSpecial)
                bst_report('Warning', sProcess, sInputs(iSpecial), 'Some studies cannot be remove (inter-subject, intra-subject, default study).');
                sInputs(iSpecial) = [];
            end
            % Nothing to delete
            if isempty(sInputs)
                return;
            end
            % Delete the studies
            db_delete_studies(unique([sInputs.iStudy]));
            % Update whole tree
            panel_protocols('UpdateTree');

        % === SUBJECTS ===
        case 3
            % Process each subject independently
            uniqueSubj = unique({sInputs.SubjectName});
            % Get subject indices
            iSubjects = [];
            for iSubj = 1:length(uniqueSubj)
                [sSubject, iSubjects(iSubj)] = bst_get('Subject', uniqueSubj{iSubj});
            end
            % Remove default subject
            if any(iSubjects == 0)
                bst_report('Warning', sProcess, sInputs, 'Cannot delete default subject.');
                iSubjects(iSubjects == 0) = [];
            end
            % No subject
            if isempty(iSubjects)
                return
            end
            % Delete subjects
            db_delete_subjects( iSubjects );

%         % === RAW FILES ===
%         case 4
%             % Continuous raw file
%             if strcmpi(sInputs(1).FileType, 'raw')
%                 for i = 1:length(sInputs)
%                     panel_record('DeleteRawFile', sInputs(i).FileName, 1);
%                 end
%             % Other files
%             else
%                 bst_report('Error', sProcess, sInputs, 'The option "Delete raw file" can be used only to remove raw files.');
%             end
    end
end

author: brainstorm-tools


text:
function [Hz] = rpm2Hz(rpm)
% Convert frequency from revolutions per minute to hertz.
% Chad A. Greene 2012
Hz = rpm/60;


text:
function Obj = SOFAresample(Obj, Fs, scale)
%SOFAresample - Modify the sampling rate of an input SOFA object.
%   Usage: Obj = SOFAresample(Obj, Fs)
%          Obj = SOFAresample(Obj, Fs, scale)
%
%   SOFAresample resamples the data in Obj to a specified sampling rate Fs.
%   If Obj is not in SingleFreeFieldHRIR, it will be converted with the help
%   of SOFAconvertConventions.
%
%   Per default, SOFAresample scales the output data to match its energy
%   to the input data. This corresponds to sampling as a process of energy
%   preservation in the continuous function and creates amplitude spectra
%   of the input and output at the same level.
%
%   By calling SOFAresample with scale 'noscale', the output data are not
%   scaled such that the output amplitudes match the input amplitudes. This
%   corresponds to sampling as a process of amplitude snapshots of the
%   continuous function and creates impulse responses of the input and
%   output at the same level.
%
%   Input parameters:
%     Obj   : SOFA object
%     Fs    : Output sampling rate (Hz)
%     scale : 'scale' for scaling the output (default) or 'noscale' otherwise.
%
%   Output parameters:
%     Obj  : New SOFA object (SimpleFreeFieldHRIR convention)

% #Author: Davi R. Carvalho (09.2022)
% #Author: Michael Mihocic: adapt header for consistency (27.09.2022)
% #Author: Piotr Majdak: clean up and parameter adaptations (4.10.2022)

% SOFA Toolbox
% Copyright (C) Acoustics Research Institute - Austrian Academy of Sciences
% Licensed under the EUPL, Version 1.2 or - as soon they will be approved by the European Commission - subsequent versions of the EUPL (the "License")
% You may not use this work except in compliance with the License.
% You may obtain a copy of the License at: https://joinup.ec.europa.eu/software/page/eupl
% Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" basis, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
% See the License for the specific language governing  permissions and limitations under the License.

%% parse inputs
if nargin < 3
    do_scale = true;
elseif nargin == 3
    if strcmpi(scale,'scale')
      do_scale=true;
    else
      do_scale=false;
    end
end
%% load packages if required and prepare conventions
if exist('OCTAVE_VERSION','builtin'); pkg load signal; end
Obj = SOFAconvertConventions(Obj);
%% Resample
  % prepare data
Fs_in = Obj.Data.SamplingRate;
X = Obj.Data.IR;
M = Obj.API.M;
N = Obj.API.N;
R = Obj.API.R;
[p,q] = rat(Fs / Fs_in, 0.0001); % calculate the rational fraction between the input and output sampling rate
  % pre-calculate the resample filter coefficients for speeding up the resample function
normFc = .965 / max(p,q);
order = 256 * max(p,q);
beta = 12;
lpFilt = firls(order, [0 normFc normFc 1],[1 1 0 0]);
lpFilt = lpFilt .* kaiser(order+1,beta)';
lpFilt = lpFilt / sum(lpFilt);
lpFilt = p * lpFilt;
  % Initialize the output matrix
M = size(X, 1);
R = size(X, 2);
N = ceil((Fs/Fs_in) * size(X, 3)); % length after resample
IR=zeros(M, R, N);
  % Resample
for m = 1:M
    for r = 1:R
        IR(m,r,:) = resample(squeeze(X(m,r,:)), p, q, lpFilt);
    end
end
% do scaling
if do_scale, IR = IR.* q/p; end

%% Compile the output structure
Obj.Data.IR = IR;
Obj.Data.SamplingRate = Fs;
Obj = SOFAupdateDimensions(Obj);
end

author: sofacoustics

repo: SOFAtoolbox


text:
function [K, P] = linard2VardistPsi1Compute(linard2kern, vardist, Z)

% LINARD2VARDISTPSI1COMPUTE description.

% VARGPLVM

K = kernCompute(linard2kern,vardist.means,Z);

P = [];

author: SheffieldML

repo: vargplvm


text:
clear all, close all

% Simulate Lorenz system
dt=0.01; T=8; t=0:dt:T;
b=8/3; sig=10; r=28;

Lorenz = @(t,x)([ sig * (x(2) - x(1))       ; ...
                  r * x(1)-x(1) * x(3) - x(2) ; ...
                  x(1) * x(2) - b*x(3)         ]);
ode_options = odeset('RelTol',1e-10, 'AbsTol',1e-11);

input=[]; output=[];
for j=1:100  % training trajectories
    x0=30*(rand(3,1)-0.5);
    [t,y] = ode45(Lorenz,t,x0);
    input=[input; y(1:end-1,:)];
    output=[output; y(2:end,:)];
    plot3(y(:,1),y(:,2),y(:,3)), hold on
    plot3(x0(1),x0(2),x0(3),'ro')
end
grid on, view(-23,18)

%%
net = feedforwardnet([10 10 10]);
net.layers{1}.transferFcn = 'logsig';
net.layers{2}.transferFcn = 'radbas';
net.layers{3}.transferFcn = 'purelin';
net = train(net,input.',output.');

%%
figure(2)
x0=20*(rand(3,1)-0.5);
[t,y] = ode45(Lorenz,t,x0);
plot3(y(:,1),y(:,2),y(:,3)), hold on
plot3(x0(1),x0(2),x0(3),'ro','Linewidth',[2])
grid on

ynn(1,:)=x0;
for jj=2:length(t)
    y0=net(x0);
    ynn(jj,:)=y0.'; x0=y0;
end
plot3(ynn(:,1),ynn(:,2),ynn(:,3),':','Linewidth',[2])

figure(3)
subplot(3,2,1), plot(t,y(:,1),t,ynn(:,1),'Linewidth',[2])
subplot(3,2,3), plot(t,y(:,2),t,ynn(:,2),'Linewidth',[2])
subplot(3,2,5), plot(t,y(:,3),t,ynn(:,3),'Linewidth',[2])

figure(2)
x0=20*(rand(3,1)-0.5);
[t,y] = ode45(Lorenz,t,x0);
plot3(y(:,1),y(:,2),y(:,3)), hold on
plot3(x0(1),x0(2),x0(3),'ro','Linewidth',[2])
grid on

ynn(1,:)=x0;
for jj=2:length(t)
    y0=net(x0);
    ynn(jj,:)=y0.'; x0=y0;
end
plot3(ynn(:,1),ynn(:,2),ynn(:,3),':','Linewidth',[2])

figure(3)
subplot(3,2,2), plot(t,y(:,1),t,ynn(:,1),'Linewidth',[2])
subplot(3,2,4), plot(t,y(:,2),t,ynn(:,2),'Linewidth',[2])
subplot(3,2,6), plot(t,y(:,3),t,ynn(:,3),'Linewidth',[2])

%%
figure(2), view(-75,15)
figure(3)
subplot(3,2,1), set(gca,'Fontsize',[15],'Xlim',[0 8])
subplot(3,2,2), set(gca,'Fontsize',[15],'Xlim',[0 8])
subplot(3,2,3), set(gca,'Fontsize',[15],'Xlim',[0 8])
subplot(3,2,4), set(gca,'Fontsize',[15],'Xlim',[0 8])
subplot(3,2,5), set(gca,'Fontsize',[15],'Xlim',[0 8])
subplot(3,2,6), set(gca,'Fontsize',[15],'Xlim',[0 8])
legend('Lorenz','NN')

author: dynamicslab


text:
function [varargout]=triSurfSplitBoundary(varargin)

% function [F,V,E,CF,CV]=triSurfSplitBoundary(F,V,E,n,CF,CV)
% ------------------------------------------------------------------------
%
%
% To do: Use for loop and recursion instead
% ------------------------------------------------------------------------

%% Parse input

switch nargin
    case 3
        F=varargin{1};
        V=varargin{2};
        E=varargin{3};
        n=size(E,1)+1;
        CF=[];
        CV=[];
    case 4
        F=varargin{1};
        V=varargin{2};
        E=varargin{3};
        n=varargin{4};
        CF=[];
        CV=[];
    case 5
        F=varargin{1};
        V=varargin{2};
        E=varargin{3};
        n=varargin{4};
        CF=varargin{5};
        CV=[];
    case 6
        F=varargin{1};
        V=varargin{2};
        E=varargin{3};
        n=varargin{4};
        CF=varargin{5};
        CV=varargin{6};
end

if isempty(CF)
    CF=ones(size(F,1),1);
end

%%
if isempty(E)
    error('Edge set is empty');
end

if size(E,1)==n
    %Do nothing: Input = output
elseif size(E,1)>n
    error('Current number of edges exceeds number of desired edges');
else
    while size(E,1)~=n
        [~,indMax]=max(edgeLengths(E,V));
        numNodesPre=size(V,1); %Number of nodes before split
        [F,V,~,CF,CV]=triEdgeSplit(F,V,E(indMax,:),CF,CV);
        numNodesPost=size(V,1); %Number of nodes after split
        [Eb]=patchBoundaryLabelEdges(F,V,CF); %All boundary edges
        logicKeep=all(ismember(Eb,E) |...
            ismember(Eb,(numNodesPre+1):numNodesPost),2);
        E=Eb(logicKeep,:);
    end
end

%% Collect output
varargout{1}=F;
varargout{2}=V;
varargout{3}=E;
varargout{4}=CF;
varargout{5}=CV;

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [fProbability, fMc] = calc_McGrid(mCatalog, fBinning)
% [fProbability, fMc] = function calc_McGrid(mCatalog, fBinning);
% --------------------------------------------
% Determine Mc using maximum likelihood score fitting the
% non-cumulative frequency magnitude distribution above and below Mc:
% below: Exponetial function
% above: Gutenberg-Richter law
% !! Data not mormalized in time !!!
%
% Incoming variables:
% mCatalog   : EQ catalog
% fBinning   : Binning interval, usually 0.1
%
% Outgoing variables:
%
% J. Woessner: woessner@seismo.ifg.ethz.ch
% last update: 05.11.02

% Initialize
vProbability = [];
vMc = [];
vABValue =[];
vX_res = [];
vNCumTmp = [];

% Determine exact time period
fPeriod1 = max(mCatalog(:,3)) - min(mCatalog(:,3));

% Determine max. and min. magnitude
fMinMag = floor(min(mCatalog(:,6)));
if fMinMag > 0
  fMinMag = 0;
end
fMaxMag = ceil(10 * max(mCatalog(:,6))) / 10;

% Calculate FMD for original catalog
[vFMD, vNonCFMD] = calc_FMD(mCatalog);
% % Time normalization
% vFMD(2,:) = ceil(vFMD(2,:)./fPeriod1);
% vNonCFMD(2,:) = ceil(vNonCFMD(2,:)./fPeriod1);

%fMc = calc_Mc(mCatalog, nCalculateMC)
vNonCFMD = fliplr(vNonCFMD);
for fMc = 0.8:0.1:2
    [nIndexLo, fMagHi, vSel, vMagnitudes] = fMagToFitBValue(mCatalog, vFMD, fMc);
    [fMeanMag, fBValue, fStdDev, fAValue] =  calc_bmemag(mCatalog(vSel,:), fBinning);
    % Compute quantity of earthquakes by power law
    vMstep = [fMinMag:0.1:fMaxMag];
    vNCum = 10.^(fAValue-fBValue.*vMstep); % Cumulative number
    % Compute non-cumulative numbers vN
    fNCumTmp = 10^(fAValue-fBValue*(fMaxMag+0.1));
    vNCumTmp  = [vNCum fNCumTmp];
    %vNonCFMD(2,:) = ceil(vNonCFMD(2,:)./fPeriod1);
    vN = abs(diff(vNCumTmp));
    %vN = vN./fPeriod1;
    mData = [vN' vNonCFMD'];
    vSel = (mData(:,2) >= fMc);
    mDataTest = mData(~vSel,:);
    mDataTmp = mData.subset(vSel);
    % Curve fitting: Non cumulative part below Mc
    options = optimset;
    options = optimset('Display','off','Tolfun',1e-6,'TolX',0.0001);
    [vX, resnorm, resid, exitflag]=lsqcurvefit(@calc_expdecay2,[0 1], mDataTest(:,2), mDataTest(:,3));
    mDataTest(:,1) = vX(1).*exp(vX(2).*mDataTest(:,2))-1;
    vX_res = [vX_res; vX resnorm exitflag];
    %% Set data together
    mDataTest = [mDataTest; mDataTmp];
    vProb_ = calc_log10poisspdf(vNonCFMD(2,:)', ceil(mDataTest(:,1)));
% %     vProb_ = calc_log10poisspdf(mData(vSel,2), mData(vSel,1));
% %     vProb_ = calc_log10poisspdf(vNonCFMD(2,:)', vN');
%     figure_w_normalized_uicontrolunits(300)
%     semilogy(vNonCFMD(1,:)', vNonCFMD(2,:)', '^', vNonCFMD(1,:)', vN, '*', vNonCFMD(1,:)',mDataTest(:,1),'o')
    % Sum the probabilities
    fProbability = (-1) * sum(vProb_);
    vProbability = [vProbability; fProbability];
    vMc = [vMc; fMc];
    vABValue = [vABValue; fAValue fBValue];
    % Clear variables
    vNCumTmp = [];
    mModelDat = [];
    vNCum = [];
    vSel = [];
    mDataTest = [];
end
mResult = [vProbability vMc vABValue vX_res];
vSel = (mResult == min(mResult(:,1)));
mResult = mResult(vSel,:);
if length(mResult(:,1)) > 1
    fProbability = max(mResult(:,1));
    fMc = max(mResult(:,2));
else
    fProbability = mResult(1,1);
    fMc = mResult(1,2);
end

% figure_w_normalized_uicontrolunits(400)
% plot(vMc, vProbability,'*');

author: CelsoReyes


text:
function [ results ] = run_DAT_RANSAC( seq, res_path, bSaveImage )

project_dir = fileparts(fileparts(fileparts(mfilename('fullpath'))));
addpath(fullfile(project_dir, 'trackers', 'DAT', 'src'));

visualization = 1;
use_normal_size = 1;
normal_height = 360;
normal_width = 640;
I = imread(seq.s_frames{1});
if use_normal_size
  scale_height = size(I, 1) / normal_height;
  scale_width = size(I, 2) / normal_width;
else
  scale_height = 1;
  scale_width = 1;
end
% Load annotations
groundtruth = seq.init_rect;
groundtruth = groundtruth ./ [scale_width, scale_height, scale_width, scale_height];

% Load default settings
cfg = ransac_parameters_dat();

% Tracking
frames = 1:length(seq.s_frames);
times = zeros(size(frames));
rects = zeros(length(frames), 4);
spf = tic;
for frame = frames
  I = imread(seq.s_frames{frame});
  % I = imresize(I, [360, 640]);
  if use_normal_size
    I = imresize(I, [normal_height, normal_width]);
  end

  ttrack = tic;
  if frame == 1
    [state, location] = tracker_dat_initialize(I, groundtruth, cfg);
  else
    [state, location] = tracker_dat_update(state, I, cfg);
  end
  times(frame) = toc(ttrack);

  % Visualization
  % rects(frame, :) = location*2;
  rects(frame, :) = location .* [scale_width, scale_height, scale_width, scale_height];
  if visualization
    figure(1), clf
    imshow(I)
    rectangle('Position', location, 'EdgeColor', 'b', 'LineWidth', 2);
    drawnow
  end
end
spf = toc(spf);

results.res = rects;
results.type='rect';
results.fps = length(frames)/spf;

end

author: flyers

repo: drone-tracking


text:
%--- help for generic/setup_nonlinear_restrictions ---
%
%  INTERNAL FUNCTION: Sets nonlinear restrictions
%
%  Note:
%
%     - uses estim_nonlinear_restrictions, which should be a cell array. Each
%     item of the array is a string of the form
%       - 'f(p1,p2,...,pn)>=h(p1,p2,...,pn)'
%       - 'f(p1,p2,...,pn)>h(p1,p2,...,pn)'
%       - 'f(p1,p2,...,pn)<=h(p1,p2,...,pn)'
%       - 'f(p1,p2,...,pn)<h(p1,p2,...,pn)'
%       - 'pj=h(p1,p2,...,pn)'
%
%     - In the statements above,
%       - eqtn [digits|variable name]
%       - vbl [digits|variable name]
%       - lag [digits]
%       - chain [char]
%       - state [digits]
%
%
%    Other functions named setup_nonlinear_restrictions
%
%       abstvar/setup_nonlinear_restrictions
%

author: jmaih


text:
function m=moda(u,N)
% MODA   Mode of a distribution
%
%    m=MODE(u,N) calculates the mode of the set of data "u" using the histogram.
%    To avoid outliers, for the calculation are only taken into account those
%    values in [mean-2sigma, mean+2sigma];
%
%    INPUT:
%
%	- u (set of data)
%       - N: Number of points for the histogram. If N=0 then 5000 points are
%            considered
%
%   Author: Santiago Aja Fernandez
%   LOCAL STATISTICS TOOLBOX
%
%   Modified: Feb 01 2008
%
u=double(u);
if N==0
	N=5000;
end
M1=mean(u(:));
V1=std(u(:));
C2=u( (u(:)>=(M1-2*V1)) & (u(:)<=(M1+2*V1)) ) ;
%C2=u;
[h,x]=hist(C2,N);
[M,M2]=max(h);
m=x(M2);
%

author: Sable

repo: mcbench-benchmarks