text:
function [f] = Force(nnel,shape,P)

%--------------------------------------------------------------------------
% Purpose :
%         Determines the force vector for the element 'iel'
% Synopsis :
%          [f] = Force(nnel,shape,P)
% Variable Description:
%           f - element force vector
%           nnel - number of nodes per element
%           shape - Shape functions for the element
%           P - applied transverse pressure
%--------------------------------------------------------------------------

fef = shape*P ;

for i = 1:nnel
    i1=(i-1)*3+1;
    i2=i1+1;
    i3=i2+1;
    f(i1,1) = fef(1) ;
    f(i2,1) = 0 ;
    f(i3,1) = 0 ;
end


text:
function value = i4_factorial2 ( n )

%*****************************************************************************80
%
%% I4_FACTORIAL2 computes the double factorial N!!
%
%  Formula:
%
%    FACTORIAL2( N ) = Product ( N * (N-2) * (N-4) * ... * 2 )  (N even)
%                    = Product ( N * (N-2) * (N-4) * ... * 1 )  (N odd)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 January 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the argument of the double factorial function.
%    If N is less than 1, the value is returned as 1.
%
%    Output, integer VALUE, the value of N!!.
%
  if ( n < 1 )
    value = 1;
    return
  end

  value = 1;

  while ( 1 < n )
    value = value * n;
    n = n - 2;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [data,units] = compute_wing_lengthr_mm(trx,n)

flies = trx.exp2flies{n};
nflies = numel(flies);
data = cell(1,nflies);
for i = 1:nflies,
  fly = flies(i);

  data{i} = trx(fly).wing_lengthr ./ trx.pxpermm;

end
units = parseunits('mm');

author: kristinbranson

repo: JAABA


text:
function g = endpoints(f)

persistent lut
if isempty(lut)
    lut = makelut(@endpoint_fcn, 3);
end

g = applylut(f,lut);

function is_end_point = endpoint_fcn(nhood)
is_end_point = nhood(2,2) & (sum(nhood(:)) == 2);

author: decouples


text:
function similarity_output = canlab_pattern_similarity(dat, pattern_weights, varargin)
% Calculate similarity between each column in a data matrix dat and a vector of pattern weights
%
% - Similarity options: dot product, cosine similarity, and correlation
% - Columns are often images, e.g., from fmri_data.dat for fmri_data objects
% - weights are pattern weights from one or more 'signature' patterns (each pattern is a column)
%
% - Assumes dat and pattern_weights matrices have equal rows (voxels) and include valid voxels
% - Removes empty column data column-wise (image-wise) in case some images have uneven voxel coverage
%
% - Used in apply_mask, extract_roi_averages, image_similarity_plot***
%
% similarity_output = canlab_pattern_similarity(dat, pattern_weights)
%
% :Inputs:
%   **dat**  voxels x images matrix of data to compare to patterns
%
%   **pattern_weights** voxels x images matrix of patterns to compare to data
%
% :Optional Inputs:
%
%   **dot_product**
%       [Default] Use dot product
%
%   **cosine_similarity**
%       Use cosine similarity metric for pattern expression instead of dot product
%
%   **correlation, corr**
%       Use correlation metric for pattern expression instead of dot product
%
%   **binary_overlap**
%       Use percent overlap of binary masks instead of dot product. needs
%       binary masks and pattern weights.
%
%   **posterior_overlap**
%       Using percent overlap of binary masks, this option calculates the
%       posterior probability of observing non-zero pattern weights (binary)
%       given binary masks
%
%   **ignore_missing**
%       Suppress printing of warnings when thre are missing values
%       (zeros/NaN) in data images. This function always removes voxels
%       from analysis with missing data values. This does not affect the
%       dot product metric, but does affect cosine similarity and
%       correlation.
%
%   **no_warnings**
%       Suppress output on missing values.
%
%
%   **treat_zero_as_data**
%       In some certain situations, zero value within data.obj could be
%       meaningful. e.g, data.obj is a thresholded map (0 means value underthrethold
%       rather than missing value) or binary map.
%
% :Outputs:
%   **similarity_output**
%       Matrix of similarity measures.
%
%
% :Examples:
% ::
% dat = rand(100, 5);  % 100 voxels, 5 images
% weights = rand(100, 2); % 100 voxels, 2 patterns
% similarity_output = canlab_pattern_similarity(dat, weights, 'cosine_similarity');
%
% :See also:
%
% apply_mask
% image_similarity_plot
% extract_roi_averages, to get individual region averages / local pattern expression
% apply_nps, which does whole-pattern and local regional expression
%
% ..
%    Notes:
%
% Dealing with missing data and partial coverage
% ---------------------------------------------------
% Pattern pattern_weights can have zeros, which may be valid values in voxels,
% i.e., with binary masks
% Data images with values of zero or NaN are considered out-of-mask, as they
% are not valid values. That is, 0 is often treated as a missing data value in images,
% with the exception of "signatures" and binary pattern masks.
% Thus, this function treats values of 0 in DATA images, not pattern masks,
% as missing values, and excludes these voxels from both image and
% mask when calculating similarity.
% Thus, there is an asymmetry between pattern mask and image data
% in considering which voxels to use.
% Otherwise, all dot product, correlation, and other similarity metrics are standard.
%
% When comparing two sets of binary images (e.g. k-means clusters) and the
% cluster of the input solution does not overlap with any of the target
% patterns/clusters, cosine similarity is attempting division by 0. Instead
% of returning NaN, we set those cosine values to 0.
%
% Effects of zeros/missing values on similarity metrics
% ---------------------------------------------------
% Dot product is affected by voxel size, coverage in image and mask, scale
% in general.
% Cosine similarity is affected by coverage in image
% If we remove voxels from analysis that are not in image (NaN or zero) first,
% then cosine similarity is unaffected by coverage, in the sense that the
% upper bound is 1.
% Correlation is affected by coverage in image and mean level in image. If
% there are zeros in image, mean-centering will be affected.
% If we remove NaN/zero values from analysis on an image-by-image basis, we
% are computing partial similarity for the areas of the image that exist.
% This is good in terms of preserving scale and analyzing the data we have,
% but changes the measurement properties of the pattern measures, because
% we're looking only at a subset of the model.
%
% We also need to calculate bad values on an image-by-image basis, not
% relying on remove_empty to exclude voxels with ineligible values
% across the entire set of images.
%
%  Created by Tor Wager - 3/7/17

