text:
function H=chi2test(x, alpha)
%   CHI2TEST: Single sample Pearson Chi Square goodness-of-fit hypothesis test.
%   H=CHI2TEST(X,ALPHA) performs the particular case of Pearson Chi Square
%   test to determine whether the null hypothesis of composite normality PDF is
%   a reasonable assumption regarding the population distribution of a random sample X
%   with the desired significance level ALPHA.
%
%   H indicates the result of the hypothesis test according to the MATLAB rules
%   of conditional statements:
%   H=1 => Do not reject the null hypothesis at significance level ALPHA.
%   H=0 => Reject the null hypothesis at significance level ALPHA.
%
%   The Chi Square hypotheses and test statistic in this particular case are:
%
%   Null Hypothesis:        X is normal with unknown mean and variance.
%   Alternative Hypothesis: X is not normal.
%
%   The random sample X is shifted by its estimated mean and normalized by its
%   estimated standard deviation. The tested bins XP of the assumed normal distribution
%   are chosen [-inf, -1.6:0.4:1.6, inf] to avoid unsufficient statistics.
%
%   Let E(x) be the expected frequency X falls within XP according to the normal
%   distribution and O(x) be the observed frequency. The Pearson statistic,
%   X2=SUM((E(x)-O(x))^2/E(x)) distributes Chi Square with length(XP)-3 degrees
%   of freedom.
%
%   The decision to reject the null hypothesis is taken when the P value (probability that Chi2
%   random value with length(XP)-3 degrees of freedomd is greater than X2) is less than
%   significance level ALPHA.
%
%   X must be a row vector representing a random sample. ALPHA must be a scalar.
%   The function doesn't check the formats of X and ALPHA, as well as a number of the
%   input and output parameters.
%
%   The asymptotic limit of the Chi Square test presented is reached when LENGTH(X)>90.
%
% Acknowledge: Dr. S. Loyka
%
% Author: G. Levin, May, 2003.
%
% References:
%   W. T. Eadie, D. Drijard, F. E. James, M Roos and B. Sadoulet, "Statistical Methods
%   in Experimental Physics", North-Holland, Sec. Reprint, 1982.

%Normalize x
N=length(x);
x=(x-mean(x))/std(x); %standardization

xp=[-inf, -1.6:.4:1.6, inf]; %tested bins
E=0.5*N*diff(erfc(-xp/sqrt(2))); %expected frequency
S=histc(x, xp);
O=S(1:end-1); %%observed frequency
%plot(xp(2:end),E,'k-',xp(2:end),O,'k.');
x2=sum((E-O).^2./E); %statistics

pval=1-gammainc(x2/2,(length(O)-3)/2); %p value

H=(pval>=alpha);


text:
%====================================================
%Practise problem for the course of CAO
%Prof P.Beckers
%====================================================
%Porpose :Drawing the orthotomic surface of a Bezier surface  %
% Student: BUI QUOC TINH
% European Master Mechanices of Contructions	(EMMC)
% University the Liege. Belgium										    		     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%====================================================

%The function calculation BernsTein polynomial derivation
function Bepoly=B_daoham(i,n,u)
% Calculation combination
if i==n|i==0
   c=1;
elseif i<n & i>=0
   c=factorial(n)/(factorial(i)*factorial(n-i));
else
   c=0;
end
%BernsTein polynomial
Bepoly=c*u^i*(1-u)^(n-i)*(i-n*u)/((1-u)*u);

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

author: Sable

repo: mcbench-benchmarks


text:
function imsegs = APPgetSpStats(imsegs)
% imsegs = APPgetSpStats(imsegs)
% Gets basic information about the superpixels
%
% Copyright(C) Derek Hoiem, Carnegie Mellon University, 2005
% Current Version: 1.0  09/30/2005

for ii = 1:length(imsegs)

	nseg = imsegs(ii).nseg;
	segimage = double( imsegs(ii).segimage );
%     segimage = imsegs(ii).segimage;

    imh = size(segimage, 1);

	adjmat = eye([nseg nseg]);

    % get adjacency
    dx = segimage ~= segimage(:,[2:end end]);
    dy = segimage ~= segimage([2:end end], :);

    ind1 = find(dy);
    ind2 = ind1 + 1;
    s1 = segimage(ind1);
    s2 = segimage(ind2);
%     adjmat(s1 + nseg*(s2-1)) = 1;
%     adjmat(s2 + nseg*(s1-1)) = 1;
    adjmat(sub2ind([nseg, nseg], s1, s2)) = 1;
    adjmat(sub2ind([nseg, nseg], s2, s1)) = 1;

    ind3 = find(dx);
    ind4 = ind3 + imh;
    s3 = segimage(ind3);
    s4 = segimage(ind4);
%     adjmat(s3 + nseg*(s4-1)) = 1;
%     adjmat(s4 + nseg*(s3-1)) = 1;
    adjmat(sub2ind([nseg, nseg], s3, s4)) = 1;
    adjmat(sub2ind([nseg, nseg], s4, s3)) = 1;

%   slower code
% 	[height, width] = size(segimage);
%
% 	for y = 1:height-1
%         for x = 1:width-1
%             s1 = segimage(y, x);
%             s2 = segimage(y+1, x);
%             s3 = segimage(y, x+1);
%             if s1 > 0
%                 npixels(s1) = npixels(s1) + 1;
%                 if s2 > 0
%                     adjmat(s1, s2) = 1;
%                     adjmat(s2, s1) = 1;
%                 end
%                 if s3 > 0
%                     adjmat(s1, s3) = 1;
%                     adjmat(s3, s1) = 1;
%                 end
%             end
%         end
% 	end
%
% 	x = width;
% 	for y = 1:height
%         s1 = segimage(y, x);
%         if s1 > 0
%             npixels(s1) = npixels(s1) + 1;
%         end
% 	end
%
% 	y = height;
% 	for x = 1:width-1
%         s1 = segimage(y, x);
%         if s1 > 0
%             npixels(s1) = npixels(s1) + 1;
%         end
% 	end

    stats = regionprops(segimage, 'Area');
    imsegs(ii).npixels = vertcat(stats(:).Area);
	imsegs(ii).adjmat = logical(adjmat);

end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
%ELIIPSOID  Plots the upper half of a parametric ellipsoid of rotation
%   Figure 13.10.For explanations see Subsection 13.26 in the book.

% define axes
a = 5;
b = 3;
c = 2;
% plot the boundary of the patch and define frame
w = 0: 0.01: 1;
x = a*cos(2*pi*w);
y = b*sin(2*pi*w);
z = zeros(size(x));
plot3(x, y, z, 'k-')
axis off, axis equal
axis([ -1.1*a 1.1*a -1.1*b 1.1*b -0.3 1.2*c ])
t = [ 'Ellipsoid x^2/a^2 + y^2/b^2 + z^2/c^2 = 1, a = ' num2str(a) ];
t = [ t ', b = ' num2str(b) ' c = ' num2str(c) ]
Ht = title(t);
set(Ht, 'FontSize', 14)
hold on
% plot now family of constant-w curves
u = 0: 0.01: 1;
for w = 0: 0.1: 1
    x = a*cos(pi*u).*cos(2*pi*w);
    y = b*cos(pi*u).*sin(2*pi*w);
    z = c*sin(pi*u);
    plot3(x, y, z, 'k-')
end
for w1 = 0.5: 0.1: 0.9
        t = [ 'w = ' num2str(w1) ];
        xt = a*cos(2*pi*w1);
        yt = b*sin(2*pi*w1);
        zt = -0.3;
        Ht = text(xt, yt, zt, t);
        set(Ht, 'FontSIze', 12)
end

% plot now family of constant-u curves
w2 = 0: 0.01: 1;
for u2 = 0: 0.1: 1
    x = a*cos(pi/2*u2).*cos(2*pi*w2);
    y = b*cos(pi/2*u2).*sin(2*pi*w2);
    z = c*sin(pi/2*u2).*ones(size(x));
    plot3(x, y, z, 'k-')
end
w3 = 0.6;
for u3   = 0: 0.1: 1
    x3   = a*cos(pi/2*u3).*cos(2*pi*w3);
    y3   = b*cos(pi/2*u3).*sin(2*pi*w3);
    z3   = c*sin(pi/2*u3).*ones(size(x3));
    t3  = [ 'u = ' num2str(u3) ];
    % xt1 = a*cos(2*pi*w3);
    % yt1 = b*sin(2*pi*w);
    % zt1 = c*sin(pi/2*u);
    Ht = text(x3, y3, z3, t3);
    set(Ht, 'FontSIze', 12)
end
hold off


text:
function mappedModel = mapModelMets(model,mnx)
%mapModelMets  Retrieve and assign standard IDs to model metabolites.
%
% USAGE:
%
%   mappedModel = mapModelMets(model,mnx);
%
% INPUTS:
%
%   model   A genome scale model structure, containing metabolite related
%           fields (e.g., mets, metNames, etc.).
%
%   mnx     An MNX metabolite database structure generated using the
%           buildMNXmodel('met') function.
%
% OUTPUS:
%
%   mappedModel   The model returned with additional metabolite-related
%                 fields, associating various metabolite IDs to MNX IDs, as
%                 well as a metMNXID field which combines the MNX IDs
%                 obtained by mapping metabolites along each of its
%                 available ID fields.
%                 *NOTE: Some of the met-related fields in mappedModel may
%                        contain multiple columns, which contain the
%                        multiple IDs that matched to one or more of the
%                        metabolites.
%

% handle input arguments
if nargin < 2
    mnx = [];
end

% get list of metID fields
ignoreFields = {'metFormulas','metMiriams','metComps','mets', ...
                'metNamesAlt','metCharges','metSMILES','metPdMap'};
metIDfields = fields(model);
metIDfields(~startsWith(metIDfields,'met') | ismember(lower(metIDfields),lower(ignoreFields))) = [];

% load metabolite information from MNX database file
if isempty(mnx)
    mnx = buildMNXmodel('met');
end

% associate each set of IDs to MNX IDs
fprintf('Mapping metabolite external IDs to MNX IDs:\n');
for i = 1:length(metIDfields)

    % skip field if it isn't present in the MNX database model structure
    if ~isfield(mnx,metIDfields{i})
        continue
    end

    fprintf('\t%s\n',metIDfields{i});
    if strcmp(metIDfields{i},'metNames')  % the 'metNames' field is handled differently than others

        % combine names and alternative names into single cell array
        if isfield(model,'metNamesAlt')
            metNames = [model.metNames,model.metNamesAlt];
        else
            metNames = model.metNames;
        end

        % ignore case
        metNames = lower(metNames);
        mnx.mnxID2name(:,2) = lower(mnx.mnxID2name(:,2));

        % Perform the name-matching process twice. The first pass will
        % search for exact name matches (ignoring case), whereas the second
        % pass will loosen the criteria by removing all special characters
        % (e.g., hyphens, parentheses, spaces, etc.) from the met names,
        % and search again for any mets that were not matched during the
        % first pass.
        for ii = 1:2

            if ii == 2
                % remove special characters from metabolite names
                metNames = regexprep(metNames,'[^a-zA-Z0-9]','');
                mnx.mnxID2name(:,2) = regexprep(mnx.mnxID2name(:,2),'[^a-zA-Z0-9]','');
            else
                model.metName2MNX = repmat({''},size(model.mets));
            end

            % extract subset of MNXID-name pairs containing matching names (for faster processing)
            keep_ind = ismember(mnx.mnxID2name(:,2),metNames);
            mnx_ids = mnx.mnxID2name(keep_ind,1);
            mnx_names = mnx.mnxID2name(keep_ind,2);

            % convert metNames from cell array to column vector of nested cells (for next processing step)
            metNamesNest = nestCell(metNames,true);
            ignore_inds = ~cellfun(@isempty,model.metName2MNX) | cellfun(@isempty,metNamesNest);

            % retrieve all matching IDs for each met
            model.metName2MNX(~ignore_inds) = cellfun(@(x) mnx_ids(ismember(mnx_names,x)),metNamesNest(~ignore_inds),'UniformOutput',false);

        end

        % remove entries that have matched to too many IDs (>100)
        model.metName2MNX(cellfun(@numel,model.metName2MNX) > 100) = {''};

        % flatten cell array
        model.metName2MNX = flattenCell(model.metName2MNX,true);

    else
        % get model met IDs, and compress each row into nested cells
        % (this is to deal with fields that have multiple columns)
        model_ids = nestCell(model.(metIDfields{i}),true);

        % extract only subset of MNX model containing matching IDs (for faster processing)
        keep_ind = ismember(mnx.mnxID2extID(:,2),metIDfields(i)) & ismember(mnx.mnxID2extID(:,3),model.(metIDfields{i}));
        ext_ids = mnx.mnxID2extID(keep_ind,3);
        mnx_ids = mnx.mnxID2extID(keep_ind,1);

        % get empty indices for current field
        empty_inds = cellfun(@isempty,model_ids);

        % retrieve all matching IDs for each met
        newField = strcat(metIDfields{i},'2MNX');
        model.(newField) = repmat({''},size(model.mets));
        model.(newField)(~empty_inds) = cellfun(@(x) mnx_ids(ismember(ext_ids,x)),model_ids(~empty_inds),'UniformOutput',false);

        % flatten cell array
        model.(newField) = flattenCell(model.(newField),true);
    end