% Programmer's Notes:
%
% 2017/09/07 Stephan Geuter
%   - added option for percent overlap of binary masks (see also
%   image_similarity_plot.m and riverplot.m
%   - Changed metric selection to string format.
%
% 2022/03/31 Ke Bo
%   - added option for treating zero value in the map as real value rather
%   than missing data
%

% ---------------------------------
% Defaults and optional inputs
% ---------------------------------

% docorr = false;     % run correlation instead of dot-product for pattern expression
sim_metric = 'dotproduct'; % Default: Correlation. SG. docosine = false;   % run cosine sim instead
doignoremissing = false; % ignore warnings for missing voxels
doprintwarnings = true;  % print warnings regarding missing voxels, etc.
treat_zero_as_data=false; % Treat zero value as missing data.

if any(strcmp(varargin, 'ignore_missing'))
    doignoremissing = true;
end

if any(strcmp(varargin, 'no_warnings'))
    doprintwarnings = false;
end

if any(strcmp(varargin, 'cosine_similarity')) % run cosine instead of dot-product
    sim_metric = 'cosine';
end

if any(strcmp(varargin, 'correlation')) || any(strcmp(varargin, 'corr')) % run correlation instead of dot-product
    sim_metric = 'corr';
end

if any(strcmp(varargin, 'binary_overlap')) % run overlap instead of dot-product
    sim_metric = 'overlap';
end

if any(strcmp(varargin, 'dotproduct')) % default. overwrites previous selections
    sim_metric = 'dotproduct';
end

if any(strcmp(varargin, 'posterior_overlap')) % run overlap instead of dot-product
    sim_metric = 'posterior_overlap';
end

if any(strcmp(varargin, 'treat_zero_as_data'))
    treat_zero_as_data = true;
end

% if docosine && docorr, error('Choose either cosine_similarity or correlation, or no optional inputs for dot product'); end

% ---------------------------------
% Variable types and sizes
% ---------------------------------
pattern_weights = double(pattern_weights); % force double b/c of matlab instabilities
dat = double(dat); % force double b/c of matlab instabilities

[n, k] = size(dat);
[n2, npatt] = size(pattern_weights);

if n ~= n2, error('Number of observations must be equal for data and patterns'); end

similarity_output = NaN .* zeros(k, npatt);

% ---------------------------------
% Missing/excluded values image-wise
% ---------------------------------

if treat_zero_as_data==true
    badvals = isnan(dat);
else
    badvals = dat == 0 | isnan(dat);  % Matrix. not used for binary overlap (SG).
end

% ---------------------------------
% Main similarity calculation
% ---------------------------------

% if ~docorr && ~docosine
if strcmp(sim_metric,'dotproduct')
    % dot product. No need to remove missing voxels because dotproduct is
    % the same either way.

    similarity_output = dotproduct(dat, pattern_weights);

else

    % all other metrics
    for i = 1:npatt

        switch sim_metric
            case 'corr'

                    similarity_output(:, i) = image_correlation(dat, pattern_weights(:, i), badvals);
            case 'cosine'

                similarity_output(:, i) = cosine_similarity(dat, pattern_weights(:, i), badvals);

            case 'overlap'

                similarity_output(:, i) = overlap_similarity(dat, pattern_weights(:, i));

            case 'posterior_overlap'

                similarity_output(:, i) = posterior_overlap_similarity(dat, pattern_weights(:, i));

            otherwise
                error('Invalid similarity metric.');
        end

    end

end % pattern sim calculation

% ---------------------------------
% Weight cases based on distance to mean
% ---------------------------------

if any(strcmp(varargin,'weighted')) % if we want to weight data based on similarity to group mean
    mean_dat = mean(dat,2); %sample mean
    distances = squareform(pdist([dat,mean_dat]')); %use other measures?
    weights=1./distances(1:end-1,end);
    weights=weights-min(weights)+1e-12; %min0
    weights=weights./mean(weights); % scale so weights are positive with mean value one
else %
    weights=ones(size(dat,2),1); %just use a weight of one
end

similarity_output = bsxfun(@times,similarity_output,weights);

% ---------------------------------
% Warnings
% ---------------------------------

if doprintwarnings && any(badvals(:)) && ~doignoremissing && ~strcmp(sim_metric,'overlap')

    for i = 1:npatt

        inmask = ~(pattern_weights(:, i) == 0 | isnan(pattern_weights(:, i)));

        bad_in_mask = bsxfun(@times, badvals, inmask);
        bad_in_mask = sum(bad_in_mask);  % how many bad values are in mask, across images

        if any(bad_in_mask)
            fprintf('Warning: Some images have zero values in some of the %3.0f voxels in weight mask. These will be excluded from similarity analysis image-wise.\n', sum(inmask));
            disp('Number of zero or NaN values within weight mask, by input image:');

            for j = 1:length(bad_in_mask)
                fprintf('%3.0f ', bad_in_mask(j));
            end
            fprintf('\n');
        end

    end  % pattern index

end % bad val warnings

end % main function

% ---------------------------------
% *
% Sub-functions
% *
% ---------------------------------

% FUNCTIONS for dot product, cosine_similarity, correlation with missing
% values (voxels) that may vary on a data image-by-image basis.

function pexp = dotproduct(X, weights)

% X = voxels x images data matrix, e.g., dat.dat
% weights is voxels x 1 weight data
% badvals is voxels x images logical matrix of voxels to exclude image-wise

% badvals do not matter
pexp = (weights' * X)';

end % function

function [ab, a, b] = image_norms(dat, pattern_weights, badvals)
% [a, b] = image_norms(dat, pattern_weights, badvals)
% dat = voxels dat images data matrix, e.g., dat
% pattern_weights is voxels dat 1 weight data
% badvals is voxels dat images logical matrix of voxels to exclude image-wise

% Norm for image data.
% All non-zero, non-NaN values valid. zero/nan have no
% effect implictly.

a = nansum(dat .^ 2)' .^ .5;

for i = 1:size(dat, 2)    % Loop because we may have different voxel exclusions in each image

    inmask = ~badvals(:, i);

    b(i, 1) = nansum(pattern_weights(inmask) .^ 2) .^ .5;  % Norm for pattern_weights, excluding out-of-image pattern_weights image-wise

end

ab = a .* b;

end % function

function cossim = cosine_similarity(dat, pattern_weights, badvals)

pexp = dotproduct(dat, pattern_weights);

ab = image_norms(dat, pattern_weights, badvals);

cossim = pexp ./ ab;

% when division by zero (e.g. non-overlapping binary masks) return 0.
% Stephan 2017/3/8
if any(ab==0)
    cossim(isnan(cossim)) = 0;
end

end % function

function r = image_correlation(dat, pattern_weights, badvals, varargin)

for i = 1:size(dat, 2)    % Loop because we may have different voxel exclusions in each image

    inmask = ~badvals(:, i);

    r(i, 1) = corr(pattern_weights(inmask), dat(inmask, i));  % Correlation, excluding out-of-image pattern_weights image-wise

end

end % function

function r = overlap_similarity(dat, pattern_weights)

% check for binary data
if numel(unique(dat(:)))>2 || sum(unique(dat(:))-[0; 1])~=0 ...
        || numel(unique(pattern_weights(:)))>2 || sum(unique(pattern_weights(:))-[0; 1])~=0
    if numel(unique(round(pattern_weights))) == 2
        pattern_weights = round(pattern_weights); % an easy (but maybe not optimal) solution for the resampled binary mask
    else
        error('Binary overlap similarity needs binary data [0 1] input.');
    end
end

% compute overlap
for i = 1:size(dat, 2)

    inmask = isfinite(dat(:,i)) & isfinite(pattern_weights);
    nVox = sum(inmask);

    r(i,1) = sum(dat(inmask,i)==1 & pattern_weights(inmask)==1) / nVox; % overlap excluding NaNs

end

end % function

function r = posterior_overlap_similarity(dat, pattern_weights)

% check for binary data
if numel(unique(dat(:)))>2 || sum(unique(dat(:))-[0; 1])~=0 ...
        || numel(unique(pattern_weights(:)))>2 || sum(unique(pattern_weights(:))-[0; 1])~=0
    if numel(unique(round(pattern_weights))) == 2
        pattern_weights = round(pattern_weights); % an easy (but maybe not optimal) solution for the resampled binary mask
    else
        error('Binary overlap similarity needs binary data [0 1] input.');
    end
end

% compute posterior overlap
for i = 1:size(dat, 2)

    % calculate the space
    inmask = sum(dat,2)~=0;
    % calculate P(pattern | data) using overlaps
    r(i,1) = sum(dat(inmask,i)==1 & pattern_weights(inmask)==1)./sum(dat(inmask,i)==1);

end

end

function R = weightedcorrs(Y, w)
%
%   WEIGHTEDCORRS returns a symmetric matrix R of weighted correlation
%   coefficients calculated from an input T-by-N matrix Y whose rows are
%   observations and whose columns are variables and an input T-by-1 vector
%   w of weights for the observations. This function may be a valid
%   alternative to CORRCOEF if observations are not all equally relevant
%   and need to be weighted according to some theoretical hypothesis or
%   knowledge.
%
%   R = WEIGHTEDCORRS(Y, w) returns a positive semidefinite matrix R,
%   i.e. all its eigenvalues are non-negative (see Example 1).

%
% % ======================================================================
%
%   See also CORRCOEF, COV, STD, MEAN.
%
% % ======================================================================
%
%-*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-*%
%                                                                                %
%            Author: Liber Eleutherios                                             %
%            E-Mail: libereleutherios@gmail.com                             %
%            Date: 23 July 2008                                                      %
%            Updated: 6 June 2012                                                 %
%                                                                                %
%-*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-* -*-*%
%
% % ======================================================================
%
% Check input
ctrl = isvector(w) & isreal(w) & ~any(isnan(w)) & ~any(isinf(w)) & all(w > 0);
if ctrl
    w = w(:) / sum(w);                                                          % w is column vector
else
    error('Check w: it needs be a vector of real positive numbers with no infinite or nan values!')
end
ctrl = isreal(Y) & ~any(isnan(Y)) & ~any(isinf(Y)) & (size(size(Y), 2) == 2);
if ~ctrl
    error('Check Y: it needs be a 2D matrix of real numbers with no infinite or nan values!')
end
ctrl = length(w) == size(Y, 1);
if ~ctrl
    error('size(Y, 1) has to be equal to length(w)!')
end
[T, N] = size(Y);                                                             % T: number of observations; N: number of variables
temp = Y - repmat(w' * Y, T, 1);                                              % Remove mean (which is, also, weighted)
temp = temp' * (temp .* repmat(w, 1, N));                                     % Covariance Matrix (which is weighted)
temp = 0.5 * (temp + temp');                                                  % Must be exactly symmetric
R = diag(temp);                                                               % Variances
R = temp ./ sqrt(R * R');

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
%Program for creating CSV File

%Author : Athi Narayanan S
%M.E, Embedded Systems,
%K.S.R College of Engineering
%Erode, Tamil Nadu, India.
%http://sites.google.com/site/athisnarayanan/

%Program Description
%This program generates a CSV file containing the colors in the output image.
%The CSV format is as follows
%PaletteNumber  R   G   B

function writeCSV(out_map, OutCSVName)

s = size(out_map);
fid = fopen(OutCSVName,'w');
for i=1:s(1)
    fprintf(fid,'%d,',i); %Palette No.
    for j=1:s(2)
        fprintf(fid,'%d,',double(out_map(i,j)));
    end
    fprintf(fid,'\n');
end
fclose(fid);

author: Sable

repo: mcbench-benchmarks


text:
function interior_point = point_inside_polygon(V)
  % POINT_INSIDE_POLYGON  return some point inside of the given polygon
  %
  % interior_point = point_inside_polygon(V)
  %
  % Input
  %  V n by 2 vertex array
  % Output
  %  interior_point  a point inside V
  %
  % Copyright 2011, Alec Jacobson (jacobson@inf.ethz.ch)

  assert(size(V,1) >= 3);
  % start with a point above the polygon in the middle
  p = 1.0 + [mean(V(:,1)) max(V(:,2))];

  % first hits are just mid points on segments
  first_hits = ...
    [V(:,1) + V([end, 1:(end-1)],1) , ...
    V(:,2) + V([end, 1:(end-1)],2)]/2;

  distance_to_p = ...
    sum((repmat(p,size(first_hits,1),1) - first_hits).^2,2);

  first_hit = first_hits(distance_to_p == min(distance_to_p),:);
  first_hit = first_hit(1,:);

  x1 = p(1);
  y1 = p(2);
  x2 = first_hit(1);
  y2 = first_hit(2);
  x3 = V(1:end,1);
  y3 = V(1:end,2);
  x4 = V([end, 1:(end-1)],1);
  y4 = V([end, 1:(end-1)],2);
  % Find intersections along this ray with each line segment's line
  % from http://en.wikipedia.org/wiki/Line-line_intersection
  second_hits = ...
      [((x1.*y2-y1.*x2).*(x3-x4) - (x1-x2).*(x3.*y4-y3.*x4))./...
      ((x1-x2).*(y3-y4) - (y1 - y2).*(x3-x4)), ...
      ((x1.*y2-y1.*x2).*(y3-y4) - (y1-y2).*(x3.*y4-y3.*x4))./...
      ((x1-x2).*(y3-y4) - (y1 - y2).*(x3-x4))];

  % Check if intersection is really on line segment, handling vertical and
  % horizontal lines
  valid_second_hits = ...
    second_hits( ...
      ((second_hits(:,1) >= x3 & second_hits(:,1) <= x4 )| ...
      (second_hits(:,1) <= x3 & second_hits(:,1) >= x4 ) | ...
      (x3 == x4) ...
      ) & ...
      ((second_hits(:,2) >= y3 & second_hits(:,2) <= y4 )| ...
      (second_hits(:,2) <= y3 & second_hits(:,2) >= y4 ) | ...
      (y3 == y4) ...
      ) & ~((x3 == x4)&(y3 == y4))...
    ,:);

  distance_to_first_hit = ...
    sum((repmat(first_hit,size(valid_second_hits,1),1) - ...
      valid_second_hits).^2,2);

  EPSILON = 10.0^-15;
  valid_second_hits = ...
    valid_second_hits(distance_to_first_hit>EPSILON,:);
  distance_to_first_hit = ...
    distance_to_first_hit(distance_to_first_hit>EPSILON,:);

  second_hit = valid_second_hits(...
    distance_to_first_hit == min(distance_to_first_hit),:);
  % shouldn't have to do this:
  second_hit = second_hit(1,:);

  interior_point = (first_hit + second_hit)/2;

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function GUI = sessionGUI_helpMenu(GUI);
% Attach a Help Menu to the session GUI figure.
%
% GUI = sessionGUI_helpMenu(GUI);
%
%
% ras, 07/06.
GUI.menus.help = uimenu('Label', 'Help', 'Separator', 'on');

cb = 'web http://white.stanford.edu/newlm';
uimenu(GUI.menus.help, 'Label', 'mrVista Home / Wiki', ...
       'Separator', 'off', 'Callback', cb);

cb = 'web http://white.stanford.edu/newlm -browser';
uimenu(GUI.menus.help, 'Label', 'mrVista Home / Wiki (external browser)', ...
       'Separator', 'off', 'Callback', cb);

uimenu(GUI.menus.help, 'Label', 'Identify Callback for a Menu Item', ...
        'Separator', 'on', 'Callback', 'helpFindCallback; ');

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = plog(varargin)
%PLOG
%
% y = PLOG(x)
%
% Computes concave perspective log, x(1)*log(x(2)/x(1)) on x>0
%
% Implemented as evalutation based nonlinear operator. Hence, the concavity
% of this function is exploited to perform convexity analysis and rigorous
% modelling.

switch class(varargin{1})

    case 'double'

        if ~isequal(prod(size(varargin{1})),2)
            error('PLOG only defined for 2x1 arguments');
        end
        x = varargin{1};
        % Safe version with defined negative values (helps fmincon when
        % outside feasible region)

        if isequal(x(1),[0])
            varargout{1} = 0;
        else
            varargout{1} = x(1)*log(x(2)/x(1));
        end

    case 'sdpvar'

        if ~isequal(prod(size(varargin{1})),2)
            error('PLOG only defined for 2x1 arguments');
        else
            varargout{1} = yalmip('define',mfilename,varargin{1});
        end

    case 'char'

        operator = CreateBasicOperator('concave','callback');
        operator.range = [-inf inf];
        operator.domain = [0 inf];
        operator.bounds = @bounds;
        operator.convexhull = @convexhull;
        operator.derivative = @derivative;

        varargout{1} = [];
        varargout{2} = operator;
        varargout{3} = varargin{3};

    otherwise
        error([upper(mfilename) ' called with weird argument']);
end

function dp = derivative(x)
z = x(2)/x(1);
dp = [log(z)-1;1./z];

function [L,U] = bounds(xL,xU)
xU(isinf(xU)) = 1e12;
x1 = xL(1)*log(xL(2)/xL(1));
x2 = xU(1)*log(xU(2)/xU(1));
x3 = xL(1)*log(xU(2)/xL(1));
x4 = xU(1)*log(xL(2)/xU(1));
L = min([x1 x2 x3 x4]);

% Stationary in x1 when x1 = x2*exp(-1)
% Increasing in x2, so max at border
p1 = [exp(-1)*xU(2);xU(2)];
x5 = p1(1)*log(p1(2)/p1(1));
if x5(1) >= xL(1) && x5(1) <= xU(1)
    U = max([x1 x2 x3 x4 x5]);
else
    U = max([x1 x2 x3 x4]);
end

function [Ax,Ay,b,K] = convexhull(xL,xU)

x1 = [xL(1);xL(2)];
x2 = [xU(1);xL(2)];
x3 = [xL(1);xU(2)];
x4 = [xU(1);xU(2)];
x5 = (xL+xU)/2;

f1 = plog(x1);
f2 = plog(x2);
f3 = plog(x3);
f4 = plog(x4);
f5 = plog(x5);

df1 = derivative(x1);
df2 = derivative(x2);
df3 = derivative(x3);
df4 = derivative(x4);
df5 = derivative(x5);

[Ax,Ay,b,K] = convexhullConcave2D(x1,f1,df1,x2,f2,df2,x3,f3,df3,x4,f4,df4,x5,f5,df5);

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
% LTFAT - Quadratic time-frequency distributions
%
%  Jordy van Velthoven, 2014 - 2018.
%
%  Quadratic distributions
%    AMBIGUITYFUNCTION   -  Ambiguity function
%    WIGNERVILLEDIST     -  Wigner-Ville distribution
%    DRIHACZEKDIST       -  Discrete Rihaczek distribution
%    QUADTFDIST          -  Generic Quadratic distribution
%
%  Plots
%    PLOTQUADTFDIST      - Plot quadratic time-frequency distribution
%
%  For help, bug reports, suggestions etc. please visit
%  http://github.com/ltfat/ltfat/issues

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [P,Uinit,output] = cp_als(X,R,varargin)
%CP_ALS Compute a CP decomposition of any type of tensor.
%
%   P = CP_ALS(X,R) computes an estimate of the best rank-R
%   CP model of a tensor X using an alternating least-squares
%   algorithm.  The input X can be a tensor, sptensor, ktensor, or
%   ttensor. The result P is a ktensor.
%
%   P = CP_ALS(X,R,'param',value,...) specifies optional parameters and
%   values. Valid parameters and their default values are:
%      'tol' - Tolerance on difference in fit {1.0e-4}
%      'maxiters' - Maximum number of iterations {50}
%      'dimorder' - Order to loop through dimensions {1:ndims(A)}
%      'init' - Initial guess [{'random'}|'nvecs'|cell array]
%      'printitn' - Print fit every n iterations; 0 for no printing {1}
%
%   [P,U0] = CP_ALS(...) also returns the initial guess.
%
%   [P,U0,out] = CP_ALS(...) also returns additional output that contains
%   the input parameters.
%
%   Note: The "fit" is defined as 1 - norm(X-full(P))/norm(X) and is
%   loosely the proportion of the data described by the CP model, i.e., a
%   fit of 1 is perfect.
%
%   NOTE: Updated in various minor ways per work of Phan Anh Huy. See Anh
%   Huy Phan, Petr Tichavsk? Andrzej Cichocki, On Fast Computation of
%   Gradients for CANDECOMP/PARAFAC Algorithms, arXiv:1204.1586, 2012.
%
%   Examples:
%   X = sptenrand([5 4 3], 10);
%   P = cp_als(X,2);
%   P = cp_als(X,2,'dimorder',[3 2 1]);
%   P = cp_als(X,2,'dimorder',[3 2 1],'init','nvecs');
%   U0 = {rand(5,2),rand(4,2),[]}; %<-- Initial guess for factors of P
%   [P,U0,out] = cp_als(X,2,'dimorder',[3 2 1],'init',U0);
%   P = cp_als(X,2,out.params); %<-- Same params as previous run
%
%   See also KTENSOR, TENSOR, SPTENSOR, TTENSOR.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) 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

%% Extract number of dimensions and norm of X.
N = ndims(X);
normX = norm(X);

%% Set algorithm parameters from input or by using defaults
params = inputParser;
params.addParamValue('tol',1e-4,@isscalar);
params.addParamValue('maxiters',50,@(x) isscalar(x) & x > 0);
params.addParamValue('dimorder',1:N,@(x) isequal(sort(x),1:N));
params.addParamValue('init', 'random', @(x) (iscell(x) || ismember(x,{'random','nvecs'})));
params.addParamValue('printitn',1,@isscalar);
params.parse(varargin{:});

%% Copy from params object
fitchangetol = params.Results.tol;
maxiters = params.Results.maxiters;
dimorder = params.Results.dimorder;
init = params.Results.init;
printitn = params.Results.printitn;

%% Error checking

%% Set up and error checking on initial guess for U.
if iscell(init)
    Uinit = init;
    if numel(Uinit) ~= N
        error('OPTS.init does not have %d cells',N);
    end
    for n = dimorder(2:end);
        if ~isequal(size(Uinit{n}),[size(X,n) R])
            error('OPTS.init{%d} is the wrong size',n);
        end
    end
else
    % Observe that we don't need to calculate an initial guess for the
    % first index in dimorder because that will be solved for in the first
    % inner iteration.
    if strcmp(init,'random')
        Uinit = cell(N,1);
        for n = dimorder(2:end)
            Uinit{n} = rand(size(X,n),R);
        end
    elseif strcmp(init,'nvecs') || strcmp(init,'eigs')
        Uinit = cell(N,1);
        for n = dimorder(2:end)
            Uinit{n} = nvecs(X,n,R);
        end
    else
        error('The selected initialization method is not supported');
    end
end

%% Set up for iterations - initializing U and the fit.
U = Uinit;
fit = 0;

if printitn>0
  fprintf('\nCP_ALS:\n');
end

%% Main Loop: Iterate until convergence

if (isa(X,'sptensor') || isa(X,'tensor')) && (exist('cpals_core','file') == 3)

    %fprintf('Using C++ code\n');
    [lambda,U] = cpals_core(X, Uinit, fitchangetol, maxiters, dimorder);
    P = ktensor(lambda,U);

else

    UtU = zeros(R,R,N);
    for n = 1:N
        if ~isempty(U{n})
            UtU(:,:,n) = U{n}'*U{n};
        end
    end

    for iter = 1:maxiters

        fitold = fit;

        % Iterate over all N modes of the tensor
        for n = dimorder(1:end)

            % Calculate Unew = X_(n) * khatrirao(all U except n, 'r').
            Unew = mttkrp(X,U,n);

            % Compute the matrix of coefficients for linear system
            Y = prod(UtU(:,:,[1:n-1 n+1:N]),3);
            Unew = Unew / Y;
            if issparse(Unew)
                Unew = full(Unew);   % for the case R=1
            end

            % Normalize each vector to prevent singularities in coefmatrix
            if iter == 1
                lambda = sqrt(sum(Unew.^2,1))'; %2-norm
            else
                lambda = max( max(abs(Unew),[],1), 1 )'; %max-norm
            end

            Unew = bsxfun(@rdivide, Unew, lambda');

            U{n} = Unew;
            UtU(:,:,n) = U{n}'*U{n};
        end

        P = ktensor(lambda,U);
        if normX == 0
            fit = norm(P)^2 - 2 * innerprod(X,P);
        else
            normresidual = sqrt( normX^2 + norm(P)^2 - 2 * innerprod(X,P) );
            fit = 1 - (normresidual / normX); %fraction explained by model
        end
        fitchange = abs(fitold - fit);

        % Check for convergence
        if (iter > 1) && (fitchange < fitchangetol)
            flag = 0;
        else
            flag = 1;
        end

        if (mod(iter,printitn)==0) || ((printitn>0) && (flag==0))
            fprintf(' Iter %2d: f = %e f-delta = %7.1e\n', iter, fit, fitchange);
        end

        % Check for convergence
        if (flag == 0)
            break;
        end
    end
end

%% Clean up final result
% Arrange the final tensor so that the columns are normalized.
P = arrange(P);
% Fix the signs
P = fixsigns(P);

if printitn>0
    if normX == 0
        fit = norm(P)^2 - 2 * innerprod(X,P);
    else
        normresidual = sqrt( normX^2 + norm(P)^2 - 2 * innerprod(X,P) );
        fit = 1 - (normresidual / normX); %fraction explained by model
    end
  fprintf(' Final f = %e \n', fit);
end

output = struct;
output.params = params.Results;
output.iters = iter;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function ok = all_set_item(item)

% function ok = all_set_item(item)
% Perform within-item all_set check. For generic items, this is the same
% as all_set.
%
% This code is part of a batch job configuration system for MATLAB. See
%      help matlabbatch
% for a general overview.
%_______________________________________________________________________
% Copyright (C) 2007 Freiburg Brain Imaging

% Volkmar Glauche
% $Id: all_set_item.m 1716 2008-05-23 08:18:45Z volkmar $

rev = '$Rev: 1716 $'; %#ok
ok = all_set(item);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% B-A Scale-Free Network Generation and Visualization
% *By Mathew Neil George*

%% Description and Cautions
% The *SFNG* m-file is used to simulate the B-A algorithm and returns scale-free
% networks of given node sizes. Understanding the B-A algorithm is key
% to using this code to its fullest. Due to Matlab resource limitations, it may not be
% possible to generate networks much larger than 15000 nodes, and increasing the
% *mlinks* variable increases processing time severely. This code was
% developed so that one could generate a network of small size, and then
% use that network as a seed to build a greater sized network, continuing
% this process until the actual desired network size is reached. This is for
% processor and time management purposes. However, realize that the initial
% seed does not have to have scale-free properties, while the later seeds
% may happen to have these properties. Therefore, it is prudent not to make the
% initial seed size much larger than a few nodes (most commonly 5
% interconnected nodes). In addition, the *mlinks* should be kept constant
% throughout the creation of the scale-free network.
%
% The *PLplot* m-file takes a scale-free network in adjacency matrix format
% and draws a best fit line to the frequency of degrees distribution of the
% nodes. Degree is the number of links that connect to and from a single node
% For scale-free networks, the frequency of degrees distribution forms a
% power-law curve, with an exponent usually between -2 and -3. This code is
% designed to allow only non-zero frequencies to be graphed in log-log format.
% The function returns the equation of the power-law fit in a cfit variable.
%
% The *CNet* m-file function creats a network graph using the *gplot*
% function with circular coordinates. It allows for a simple, yet
% intuitive, visualization of a given network.

%% Parameters
% *SFNG*
%
% * *Nodes* is the desired network size, including the seed network size
%   (i.e. Nodes minus seed network size equals the number of nodes to be
%   added).
%
% * *mlinks* controls the number of links a new node can make to the existing
%   network nodes.
%
% * *seed* is the original network to which the B-A algorithm links
%   additional nodes with a specific preferential attachment procedure.
%   This undirected adjacency matrix can be created manually, or one could
%   use the *Adjacency Matrix GUI*. Each node must have at least one link.
%   The *seed* variable can be replaced with a developed scale-free network
%   to generate a larger one. Make sure the new *Nodes* variable is greater
%   than the size of the *seed* network.
%
% *PLplot*
%
% * *Net* is the input network which is to be graphed.
%
% *CNet*
%
% * *Net* is the input network which is to be graphed.
%
% Note that variables *Nodes*, *mlinks*, and *size* must be whole numbers and
% variables *seed* and *Net* must be undirected adjacency matrices. The
% diagonol elements of any adjacency matrix used with these functions must
% all be zero.

%% Sample Output
% Here is a small example to demonstrate how to use the code. This code
% creates a seed network of 5 nodes, generates a scale-free network of 300 nodes from
% the seed network, and then performs the two graphing procedures.
seed =[0 1 0 0 1;1 0 0 1 0;0 0 0 1 0;0 1 1 0 0;1 0 0 0 0]
Net = SFNG(300, 1, seed);
PL_Equation = PLplot(Net)
CNet(Net)

%% References
% One explanation of the B-A Algorithm can be found on this PDF website
% http://arxiv.org/PS_cache/cond-mat/pdf/0107/0107420.pdf
%
% Undirected Adjecency Matrices are defined on Wikipedia.org
% http://en.wikipedia.org/wiki/Adjacency_matrix
%
% The *Adjacency Matrix GUI* file by Steve Chuang can be found on the Matlab File Exchange
% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=6937&objectType=file

%% Acknowledgements
% Special thanks to Mark Ballerini with the Massapequa High School Science
% Research Program and Zoltan Dezso at the University of Notre Dame for
% their invaluable help in researching network theory as well as to my
% family for providing motivation and encouragement in pursuing science.


text:
function results = vl_test_cummax(varargin)
% VL_TEST_CUMMAX
vl_test_init ;

function test_basic()
vl_assert_almost_equal(...
  vl_cummax(1), 1) ;
vl_assert_almost_equal(...
  vl_cummax([1 2 3 4], 2), [1 2 3 4]) ;

function test_multidim()
a = [1 2 3 4 3 2 1] ;
b = [1 2 3 4 4 4 4] ;
for k=1:6
  dims = ones(1,6) ;
  dims(k) = numel(a) ;
  a = reshape(a, dims) ;
  b = reshape(b, dims) ;
  vl_assert_almost_equal(...
    vl_cummax(a, k), b) ;
end

function test_storage_classes()
types = {@double, @single, ...
         @int32, @uint32, ...
         @int16, @uint16, ...
         @int8, @uint8} ;
if vl_matlabversion() > 71000
  types = horzcat(types, {@int64, @uint64}) ;
end
for a = types
  a = a{1} ;
  for b = types
    b = b{1} ;
    vl_assert_almost_equal(...
      vl_cummax(a(eye(3))), a(toeplitz([1 1 1], [1 0 0 ]))) ;
  end
end

repo: StructCompletion


text:
function [ Y ] = TMAV( X, h )
% TMED  - The fuzzy filter with Triangular function with MEDian center
% x     - input (noisy image)
% y     - output (de-noised image)
% F     - window function, weighting function

if nargin < 2
    h=3;    % local search window size: h=3; h=5; h=7; h=9; h=11;
    f=1;    % padding value for window: f=1; f=2; f=3; f=4; f=5;
else
    f=(h-1)/2;    % padding value for window: f=1; f=2; f=3; f=4; f=5;
end
F = zeros(h*h,1);
Y  = zeros(size(X));

X = double(padarray(X,[f f],'symmetric'));
[m n ~] = size(X);

for i=1+f:1:m-f
   for j=1+f:1:n-f
       x=reshape(X(i-f:i+f, j-f:j+f),[],1);
       xmax = max(x);
       xmin = min(x);
       xmav = mean(x);
       xmv = max(xmax-xmav, xmav-xmin);

       F(:,:) = 0;
       if xmv==0
           F(:,:) = 1;
       else
           F = (abs(x-xmav)/xmv);
       end

       Y(i-f,j-f) = sum(sum(F.*x))/sum(sum(F));
       clear xmax xmin xmed xmm;

   end
end
end


text:
function i4vec_index_delete_all_test ( )

%*****************************************************************************80
%
%% I4VEC_INDEX_DELETE_ALL_TEST tests I4VEC_INDEX_DELETE_ALL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  n_max = 25;
  n = 0;
  x = [];
  indx = [];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_INDEX_DELETE_ALL_TEST\n' );
  fprintf ( 1, '  I4VEC_INDEX_DELETE_ALL deletes all copies of a\n' );
  fprintf ( 1, '  particular value.\n' );

  xval = 8;
  [ n, x, indx ] = i4vec_index_insert ( n, x, indx, xval );

  xval = 7;
  [ n, x, indx ] = i4vec_index_insert ( n, x, indx, xval );

  b = 0;
  c = 20;
  seed = 123456789;

  for i = 1 : 20
    [ xval, seed ] = i4_uniform_ab ( b, c, seed );
    fprintf ( 1, '  %6d\n', xval );
    [ n, x, indx ] = i4vec_index_insert ( n, x, indx, xval );
  end

  xval = 7;
  [ n, x, indx ] = i4vec_index_insert ( n, x, indx, xval );

  xval = 8;
  [ n, x, indx ] = i4vec_index_insert ( n, x, indx, xval );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Indexed list of entries:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I  INDX(I)  X(I)  X(INDX(I))\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %6d  %6d  %6d  %6d\n', i, indx(i), x(i), x(indx(i)) );
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Call I4VEC_INDEX_DELETE_ALL to delete values of 7:\n' );

  xval = 7;
  [ n, x, indx ] = i4vec_index_delete_all ( n, x, indx, xval );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Indexed list of entries:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  I  INDX(I)  X(I)  X(INDX(I))\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %6d  %6d  %6d  %6d\n', i, indx(i), x(i), x(indx(i)) );
  end

  return
end


text:
function [data_float, fs]=aurora2read(fname)

% Read data from the Aurora2 database

fid=fopen(fname,'r','b');
data=fread(fid,'int16');
fclose(fid);
fs=8000;

%str1=strread(fname,'%s','delimiter','.');

data_float = double(data)/2^15; %% Normalize int16(y) by 2^15

% wavwrite(data_float,fs,strcat(str1{1},'.wav'));

author: zhenghuatan

repo: rVAD

save_path: github-repos/MATLAB/zhenghuatan-rVAD


text:
function [model] = rvm_train(data, labels, options)
% RVM_TRAIN Trains an RVM Model using SB1 Tipping Toolbox
%
%   input ----------------------------------------------------------------
%
%       o data        : (N x D), N data points of D dimensionality.
%
%       o labels      : (N x 1), Either 1, -1 for binary
%
%       o options     : struct
%
%
%   output ----------------------------------------------------------------
%
%       o model       : struct.
%
%
%% %    RVM OPTIONS
%       ALPHA   Scalar initial value for hyperparameters
%       BETA    Initial value for inverse noise variance (in regression)
%               Set this negative to fix the value, rather than estimate
%       KERNEL  Kernel type: see SB1_KERNELFUNCTION for options
%       LEN     Kernel length scale
%       USEBIAS Set to non-zero to utilise a "bias" offset
%       MAXITS  Maximum iterations to run for.

%% Parse RVM Options
% Parsing Parameter for RVM
N	= length(data);
useBias = options.useBias;
kernel_	= options.kernel_;
width   = options.width;
maxIts	= options.maxIts;
monIts	= round(maxIts/20);

% Set verbosity of output (0 to 4)
setEnvironment('Diagnostic','verbosity',3);
% Set file ID to write to (1 = stdout)
setEnvironment('Diagnostic','fid',1);

%% Set up initial hyperparameters - precise settings should not be critical
initAlpha	= (1/N)^2;
% Set beta to zero for classification
initBeta	= 0;
% Check that labels are 0 for negative class
labels(find(labels==-1)) = 0;

%% Train RVM Model

% "Train" a sparse Bayes kernel-based model (relevance vector machine)
[weights, used, bias, marginal, alpha, beta, gamma] = ...
    SB1_RVM(data,labels,initAlpha,initBeta,kernel_,width,useBias,maxIts,monIts);

%% Model Parameters
%       WEIGHTS Parameter values of estimated model (sparse)
%       USED    Index vector of "relevant" kernels (data points)
%       BIAS    Value of bias or offset parameter
%       ML      Log marginal likelihood of model
%       ALPHA   Estimated hyperparameter values (sparse)
%       BETA    Estimated inverse noise variance for regression
%       GAMMA   "Well-determinedness" factors for relevant kernels

model.weights  = weights;
model.kernel_  = kernel_;
model.width    = width;
model.RVs      = data(used,:);
model.bias     = bias;
model.marginal = marginal;
model.alpha    = alpha;
model.beta     = beta;
model.gamma    = gamma;

end

author: epfl-lasa


text:
% OtMatrixCubeSet Put a matrix into a matrix cube.
%
% Usage:
%    OtMatrixCubeSet(hmcube, pmtu, matrix)
% where:
%    hmcube matrix cube handle as created by OtMatrixCubeNew()
%    pmtu   is a PMTU struct as created by OtPmtu(),
%    matrix is the matrix to store in the cube
%

% OtMatlab: Matlab interface to the Omnitrans transport planning software

% Copyright (c) 2008, VORtech Computing
% Copyright (c) 2010, Omnitrans International

% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:

%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%     * Neither the name of VORtech Computing and Omnitrans International nor the names
%       of its contributors may be used to endorse or promote products derived
%       from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.

% 2006-04-19 Created (Johan Meijdam, VORtech)

function iret = OtMatrixCubeSet(hmcube, pmtu, matrix)
	iret = Ot('matrixcube.set', hmcube, pmtu, matrix);
end


text:
function response = Waffely(n,p)

persistent oscyflag
if n==1; oscyflag=1; end
if oscyflag
    response = 'cooperate';
    oscyflag=0;
else
    response = 'defect';
    oscyflag=1;
end


text:
%  演示调用 平方根法 分解
clear;
A=[6 7 5; 7 13 8; 5 8 6];   %对称正定矩阵
b=[9; 10; 9];

% 平方根分解法
[L, x] = my_cholesky(A, b)

% 直接调用matlab左除运算符
x_ = A\b

repo: NumericalAnalysis


text:
classdef TestMatchGMS
    %TestMatchGMS

    methods (Static)
        function test_1
            im1 = imread(fullfile(mexopencv.root(),'test','left01.jpg'));
            im2 = imread(fullfile(mexopencv.root(),'test','right01.jpg'));
            sz1 = [size(im1,2), size(im1,1)];
            sz2 = [size(im2,2), size(im2,1)];
            obj = cv.ORB('MaxFeatures',5000, 'FastThreshold',0);
            [kpts1, feat1] = obj.detectAndCompute(im1);
            [kpts2, feat2] = obj.detectAndCompute(im2);
            matcher = cv.DescriptorMatcher('BruteForce-Hamming');
            m = matcher.match(feat1, feat2);

            mGMS = cv.matchGMS(sz1, kpts1, sz2, kpts2, m);
            validateattributes(mGMS, {'struct'}, {'vector'});
            assert(all(ismember(TestDescriptorMatcher.fields, ...
                fieldnames(mGMS))));
        end

        function test_error_argnum
            try
                cv.matchGMS();
                throw('UnitTest:Fail');
            catch e
                assert(strcmp(e.identifier,'mexopencv:error'));
            end
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function V = variablesPack(X, I, U, P, H)

N = length(H);
nx = size(X,1);
ni = size(I,1);
np = size(P,1);

nv = nvars(N, nx, ni, nu, np);

[X_indizes, I_indizes, U_indizes, P_indizes, H_indizes] = ocl.simultaneous.indizes(N, nx, ni, nu, np);

V = zeros(nv, 1);
V(X_indizes) = X;
V(I_indizes) = I;
V(U_indizes) = U;
V(P_indizes) = P;
V(H_indizes) = H;

author: OpenOCL

repo: OpenOCL

save_path: github-repos/MATLAB/OpenOCL-OpenOCL


text:
function K = simwhiteXsimwhiteKernCompute(simKern1, simKern2, t1, t2)

% SIMWHITEXSIMWHITEKERNCOMPUTE Compute a cross kernel between two SIM-WHITE
% kernels.
% FORMAT
% DESC computes cross kernel terms between two SIM-WHITE kernels for
% the multiple output kernel.
% ARG simKern1 : the kernel structure associated with the first SIM-WHITE
% kernel.
% ARG simKern2 : the kernel structure associated with the second SIM-WHITE
% kernel.
% ARG t1 : inputs for which kernel is to be computed.
% RETURN K : block of values from kernel matrix.
%
% FORMAT
% DESC computes cross kernel terms between two SIM-WHITE kernels for
% the multiple output kernel.
% ARG simKern1 : the kernel structure associated with the first SIM-WHITE
% kernel.
% ARG simKern2 : the kernel structure associated with the second SIM-WHITE
% kernel.
% ARG t1 : row inputs for which kernel is to be computed.
% ARG t2 : column inputs for which kernel is to be computed.
% RETURN K : block of values from kernel matrix.
%
% SEEALSO : multiKernParamInit, multiKernCompute, simwhiteKernParamInit
%
% COPYRIGHT : David Luengo, 2009

% KERN

if nargin < 4
  t2 = t1;
end
if size(t1, 2) > 1 | size(t2, 2) > 1
  error('Input can only have one column');
end
if simKern1.variance ~= simKern2.variance
  error('Kernels cannot be cross combined if they have different variances.')
end

% Parameters of the kernels required in the computation
variance = simKern1.variance;
sensitivity1 = simKern1.sensitivity;
sensitivity2 = simKern2.sensitivity;
decay1 = simKern1.decay;
decay2 = simKern2.decay;

isStationary = (simKern1.isStationary == true) & (simKern2.isStationary == true);

% Auxiliary constants and matrices
c = variance * sensitivity1 * sensitivity2 / (decay1 + decay2);
T1 = repmat(t1, 1, size(t2, 1));
T2 = repmat(t2.', size(t1, 1), 1);
ind = (T1 < T2);
Dv = decay2 .* ind + decay1 .* (~ind);
K = exp(-Dv .* abs(T1-T2));
if (isStationary == false)
    K = K - exp(-(decay1 * T1 + decay2 * T2));
end
K = c*K;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function out=mp(x,y)
%MP multiple precision class constructor.
%   p = mp(x,y) creates a mp object from the matrices x and y,
%   where x contains the double to be converted into an mp object
%         y contains the precision
% Special calls to mp include:
%   mp('pi',precision) => returns pi to precision (precision is optional)

maxDoublePrec=16; %digits
if nargin == 0
 out.rval=[];
 out.ival=[];
 out.precision=[];
 out=class(out,'mp');
elseif isa(x,'mp')
 %out=x;
 if nargin==1
  out=x;
 else
  for ii=1:numel(x)
   out(ii)=mp(x(ii).rval,y(min(numel(y),ii)));
   if ~isreal(x(ii))
    out(ii)=out(ii)+mp(x(ii).ival,y(min(numel(y),ii)))*i;
   end
  end
 end % if nargin==2
else
 mp_defaults
 precision=default_precision;
 out_rval=cell(size(x));
 out_ival=cell(size(x));
 if nargin==2
  precision=double(y(1));
 end % if nargin==2
 if isa(x,'double')
  for ii=1:numel(x)
   [str,exponent]=mpfr_construct_dd(real(x(ii)),precision);
   % throw away anything past maxDoublePrec, set to 0's
   if length(str)>maxDoublePrec
    str(maxDoublePrec+1:end)='0';
   end
   out_rval{ii}=mpExpForm(mpAddDecimal(str),exponent);
   if ~isreal(x(ii))
    [str,exponent]=mpfr_construct_dd(imag(x(ii)),precision);
    % throw away anything past maxDoublePrec, set to 0's
    if length(str)>maxDoublePrec
     str(maxDoublePrec+1:end)='0';
    end
    out_ival{ii}=mpExpForm(mpAddDecimal(str),exponent);
   end
  end % for ii=1:size(x,
 elseif isa(x,'cell')
  for ii=1:numel(x)
   [str,exponent]=mpfr_construct_cd(x{ii},precision);
   out_rval{ii}=mpExpForm(mpAddDecimal(str),exponent);
  end % for ii=1:size(x,
 elseif isa(x,'char')
  out_rval=cell(1);
  out_ival=cell(1);
  if any(strfind(lower(x),'pi'))
   out_rval=mpfr_pi(precision);
%%%   [str,exponent]=mpfr_pi(precision);
%%%   out_rval=mpExpForm(mpAddDecimal(str),exponent);
  else
   [str,exponent]=mpfr_construct_cd(x,precision);
   out_rval{1,1}=mpExpForm(mpAddDecimal(str),exponent);
  end
 end
 out=class(struct('rval',out_rval,...
                  'ival',out_ival,...
                  'precision',precision),'mp');
end

%%%out_rval
%%%out_rexp
%%%out_ival
%%%out_iexp

%'rrrrrrr',kb
% x=magic(3),s1='023e1',s2='-.002e-1'
% mp(x)

author: opencobra

repo: cobratoolbox


text:
function f = crowding_distance(x,problem)
% This function calculates the crowding distance

%
%  Copyright (c) 2009, Aravind Seshadri
%  All rights reserved.
%
%  Redistribution and use in source and binary forms, with or without
%  modification, are permitted provided that the following conditions are
%  met:
%
%     * Redistributions of source code must retain the above copyright
%       notice, this list of conditions and the following disclaimer.
%     * Redistributions in binary form must reproduce the above copyright
%       notice, this list of conditions and the following disclaimer in
%       the documentation and/or other materials provided with the distribution
%
%  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%  POSSIBILITY OF SUCH DAMAGE.

[N,M] = size(x);
switch problem
    case 1
        M = 2;
        V = 6;
    case 2
        M = 3;
        V = 12;
end

% Crowding distance for each front
for i = 1 : length(F(front).f)
    y(i,:) = x(F(front).f(i),:);
end
for i = 1 : M
    [sorted(i).individual,sorted(i).index] = sort(y(:,V + i));
    distance(sorted(i).index(1)).individual = Inf;
    distance(sorted(i).index(length(sorted(i).index))).individual = Inf;
end

[num,len] = size(y);
% Initialize all the distance of individuals as zero.
for i = 1 : M
    for j = 2 : num - 1
        distance(j).individual = 0;
    end
    objective(i).range = ...
                sorted(i).individual(length(sorted(i).individual)) - ...
                sorted(i).individual(1);
        % Maximum and minimum objectives value for the ith objective
end
% Caluclate the crowding distance for front one.
for i = 1 : M
    for j = 2 : num - 1
        distance(j).individual = distance(j).individual + ...
            (sorted(i).individual(j + 1) - sorted(i).individual(j - 1))/...
            objective(i).range;
        y(sorted(i).index(j),M + V + 2) = distance(j).individual;
    end
end


text:
function [output, firstIndexPoint] = read_hs_file( filename, outfile)

%read_hs_file Reads in BTI-Headshape files
%   filename: file with the headshape informations
%   outfile: if present, a ctf ".shape" file is written
%   output: if present, a 3xN matrix containing the headshape-points
%
%   (C) 2007 by Thomas Hartmann

% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
% for the documentation and details.
%
%    FieldTrip is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    FieldTrip is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id: read_bti_hs.m 945 2010-04-21 17:41:20Z roboos $

if nargin == 1
    outfile = [];
end %if

fid = fopen(filename, 'r', 'b');
version = fread(fid, 1, '*uint32');
timestamp = fread(fid, 1, '*int32');
checksum = fread(fid, 1, '*int32');
nPoints = fread(fid, 1, '*int32');

firstIndexPoint = fread(fid, [3, 5], 'double')';

points = fread(fid, [3, double(nPoints)], 'double');

fclose(fid);

if(nargout > 0)
    output = points';
end %if

if(nargin == 2)
    fid = fopen(outfile, 'wt');
    fprintf(fid, '%d\n', nPoints);
    for i = 1:size(points, 2)
        fprintf(fid, '%.3f\t%.3f\t%.3f\n', points(1, i), points(2, i), points(3, i));
    end %for
    fclose(fid);

end %if

author: PatternRecognition

repo: OpenBMI


text:
function [Ft3]=tri6_subtri3(Ft6,Vt6)

% function [Ft3]=tri6_subtri3(Ft6,Vt6)
% ------------------------------------------------------------------------
% The input 6-node triangles are either oriented such that the first node
% is a corner node or that the first node is a mid-edge node. The two cases
% are highlighted below.
%
% Case I
%                     1
%                    / \
%                   2   6
%                  /     \
%                 3___4___5
% Case II
%                     6
%                    / \
%                   5   1
%                  /     \
%                 4___3___2
%
% The output will follow one of these rules and maintains face normals
% directions:
%
% Case I
%                     1
%                    / \
%                   2___6
%                  / \ / \
%                 3___4___5
%
% Case II
%                     6
%                    / \
%                   5___1
%                  / \ / \
%                 4___3___2
%
% ------------------------------------------------------------------------
%%
%Check face order based on first
a=Vt6(Ft6(1,2),:)-Vt6(Ft6(1,1),:); %2 1 edge
b=Vt6(Ft6(1,6),:)-Vt6(Ft6(1,1),:); %6 1 edge
c=Vt6(Ft6(1,3),:)-Vt6(Ft6(1,2),:); %3 2 edge

am=sqrt(sum(a.^2));
bm=sqrt(sum(b.^2));
cm=sqrt(sum(c.^2));

phi=abs(acos(dot(a,b)./(am*bm)));
phi=mod(phi,pi); %Angle between 2 1 and 6 1 edge
theta=abs(acos(dot(a,c)./(am*cm)));
theta=mod(theta,pi); %Angle between 2 1 and 3 2 edge

if phi>theta
    Ft3=[Ft6(:,[1 2 6]); Ft6(:,[2 3 4]); Ft6(:,[4 5 6]); Ft6(:,[2 4 6])]; %Clockwise
else
    Ft3=[Ft6(:,[1 2 3]); Ft6(:,[3 4 5]); Ft6(:,[5 6 1]); Ft6(:,[1 3 5])]; %Anti-clockwise
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function X = redo_scaling(X_scal,param)

% redo scaling (from scaled data to original data)
%
% INPUT
% X_scal:   pretreated data matrix (samples x variables)
% param:    output data structure from data_pretreatment routine
%
% OUTPUT
% X:        data matrix (samples x variables)
%
% version 1.0 - september 2009
% Davide Ballabio
% Milano Chemometrics and QSAR Research Group
% www.disat.unimib.it/chm

a = param.a;
s = param.s;
m = param.m;
M = param.M;
pret_type = param.pret_type;

if strcmp(pret_type,'cent')
    for j=1:size(X_scal,2)
        X(:,j) = X_scal(:,j) + a(j);
    end
elseif strcmp(pret_type,'scal')
    for j=1:size(X_scal,2)
        X(:,j) = X_scal(:,j)*s(j);
    end
elseif strcmp(pret_type,'auto')
    for j=1:size(X_scal,2)
        X(:,j) = X_scal(:,j)*s(j) + a(j);
    end
elseif strcmp(pret_type,'rang')
    for j=1:size(X_scal,2)
        X(:,j) = X_scal(:,j)*(M(j) - m(j)) + m(j);
    end
else
    X = X_scal;
end

author: kmansouri

repo: OPERA

save_path: github-repos/MATLAB/kmansouri-OPERA


text:
function tagS = gsps_module_tags
%"gsps_module_tags"
%   Return the tags used to represent a ROI contour specified by
%   A.33 in PS3.3 of 2006 DICOM specification. Tags are returned in a
%   struct array with 3 fields:
%
%   Tag: String containing hex DICOM tag of a field.
%  Type: String describing type of field, with 5 options:
%         '1' Field must exist, data must exist and be valid.
%         '2' Field must exist, data can exist and be valid, or be NULL.
%         '3' Field is optional, if the field exists data can exist and be
%             valid or be NULL.
%         '1C' Field must exist under certain conditions and contain valid
%             data.
%         '2C' Field must exist under certain conditions and can contain
%             valid data or be NULL.
%Children: For sequences, a tagS with the same format as this struct array
%          containing the tags of child fields.
%
%APA 07/12/2019
%
%Usage:
%   tagS = ROI_contour_module_tags
%
% 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 tagS structure.
tagS = struct('tag', {}, 'tagdec', {}, 'type', {}, 'children', {});

%Create an empty tagS template for sequence creation.
template = tagS;

%Add tags based on PS3.3 attribute lists

% Presentation Creation Date
tagS(end+1) = struct('tag', '00700082', 'tagdec', 7340162, 'type', '1', 'children', []);

% Presentation Creation Time
tagS(end+1) = struct('tag', '00700083', 'tagdec', 7340163, 'type', '1', 'children', []);

% Presentation LUT Shape
tagS(end+1) = struct('tag', '20500020', 'tagdec', 542113824, 'type', '1C', 'children', []);

% Displayed Area Selection Sequence
tagS(end+1) = struct('tag', '0070005A', 'tagdec', 7340122, 'type', '1', 'children', []);

    child_0 = template;
    %Referenced Image Sequence
    child_0(end+1) = struct('tag', '00081140', 'tagdec', 528704, 'type', '1C', 'children', image_SOP_instance_reference_macro_tags);
    % Displayed Area Top Left Hand Corner
    child_0(end+1) = struct('tag', '00700052', 'tagdec', 7340114, 'type', '1', 'children', []);
    % Displayed Area Bottom Right Hand Corner
    child_0(end+1) = struct('tag', '00700053', 'tagdec', 7340115, 'type', '1', 'children', []);
    % Presentation Size Mode
    child_0(end+1) = struct('tag', '00700100', 'tagdec', 7340288, 'type', '1', 'children', []);
    % Presentation Pixel Spacing
    child_0(end+1) = struct('tag', '00700101', 'tagdec', 7340289, 'type', '1C', 'children', []);
    % Presentation Pixel Aspect Ratio
    child_0(end+1) = struct('tag', '00700102', 'tagdec', 7340290, 'type', '1C', 'children', []);
    % Presentation Pixel Magnification Ratio
    child_0(end+1) = struct('tag', '00700103', 'tagdec', 7340291, 'type', '1C', 'children', []);

tagS(end).children = child_0;

% Referenced Series Sequence
tagS(end+1) = struct('tag', '00081115', 'tagdec', 528661, 'type', '1', 'children', []);
    child_0 = template;
    % Series Instance UID
    child_0(end+1) = struct('tag', '0020000E', 'tagdec', 2097166, 'type', '1', 'children', []);
    % Referenced Image Sequence
    child_0(end+1) = struct('tag', '00081140', 'tagdec', 528704, 'type', '1', 'children', image_SOP_instance_reference_macro_tags);

tagS(end).children = child_0;

% Referenced Study Sequence

%Graphic Annotation Sequence
tagS(end+1) = struct('tag', '00700001', 'tagdec', 7340033, 'type', '1', 'children', []);

child_0 = template;

%Referenced Image Sequence
child_0(end+1) = struct('tag', '00081140', 'tagdec', 528704, 'type', '1C', 'children', image_SOP_instance_reference_macro_tags);

%Graphic layer
child_0(end+1) = struct('tag', '00700002', 'tagdec', 7340034, 'type', '1', 'children', []);

%Text object sequence
child_0(end+1) = struct('tag', '00700008', 'tagdec', 7340040, 'type', '1C', 'children', []);

        child_1 = template;

        %Bounding Box Annotation Units
        %child_1(end+1) = struct('tag', '00700003', 'tagdec',  , 'type', '1C', 'children', []);

        %Anchor Point Annotation Units
        child_1(end+1) = struct('tag', '00700004', 'tagdec', 7340036, 'type', '1C', 'children', []);

        %Unformatted Text Value
        child_1(end+1) = struct('tag', '00700006', 'tagdec', 7340038 , 'type', '1', 'children', []);

        %Bounding Box Top Left Hand Corner
        %child_1(end+1) = struct('tag', '00700010', 'tagdec', 7340048, 'type', '1C', 'children', []);

        %Bounding Box Bottom Right Hand Corner
        %child_1(end+1) = struct('tag', '00700011', 'tagdec', 7340049, 'type', '1C', 'children', []);

        %Bounding Box Text Horizontal Justification
        %child_1(end+1) = struct('tag', '00700012', 'tagdec', 7340050, 'type', '1C', 'children', []);

        %Anchor Point
        child_1(end+1) = struct('tag', '00700014', 'tagdec', 7340052, 'type', '1C', 'children', []);

        %Anchor Point Visibility
        child_1(end+1) = struct('tag', '00700015', 'tagdec', 7340053, 'type', '1C', 'children', []);

child_0(end).children = child_1;

%Graphic object sequence
child_0(end+1) = struct('tag', '00700009', 'tagdec', 7340041, 'type', '1', 'children', []);

        child_2 = template;

        %Graphic Annotation Units
        child_2(end+1) = struct('tag', '00700005', 'tagdec', 7340037, 'type', '1', 'children', []);

        %Graphic Dimensions
        child_2(end+1) = struct('tag', '00700020', 'tagdec', 7340064, 'type', '1', 'children', []);

        %Number of Graphic Points
        child_2(end+1) = struct('tag', '00700021', 'tagdec', 7340065, 'type', '1', 'children', []);

        %Graphic Data
        child_2(end+1) = struct('tag', '00700022', 'tagdec', 7340066, 'type', '1', 'children', []);

        %Graphic Type
        child_2(end+1) = struct('tag', '00700023', 'tagdec', 7340067, 'type', '1', 'children', []);

        %Graphic Filled
        child_2(end+1) = struct('tag', '00700024', 'tagdec', 7340068, 'type', '1C', 'children', []);

child_0(end).children = child_2;

tagS(end).children = child_0;

%Graphic Layer Sequence
tagS(end+1) = struct('tag', '00700060', 'tagdec', 7340128, 'type', '1', 'children', []);
        child_1 = template;

        %0070,0002  Graphic Layer
        child_1(end+1) = struct('tag', '00700002', 'tagdec', 7340034, 'type', '1', 'children', []);

        %0070,0062  Graphic Layer Order
        child_1(end+1) = struct('tag', '00700062', 'tagdec', 7340130, 'type', '1', 'children', []);

        %0070,0066  Graphic Layer Recommended Display Grayscale Value
        child_1(end+1) = struct('tag', '00700066', 'tagdec', 7340134, 'type', '3', 'children', []);

        %0070,0401  Graphic Layer Recommended Display CIELab Value
        child_1(end+1) = struct('tag', '00700401', 'tagdec', 7341057, 'type', '3', 'children', []);

        %0070,0068  Graphic Layer Description
        child_1(end+1) = struct('tag', '00700068', 'tagdec', 7340136, 'type', '3', 'children', []);

tagS(end).children = child_1;

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function tcell_flow_png ( )

%*****************************************************************************80
%
%% TCELL_FLOW_PNG plots TCELL data.
%
%  Discussion:
%
%    This MATLAB script file reads the T-Cell flow data:
%
%      geometry (XY values at nodes, assumed to be in 'xy.txt')
%      flow (UV values at nodes, in a sequence of files starting with 'up001.txt')
%
%    and plots the velocity vectors (U,V)(X,Y), and saves each plot as a
%    PNG file, presumably so the PNG files can be gathered into an
%    animation.
%
%    The file plots either the velocity vector field, or the velocity
%    direction field, depending on the value of the internal logical
%    parameter "normalized".
%
%    The MATLAB routine quiver internally scales the vectors, but this can be
%    adjusted by using a value of SCALE that is not 1.
%
%    For the unnormalized case, the velocity field is computed to a fixed
%    scale determined by finding the largest velocity vector over time.
%
%    This routine requires some auxiliary routines in order to "increment"
%    the name of the current velocity data file to get the name of the next
%    one, with the main such routine being called FILE_NAME_INC.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 July 2003
%
%  Author:
%
%    John Burkardt
%
  FALSE = 0;
  TRUE = 1;
  scale = 1.0;
  normalized = TRUE;

  nframes = 500;
%
%  Get the XY coordinates of the nodes.
%
  load xy.txt;

  x = xy(:,1);
  y = xy(:,2);
%
%  Thin the data.
%
  thin_factor = 4;
  thin_dex = thin_index ( x, y, thin_factor );
  thin_num = length ( thin_dex );
  x = x(thin_dex);
  y = y(thin_dex);
%
%  For unnormalized plots, you need to make sure that a fixed scale is preserved from
%  step to step.  The only way I can see to do this requires that I determine the
%  maximum velocity over all time, and then append one extract node and velocity
%  to the data structure.
%
  if ( normalized == FALSE )

    upfile = 'up000.txt';
    vnorm_max = 0.0;

    for i = 1 : nframes
      fprintf ( 1, '  Checking file %d.\n', i );
      upfile = file_name_inc ( upfile );
      uv = load ( upfile );
      u = uv(thin_dex,1);
      v = uv(thin_dex,2);
      norm = sqrt ( u.^2 + v.^2 );
      vnorm_max = max ( vnorm_max, max ( norm ) );
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Maximum visible velocity magnitude is %f\n', vnorm_max );

  end
%
%  Set the coordinates of the boundary lines, and an extra set of
%  lines that will be invisible, but which make some space for the
%  file name to be displayed within the plot area.
%
  bx1 = [ 0.00, 1.00, 1.00, 0.00, 0.00, 0.00 ];
  by1 = [ 1.00, 1.00, 0.99, 0.99, 0.99, 1.00 ];

  bx2 = [ 0.00, 0.25, 0.25, 0.75, 0.75, 1.00, 1.00, 0.74, 0.74, 0.26, 0.26, 0.00, 0.00];
  by2 = [ 0.50, 0.50, 0.00, 0.00, 0.50, 0.50, 0.51, 0.51, 0.01, 0.01, 0.51, 0.51, 0.50];

  bx3 = [ -0.10,  1.10, 1.10, -0.10, -0.10 ];
  by3 = [ -0.10, -0.10, 1.10,  1.10, -0.10 ];
%
%  Set the name of the 0th (nonexistent) velocity file.
%  All the velocity files with have this format, with the
%  numeric part of the name incremented to get the next one.
%  In particular, the first velocity file is called UP001.TXT.
%
  upfile = 'up000.txt';
  pngfile = 'up000.png';

  for i = 1 : nframes

    fprintf ( 1, '  Converting file %d.\n', i );

    upfile = file_name_inc ( upfile );
    pngfile = file_name_inc ( pngfile );
    uv = load ( upfile );
    u = uv(thin_dex,1);
    v = uv(thin_dex,2);

    if ( normalized == TRUE )
      norm = sqrt ( u.^2 + v.^2 );
      nonzero = find ( norm ~= 0.0 );
      u(nonzero) = u(nonzero) ./ norm(nonzero);
      v(nonzero) = v(nonzero) ./ norm(nonzero);
      vnorm_max = 1.0;
    end

    x(thin_num+1) = 0.00;
    y(thin_num+1) = 1.05;
    u(thin_num+1) = vnorm_max;
    v(thin_num+1) = 0.0;

    quiver ( x, y, u, v, scale )
%
%  Draw the boundary lines, and the invisible bounding line.
%
    line ( bx1, by1, 'color', 'r' )
    line ( bx2, by2, 'color', 'r' )
    line ( bx3, by3, 'color', 'w' )

    axis equal
    if ( normalized == TRUE )
      title ( 'T-Cell Direction Field' )
    else
      title ( 'T-Cell Flow Field' )
    end
    text ( 0.420, 1.03, upfile )
%
%  Here's where we take a snapshot of the current image, and save it to
%  a PNG file.
%
    F = getframe;
    [X,map] = frame2im ( F );
    imwrite ( X, pngfile, 'PNG' );

  end

  return
end
function file_name = file_name_inc ( file_name )

%*****************************************************************************80
%
%% FILE_NAME_INC generates the next filename in a series.
%
%  Discussion:
%
%    It is assumed that the digits in the name, whether scattered or
%    connected, represent a number that is to be increased by 1 on
%    each call.  If this number is all 9's on input, the output number
%    is all 0's.  Non-numeric letters of the name are unaffected..
%
%    If the name is empty, then the routine stops.
%
%    If the name contains no digits, the empty string is returned.
%
%  Example:
%
%      Input            Output
%      -----            ------
%      'a7to11.txt'     'a7to12.txt'  (typical case.  Last digit incremented)
%      'a7to99.txt'     'a8to00.txt'  (last digit incremented, with carry.)
%      'a9to99.txt'     'a0to00.txt'  (wrap around)
%      'cat.txt'        ' '           (no digits in input name.)
%      ' '              STOP!         (error.)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 September 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string FILE_NAME, the string to be incremented.
%
%    Output, string FILE_NAME, the incremented string.
%
  lens = length ( file_name );

  if ( lens <= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'FILE_NAME_INC - Fatal error!\n' );
    fprintf ( 1, '  The input filename is empty.\n' );
    error ( 'FILE_NAME_INC - Fatal error!' );
  end

  change = 0;

  for i = lens : -1 : 1

    c = file_name(i);

    if ( '0' <= c & c <= '8' )

      change = change + 1;

      c = c + 1;

      file_name(i) = c;

      return

    elseif ( c == '9' )

      change = change + 1;

      c = '0';

      file_name(i) = c;

    end

  end

  if ( change == 0 )
    file_name = ' ';
  end

  return
end
function thin_dex = thin_index ( x, y, thin_factor )

%*****************************************************************************80
%
%%  THIN_INDEX determines thinning indices for a X, Y data.
%
%  Discussion:
%
%    A set of X, Y data is given, that is presumably, not too far off
%    from being on a rectangular grid.
%
%    The input value of THIN_FACTOR indicates by how much the data should
%    be thinned.
%
%    The X and Y ranges are computed, and only those data points are
%    retained for which both X and Y lie in an appropriate subrange.
%
%    For instance, a THIN_FACTOR of 2 would essentially save data
%    that lay in the black squares of a checkerboard.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 June
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(NODE_NUM), Y(NODE_NUM), the X and Y coordinates
%    of the nodes.
%
%    Input, integer THIN_FACTOR, the thinning factor.
%
%    Output, integer THIN_DEX(NODE_NUM), contains in (1:THIN_NUM) the
%    indices into X and Y of the vectors to be retained after thinning.
%
  TRUE = 1;
  FALSE = 0;

  node_num = length ( x );

  x_unique_num = 0;

  for i = 1 : node_num

    unique = TRUE;

    for j = 1 : x_unique_num
      if ( x(i) == x_unique(j) )
        unique = FALSE;
        break;
      end
    end

    if ( unique )
      x_unique_num = x_unique_num + 1;
      x_unique(x_unique_num) = x(i);
    end

  end

  sort ( x_unique );

  y_unique_num = 0;

  for i = 1 : node_num

    unique = TRUE;

    for j = 1 : y_unique_num
      if ( y(i) == y_unique(j) )
        unique = FALSE;
        break;
      end
    end

    if ( unique )
      y_unique_num = y_unique_num + 1;
      y_unique(y_unique_num) = y(i);
    end

  end

  sort ( y_unique );

  thin_num = 0;

  for i = 1 : node_num

    for j = 1 : x_unique_num-1
      if ( x_unique(j) <= x(i) & x(i) <= x_unique(j+1) )
        x_bin = j;
        break;
      end
    end

    for j = 1 : y_unique_num-1
      if ( y_unique(j) <= y(i) & y(i) <= y_unique(j+1) )
        y_bin = j;
        break;
      end
    end

    if ( mod ( y_bin, thin_factor ) == thin_factor / 2 & ...
         mod ( x_bin, thin_factor ) == thin_factor / 2 )

      thin_num = thin_num + 1;
      thin_dex(thin_num) = i;

    end

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function a = c4mat_test ( n )

%*****************************************************************************80
%
%% C4MAT_TEST sets up a test matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 April 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, complex A(N,N), the Fourier matrix.
%
  complex_i = complex ( 0.0, 1.0 );

  for i = 1 : n
    for j = 1 : n

      angle = 2.0 * pi * ( i - 1 ) * ( j - 1 ) / n;

      a(i,j) = exp ( complex_i * angle ) / sqrt ( n );

    end
  end

  return
end


text:
function [curv] = mne_read_curvature(fname)
%
% [curf] = mne_read_surface(fname)
%
% Reads a FreeSurfer curvature file
%
% fname       - The file to read
% curv        - The curvature values
%

%
%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%
%   Revision 1.4  2006/07/31 05:12:39  msh
%   fread3 was still called instead of mne_fread3
%
%   Revision 1.3  2006/05/22 11:01:47  msh
%   Deleted superfluous text from the comments.
%
%   Revision 1.2  2006/05/22 10:55:02  msh
%   Fixed help text.
%
%   Revision 1.1  2006/05/22 10:44:44  msh
%   Added surface and curvature reading routines.
%   Fixed error in mne_read_source_spaces: triangle normals were not normalized
%

me='MNE:mne_read_curvature';

%
%   The input file is big endian
%
fid = fopen(fname,'rb','ieee-be');

if (fid < 0)
   error(me,'Cannot open file %s', fname);
end;

magic = mne_fread3(fid) ;
NEW_VERSION_MAGIC_NUMBER = 16777215;
if (magic == NEW_VERSION_MAGIC_NUMBER)
    nvert = fread(fid, 1, 'int32');
    nface = fread(fid, 1, 'int32');
    val_per_vertex = fread(fid, 1, 'int32');
    if val_per_vertex ~= 1
        fclose(fid);
        error(me,'Values per vertex not equal to one');
    end
    curv = fread(fid, nvert, 'float') ;
    fprintf(1,'\t%d values read from a new style curvature file.\n',nvert);
else
    nvert = magic;
    nface = mne_fread3(fid);
    curv = fread(fid, nvert, 'int16') ./ 100 ;
    fprintf(1,'\t%d values read from an old style curvature file.\n',nvert);
end

fclose(fid);

return;


text:
function [X, maxdot] = packing_on_the_sphere(d, n, epsilon, X0)
% Return a set of points spread out on the sphere.
%
% function [X, maxdot] = packing_on_the_sphere(d, n, epsilon, X0)
%
% Using optimization on the oblique manifold, that is, the product of
% spheres, this function returns a set of n points with unit norm in R^d in
% the form of a matrix X of size nxd, such that the points are spread out
% on the sphere. Ideally, we would minimize the maximum inner product
% between any two points X(i, :) and X(j, :), i~=j, but that is a nonsmooth
% cost function. Instead, we replace the max function by a classical
% log-sum-exp approximation and (attempt to) solve:
%
% min_{X in OB(d, n)} log( .5*sum_{i~=j} exp( xi'*xj/epsilon ) ),
%
% with xi = X(:, i) and epsilon is some "diffusion constant". As epsilon
% goes to zero, the cost function is a sharper approximation of the max
% function (under some assumptions), but the cost function becomes stiffer
% and hence harder to optimize.
%
% The second output, maxdot, is the maximum inner product between any two
% points in the returned X. This number is the one we truly are trying to
% minimize.
%
% Notice that this cost function is invariant under rotation of X:
% f(X) = f(XQ) for all orthogonal Q in O(d).
% This calls for optimization over the set of symmetric positive
% semidefinite matrices of size n and rank d with unit diagonal, which can
% be thought of as the quotient of the oblique manifold OB(d, n) by O(d):
% See elliptopefactory.
%
% This is known as the Thomson or, more specifically, the Tammes problem:
% http://en.wikipedia.org/wiki/Tammes_problem
% An interesting page by Neil Sloane collecting best known packings is
% available here http://neilsloane.com/packings/

% This file is part of Manopt and is copyrighted. See the license file.
%
% Main author: Nicolas Boumal, July 2, 2013
% Contributors:
%
% Change log:
%   Aug. 14, 2013 (NB) : Code now compatible to experiment with both the
%                        obliquefactory and the elliptopefactory.
%
%   Jan.  7, 2014 (NB) : Added reference to Neil Sloane's page and the
%                        maxdot output.
%
%   June 24, 2014 (NB) : Now shifting exponentials to alleviate numerical
%                        trouble when epsilon is too small.
%
%   Aug. 31, 2021 (XJ) : Added AD to compute the gradient

    if ~exist('d', 'var') || isempty(d)
        % Dimension of the embedding space: R^d
        d = 3;
    end
    if ~exist('n', 'var') || isempty(n)
        % Number n of points to place of the sphere in R^d.
        % For example, n=12 yields an icosahedron:
        % https://en.wikipedia.org/wiki/Icosahedron
        % Notice though that platonic solids are not always optimal.
        % Try for example n = 8: you don't get a cube.
        n = 24;
    end
    if ~exist('epsilon', 'var') || isempty(epsilon)
        % This value should be as close to 0 as affordable.
        % If it is too close to zero, optimization first becomes much
        % slower, than simply doesn't work anymore becomes of floating
        % point overflow errors (NaN's and Inf's start to appear).
        % If it is too large, then log-sum-exp is a poor approximation of
        % the max function, and the spread will be less uniform.
        % An okay value seems to be 0.01 or 0.001 for example. Note that a
        % better strategy than using a small epsilon straightaway is to
        % reduce epsilon bit by bit and to warm-start subsequent
        % optimization in that way. Trustregions will be more appropriate
        % for these fine tunings.
        epsilon = 0.0015;
    end

    % Pick your manifold (the elliptope factory quotients out the global
    % rotation invariance of the problem, which is more natural but
    % conceptually a bit more complicated --- for usage with the toolbox it
    % is the same though: just uncomment the appropriate line).
    manifold = obliquefactory(d, n, true);
    % manifold = elliptopefactory(n, d);

    % Generate a random initial guess if none was given.
    if ~exist('X0', 'var') || isempty(X0)
        X0 = manifold.rand();
    end

    % Define the cost function with caching system used: the store
    % structure we receive as input is tied to the input point X. Everytime
    % this cost function is called at this point X, we will receive the
    % same store structure back. We may modify the store structure inside
    % the function and return it: the changes will be remembered for next
    % time.
    function [f, store] = cost(X, store)
        if ~isfield(store, 'ready')
            XXt = X*X';
            % Shift the exponentials by the maximum value to reduce
            % numerical trouble due to possible overflows.
            s = max(max(triu(XXt, 1)));
            expXXt = exp((XXt-s)/epsilon);
            % Zero out the diagonal
            expXXt(1:(n+1):end) = 0;
            u = sum(sum(triu(expXXt, 1)));
            store.XXt = XXt;
            store.s = s;
            store.expXXt = expXXt;
            store.u = u;
            store.ready = true;
        end
        u = store.u;
        s = store.s;
        f = s + epsilon*log(u);
    end

    % Define the gradient of the cost. When the gradient is called at a
    % point X for which the cost was already called, the store structure we
    % receive remember everything that the cost function stored in it, so
    % we can reuse previously computed elements.
    function [g, store] = grad(X, store)
        if ~isfield(store, 'ready')
            [~, store] = cost(X, store);
        end
        % Compute the Euclidean gradient
        eg = store.expXXt*X / store.u;
        % Convert to the Riemannian gradient (by projection)
        g = manifold.egrad2rgrad(X, eg);
    end

    % Setup the problem structure with its manifold M and cost+grad
    % functions.
    problem.M = manifold;
    problem.cost = @cost;
    problem.grad = @grad;

    % An alternative way to compute the grad is to use automatic
    % differentiation provided in the deep learning toolbox (slower)
    % Notice that the function triu is not supported for AD so far.
    % Replace it with ctriu described in the file manoptADhelp.m
    % problem.cost = @cost_AD;
    %    function f = cost_AD(X)
    %        XXt = X*X';
    %        s = max(max(ctriu(XXt, 1)));
    %        expXXt = exp((XXt-s)/epsilon);
    %        expXXt(1:(n+1):end) = 0;
    %        u = sum(sum(ctriu(expXXt, 1)));
    %        f = s + epsilon*log(u);
    %    end
    % Call manoptAD to prepare AD for the problem structure
    % problem = manoptAD(problem,'egrad');

    % For debugging, it's always nice to check the gradient a few times.
    % checkgradient(problem);
    % pause;

    % Call a solver on our problem with a few options defined. We did not
    % specify the Hessian but it is still okay to call trustregion: Manopt
    % will approximate the Hessian with finite differences of the gradient.
    opts.tolgradnorm = 1e-8;
    opts.maxtime = 1200;
    opts.maxiter = 1e5;
    % X = trustregions(problem, X0, opts);
    X = conjugategradient(problem, X0, opts);

    % Evaluate the maximum inner product between any two points of X.
    XXt = X*X';
    dots = XXt(find(triu(ones(n), 1))); %#ok<FNDSB>
    maxdot = max(dots);

    % Similarly, even though we did not specify the Hessian, we may still
    % estimate its spectrum at the solution. It should reflect the
    % invariance of the cost function under a global rotatioon of the
    % sphere, which is an invariance under the group O(d) of dimension
    % d(d-1)/2 : this translates into d(d-1)/2 zero eigenvalues in the
    % spectrum of the Hessian.
    % The approximate Hessian is not a linear operator, and is it a
    % fortiori not symmetric. The result of this computation is thus not
    % reliable. It does display the zero eigenvalues as expected though.
    if manifold.dim() < 300
        evs = real(hessianspectrum(problem, X));
        figure;
        stem(1:length(evs), sort(evs), '.');
        title(['Eigenvalues of the approximate Hessian of the cost ' ...
               'function at the solution']);
    end

    % Show how the inner products X(:, i)'*X(:, j) are distributed.
    figure;
    histogram(real(acos(dots)), 20);
    title('Histogram of the geodesic distances');

    % This is the quantity we actually want to minimize.
    fprintf('Maximum inner product between two points: %g\n', maxdot);

    % Give some visualization if the dimension allows
    if d == 2
        % For the circle, the optimal solution consists in spreading the
        % points with angles uniformly sampled in (0, 2pi). This
        % corresponds to the following value for the max inner product:
        fprintf('Optimal value for the max inner product: %g\n', cos(2*pi/n));
        figure;
        t = linspace(-pi, pi, 201);
        plot(cos(t), sin(t), '-', 'LineWidth', 3, 'Color', [152,186,220]/255);
        daspect([1 1 1]);
        box off;
        axis off;
        hold on;
        plot(X(:, 1), X(:, 2), 'r.', 'MarkerSize', 25);
        hold off;
    end
    if d == 3
        figure;
        set(gcf, 'Color', 'w');
        % Plot the sphere
        [sphere_x, sphere_y, sphere_z] = sphere(50);
        handle = surf(sphere_x, sphere_y, sphere_z);
        set(handle, 'FaceColor', [152,186,220]/255);
        set(handle, 'FaceAlpha', .5);
        set(handle, 'EdgeColor', [152,186,220]/255);
        set(handle, 'EdgeAlpha', .5);
        daspect([1 1 1]);
        box off;
        axis off;
        hold on;
        % Add the chosen points
        Y = 1.02*X';
        plot3(Y(1, :), Y(2, :), Y(3, :), 'r.', 'MarkerSize', 25);
        % And connect the points which are at minimal distance,
        % within some tolerance.
        min_distance = real(acos(maxdot));
        connected = real(acos(XXt)) <= 1.20*min_distance;
        [Ic, Jc] = find(triu(connected, 1));
        for k = 1 : length(Ic)
            i = Ic(k); j = Jc(k);
            plot3(Y(1, [i j]), Y(2, [i j]), Y(3, [i j]), 'k-');
        end
        hold off;
    end

end

author: NicolasBoumal

repo: manopt


text:
function jac = p19_jac ( neqn, t, y )

%*****************************************************************************80
%
%% P19_JAC evaluates the jacobian for problem p19.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NEQN, the number of equations.
%
%    Input, real T, Y(NEQN), the arguments of the jacobian.
%
%    Output, real JAC(NEQN,NEQN), the jacobian matrix.
%
  jac = zeros ( neqn, neqn );

  d = ( sqrt ( ( y(1).^2 + y(2).^2 ) ) ).^5;

  jac(1,3) = 1.0;
  jac(2,4) = 1.0;
  jac(3,1) = ( 2.0 * y(1).^2 - y(2).^2 ) / d;
  jac(3,2) = 3.0 * y(1) * y(2) / d;
  jac(4,1) = 3.0 * y(1) * y(2) / d;
  jac(4,2) = ( - y(1).^2 + 2.0 * y(2).^2 ) / d;

  return
end