end
fprintf('Done.\n');

% now combine all the MNX IDs for each metabolite
mnxIDfields = fields(model);
mnxIDfields(~(startsWith(mnxIDfields,'met') & endsWith(mnxIDfields,'2MNX'))) = [];
metMNXIDs = {};  % intialize cell array of met MNX IDs
for i = 1:length(mnxIDfields)
    % append each field as new column(s)
    metMNXIDs = [metMNXIDs,model.(mnxIDfields{i})];
end
empty_inds = cellfun(@isempty,metMNXIDs);
% metMNXIDs(empty_inds) = {''};

% obtain unique set of MNX IDs for each metabolite
met_index = transpose(1:length(model.mets));
model.metMNXID = arrayfun(@(i) unique(metMNXIDs(i,~empty_inds(i,:))),met_index,'UniformOutput',false);
model.metMNXID = flattenCell(model.metMNXID,true);

% assign output
mappedModel = model;

end

author: SysBioChalmers

repo: Human-GEM


text:
function eventLog = ma_executeBasicMassDump(massDumpEvent, initialState, eventNum)
%ma_executeMassDump Summary of this function goes here
%   Detailed explanation goes here

    ut = initialState(1);
    rVectUT = initialState(2:4);
    vVectUT = initialState(5:7);
    bodyID = initialState(8);

    dmasses = massDumpEvent.dumpValue;
    masses = initialState(9:12) - dmasses;

    eventLog(1,:) = [ut, rVectUT, vVectUT, bodyID, masses, eventNum];
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function sourceLoc = dba_anatmodel(ind, sScout, CortexMat, srcType)
%
% Compute sources locations for deep structures
%
% Yohan Attal - HM-TC project 2013

disp = 0;

switch srcType
    case 'surf'
        sourceLoc = CortexMat.Vertices(ind,:);

    case 'vol'
        [ind_s, I_ind] = sort(ind);
        vertTmp = CortexMat.Vertices(sScout.Vertices,:);
        tessTmp.Vertices = vertTmp(I_ind,:);
        iFaces = sum(ismember(CortexMat.Faces, ind_s),2)==3;
        tessTmp.Faces = CortexMat.Faces(iFaces,:);
        for iv=1:numel(tessTmp.Faces)
            tessTmp.Faces(iv) = find(ind_s==tessTmp.Faces(iv));
        end
        sourceLoc = dba_vol_grids(tessTmp);
end
if disp && isequal(srcType,'vol')
    figure('color','w'); hold on; pp = patch(tessTmp);
    set(pp, 'FaceColor', [0 0 1] , 'FaceAlpha', .5, 'EdgeAlpha', 0)
    plot3(sourceLoc(:,1), sourceLoc(:,2),sourceLoc(:,3),'ro')
end
end

author: brainstorm-tools


text:
function [ pn, dist ] = shape_point_near_2d ( center, p1, nside, p, pn, dist )

%*****************************************************************************80
%
%% SHAPE_POINT_NEAR_2D: nearest point ( regular shape, point ) in 2D.
%
%  Discussion:
%
%    The "regular shape" is assumed to be an equilateral and equiangular
%    polygon, such as the standard square, pentagon, hexagon, and so on.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real CENTER(2,1), the center of the shape.
%
%    Input, real P1(2,1), the first vertex of the shape.
%
%    Input, integer NSIDE, the number of sides in the shape.
%
%    Input, real P(2,1), the point to be checked.
%
%    Output, real PN(2,1), the point on the shape that is nearest
%    to the given point.
%
%    Output, real DIST, the distance between the points.
%

%
%  Determine the angle subtended by a single side.
%
  sector_angle = 360.0 / nside;
%
%  How long is the half-diagonal?
%
  radius = sqrt ( sum ( ( p1(1:2,1) - center(1:2,1) ).^2 ) );
%
%  If the radius is zero, then the shape is a point and the computation is easy.
%
  if ( radius == 0.0 )
    pn(1:2,1) = center(1:2,1);
    dist = sqrt ( sum ( ( p(1:2,1) - pn(1:2,1) ).^2 ) );
    return
  end
%
%  If the test point is at the center, then the computation is easy.
%  The angle subtended by any side is ( 2 * PI / NSIDE ) and the
%  nearest distance is the midpoint of any such side.
%
  if ( p(1:2,1) == center(1:2,1) )
    angle = pi / nside;
    pd(1,1) =   ( p(1,1) - center(1,1) ) * cos ( angle ) ...
              + ( p(2,1) - center(2,1) ) * sin ( angle );
    pd(2,1) = - ( p(1,1) - center(1,1) ) * sin ( angle ) ...
              + ( p(2,1) - center(2,1) ) * cos ( angle );
    pn(1,1) = center(1,1) + pd(1,1) * cos ( angle );
    pn(2,1) = center(2,1) + pd(2,1) * sin ( angle );
    dist = radius * cos ( angle );
    return
  end
%
%  Determine the angle between the ray to the first corner,
%  and the ray to the test point.
%
  angle = angle_deg_2d ( p1(1:2,1), center(1:2,1), p(1:2,1) );
%
%  Determine the sector of the point.
%
  sector_index = floor ( angle / sector_angle ) + 1;
%
%  Generate the two corner points that terminate the SECTOR-th side.
%
  angle2 = ( sector_index - 1 ) * sector_angle;
  angle2 = degrees_to_radians ( angle2 );

  pa = vector_rotate_base_2d ( p1, center, angle2 );

  angle2 = ( sector_index ) * sector_angle;
  angle2 = degrees_to_radians ( angle2 );

  pb = vector_rotate_base_2d ( p1, center, angle2 );
%
%  Determine the point on the SECTOR-th side of the shape which is
%  nearest.
%
  [ pn, dist, t ] = segment_point_near_2d ( pa, pb, p );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% L2Norm is the dagnn wrapper which normalizes the features to be norm 1 in
% 2-norm at each location

% Copyright (C) 2015 Tsung-Yu Lin, Aruni RoyChowdhury, Subhransu Maji.
% All rights reserved.
%
% This file is part of the BCNN. It uses MATCONVNET package and is made
% available under the terms of the BSD license (see the COPYING file).

classdef L2Norm < dagnn.ElementWise
  properties
    param = 1e-10;
  end

  methods
    function outputs = forward(obj, inputs, params)
      outputs{1} = vl_nnl2norm(inputs{1}, obj.param) ;
    end

    function [derInputs, derParams] = backward(obj, inputs, params, derOutputs)
      derInputs{1} = vl_nnl2norm(inputs{1}, obj.param, derOutputs{1}) ;
      derParams = {} ;
    end

    function rfs = getReceptiveFields(obj)
      rfs.size = [1 1] ;
      rfs.stride = [1 1] ;
      rfs.offset = [1 1] ;
    end

    function obj = SpatialNorm(varargin)
      obj.load(varargin) ;
    end
  end
end

repo: DBCNN


text:
function [nn,ctrs] = myhist3(varargin)
% almost the same as Matlab's hist3, except it can take weights.

i = find(strcmpi(varargin,'weights'));
if ~isempty(i),
  weights = varargin{i(end)+1};
  varargin([i,i+1]) = [];
end

[cax,args,nargs] = axescheck(varargin{:});

if nargs < 1
    error('stats:hist3:TooFewInputs', 'Requires X.')
end
x = args{1};

% See if nbins/ctrs was given as the second argument, or only name/value
% pairs.
if nargs > 1 && ~ischar(args{2})
    binSpec = args{2};
    args = args(3:end); % strip off x and nbins/ctrs
else
    binSpec = [];
    args = args(2:end); % strip off x
end

% Process input parameter name/value pairs, assume unrecognized ones are
% graphics properties.
pnames = {'nbins','ctrs','edges'};
dflts =  { [],     [],       []};
[errid,errmsg,nbins,ctrs,edges,plotArgs] = statgetargs(pnames, dflts, args{:});
if ~isempty(errmsg)
    error(['stats:hist3:' errid], errmsg);
end

% Make sure they haven't mixed 'nbins'/'ctrs'/'edges' name/value pairs with
% the CTRS or NBINS unnamed second arg syntax, or used more than one of
% those parameter name/value pairs.
if (isempty(nbins)+isempty(ctrs)+isempty(edges)+isempty(binSpec)) < 3
    error('stats:hist3:AmbiguousBinSpec', 'Ambiguous specification of bins.');
elseif ~isempty(binSpec)
    if iscell(binSpec)  % hist3(x,ctrs)
        ctrs = binSpec;
    else                % hist3(x,nbins)
        nbins = binSpec;
    end
end

if ~isempty(nbins)
    % Use the specified number of bars in each direction, centers and edges
    % to be determined.
    histBehavior = true;
    if ~(isnumeric(nbins) && numel(nbins)==2)
        error('stats:hist3:BadNbins', ...
              'The number of bins must be specified with a 2-element numeric vector.');
    end
    autobins = true;

elseif ~isempty(ctrs)
    % Use the specified bin centers.
    histBehavior = true;
    if ~(iscell(ctrs) && numel(ctrs)==2 && isnumeric(ctrs{1}) && isnumeric(ctrs{2}))
        error('stats:hist3:BadCtrs', ...
              'Bin centers must be specified with a cell array containing two numeric vectors.');
    end
    ctrs = {ctrs{1}(:)' ctrs{2}(:)'};
    autobins = false;
    nbins = [length(ctrs{1}) length(ctrs{2})];

elseif ~isempty(edges)
    % Use the specified bin edges.
    histBehavior = false;
    if ~(iscell(edges) && numel(edges)==2 && isnumeric(edges{1}) && isnumeric(edges{2}))
        error('stats:hist3:BadEdges', ...
              'Bin edges must be specified with a cell array containing two numeric vectors.');
    end
    edges = {edges{1}(:)' edges{2}(:)'};
    autobins = false;
    % Just as with histc, there will be #edges bins
    nbins = [length(edges{1}) length(edges{2})];

else
    % Assume a 10x10 grid of bars, centers and edges to be determined.
    histBehavior = true;
    autobins = true;
    nbins = [10 10];
end

[nrows,ncols] = size(x);
if ncols ~= 2
    error('stats:hist3:WrongNumCols', 'X must be a matrix with two columns.');
end

% Special case for empty data (follows what HIST does).
if isempty(x)
    if autobins
       ctrs = {1:nbins(1) 1:nbins(2)};
    end
    n = zeros(nbins); % Nothing to count, return nbins(1) by nbins(2) zeros

else
    % Bin each observation in the x-direction, and in the y-direction.
    bin = zeros(nrows,2);
    for i = 1:2
        minx = min(x(:,i));
        maxx = max(x(:,i));

        % If only the number of bins was given, compute edges and centers
        % for equal-sized bins spanning the data.
        if autobins
            if isinf(minx) || isinf(maxx)
                error('stats:hist3:InfData', ...
                      'Bin centers or edges must be specified when data contain infinite values.');
            elseif minx == maxx
                minx = minx - floor(nbins(i)/2) - 0.5;
                maxx = maxx + ceil(nbins(i)/2) - 0.5;
            end
            binwidth{i} = (maxx - minx) / nbins(i);
            edges{i} = minx + binwidth{i}*(0:nbins(i));
            ctrs{i} = edges{i}(1:nbins(i)) + binwidth{i}/2;
            % Make histc mimic hist behavior:  everything < ctrs(1) gets
            % counted in first bin, everything > ctrs(end) gets counted in
            % last bin.  ctrs, edges, and binwidth do not reflect that, but
            % histcEdges does.
            histcEdges = [-Inf edges{i}(2:end-1) Inf];

        % If the bin centers were given, compute their edges and widths.
        elseif histBehavior
            c = ctrs{i};
            dc = diff(c);
            edges{i} = [c(1) c] + [-dc(1) dc dc(end)]/2;
            binwidth{i} = diff(edges{i});
            % Make histc mimic hist behavior:  everything < ctrs(1) gets
            % counted in first bin, everything > ctrs(end) gets counted in
            % last bin.  ctrs, edges, and binwidth do not reflect that, but
            % histcEdges does.
            histcEdges = [-Inf edges{i}(2:end-1) Inf];

        % If the bin edges were given, compute their widths and centers (if
        % asked for).
        else % if ~histBehavior
            e = edges{i};
            de = diff(e);
            histcEdges = e;
            % Make the display mimic bar's histc behavior: an implied bin
            % above edges(end), the same width as the last explicit one.
            % ctrs, edges, and binwidth need that explicitly, histcEdges
            % doesn't.
            edges{i} = [e e(end)+de(end)];
            binwidth{i} = [de de(end)];
            if nargout > 1
                c = zeros(size(de));
                c(1) = e(1) + de(1)/2;
                for j = 2:length(c)
                    c(j) = 2*e(j) - c(j-1);
                end
                % When edges are specified, it may not be possible to return
                % centers for which the edges are midpoints.  Warn if that's
                % the case.
                if any(c <= e(1:end-1)) || ...
                   abs(c(end) - (e(end)-de(end)/2)) > 1000*eps(de(end));
                    warning('stats:hist3:InconsistentEdges', ...
                            'Cannot compute centers that are consistent with EDGES.');
                    c = e(1:end-1) + de/2;
                end
                ctrs{i} = [c e(end)+de(end)/2];
            end
        end

        % Get the 1D bin numbers for this column of x.  Make sure +Inf
        % goes into the nth bin, not the (n+1)th.

        % bin(j,1) holds the x bin idx for data point j
        % bin(j,2) holds the y bin idx for data point j
        [dum,bin(:,i)] = histc(x(:,i),histcEdges,1);
        bin(:,i) = min(bin(:,i),nbins(i));
    end

    % Combine the two vectors of 1D bin counts into a grid of 2D bin
    % counts
    if exist('weights','var'),
      nzidx = all(bin>0,2);
      n = accumarray(bin(nzidx,:),weights(nzidx),nbins);
    else
      n = accumarray(bin(all(bin>0,2),:),1,nbins);
    end

end

if 0 < nargout
    nn = n;
    return
end

del = .001; % space between bars, relative to bar size

% Build x-coords for the eight corners of each bar.
xx = edges{1};
xx = [xx(1:nbins(1))+del*binwidth{1}; xx(2:nbins(1)+1)-del*binwidth{1}];
xx = [reshape(repmat(xx(:)',2,1),4,nbins(1)); repmat(NaN,1,nbins(1))];
xx = [repmat(xx(:),1,4) repmat(NaN,5*nbins(1),1)];
xx = repmat(xx,1,nbins(2));

% Build y-coords for the eight corners of each bar.
yy = edges{2};
yy = [yy(1:nbins(2))+del*binwidth{2}; yy(2:nbins(2)+1)-del*binwidth{2}];
yy = [reshape(repmat(yy(:)',2,1),4,nbins(2)); repmat(NaN,1,nbins(2))];
yy = [repmat(yy(:),1,4) repmat(NaN,5*nbins(2),1)];
yy = repmat(yy',nbins(1),1);

% Build z-coords for the eight corners of each bar.
zz = zeros(5*nbins(1), 5*nbins(2));
zz(5*(1:nbins(1))-3, 5*(1:nbins(2))-3) = n;
zz(5*(1:nbins(1))-3, 5*(1:nbins(2))-2) = n;
zz(5*(1:nbins(1))-2, 5*(1:nbins(2))-3) = n;
zz(5*(1:nbins(1))-2, 5*(1:nbins(2))-2) = n;

cax = newplot(cax);
holdState = ishold(cax);

% Plot the bars in a light steel blue.
cc = repmat(cat(3,.75,.85,.95), [size(zz) 1]);

% Plot the surface, using any specified graphics properties to override
% defaults.
h = surf(cax, xx, yy, zz, cc, 'tag','hist3', plotArgs{:});

if ~holdState
    % Set ticks for each bar if fewer than 16 and the centers/edges are
    % integers.  Otherwise, leave the default ticks alone.
    if (nbins(1)<16)
        if histBehavior && all(floor(ctrs{1})==ctrs{1})
            set(cax,'xtick',ctrs{1});
        elseif ~histBehavior && all(floor(edges{1})==edges{1})
            set(cax,'xtick',edges{1});
        end
    end
    if (nbins(2)<16)
        if histBehavior && all(floor(ctrs{2})==ctrs{2})
            set(cax,'ytick',ctrs{2});
        elseif ~histBehavior && all(floor(edges{2})==edges{2})
            set(cax,'ytick',edges{2});
        end
    end

    % Set the axis limits to have some space at the edges of the bars.
    dx = range(edges{1})*.05;
    dy = range(edges{2})*.05;
    set(cax,'xlim',[edges{1}(1)-dx edges{1}(end)+dx]);
    set(cax,'ylim',[edges{2}(1)-dy edges{2}(end)+dy]);

    view(cax,3);
    grid(cax,'on');
    set(get(cax,'parent'),'renderer','zbuffer');
end

if nargout > 0
    nn = n;
end

author: kristinbranson

repo: JAABA


text:
function g = sigmoidGradient(z)
%SIGMOIDGRADIENT returns the gradient of the sigmoid function
%evaluated at z
%   g = SIGMOIDGRADIENT(z) computes the gradient of the sigmoid function
%   evaluated at z. This should work regardless if z is a matrix or a
%   vector. In particular, if z is a vector or matrix, you should return
%   the gradient for each element.

g = zeros(size(z));

% ====================== YOUR CODE HERE ======================
% Instructions: Compute the gradient of the sigmoid function evaluated at
%               each value of z (z can be a matrix, vector or scalar).

g = sigmoid(z).*(1 - sigmoid(z));

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

end

author: fewtime

repo: ML

save_path: github-repos/MATLAB/fewtime-ML


text:
function name = p13_name ( )

%*****************************************************************************80
%
%% P13_NAME returns the name of problem 13.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string NAME, the name of the problem.
%
  name = 'BFN03';

  return
end


text:
function cut = knapsack_create_cover_cut(a,b,x,alg,gubs)
% Derive cover inequalities cut*[1;x]>=0 for a*x <= b

if (all(a>=0) && sum(a)<=b) || (all(a<=0) && (b>=0))
    % quick exit sum x <= large or sum x >= neg
    cut = [];
    return
end

if nargin < 3
    x = ones(length(a),1);
    alg = 'crowder';
    gubs = [];
elseif nargin < 4
    alg = 'crowder';
    gubs = [];
elseif nargin < 5
    gubs = spalloc(1,length(a),0);
end

if ~all(a)
    nz = find(a);
    a_ = a(nz);
    cut_ = knapsack_create_cover_cut(a_,b,x(nz),alg,gubs(nz));
    if ~isempty(cut_)
        cut = spalloc(1,length(a)+1,0);
        cut(1) = cut_(1);
        cut(1+nz) = cut_(2:end);
    else
        cut = [];
    end
    return
elseif any(a<0)
    % Negative values in a*x <= b
    % models -a*x >= -b i.e. typically something like sum stuff >= bound
    %        -a*(1-y) >= -b
    %        -a*y <= b-sum(a)
    neg = find(a < 0);
    a_ = a;
    a_(neg) = -a(neg);
    x_ = x;
    x_(neg) = 1-x(neg);
    cut = knapsack_create_cover_cut(a_,b-sum(a(neg)),x_,alg,gubs);
    % we know have a cut, cut(1) + cut(2:end)*y >=0
    %                     cut(1) + cut(2:end)(1-x) >= 0
    %                     cut(1) + sum(cut(2:end)) + (-cut(2))*x >= 0
    if ~isempty(cut)
        cut(1) = cut(1)+sum(cut(1 + neg));
        cut(1 + neg) = -cut(1 + neg);
    end
    return
end

switch alg
    case {'crowder',''}
        % Heuristics from H. Crowder, E. Johnson, M. Padberg
        % Solving large-scale 0–1 linear programming programs
        [val,loc] = sort((1-x)./a(:),'ascend');
    case 'gu'
        [val,loc] = sort((x),'descend');
    case 'gu-reverse'
        [val,loc] = sort((x),'ascend');

    otherwise
        error('Unsupport cover cut separation')
end
% Initial cover
C = min(find(cumsum(a(loc))>b));

% Apply Balas lifting
q = knapsack_cover_lift_balas(a,loc(1:C));
% qL = knapsack_cover_lift_letchford(a,b,loc(1:C));
% if q~=qL
%    'HEJ'
% end
% Return row where row*[1;x] hopefully is violated
cut = spalloc(1,length(a)+1,0);
cut(1) = C-1;
cut(2:end)=-q;

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function data = clusterincluster(N, r1, r2, w1, w2, arms)

    if nargin < 1
        N = 1000;
    end
    if nargin < 2
        r1 = 1;
    end
    if nargin < 3
        r2 = 5*r1;
    end
    if nargin < 4
        w1 = 0.8;
    end
    if nargin < 5
        w2 = 1/3;
    end
    if nargin < 6
        arms = 64;
    end

    data = [];

    N1 = floor(N/2);
    N2 = N-N1;

    phi1 = rand(N1,1) * 2 * pi;
    dist1 = r1 + randint(N1,1,3)/3 * r1 * w1;
    d1 = [dist1 .* cos(phi1) dist1 .* sin(phi1) zeros(N1,1)];

    perarm = round(N2/arms);
    N2 = perarm * arms;
    radperarm = (2*pi)/arms;
    phi2 = ((1:N2) - mod(1:N2, perarm))/perarm * (radperarm);
    phi2 = phi2';
    dist2 = r2 * (1 - w2/2) + r2 * w2 * mod(1:N2, perarm)'/perarm;
    d2 = [dist2 .* cos(phi2) dist2 .* sin(phi2) ones(N2,1)];

    data = [d1;d2];

    scatter(data(:,1), data(:,2), 20, data(:,3)); axis square;
end

repo: gspbox


text:
function normalise = spm_cfg_normalise
% SPM Configuration file for toolbox 'Old Normalise'
%__________________________________________________________________________
% Copyright (C) 2005-2012 Wellcome Trust Centre for Neuroimaging

% $Id: spm_cfg_normalise.m 7155 2017-08-17 10:55:05Z john $

if ~isdeployed, addpath(fullfile(spm('dir'),'toolbox','OldNorm')); end

%--------------------------------------------------------------------------
% source Source Image
%--------------------------------------------------------------------------
source         = cfg_files;
source.tag     = 'source';
source.name    = 'Source Image';
source.help    = {'The image that is warped to match the template(s).  The result is a set of warps, which can be applied to this image, or any other image that is in register with it.'};
source.filter  = 'image';
source.ufilter = '.*';
source.num     = [1 1];

%--------------------------------------------------------------------------
% wtsrc Source Weighting Image
%--------------------------------------------------------------------------
wtsrc         = cfg_files;
wtsrc.tag     = 'wtsrc';
wtsrc.name    = 'Source Weighting Image';
wtsrc.val     = {''};
wtsrc.help    = {'Optional weighting images (consisting of pixel values between the range of zero to one) to be used for registering abnormal or lesioned brains.  These images should match the dimensions of the image from which the parameters are estimated, and should contain zeros corresponding to regions of abnormal tissue.'};
wtsrc.filter  = 'image';
wtsrc.ufilter = '.*';
wtsrc.num     = [0 1];

%--------------------------------------------------------------------------
% subj Subject
%--------------------------------------------------------------------------
subj         = cfg_branch;
subj.tag     = 'subj';
subj.name    = 'Subject';
subj.val     = {source wtsrc };
subj.help    = {'Data for this subject.  The same parameters are used within subject.'};

%--------------------------------------------------------------------------
% esubjs Data
%--------------------------------------------------------------------------
esubjs         = cfg_repeat;
esubjs.tag     = 'esubjs';
esubjs.name    = 'Data';
esubjs.help    = {'List of subjects. Images of each subject should be warped differently.'};
esubjs.values  = {subj };
esubjs.num     = [1 Inf];

%--------------------------------------------------------------------------
% template Template Image
%--------------------------------------------------------------------------
template         = cfg_files;
template.tag     = 'template';
template.name    = 'Template Image';
template.help    = {'Specify a template image to match the source image with. The contrast in the template must be similar to that of the source image in order to achieve a good registration.  It is also possible to select more than one template, in which case the registration algorithm will try to find the best linear combination of these images in order to best model the intensities in the source image.'};
template.filter  = 'image';
template.ufilter = '.*';
template.dir     = fileparts(mfilename('fullpath'));
template.num     = [1 Inf];

%--------------------------------------------------------------------------
% weight Template Weighting Image
%--------------------------------------------------------------------------
weight         = cfg_files;
weight.tag     = 'weight';
weight.name    = 'Template Weighting Image';
weight.val     = {''};
weight.help    = {
                  'Applies a weighting mask to the template(s) during the parameter estimation.  With the default brain mask, weights in and around the brain have values of one whereas those clearly outside the brain are zero.  This is an attempt to base the normalisation purely upon the shape of the brain, rather than the shape of the head (since low frequency basis functions can not really cope with variations in skull thickness).'
                  ''
                  'The option is now available for a user specified weighting image. This should have the same dimensions and mat file as the template images, with values in the range of zero to one.'
}';
weight.filter  = 'image';
weight.ufilter = '.*';
weight.num     = [0 1];

%--------------------------------------------------------------------------
% smosrc Source Image Smoothing
%--------------------------------------------------------------------------
smosrc         = cfg_entry;
smosrc.tag     = 'smosrc';
smosrc.name    = 'Source Image Smoothing';
smosrc.help    = {'Smoothing to apply to a copy of the source image. The template and source images should have approximately the same smoothness. Remember that the templates supplied with SPM have been smoothed by 8mm, and that smoothnesses combine by Pythagoras'' rule.'};
smosrc.strtype = 'e';
smosrc.num     = [1 1];
smosrc.def     = @(val)spm_get_defaults('old.normalise.estimate.smosrc', val{:});

%--------------------------------------------------------------------------
% smoref Template Image Smoothing
%--------------------------------------------------------------------------
smoref         = cfg_entry;
smoref.tag     = 'smoref';
smoref.name    = 'Template Image Smoothing';
smoref.help    = {'Smoothing to apply to a copy of the template image. The template and source images should have approximately the same smoothness. Remember that the templates supplied with SPM have been smoothed by 8mm, and that smoothnesses combine by Pythagoras'' rule.'};
smoref.strtype = 'e';
smoref.num     = [1 1];
smoref.def     = @(val)spm_get_defaults('old.normalise.estimate.smoref', val{:});

%--------------------------------------------------------------------------
% regtype Affine Regularisation
%--------------------------------------------------------------------------
regtype         = cfg_menu;
regtype.tag     = 'regtype';
regtype.name    = 'Affine Regularisation';
regtype.help    = {'Affine registration into a standard space can be made more robust by regularisation (penalising excessive stretching or shrinking).  The best solutions can be obtained by knowing the approximate amount of stretching that is needed (e.g. ICBM templates are slightly bigger than typical brains, so greater zooms are likely to be needed). If registering to an image in ICBM/MNI space, then choose the first option.  If registering to a template that is close in size, then select the second option.  If you do not want to regularise, then choose the third.'};
regtype.labels  = {
                  'ICBM space template'
                  'Average sized template'
                  'No regularisation'
}';
regtype.values  = {
                  'mni'
                  'subj'
                  'none'
}';
regtype.def     = @(val)spm_get_defaults('old.normalise.estimate.regtype', val{:});

%--------------------------------------------------------------------------
% cutoff Nonlinear Frequency Cutoff
%--------------------------------------------------------------------------
cutoff         = cfg_entry;
cutoff.tag     = 'cutoff';
cutoff.name    = 'Nonlinear Frequency Cutoff';
cutoff.help    = {'Cutoff of DCT bases.  Only DCT bases of periods longer than the cutoff are used to describe the warps. The number used will depend on the cutoff and the field of view of the template image(s).'};
cutoff.strtype = 'e';
cutoff.num     = [1 1];
cutoff.def     = @(val)spm_get_defaults('old.normalise.estimate.cutoff', val{:});

%--------------------------------------------------------------------------
% nits Nonlinear Iterations
%--------------------------------------------------------------------------
nits         = cfg_entry;
nits.tag     = 'nits';
nits.name    = 'Nonlinear Iterations';
nits.help    = {'Number of iterations of nonlinear warping performed.'};
nits.strtype = 'w';
nits.num     = [1 1];
nits.def     = @(val)spm_get_defaults('old.normalise.estimate.nits', val{:});

%--------------------------------------------------------------------------
% reg Nonlinear Regularisation
%--------------------------------------------------------------------------
reg         = cfg_entry;
reg.tag     = 'reg';
reg.name    = 'Nonlinear Regularisation';
reg.help    = {'The amount of regularisation for the nonlinear part of the spatial normalisation. Pick a value around one.  However, if your normalised images appear distorted, then it may be an idea to increase the amount of regularisation (by an order of magnitude) - or even just use an affine normalisation. The regularisation influences the smoothness of the deformation fields.'};
reg.strtype = 'e';
reg.num     = [1 1];
reg.def     = @(val)spm_get_defaults('old.normalise.estimate.reg', val{:});

%--------------------------------------------------------------------------
% eoptions Estimation Options
%--------------------------------------------------------------------------
eoptions      = cfg_branch;
eoptions.tag  = 'eoptions';
eoptions.name = 'Estimation Options';
eoptions.val  = {template weight smosrc smoref regtype cutoff nits reg };
eoptions.help = {'Various settings for estimating warps.'};

%--------------------------------------------------------------------------
% est Old Normalise: Estimate
%--------------------------------------------------------------------------
est         = cfg_exbranch;
est.tag     = 'est';
est.name    = 'Old Normalise: Estimate';
est.val     = {esubjs eoptions };
est.help    = {'Computes the warp that best registers a source image (or series of source images) to match a template, saving it to a file imagename''_sn.mat''.'};
est.prog    = @spm_run_normalise;
est.vout    = @vout_estimate;

%--------------------------------------------------------------------------
% matname Parameter File
%--------------------------------------------------------------------------
matname         = cfg_files;
matname.tag     = 'matname';
matname.name    = 'Parameter File';
matname.help    = {'Select the ''_sn.mat'' file containing the spatial normalisation parameters for that subject.'};
matname.filter  = 'mat';
matname.ufilter = '.*_sn\.mat$';
matname.num     = [1 1];

%--------------------------------------------------------------------------
% resample Images to Write
%--------------------------------------------------------------------------
resample         = cfg_files;
resample.tag     = 'resample';
resample.name    = 'Images to Write';
resample.help    = {'These are the images for warping according to the estimated parameters. They can be any images that are in register with the "source" image used to generate the parameters.'};
resample.filter  = 'image';
resample.ufilter = '.*';
resample.num     = [1 Inf];

%--------------------------------------------------------------------------
% subj Subject
%--------------------------------------------------------------------------
subj         = cfg_branch;
subj.tag     = 'subj';
subj.name    = 'Subject';
subj.val     = {matname resample };
subj.help    = {'Data for this subject.  The same parameters are used within subject.'};

%--------------------------------------------------------------------------
% wsubjs Data
%--------------------------------------------------------------------------
wsubjs         = cfg_repeat;
wsubjs.tag     = 'wsubjs';
wsubjs.name    = 'Data';
wsubjs.help    = {'List of subjects. Images of each subject should be warped differently.'};
wsubjs.values  = {subj };
wsubjs.num     = [1 Inf];

%--------------------------------------------------------------------------
% preserve Preserve
%--------------------------------------------------------------------------
preserve         = cfg_menu;
preserve.tag     = 'preserve';
preserve.name    = 'Preserve';
preserve.help    = {
                    'Preserve Concentrations: Spatially normalised images are not "modulated". The warped images preserve the intensities of the original images.'
                    ''
                    'Preserve Total: Spatially normalised images are "modulated" in order to preserve the total amount of signal in the images. Areas that are expanded during warping are correspondingly reduced in intensity.'
}';
preserve.labels = {
                   'Preserve Concentrations'
                   'Preserve Amount'
}';
preserve.values = {0 1};
preserve.def     = @(val)spm_get_defaults('old.normalise.write.preserve', val{:});

%--------------------------------------------------------------------------
% bb Bounding box
%--------------------------------------------------------------------------
bb         = cfg_entry;
bb.tag     = 'bb';
bb.name    = 'Bounding box';
bb.help    = {'The bounding box (in mm) of the volume which is to be written (relative to the anterior commissure).'};
bb.strtype = 'e';
bb.num     = [2 3];
bb.def     = @(val)spm_get_defaults('old.normalise.write.bb', val{:});

%--------------------------------------------------------------------------
% vox Voxel sizes
%--------------------------------------------------------------------------
vox         = cfg_entry;
vox.tag     = 'vox';
vox.name    = 'Voxel sizes';
vox.help    = {'The voxel sizes (x, y & z, in mm) of the written normalised images.'};
vox.strtype = 'e';
vox.num     = [1 3];
vox.def     = @(val)spm_get_defaults('old.normalise.write.vox', val{:});

%--------------------------------------------------------------------------
% interp Interpolation
%--------------------------------------------------------------------------
interp         = cfg_menu;
interp.tag     = 'interp';
interp.name    = 'Interpolation';
interp.help    = {
                  ['The method by which the images are sampled when ' ...
                  'being written in a different space. ' ...
                  '(Note that Inf or NaN values are treated as zero, ' ...
                  'rather than as missing data)']
                  '    Nearest Neighbour:'
                  '      - Fastest, but not normally recommended.'
                  '    Trilinear Interpolation:'
                  '      - OK for PET, realigned fMRI, or segmentations'
                  '    B-spline Interpolation:'
                  ['      - Better quality (but slower) interpolation' ...
                  '/* \cite{thevenaz00a}*/, especially with higher ' ...
                  'degree splines. Can produce values outside the ' ...
                  'original range (e.g. small negative values from an ' ...
                  'originally all positive image).']
}';
interp.labels = {
                 'Nearest neighbour'
                 'Trilinear'
                 '2nd Degree B-spline'
                 '3rd Degree B-Spline '
                 '4th Degree B-Spline '
                 '5th Degree B-Spline'
                 '6th Degree B-Spline'
                 '7th Degree B-Spline'
}';
interp.values = {0 1 2 3 4 5 6 7};
interp.def     = @(val)spm_get_defaults('old.normalise.write.interp', val{:});

%--------------------------------------------------------------------------
% wrap Wrapping
%--------------------------------------------------------------------------
wrap         = cfg_menu;
wrap.tag     = 'wrap';
wrap.name    = 'Wrapping';
wrap.help    = {
                'These are typically:'
                '    No wrapping: for PET or images that have already been spatially transformed. '
                '    Wrap in  Y: for (un-resliced) MRI where phase encoding is in the Y direction (voxel space).'
}';
wrap.labels  = {
               'No wrap'
               'Wrap X'
               'Wrap Y'
               'Wrap X & Y'
               'Wrap Z'
               'Wrap X & Z'
               'Wrap Y & Z'
               'Wrap X, Y & Z'
}';
wrap.values  = {[0 0 0] [1 0 0] [0 1 0] [1 1 0] [0 0 1] [1 0 1] [0 1 1]...
               [1 1 1]};
wrap.def     = @(val)spm_get_defaults('old.normalise.write.wrap', val{:});

%--------------------------------------------------------------------------
% prefix Filename Prefix
%--------------------------------------------------------------------------
prefix         = cfg_entry;
prefix.tag     = 'prefix';
prefix.name    = 'Filename Prefix';
prefix.help    = {'Specify the string to be prepended to the filenames of the normalised image file(s). Default prefix is ''w''.'};
prefix.strtype = 's';
prefix.num     = [1 Inf];
prefix.def     = @(val)spm_get_defaults('old.normalise.write.prefix', val{:});

%--------------------------------------------------------------------------
% roptions Writing Options
%--------------------------------------------------------------------------
roptions      = cfg_branch;
roptions.tag  = 'roptions';
roptions.name = 'Writing Options';
roptions.val  = {preserve bb vox interp wrap prefix };
roptions.help = {'Various options for writing normalised images.'};

%--------------------------------------------------------------------------
% write Old Normalise: Write
%--------------------------------------------------------------------------
write         = cfg_exbranch;
write.tag     = 'write';
write.name    = 'Old Normalise: Write';
write.val     = {wsubjs roptions };
write.help    = {'Allows previously estimated warps (stored in imagename''_sn.mat'' files) to be applied to series of images.'};
write.prog    = @spm_run_normalise;
write.vout    = @vout_write;

%--------------------------------------------------------------------------
% source Source Image
%--------------------------------------------------------------------------
source         = cfg_files;
source.tag     = 'source';
source.name    = 'Source Image';
source.help    = {'The image that is warped to match the template(s).  The result is a set of warps, which can be applied to this image, or any other image that is in register with it.'};
source.filter  = 'image';
source.ufilter = '.*';
source.num     = [1 1];

%--------------------------------------------------------------------------
% wtsrc Source Weighting Image
%--------------------------------------------------------------------------
wtsrc         = cfg_files;
wtsrc.tag     = 'wtsrc';
wtsrc.name    = 'Source Weighting Image';
wtsrc.val     = {''};
wtsrc.help    = {'Optional weighting images (consisting of pixel values between the range of zero to one) to be used for registering abnormal or lesioned brains.  These images should match the dimensions of the image from which the parameters are estimated, and should contain zeros corresponding to regions of abnormal tissue.'};
wtsrc.filter  = 'image';
wtsrc.ufilter = '.*';
wtsrc.num     = [0 1];

%--------------------------------------------------------------------------
% resample Images to Write
%--------------------------------------------------------------------------
resample         = cfg_files;
resample.tag     = 'resample';
resample.name    = 'Images to Write';
resample.help    = {'These are the images for warping according to the estimated parameters. They can be any images that are in register with the "source" image used to generate the parameters.'};
resample.filter  = 'image';
resample.ufilter = '.*';
resample.num     = [1 Inf];

%--------------------------------------------------------------------------
% subj Subject
%--------------------------------------------------------------------------
subj         = cfg_branch;
subj.tag     = 'subj';
subj.name    = 'Subject';
subj.val     = {source wtsrc resample };
subj.help    = {'Data for this subject.  The same parameters are used within subject.'};

%--------------------------------------------------------------------------
% ewsubjs Data
%--------------------------------------------------------------------------
ewsubjs         = cfg_repeat;
ewsubjs.tag     = 'ewsubjs';
ewsubjs.name    = 'Data';
ewsubjs.help    = {'List of subjects. Images of each subject should be warped differently.'};
ewsubjs.values  = {subj };
ewsubjs.num     = [1 Inf];

%--------------------------------------------------------------------------
% estwrite Old Normalise: Estimate & Write
%--------------------------------------------------------------------------
estwrite      = cfg_exbranch;
estwrite.tag  = 'estwrite';
estwrite.name = 'Old Normalise: Estimate & Write';
estwrite.val  = {ewsubjs eoptions roptions };
estwrite.help = {'Computes the warp that best registers a source image (or series of source images) to match a template, saving it to the file imagename''_sn.mat''. This option also allows the contents of the imagename''_sn.mat'' files to be applied to a series of images.'};
estwrite.prog = @spm_run_normalise;
estwrite.vout = @vout_estwrite;

%--------------------------------------------------------------------------
% oldnorm Old Normalise
%--------------------------------------------------------------------------
normalise         = cfg_choice;
normalise.tag     = 'oldnorm';
normalise.name    = 'Old Normalise';
normalise.help    = {
                     'This very ancient module /* \cite{ashburner97b,ashburner99a} */ spatially (stereotactically) normalises MRI, PET or SPECT images into a standard space defined by some ideal model or template image[s].  The template images supplied with SPM conform to the space defined by the ICBM, NIH P-20 project, and approximate that of the the space described in the atlas of Talairach and Tournoux (1988). The transformation can also be applied to any other image that has been coregistered with these scans. A few researchers may wish to continue using this strategy, but (when good quality anatomical MRI scans are available) the DARTEL approach is now generally recommended instead.'
                     ''
                     'Generally, the algorithms work by minimising the sum of squares difference between the image which is to be normalised, and a linear combination of one or more template images.  For the least squares registration to produce an unbiased estimate of the spatial transformation, the image contrast in the templates (or linear combination of templates) should be similar to that of the image from which the spatial normalisation is derived.  The registration simply searches for an optimum solution.  If the starting estimates are not good, then the optimum it finds may not find the global optimum.'
                     ''
                     'The first step of the normalisation is to determine the optimum 12-parameter affine transformation.  Initially, the registration is performed by matching the whole of the head (including the scalp) to the template.  Following this, the registration proceeded by only matching the brains together, by appropriate weighting of the template voxels.  This is a completely automated procedure (that does not require ``scalp editing'') that discounts the confounding effects of skull and scalp differences.   A Bayesian framework is used, such that the registration searches for the solution that maximises the a posteriori probability of it being correct /* \cite{ashburner97b} */.  i.e., it maximises the product of the likelihood function (derived from the residual squared difference) and the prior function (which is based on the probability of obtaining a particular set of zooms and shears).'
                     ''
                     'The affine registration is followed by estimating nonlinear deformations, whereby the deformations are defined by a linear combination of three dimensional discrete cosine transform (DCT) basis functions /* \cite{ashburner99a} */.  The default options result in each of the deformation fields being described by 1176parameters, where these represent the coefficients of the deformations in three orthogonal directions.  The matching involved simultaneously minimising the membrane energies of the deformation fields and the residual squared difference between the images and template(s).'
                     ''
                     'The primarily use is for stereotactic normalisation to facilitate inter-subject averaging and precise characterisation of functional anatomy /* \cite{ashburner97bir} */.  It is not necessary to spatially normalise the data (this is only a pre-requisite  for  inter-subject averaging or reporting in the Talairach space).  If you wish to circumnavigate this step  (e.g. if you have single slice data or do not have an appropriate high resolution MRI scan) simply specify where you think the  anterior commissure  is  with  the  ORIGIN in the header of the first scan (using the ''Display'' facility) and proceed directly  to ''Smoothing''or ''Statistics''.'
                     ''
                     'All normalised images are written to the same subdirectory as the original images, prefixed with a ''w''.  The details of the transformations are displayed in the results window, and the parameters are saved in the "*_sn.mat" file.'
}';
normalise.values   = {est write estwrite};

%==========================================================================
function dep = vout_estimate(job)
for k=1:numel(job.subj)
    dep(k)            = cfg_dep;
    dep(k).sname      = sprintf('Norm Params File (Subj %d)',k);
    dep(k).src_output = substruct('()',{k},'.','params');
    dep(k).tgt_spec   = cfg_findspec({{'filter','mat','strtype','e'}});
end

%==========================================================================
function dep = vout_write(job)
for k=1:numel(job.subj)
    dep(k)            = cfg_dep;
    dep(k).sname      = sprintf('Normalised Images (Subj %d)',k);
    dep(k).src_output = substruct('()',{k},'.','files');
    dep(k).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
end

%==========================================================================
function dep = vout_estwrite(job)
depe = vout_estimate(job);
depw = vout_write(job);
dep = [depe depw];

author: spm


text:
function q = p06_q ( m, c, w )

%*****************************************************************************80
%
%% P06_Q evaluates the integral for problem p06.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Alan Genz,
%    A Package for Testing Multiple Integration Subroutines,
%    in Numerical Integration: Recent Developments, Software
%    and Applications,
%    edited by Patrick Keast and Graeme Fairweather,
%    Reidel, 1987, pages 337-340,
%    ISBN: 9027725144,
%    LC: QA299.3.N38.
%
%  Parameters:
%
%    Input, integer M, the dimension of the argument.
%
%    Input, real C(M), W(M), the problem parameters.
%
%    Output, real Q, the integral.
%

%
%  To simplify the calculation, force W(3:M) to be at least 1.0.
%
  w(3:m) = 1.0;

  q = 1.0;

  for i = 1 : m

    if ( w(i) <= 0.0 )

      q = q * 0.0;

    elseif ( w(i) <= 1.0 )

      if ( c(i) == 0.0 )
        q = q * w(i);
      else
        q = q * ( exp ( c(i) * w(i) ) - 1.0 ) / c(i);
      end

    else

      if ( c(i) ~= 0.0 )
        q = q * ( exp ( c(i) * w(i) ) - 1.0 ) / c(i);
      end

    end

  end

  return
end


text:
%% #######################     HEADER START    ############################
%*************************************************************************
%
% Filename:				Retrieve_Directory_Structure.m
%
% Author:				A. Mering
% Created:				22-Feb-2013
%
% Changed on:			XX-XX-XXXX  by USERNAME		SHORT CHANGE DESCRIPTION
%						XX-XX-XXXX  by USERNAME		SHORT CHANGE DESCRIPTION
%
%*************************************************************************
%
% Description:
%		Recursively analyze the directory structure and get all files according to defined file type(s).
%
%
% Input parameter:
%		- root_dir:			Root directory from where to start the analysis
%		- file_types:		Cell array of file types to be take into account
%
% Output parameter:
%		- file_structure:	Structure containing information about all the files occuring
%		- num_files_below:  Number of files matching search expression within lower lying directories
%		- exclusion_list:	List of file names to be excluded
%
%*************************************************************************
%
% Intrinsic Subfunctions
%		- none
%
% Intrinsic Callbacks
%		- none
%
% #######################      HEADER END     ############################

%% #######################    FUNCTION START   ############################
function [directory_structure, num_files] = Retrieve_Directory_Structure(root_dir, file_types, exclusion_list)

%% Input Validation
% start with input checks

% Is the root_dir of type char and does the directory exist?
if ~ischar(root_dir) || ~exist(root_dir, 'dir')
	error('Input directory is not an existing directory!')
end

% Is filetype of type cell?
if ~iscell(file_types)
	% no? then make a cell out of it
	file_types = {file_types};
end

% Are all cell elements of type char?
if ~all(cellfun(@ischar, file_types))
	error('Invalid input for file types')
end

% Are all cell elements of the form '*....'?
if any(cellfun('isempty', regexp(file_types, '^\*\..+$')))
	error('File types should be of the form ''*.xyz''!')
end

%% PostProcess input parameters
% append '$' to identify files via regexp
file_types = cellfun(@(x) sprintf('%s$',x), file_types, 'UniformOutput', false);
file_types = regexprep(file_types, '\$*$', '\$');

%% Initialize
root_content = dir(root_dir);

max_num_found_files_string_length = ceil(log10(sum(~cell2mat({root_content.isdir})+1)));
max_num_found_subdirs_string_length = ceil(log10(sum(cell2mat({root_content.isdir})+1)));

num_found_files = 0;
num_found_subdirs = 0;
num_dir_level_below = 0;
num_files_below = 0;
num_files = 0;

directory_structure = struct();

%% Get directory content
for n = 1 : size(root_content)
	current_item = root_content(n);

	% Check for directory -> Recursion
	if current_item.isdir

		if isempty(regexp(current_item.name, '^\.{1,2}$', 'once'))
			[directory_content, num_files_within] = Retrieve_Directory_Structure([root_dir,'\',current_item.name], file_types, exclusion_list);

			num_files_below = num_files_below + num_files_within;

			% Only append directories with found items
			if ~isempty(fieldnames(directory_content))

				num_dir_level_below = max([num_dir_level_below, directory_content.Statistics_Num_Dir_Level_Below + 1]);

				num_found_subdirs = num_found_subdirs + 1;

				structure_label = sprintf(['Directory%0',num2str(max_num_found_subdirs_string_length),'i'],num_found_subdirs);

				directory_structure.(structure_label) = directory_content;
				directory_structure.(structure_label).DirName = current_item.name;

				directory_structure.(structure_label) = orderfields(directory_structure.(structure_label));
			end
		end

	else

		% Check filetypes

		if ~ismember(current_item.name, exclusion_list)
			found_type_pos = find(~cellfun('isempty', regexp(current_item.name, regexprep(file_types, '^\*', '\'))),1);

			if any(strcmp(file_types, '*.*$')) || ~isempty(found_type_pos)
				num_found_files = num_found_files + 1;

				structure_label = sprintf(['File%0',num2str(max_num_found_files_string_length),'i'],num_found_files);

				directory_structure.(structure_label) = current_item.name;

			end
		end
	end

end

num_files = num_files_below + num_found_files;

directory_structure = orderfields(directory_structure);
if num_found_subdirs  ~= 0 || num_found_files ~= 0
	directory_structure.Statistics_Num_Dirs = num_found_subdirs;
	directory_structure.Statistics_Num_Dir_Level_Below = num_dir_level_below;
	directory_structure.Statistics_Num_Files_Current = num_found_files;
	directory_structure.Statistics_Num_Files_Below = num_files_below;

	[~,directory_structure.DirName ] = fileparts(root_dir);
end

% #######################     FUNCTION END    ############################

%% #######################  SUBFUNCTION START  ############################

% #######################   SUBFUNCTION END   ############################

%% #######################    CALLBACK START   ############################

% #######################     CALLBACK END    ############################

author: opencobra

repo: cobratoolbox


text:
function [TruePositives, FalseNegatives] = testAromaticAADegradation(model, microbeID, biomassReaction, database)
% Performs an FVA and reports those AromaticAA pathway end reactions (exchange reactions)
% that can carry flux in the model and should carry flux according to
% data (true positives) and those AromaticAA pathway end reactions that
% cannot carry flux in the model but should be secreted according to in
% vitro data (false negatives).
%
% INPUT
% model             COBRA model structure
% microbeID         Microbe ID in carbon source data file
% biomassReaction   Biomass objective functions (low flux through BOF
%                   required in analysis)
% database          Structure containing rBioNet reaction and metabolite
%                   database
%
% OUTPUT
% TruePositives     Cell array of strings listing all aromatic amino acid
% degradation products
% that can be secreted by the model and in comnparative genomic data.
% FalseNegatives    Cell array of strings listing all aromatic amino acid
% degradation products
% that cannot be secreted by the model but should be secreted according to comparative genomic data.
%
% .. Author:
%      Almut Heinken, Dec 2017
%                     March  2022 - changed code to string-matching to make
%                     it more robust

global CBT_LP_SOLVER
if isempty(CBT_LP_SOLVER)
    initCobraToolbox
end

% read aromatic amino acid degradation product table
dataTable = readInputTableForPipeline('AromaticAATable.txt');

corrRxns = {'Phenylpropanoate','EX_pppn(e)';'4-Hydroxyphenylpropanoate','EX_r34hpp(e)';'Indolepropionate','EX_ind3ppa(e)';'Isocaproate','EX_isocapr(e)'};

TruePositives = {};  % true positives (uptake in vitro and in silico)
FalseNegatives = {};  % false negatives (uptake in vitro not in silico)

% find microbe index in data table
mInd = find(strcmp(dataTable(:,1), microbeID));
if isempty(mInd)
    warning(['Microbe "', microbeID, '" not found in aromatic amino acid degradation data file.'])
else
    % perform FVA to identify uptake metabolites
    % set BOF
    if ~any(ismember(model.rxns, biomassReaction)) || nargin < 3
        error(['Biomass reaction "', biomassReaction, '" not found in model.'])
    end
    model = changeObjective(model, biomassReaction);
    % set a low lower bound for biomass
    %     model = changeRxnBounds(model, biomassReaction, 1e-3, 'l');
    % list exchange reactions
    exchanges = model.rxns(strncmp('EX_', model.rxns, 3));
    % open all exchanges
    model = changeRxnBounds(model, exchanges, -1000, 'l');
    model = changeRxnBounds(model, exchanges, 1000, 'u');

    % get the reactions to test
    rxns = {};
    for i=2:size(dataTable,2)
        if contains(version,'(R202') % for Matlab R2020a and newer
            if dataTable{mInd,i}==1
                findCorrRxns = find(strcmp(corrRxns(:,1),dataTable{1,i}));
                rxns = union(rxns,corrRxns(findCorrRxns,2:end));
            end
        else
            if strcmp(dataTable{mInd,i},'1')
                findCorrRxns = find(strcmp(corrRxns(:,1),dataTable{1,i}));
                rxns = union(rxns,corrRxns(findCorrRxns,2:end));
            end
        end
    end

    % flux variability analysis on reactions of interest
    rxns = unique(rxns);
    rxns = rxns(~cellfun('isempty', rxns));
    rxnsInModel=intersect(rxns,model.rxns);
    if ~isempty(rxnsInModel)
        currentDir=pwd;
        try
            [minFlux, maxFlux, ~, ~] = fastFVA(model, 0, 'max', 'ibm_cplex', ...
                rxnsInModel, 'S');
        catch
            warning('fastFVA could not run, so fluxVariability is instead used. Consider installing fastFVA for shorter computation times.');
            cd(currentDir)
            [minFlux, maxFlux] = fluxVariability(model, 0, 'max', rxnsInModel);
        end

        % active flux
        flux = rxnsInModel(maxFlux > 1e-6);
    else
        flux = {};
    end

    % which reaction should carry flux according to in vitro data
    for i=2:size(dataTable,2)
        rxn={};
        if contains(version,'(R202') % for Matlab R2020a and newer
            if dataTable{mInd,i}==1
                rxn = corrRxns{find(strcmp(corrRxns(:,1),dataTable{1,i})),2};
            end
        else
            if strcmp(dataTable{mInd,i},'1')
                rxn = corrRxns{find(strcmp(corrRxns(:,1),dataTable{1,i})),2};
            end
        end
        if ~isempty(rxn)
            % add any that are not in model/not carrying flux to the false negatives
            if ~isempty(intersect(rxn,flux))
                TruePositives = union(TruePositives,rxn);
            else
                FalseNegatives=union(FalseNegatives,rxn);
            end
        end
    end
end

% replace reaction IDs with metabolite names
if ~isempty(TruePositives)
    TruePositives = TruePositives(~cellfun(@isempty, TruePositives));
    TruePositives=strrep(TruePositives,'EX_','');
    TruePositives=strrep(TruePositives,'(e)','');

    for i=1:length(TruePositives)
        TruePositives{i}=database.metabolites{find(strcmp(database.metabolites(:,1),TruePositives{i})),2};
    end
end

% warn about false negatives
if ~isempty(FalseNegatives)
    FalseNegatives = FalseNegatives(~cellfun(@isempty, FalseNegatives));
    FalseNegatives=strrep(FalseNegatives,'EX_','');
    FalseNegatives=strrep(FalseNegatives,'(e)','');
    for i = 1:length(FalseNegatives)
        FalseNegatives{i}=database.metabolites{find(strcmp(database.metabolites(:,1),FalseNegatives{i})),2};
    end
end

end


text:
%% Display Start message
disp('Running QAM BER/SER Simulation set')
%% Plot theoretical curves
SNRs = -4:28;
disp('Plot theoretical curves')
[h_fig, h_lines] = QAM_BER_Curves(SNRs);
%% Run Monte Carlo simulations
disp('Run Monte Carlo Simulations')

% Create place-holder plots
QAM_BER = zeros(9,length(SNRs));
hold on
simLines = semilogy(SNRs, QAM_BER,'*');

tic
[QAM_BER, QAM_SER] = QAM_Simulate(SNRs, simLines);
toc
%% Plot simulation results
% The simulation results are plotted inside the QAM_Simulate function.

%% Add Legend
% Do MATLAB graphics magic to create concise legend
% See "Controling Legends" in MATLAB doc

% group lines together
simGrp = hggroup('DisplayName','Simulation');
theoGrp = hggroup('DisplayName','Theoretical');
set(simLines,'Parent',simGrp)
set(h_lines,'Parent',theoGrp)
set(get(get(simGrp,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','on'); % Include this hggroup in the legend
set(get(get(theoGrp,'Annotation'),'LegendInformation'),...
    'IconDisplayStyle','on'); % Include this hggroup in the legend
legend show

disp('This figure is identical to [1, Fig. 5].')
disp('[1] Cho, K., and Yoon, D., "On the general BER expression of one- and')
disp('    two-dimensional amplitude modulations", IEEE Trans. Commun.,')
disp('    Vol. 50, Number 7, pp. 1074-1080, 2002.')


text:
function pass = test_ellipjODE(pref)
% Test to check that the chebfun command for Jacobi elliptic functions
% ELLIPJ produces a similar result to the solution of the nonlinear
% differential equation that generates it.

% NOTE: Taken from V4 test chebop_ellipjode.

if ( nargin == 0 )
    pref = cheboppref();
end
tol = 1e-10;

% Elliptic parameter
m = .1;
K = ellipke(m);
d = [0, K];
x = chebfun(@(x) x, d);

%% jacobi elliptic functions
[sn, cn, dn] = ellipj(x, m);

% SN is the solution of an ODE
N = chebop(d);
N.op = @(x,u) diff(u,2) + (1+m)*u - 2*m*u.^3;
N.lbc = @(u) u;% The solution here is close to being singular, and we only enjoy linear convergence  so loosen the tolerance
N.rbc = @(u) u - 1 ;
u = N\0;
err(1) = norm(u - sn, inf);

%% CN
N.op = @(x, u) diff(u, 2) + (1-2*m)*u + 2*m*u.^3;
N.lbc = @(u) u - 1 ;
N.rbc = @(u) u ;
u = mldivide(N, 0, pref);
err(2) = norm(u - cn, inf);

% %% DN
% % The solution here is close to being singular, and we only enjoy linear
% % convergence of Newton iteratin. Loosen the tolerance.
% pref.bvpTol = 5e-5;
% N.op = @(x, u) diff(u, 2) - (2-m)*u + 2*u.^3;
% N.lbc = @(u) u - 1 ;
% N.rbc = @(u) u - sqrt(1-m);
% u = mldivide(N, 0, pref);
% err(3) = norm(u - dn, inf);

%%
% pass = err < [tol, tol, 2e6*tol];
pass = err < [tol, tol];

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function y = sqrt_cordic ( x, n )

%*****************************************************************************80
%
%% SQRT_CORDIC returns the square root of a value using the CORDIC method.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, the number whose square root is desired.
%
%    Input, integer N, the number of iterations to take.
%    This is essentially the number of binary digits of accuracy.
%
%    Output, real Y, the approximate square root of X.
%
  if ( x < 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SQRT_CORDIC - Fatal error!\n' );
    fprintf ( 1, '  X < 0.\n' );
    error ( 'SQRT_CORDIC - Fatal error!' )
  end

  if ( x == 0.0 )
    y = 0.0;
    return
  end

  if ( x == 1.0 )
    y = 1.0;
    return
  end

  poweroftwo = 1.0;

  if ( x < 1.0 )

    while ( x <= poweroftwo * poweroftwo )
      poweroftwo = poweroftwo / 2.0;
    end

    y = poweroftwo;

  elseif ( 1.0 < x )

    while ( poweroftwo * poweroftwo <= x )
      poweroftwo = 2.0 * poweroftwo;
    end

    y = poweroftwo / 2.0;

  end

  for i = 1 : n
    poweroftwo = poweroftwo / 2.0;
    if ( ( y + poweroftwo ) * ( y + poweroftwo ) <= x )
      y = y + poweroftwo;
    end
  end

  return
end


text:
% Description : This is a matlab file for video enhancement.
% The following file works with color, b&w images.
% User defined functions : 1) image_enhancement_sw - modified retinex algorithm
%                          2) gray_level_images    - contrast enhancement
% Author : Vijay Sridharan.
% Date   : 20th July, 2012.
%% Video Enhancement for color and black and white images
video = VideoReader('air.avi');
    for i = 1:video.NumberOfFrames
        img = read(video,i);
        imwrite(img,sprintf('img%d.jpg',i));
    end
    fprintf('Please wait....');
    filebase = dir('*.jpg');                      % If you are using a continous frame of images, start from here.
    num_files = numel(filebase);
    images = cell(1, num_files);
     MS=cell(1,num_files);
     for k = 1:num_files
         images{k} = imread(filebase(k).name);
         [rows columns color]=size(images{1});
         if (color==3)
             MS{k}=image_enhancement_sw(images{k});
             M(k)=im2frame(MS{k});
         else
             MS{k}=gray_level_images(images{k}); % For gray level images, movie/implay doesn't work. So, I would suggest you
                                                 % to download FIJI or ImageJ to view the sequence of images as Videos. The
         end                                     % default frame rate is 8fps in FIJI
     end
movie(M)


% sin_x.m: This m-file calculates and plots the
% function sin(x) for 0 <= x <= 6.
x = 0:0.1:6
y = sin(x)
plot(x,y)


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
% function [fh,ph,th] = checkDerivative(fctn,x0);
%
% checks the implementation of a derivative by comparing the function
% with the Taylor-poly
%
%   \| f(x0 + h ) - TP_p(x0,h) \|   !=   O( h^{p+1} )
%
% Input:
%  fctn    function handle
%  x0      expanding point
%
% Output:
%  fh      figure handle to graphical output
%  ph      plot handle
%  th      text handle
% call checkDerivative for a minimal example.
%==============================================================================

function varargout = checkDerivative(fctn,x0,varargin)

if nargin == 0, % help and minimal example
  help(mfilename);
  fctn = @xSquare;
  x0   = 1;
  checkDerivative(fctn,x0);
  return;
end;

fig = [];
for k=1:2:length(varargin),     % overwrites default parameter
  eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;

FAIRmessage(mfilename,'.')
fprintf('test derivative of function <%s>\n',func2str(fctn));

fprintf('T0 = |f0 - ft|, T1 = |f0+h*f0'' - ft|\n');
[f0,df] = feval(fctn,x0);
if ~isnumeric(df),
  try
    [f0,para,df] = feval(fctn,x0);
  catch
    if isfield(df,'Q'), df = df.Q;  end;
  end;
end;

h = logspace(-1,-10,10);
v = randn(size(x0));
if isnumeric(df),
  dvf = df*v;
elseif isa(df,'function_handle'),
  dvf = df(v);
else
  keyboard;
end;

for j=1:length(h),
  ft = feval(fctn,x0+h(j)*v);      % function value
  T0(j) = norm(f0-ft);             % TaylorPoly 0
  T1(j) = norm(f0+h(j)*dvf - ft);  % TaylorPoly 1
  fprintf('h=%12.4e     T0=%12.4e    T1=%12.4e\n',h(j),T0(j),T1(j));
end;

fh = FAIRfigure(fig);
ph = loglog(h,[T0;T1]); set(ph(2),'linestyle','--')
th = title(sprintf('%s: |f-f(h)|,|f+h*dvf -f(h)| vs. h',mfilename));

if nargout>0,
  varargout = {fh,ph,th};
end;
FAIRmessage('.');
%------------------------------------------------------------------------------
function [y,dy] = xSquare(x)
if nargin == 0, return; end;
y = x.^2; dy = reshape(2*x,1,[]);
%==============================================================================


text:
%ROC Receiver-Operator Curve
%
%   E = ROC(B,C,N)
%
% INPUT
%   A  Dataset
%   W  Trained classifier, or
%   B  Classification result, B = A*W*CLASSC
%   C  Index of desired class (default: C = 1)
%   N  Number of points on the Receiver-Operator Curve (default: 100)
%
% OUTPUT
%   E  Structure containing the error of the two classes
%
% DESCRIPTION
% Computes N points on the receiver-operator curve of the classifier W for
% class C in the labeled dataset B, which is typically the result of
% B = A*W; or for the dataset A labelled by applying the (cell array of)
% trained classifiers W.
%
% Note that a Receiver-Operator Curve is related to a specific class (class C)
% for which the errors are plotted horizontally. The total error on all other
% classes is plotted vertically. The class index C refers to its position in
% the label list of the dataset (A or B). It can be found by GETCLASSI.
%
% The curve is computed for N thresholds of the posteriori probabilities
% stored in B. The resulting error frequencies for the two classes are
% stored in the structure E. E.XVALUES contains the errors in the first
% class, E.ERROR contains the errors in the second class. In multi-class
% problems these are the mean values in a single class, respectively the
% mean values in all other classes. This may not be very useful, but not
% much more can be done as for multi-class cases the ROC is equivalent to a
% multi-dimensional surface.
%
% Use PLOTE(E) for plotting the result. In the plot the two types of error
% are annotated as 'Error I' (error of the first kind) and 'Error II' (error
% of the second kind). All error estimates are weighted according the class
% prior probabilities. Remove the priors in A or B (by setprior(A,[])) to
% produce a vanilla ROC.
%
% EXAMPLES
%	Train set A and test set T:
%	  B = T*NMC(A); E = ROC(T,50); PLOTE(E); % Plots a single curve
%	  E = ROC(T,A*{NMC,UDC,QDC});  PLOTE(E); % Plots 3 curves
%

function [roc, thr] = roc_detection(a,n)

	% Depending on the call, CLAS may the third or second argument.
	% and N the third or the fourth.

	if nargin < 2 || isempty(n), n = 100; end

	datname = getname(a);
	lablist = getlablist(a,'string');

    required_labs = {'FP' 'TP'};
    aux_lablist = intersect(cellstr(lablist), required_labs );

    cant_required_labs = length(required_labs);

    if( length(aux_lablist) ~= cant_required_labs )
        fprintf(2, ['Esta funcion esta pensada para usarse en datasets de deteccion con labels:\n' colvec([char(required_labs) repmat('\n', cant_required_labs, 1)]')' ] );
        error();
    end

    clas = find( strcmpi(cellstr(lablist), 'TP') );
	clasname = lablist(clas,:);
	%DXD: also check the class sizes:
	cs = classsizes(a);
	if any(cs == 0)
		error('Ambas clases deben contener ejemplos');
	end

	% Set up the ROC structure.

    % If a cell array of classifiers was given, apply each one here.

    a = a*normm; % make sure we have a normalised classification matrix

    [m,c] = size(a);
    nlab = getnlab(a);
    d = sort(a(:));

    % Attempt to compute a good threshold range on the first classified
    % dataset.
%     thr = [max(0, min(d)-eps) rowvec(d(round(linspace(2,length(d)-1, n-1)))) 1];
    thr = linspace(0,1,n+1);

    % NLAB_OUT will be one where B is larger than THR.
    I = matchlablist(getlablist(a),getfeatlab(a)); % Correct possible changes class orders
    nlab_out = (repmat(+a(:,I(clas)),1,n+1) >= repmat(thr,m,1));

    % aciertos will be 1 where the numeric label is unequal to NLAB_OUT
    % (i.e., where errors occur).
    bTP = nlab == clas;
    aciertos = (repmat(bTP,1,n+1) == nlab_out);

    sensitivity = mean(aciertos(bTP,:),1); % S
    pospred = sum(aciertos(bTP,:)) ./ sum(nlab_out); % +P

    roc = [ colvec(sensitivity) colvec(pospred) ];

    mod = sqrt(sum(roc.^2,2));
    [max_mod max_mod_idx] = max(mod);

    figure(100);
    h = subplot(1,2,1);
    plot(h(1), sensitivity(:), pospred(:), 'bo-' )
    hold(h(1), 'on');
    plot(h(1), sensitivity(max_mod_idx), pospred(max_mod_idx), 'rx:', 'MarkerSize',11)
    plot(h(1), sensitivity(max_mod_idx), pospred(max_mod_idx), 'mo:', 'MarkerSize',11)
    hold(h(1), 'off');
    axis(h(1), 'square');
    box(h(1), 'off')
    xlabel(h(1), 'Sensitivity')
    ylabel(h(1), 'Positive predictivity')
%     x_lim = xlim();
%     y_lim = ylim();
%     xlim([min([x_lim ylim]) 1.1]);
%     xlim([min([x_lim ylim]) 1.1]);

    h(2) = subplot(1,2,2);
    plot(h(2), thr, mod, 'bo-' )
    hold(h(2), 'on');
    plot(h(2), thr(max_mod_idx), mod(max_mod_idx), 'rx', 'MarkerSize',11)
    plot(h(2), thr(max_mod_idx), mod(max_mod_idx), 'mo', 'MarkerSize',11)
    plot(h(2), [0.5 0.5], ylim(), 'k--')
    hold(h(2), 'off');
    box(h(2), 'off')
    ylabel(h(2), 'mod(S,P+)')
    xlabel(h(2), 'Operating point')

    ConfusionMat = [ sum(aciertos(~bTP,max_mod_idx))  sum(~aciertos(~bTP,max_mod_idx)); sum(~aciertos(bTP,max_mod_idx)) sum(aciertos(bTP,max_mod_idx)) ];

    DisplayResults('dsResult', ConfusionMat, 'SupportDataset', a);

author: marianux

repo: ecg-kit

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


text:
function [ i, j ] = ij_next ( i, j, n )

%*****************************************************************************80
%
%% IJ_NEXT returns the next matrix index.
%
%  Discussion:
%
%    For N = 3, the sequence of indices returned is:
%
%      (1,1), (1,2), (1,3), (2,1), (2,2), (2,3), (3,1), (3,2), (3,3), (0,0).
%
%    Note that once the value (N,N) is returned, the next value returned
%    will be (0,0).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, J, the current pair of indices.
%
%    Input, integer N, the maximum value for I and J.
%
%    Output, integer I, J, the next pair of indices.  If either index is
%    illegal on input, the output value of (I,J) will be (1,1).
%
  if ( n < 1 )
    i = 0;
    j = 0;
    return
  end

  if ( i < 1 | n < i | j < 1 | n < j )
    i = 1;
    j = 1;
    return
  end

  if ( j < n )
    j = j + 1;
  elseif ( i < n )
    i = i + 1;
    j = 1;
  else
    i = 0;
    j = 0;
  end

  return
end


text:
classdef Blob < handle
  % Wrapper class of caffe::Blob in matlab

  properties (Access = private)
    hBlob_self
  end

  methods
    function self = Blob(hBlob_blob)
      CHECK(is_valid_handle(hBlob_blob), 'invalid Blob handle');

      % setup self handle
      self.hBlob_self = hBlob_blob;
    end
    function shape = shape(self)
      shape = caffe_('blob_get_shape', self.hBlob_self);
    end
    function reshape(self, shape)
      shape = self.check_and_preprocess_shape(shape);
      caffe_('blob_reshape', self.hBlob_self, shape);
    end
    function data = get_data(self)
      data = caffe_('blob_get_data', self.hBlob_self);
    end
    function set_data(self, data)
      data = self.check_and_preprocess_data(data);
      caffe_('blob_set_data', self.hBlob_self, data);
    end
    function diff = get_diff(self)
      diff = caffe_('blob_get_diff', self.hBlob_self);
    end
    function set_diff(self, diff)
      diff = self.check_and_preprocess_data(diff);
      caffe_('blob_set_diff', self.hBlob_self, diff);
    end
  end

  methods (Access = private)
    function shape = check_and_preprocess_shape(~, shape)
      CHECK(isempty(shape) || (isnumeric(shape) && isrow(shape)), ...
        'shape must be a integer row vector');
      shape = double(shape);
    end
    function data = check_and_preprocess_data(self, data)
      CHECK(isnumeric(data), 'data or diff must be numeric types');
      self.check_data_size_matches(data);
      if ~isa(data, 'single')
        data = single(data);
      end
    end
    function check_data_size_matches(self, data)
      % check whether size of data matches shape of this blob
      % note: matlab arrays always have at least 2 dimensions. To compare
      % shape between size of data and shape of this blob, extend shape of
      % this blob to have at least 2 dimensions
      self_shape_extended = self.shape;
      if isempty(self_shape_extended)
        % target blob is a scalar (0 dim)
        self_shape_extended = [1, 1];
      elseif isscalar(self_shape_extended)
        % target blob is a vector (1 dim)
        self_shape_extended = [self_shape_extended, 1];
      end
      % Also, matlab cannot have tailing dimension 1 for ndim > 2, so you
      % cannot create 20 x 10 x 1 x 1 array in matlab as it becomes 20 x 10
      % Extend matlab arrays to have tailing dimension 1 during shape match
      data_size_extended = ...
        [size(data), ones(1, length(self_shape_extended) - ndims(data))];
      is_matched = ...
        (length(self_shape_extended) == length(data_size_extended)) ...
        && all(self_shape_extended == data_size_extended);
      CHECK(is_matched, ...
        sprintf('%s, input data/diff size: [ %s] vs target blob shape: [ %s]', ...
        'input data/diff size does not match target blob shape', ...
        sprintf('%d ', data_size_extended), sprintf('%d ', self_shape_extended)));
    end
  end
end

author: aimerykong

repo: deepImageAestheticsAnalysis


text:
function [tri, coord] = puffball(mask)

%Weighted skeleton, where each non-zero value corresponds to the maximum
%size of a sphere that can be centered there.
wSkel = get_SkelRadius(mask);

heightFunction = puffbalIInflation(mask,wSkel);

[tri,coord] = mesh_from_height(mask, heightFunction);

end

function [tri, coord] =  mesh_from_height(mask, rec)

[M,N] = size(mask);
[x,y] = meshgrid(1:N,1:M);

tri = delaunay(x,y);

ind_tri = mask(tri);
ind_tri = ind_tri(:,1) | ind_tri(:,2) | ind_tri(:,3);

tri = tri(ind_tri,:);
aux = unique(tri);
boundary = setdiff(aux, find(mask));

points=zeros(M*N,1);
nPoints = length(aux);
points(aux) = 1:nPoints;

tri = points(tri);
boundary = points(boundary);
x = x(aux);
y = y(aux);
z = rec(aux);

for k =1:length(boundary)
    pBound = boundary(k);
    sel = sum(tri==pBound,2)>0;
    selPoints = setdiff(unique(tri(sel,:)),pBound);

    x1 = median(x(selPoints));
    y1 = median(y(selPoints));

    x(pBound) = x1;
    y(pBound) = y1;

end

points = zeros(nPoints,1);
points(boundary) = boundary;
points(points==0) = nPoints + (1:(nPoints-length(boundary)));

reflected_tri = points(tri);
tri = [tri; reflected_tri];
x = [x;x];
y = [y;y];
z = [z; -z];
x(boundary + nPoints) = [];
y(boundary + nPoints) = [];
z(boundary + nPoints) = [];

coord = [x y z];

%figure; trisurf(tri,x,y,z,'FaceColor','r');
%axis equal
%set(gca,'YDir','rev');

end

function [ h ] = puffbalIInflation( mask,wSkel )
% TAKE THE UNION (SOFT-MAX) OF MAXIMAL SPHERES %
[Y, X] = meshgrid(1:size(mask,2),1:size(mask,1));
h = ones(size(mask));

[y,x] = find(wSkel);

k = 1;

for i = 1:length(x)
    r = wSkel(y(i),x(i))^2 - (X-y(i)).^2 - (Y-x(i)).^2;
    h(r>0) = h(r>0)+exp(k*sqrt(r(r>0)));
end

h = log(h)/k;

end

function [wSkel, allRadius] = get_SkelRadius(mask)
%Weighted skeleton, where each non-zero value corresponds to the maximum
%size of a sphere that can be centered there.

% CALCULATE GRASSFIRE HEIGHT FUNCTION %
% A 3x3-tap filter to smoothly erode an anti-aliased edge

fil = [0.1218 0.4123 0.1218; 0.4123 0.9750 0.4123; ...
    0.1218 0.4123 0.1218]/1.2404;
nmask = double(mask);
allRadius = zeros(size(mask));
while ~isempty(find(nmask,1))
    allRadius = allRadius+nmask/1.67; % Each iteration erodes the edge .6 pixels
    nmaskpad = padarray(nmask,[1 1],'replicate');
    nmaskpad = conv2(nmaskpad,fil,'same')-1.4241;
    nmask = max(min(nmaskpad(2:end-1,2:end-1),1),0);
end

% LOCATE THE MEDIAL AXIS %
[dx, dy] = gradient(allRadius);
dsurf = sqrt(dx.^2+dy.^2);
% Medial axis points have a grassfire gradient measurably less than 1
radThreshold = min(max(allRadius(:)),2);
wSkel = bwmorph(dsurf<0.958&allRadius>=radThreshold,'skel',Inf).*allRadius;

end

repo: CategoryShapes


text:
%The Gabor filter is basically a Gaussian (with variances sx and sy along x and y-axes respectively)
%modulated by a complex sinusoid (with centre frequencies U and V along x and y-axes respectively)
%described by the following equation
%%
%               1                -1     x  ^    y  ^
%%% G(x,y) = ---------- * exp ([----{(----) 2+(----) 2}+2*pi*i*(Ux+Vy)])
%            2*pi*sx*sy           2    sx       sy

%% But we use one dimentional gabor filter as below

%              1              -1    x  ^
%%% G(x) = ---------- * exp ([----(----) 2+2*pi*i*Ux])
%           2*pi*sx            2    sx

%And

%              1              -1    y  ^
%%% G(y) = ---------- * exp ([----(----) 2+2*pi*i*Uy])
%           2*pi*sy            2    sy
% By using these filters,reduces the filtering operations from O(M2N2) to O(MN2).

%The output is
%%% Iu(x,y) = Conv(I(x,y),g(y));
%%% Iv(x,y) = Conv(I(x,y),g(x));
%%%Ifilt = sqrt(Iu^2+Iv^2);
%% Where I is input image.

%% Describtion :

%% I : Input image
%% Sx & Sy : Variances along x and y-axes respectively
%% U & V : Centre frequencies  along x and y-axes respectively

%% Gx & Gy : The output filters as described above
%% gabout : The output filtered image

%%  Author : Ahmad poursaberi  e-mail : a.poursaberi@ece.ut.ac.ir
%%          Faulty of Engineering, Electrical&Computer Department,Tehran
%%          University,Iran,June 2004

function [Gx,Gy,gabout] = gaborfilter1D(I,Sx,Sy,U,V);

if isa(I,'double')~=1
    I = double(I);
end
warning off
for x = -fix(Sx):fix(Sx)
    for y = -fix(Sy):fix(Sy)
        Gx(fix(Sx)+x+1) = (1/(2*pi*Sx))*exp(-.5*((x/Sx)^2)+2*pi*i*U*x);
        Gy(fix(Sy)+y+1) = (1/(2*pi*Sy))*exp(-.5*((y/Sy)^2)+2*pi*i*V*y);

    end
end

gaboutu = conv2(I,double(Gy),'same');
gaboutv = conv2(I,double(Gx),'same');
gabout = uint8(sqrt(gaboutu.*gaboutu + gaboutv.*gaboutv));


text:
function [E, V] = VBA_dirichlet_moments (a)
% // VBA toolbox //////////////////////////////////////////////////////////
%
% [E, V] = VBA_dirichlet_moments (a)
% Derives the first- and second-order moments of a Dirichlet density
%
% IN:
%   - a: parameters of the Dirichlet distribution (pseudo-counts)
% OUT:
%   - E: first-order moment (expectation)
%   - V: second-order moment (variance)
%
% /////////////////////////////////////////////////////////////////////////

% Let X = (X_1, ..., X_k) ~ Dir(a)

% shorthand
a0 = sum(a);

% First-order moment E[X]
E = a./a0;

% second order moment Var[X]
V = (a .* (a0-a)) ./ (a0.^2 * (a0+1));

author: MBB-team

repo: VBA-toolbox


text:
function beamlet = createIMBeamlet(doseV, indV, beamNum, fullLength)
%"createIMBeamlet"
%   Take a vector of dose values and a vector of indices into the scan
%   array and build an IM.beamlets element.
%
%JRA 9/20/04
%
%Usage:
%   function beamlet = createIMBeamlet(doseV, indV, beamNum, fullLength)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

%Initialize the beamlet.
% IM = initIMRTProblem;
% beamlet = IM.beamlets;
beamlet = initBeamlet;
beamlet(1).format  = 'uint8';
beamlet(1).beamNum = beamNum;
beamlet(1).fullLength = fullLength;

if isempty(indV) || isempty(doseV)
%     disp(['Warning: a pencil beam in beam ' num2str(beamNum) ' does not contribute any dose to this structure.'])
    return;
end

ind2V = find(doseV);

if isempty(ind2V)
    valsV = [];
    maxVal = 0;
    smallVals = logical([]);
else
	valsV = doseV(ind2V);

	maxVal = max(valsV(:));

	%Bool vector to note low values, helps avoid uint8 roundoff error.
	smallVals = valsV < (maxVal/(2^8 - 1));

	%Store non small values as normal uint8s.
	valsV(~smallVals) = (valsV(~smallVals)/maxVal) * (2^8 - 1);

	%Store small values as uint8s with another factor of 256.
	valsV(smallVals) = valsV(smallVals)/(maxVal)*(2^8 - 1)*(2^8 - 1);
end

%Save vector to ID low values when builidng inflM later.
%Uses a logical packer to decrease size by about 8.

beamlet(1).lowDosePoints = packLogicals(smallVals);
beamlet(1).influence = uint8(valsV);
beamlet(1).indexV = uint32(indV(ind2V));
beamlet(1).maxInfluenceVal = maxVal;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function s=get_name(a)

s=['bayesian model selection'];
eval_name                   %% print values of hyperparameters

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [Population,FrontNo,WHVLoss] = EnvironmentalSelection(Population,N,wz,AA,RA)
% The environmental selection of I-SIBEA

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

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

    %% Calculate the WHV loss of each solution front by front
    WHVLoss = CalWHVLoss(Population.objs,FrontNo,wz,AA,RA);

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

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

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function toms179_test02 ( )

%*****************************************************************************80
%
%% TOMS179_TEST02 demonstrates the use of MDBETA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2008
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TOMS179_TEST02:\n' );
  fprintf ( 1, '  MDBETA estimates the value of the modified Beta function.\n' );
  fprintf ( 1, '  Compare with tabulated values.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X         P         Q         ' );
  fprintf ( 1, 'Beta                      Beta                    Diff\n' );
  fprintf ( 1, '                                    ' );
  fprintf ( 1, '(Tabulated)               (MDBETA)\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, p, q, x, fx ] = beta_cdf_values ( n_data );

    if ( n_data == 0 )
      break
    end

    [ fx2, ifault ] = mdbeta ( x, p, q );

    fprintf ( 1, '  %8.4f  %8.4f  %8.4f  %24.16e  %24.16e  %10.4e\n', ...
    p, q, x, fx, fx2, abs ( fx - fx2 ) );

  end

  return
end


text:
function rgb = skin_dark

% This returns a predefined color as [red green blue] values
%   red               = [255   0   0]/255;
%   green             = [  0 192   0]/255;
%   blue              = [  0   0 255]/255;
%   magenta           = [255 255   0]/255;
%   cyan              = [  0 255 255]/255;
%   yellow            = [255 255   0]/255;
%   white             = [255 255 255]/255;
%   black             = [  0   0   0]/255;
%
%   skull             = [140  85  85]/255
%   cortex            = [255 213 119]/255;
%   cortex_light      = [199 194 169]/255;
%   cortex_dark       = [100  97  85]/255;
%   skin              = [249 223 192]/255;
%   skin_light        = [249 223 192]/255;
%   skin_medium_light = [225 194 158]/255;
%   skin_medium       = [188 142 106]/255;
%   skin_medium_dark  = [155 102	65]/255;
%   skin_dark         = [ 91  71  61]/255;
%
% The different skin-based colors follow the Fitzpatrick scale with type I and II
% combined, and return RGB values that approximate those used by Apple in the emoji
% skin tones. See also https://emojipedia.org/emoji-modifier-sequence/
%
% If no specific skin tone is specified, this function returns a light skin color.
% This corresponds with that of one of the developers who approximated his own skin
% color more than 15 years ago upon the first implementation of this function.

% returns a predefined color as [red green blue] values
rgb = [ 91	 71	 61]/255;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function pick = nms_matlab(boxes, overlap)
% top = nms_fast(boxes, overlap)
% Non-maximum suppression. (FAST VERSION)
% Greedily select high-scoring detections and skip detections
% that are significantly covered by a previously selected
% detection.
% NOTE: This is adapted from Pedro Felzenszwalb's version (nms.m),
% but an inner loop has been eliminated to significantly speed it
% up in the case of a large number of boxes
% Tomasz Malisiewicz (tomasz@cmu.edu)

if isempty(boxes)
    pick = [];
    return;
end

x1 = boxes(:,1);
y1 = boxes(:,2);
x2 = boxes(:,3);
y2 = boxes(:,4);
s = boxes(:,end);

area = (x2-x1+1) .* (y2-y1+1);
[vals, I] = sort(s);

pick = s*0;
counter = 1;
while ~isempty(I)

    last = length(I);
    i = I(last);
    pick(counter) = i;
    counter = counter + 1;

    xx1 = max(x1(i), x1(I(1:last-1)));
    yy1 = max(y1(i), y1(I(1:last-1)));
    xx2 = min(x2(i), x2(I(1:last-1)));
    yy2 = min(y2(i), y2(I(1:last-1)));

    w = max(0.0, xx2-xx1+1);
    h = max(0.0, yy2-yy1+1);

    o = w.*h ./ area(I(1:last-1));

    I([last; find(o>overlap)]) = [];
end

pick = pick(1:(counter-1));

author: peiyunh

repo: tiny

save_path: github-repos/MATLAB/peiyunh-tiny


text:
function a = ctranspose(a)
%CTRANSPOSE Complex conjugate transpose for tenmat.
%
%   C = CTRANSPOSE(A) swaps the row and column indices of A.
%
%   See also TENMAT.
%
%MATLAB Tensor Toolbox.
%Copyright 2012, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2012) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

tmp = a.rindices;
a.rindices = a.cindices;
a.cindices = tmp;
a.data = ctranspose(a.data);

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt