text:
% ------------------------------------------------------
% SwarmOps - Heuristic optimization for Matlab
% Copyright (C) 2003-2010 Magnus Erik Hvass Pedersen.
% Please see the file license.txt for license details.
% SwarmOps on the internet: http://www.Hvass-Labs.org/
% ------------------------------------------------------

% Behavioural parameters for Particle Swarm Optimization (PSO)
% tuned by Pedersen (1). The parameter-array consists of
% the following parameters:
% - Swarm-size (denoted s)
% - Inertia weight (denoted omega)
% - Particle's best weight (denoted phiP)
% - Swarm's best weight (denoted phiG)
%
% Select the parameters that most closely match the
% characteristics of your optimization problem.
% For example, if you want to optimize a problem where
% the search-space has 25 dimensions and you can perform
% 100000 evaluations, then you could first try using the
% parameters PSO_20DIM_40000EVALS. If that does not yield
% satisfactory results then you could try
% PSO_30DIM_60000EVALS or perhaps PSO_20DIM_400000EVALS_A.
% If that does not work then you will either need to tune
% the parameters for the problem at hand, or you should
% try using another optimizer.
%
% Literature references:
% (1) M.E.H. Pedersen.
%     Good parameters for Particle Swarm Optimization.
%     Technical Report HL1001, Hvass Laboratories, 2010.

% Parameters for non-parallel version:

PSO_2DIM_400EVALS_A       = [ 25,  0.3925,  2.5586, 1.3358];
PSO_2DIM_400EVALS_B       = [ 29, -0.4349, -0.6504, 2.2073];
PSO_2DIM_4000EVALS_A      = [156,  0.4091,  2.1304, 1.0575];
PSO_2DIM_4000EVALS_B      = [237, -0.2887,  0.4862, 2.5067];
PSO_5DIM_1000EVALS_A      = [ 63, -0.3593, -0.7238, 2.0289];
PSO_5DIM_1000EVALS_B      = [ 47, -0.1832,  0.5287, 3.1913];
PSO_5DIM_10000EVALS_A     = [223, -0.3699, -0.1207, 3.3657];
PSO_5DIM_10000EVALS_B     = [203,  0.5069,  2.5524, 1.0056];
PSO_10DIM_2000EVALS_A     = [ 63,  0.6571,  1.6319, 0.6239];
PSO_10DIM_2000EVALS_B     = [204, -0.2134, -0.3344, 2.3259];
PSO_10DIM_20000EVALS      = [ 53, -0.3488, -0.2746, 4.8976];
PSO_20DIM_40000EVALS      = [ 69, -0.4438, -0.2699, 3.3950];
PSO_20DIM_400000EVALS_A   = [149, -0.3236, -0.1136, 3.9789];
PSO_20DIM_400000EVALS_B   = [ 60, -0.4736, -0.9700, 3.7904];
PSO_20DIM_400000EVALS_C   = [256, -0.3499, -0.0513, 4.9087];
PSO_30DIM_60000EVALS      = [134, -0.1618,  1.8903, 2.1225];
PSO_30DIM_600000EVALS     = [ 95, -0.6031, -0.6485, 2.6475];
PSO_50DIM_100000EVALS     = [106, -0.2256, -0.1564, 3.8876];
PSO_100DIM_200000EVALS    = [161, -0.2089, -0.0787, 3.7637];
PSO_HANDTUNED             = [ 50,  0.7290,  1.4945, 1.4945];
PSO_DEFAULT               = PSO_20DIM_400000EVALS_A;

% Parameters for parallel version (above may also work):

PSO_PAR_5DIM_10000EVALS   = [ 72, -0.4031, -0.5631, 3.4277];
PSO_PAR_30DIM_60000EVALS  = [ 64, -0.2063, -2.7449, 2.3198];
PSO_PAR_DEFAULT           = PSO_PAR_5DIM_10000EVALS;

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


text:
function [err] = pco_uint32err(error_code)
 if(error_code<0)
  e=error_code*-1;
  err=uint32(e)-1;
  err=bitcmp(err,'uint32');
 else
  err=uint32(error_code);
 end
end

author: Shrediquette

repo: PIVlab

save_path: github-repos/MATLAB/Shrediquette-PIVlab


text:
function y = fwd(net,x)

% FWD
%
% Compute the output of a support vector classification network.
%
%    y = fwd(net, x);
%
% where x is a matrix of input patterns, where each column represents a
% variable and each row represents and observation.

%
% File        : @svc/fwd.m
%
% Date        : Tuesday 12th September 2000
%
% Author      : Dr Gavin C. Cawley
%
% Description : Part of an object-oriented implementation of Vapnik's Support
%               Vector Machine, as described in [1].
%
% References  : [1] V.N. Vapnik,
%                   "The Nature of Statistical Learning Theory",
%                   Springer-Verlag, New York, ISBN 0-387-94559-8,
%                   1995.
%
% History     : 17/08/1999 - v1.00
%               12/09/2000 - v1.01 minor improvments to comments and help
%                                  messages
%
% Copyright   : (c) Dr Gavin C. Cawley, September 2000.
%
%    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 2 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, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
%

w = repmat(net.w, size(x,1), 1);

y = sum((w.*evaluate(net.kernel,x,net.sv))') - net.bias;

y = y';

% bye bye...

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function value = ellipsoid_volume_3d ( r1, r2, r3 )

%*****************************************************************************80
%
%% ELLIPSOID_VOLUME_3D returns the volume of an ellipsoid in 3d.
%
%  Integration region:
%
%    Points (X,Y,Z) such that
%
%      ( ( X - XC ) / R1 )**2
%    + ( ( Y - YC ) / R2 )**2
%    + ( ( Z - ZC ) / R3 )**2 <= 1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 May 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real R1, R2, R3, the semi-axes of the ellipsoid.
%
%    Output, real VALUE, the volume of the ellipsoid.
%
  value = ( 4.0E+00 / 3.0E+00 ) * pi * r1 * r2 * r3;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%compute tail speed in the central-head direction
function [data,units]=compute_veltailch(trx,n)

larvae=trx.exp2flies{n};
numlarvae=numel(larvae);
veltailch=cell(1,numlarvae);
for i=1:numlarvae
    larva=larvae(i);
    veltailch{1,i}=trx(larva).velmagtail.*(cos(trx(larva).velangtail).*cos(trx(larva).centralheadang(1,1:end-1))+sin(trx(larva).velangtail).*sin(trx(larva).centralheadang(1,1:end-1)));
end

units=parseunits('mm/s');
data=veltailch;

author: kristinbranson

repo: JAABA


text:
% demo file for applying the NoRMCorre motion correction algorithm on
% 1-photon widefield imaging data
% Example file is provided from the miniscope project page
% www.miniscope.org

clear;
gcp;
%% download data and convert to single precision
name = 'msCam13.avi';
if ~exist(name,'file')  % download file if it doesn't exist in the directory
    url = 'https://caiman.flatironinstitute.org/~neuro/normcorre_datasets/msCam13.avi';
    fprintf('downloading the file...');
    outfilename = websave(name,url);
    fprintf('done.');
end

%addpath(genpath('../../NoRMCorre'));
Yf = read_file(name);
Yf = single(Yf);
[d1,d2,T] = size(Yf);

%% perform some sort of deblurring/high pass filtering

if (0)
    hLarge = fspecial('average', 40);
    hSmall = fspecial('average', 2);
    for t = 1:T
        Y(:,:,t) = filter2(hSmall,Yf(:,:,t)) - filter2(hLarge, Yf(:,:,t));
    end
    %Ypc = Yf - Y;
    bound = size(hLarge,1);
else
    gSig = 7;
    gSiz = 3*gSig;
    psf = fspecial('gaussian', round(2*gSiz), gSig);
    ind_nonzero = (psf(:)>=max(psf(:,1)));
    psf = psf-mean(psf(ind_nonzero));
    psf(~ind_nonzero) = 0;   % only use pixels within the center disk
    %Y = imfilter(Yf,psf,'same');
    %bound = 2*ceil(gSiz/2);
    Y = imfilter(Yf,psf,'symmetric');
    bound = 0;
end
%% first try out rigid motion correction
    % exclude boundaries due to high pass filtering effects
options_r = NoRMCorreSetParms('d1',d1-bound,'d2',d2-bound,'bin_width',200,'max_shift',20,'iter',1,'correct_bidir',false);

%% register using the high pass filtered data and apply shifts to original data
tic; [M1,shifts1,template1] = normcorre_batch(Y(bound/2+1:end-bound/2,bound/2+1:end-bound/2,:),options_r); toc % register filtered data
    % exclude boundaries due to high pass filtering effects
tic; Mr = apply_shifts(Yf,shifts1,options_r,bound/2,bound/2); toc % apply shifts to full dataset
    % apply shifts on the whole movie
%% compute metrics
[cY,mY,vY] = motion_metrics(Y(bound/2+1:end-bound/2,bound/2+1:end-bound/2,:),options_r.max_shift);
[cYf,mYf,vYf] = motion_metrics(Yf,options_r.max_shift);

[cM1,mM1,vM1] = motion_metrics(M1,options_r.max_shift);
[cM1f,mM1f,vM1f] = motion_metrics(Mr,options_r.max_shift);

%% plot rigid shifts and metrics
shifts_r = squeeze(cat(3,shifts1(:).shifts));
figure;
    subplot(311); plot(shifts_r);
        title('Rigid shifts','fontsize',14,'fontweight','bold');
        legend('y-shifts','x-shifts');
    subplot(312); plot(1:T,cY,1:T,cM1);
        title('Correlation coefficients on filtered movie','fontsize',14,'fontweight','bold');
        legend('raw','rigid');
    subplot(313); plot(1:T,cYf,1:T,cM1f);
        title('Correlation coefficients on full movie','fontsize',14,'fontweight','bold');
        legend('raw','rigid');

%% now apply non-rigid motion correction
% non-rigid motion correction is likely to produce very similar results
% since there is no raster scanning effect in wide field imaging

options_nr = NoRMCorreSetParms('d1',d1-bound,'d2',d2-bound,'bin_width',50, ...
    'grid_size',[128,128]*2,'mot_uf',4,'correct_bidir',false, ...
    'overlap_pre',32,'overlap_post',32,'max_shift',20);

tic; [M2,shifts2,template2] = normcorre_batch(Y(bound/2+1:end-bound/2,bound/2+1:end-bound/2,:),options_nr,template1); toc % register filtered data
tic; Mpr = apply_shifts(Yf,shifts2,options_nr,bound/2,bound/2); toc % apply the shifts to the removed percentile

%% compute metrics

[cM2,mM2,vM2] = motion_metrics(M2,options_nr.max_shift);
[cM2f,mM2f,vM2f] = motion_metrics(Mpr,options_nr.max_shift);

%% plot shifts

shifts_r = squeeze(cat(3,shifts1(:).shifts));
shifts_nr = cat(ndims(shifts2(1).shifts)+1,shifts2(:).shifts);
shifts_nr = reshape(shifts_nr,[],ndims(Y)-1,T);
shifts_x = squeeze(shifts_nr(:,2,:))';
shifts_y = squeeze(shifts_nr(:,1,:))';

patch_id = 1:size(shifts_x,2);
str = strtrim(cellstr(int2str(patch_id.')));
str = cellfun(@(x) ['patch # ',x],str,'un',0);

figure;
    ax1 = subplot(311); plot(1:T,cY,1:T,cM1,1:T,cM2); legend('raw data','rigid','non-rigid'); title('correlation coefficients for filtered data','fontsize',14,'fontweight','bold')
            set(gca,'Xtick',[],'XLim',[0,T-3])
    ax2 = subplot(312); plot(shifts_x); hold on; plot(shifts_r(:,2),'--k','linewidth',2); title('displacements along x','fontsize',14,'fontweight','bold')
            set(gca,'Xtick',[])
    ax3 = subplot(313); plot(shifts_y); hold on; plot(shifts_r(:,1),'--k','linewidth',2); title('displacements along y','fontsize',14,'fontweight','bold')
            xlabel('timestep','fontsize',14,'fontweight','bold')
    linkaxes([ax1,ax2,ax3],'x')

%% display downsampled data

tsub = 5;

Y_ds = downsample_data(Y,'time',tsub);
Yf_ds = downsample_data(Yf,'time',tsub);
M1_ds = downsample_data(M1,'time',tsub);
M1f_ds = downsample_data(Mr,'time',tsub);
M2_ds = downsample_data(M2,'time',tsub);
M2f_ds = downsample_data(Mpr,'time',tsub);
nnY_ds = quantile(Y_ds(:),0.0005);
mmY_ds = quantile(Y_ds(:),0.9995);
nnYf_ds = quantile(Yf_ds(:),0.0005);
mmYf_ds = quantile(Yf_ds(:),0.99995);
%%

make_avi = false; % save a movie
if make_avi
    vidObj = VideoWriter('filtered.avi');
    set(vidObj,'FrameRate',30);
    open(vidObj);
end
fig = figure;
    screensize = get(0,'Screensize' );
    fac = min(min((screensize(3:4)-100)./[3*d2,d1]),10);
    set(gcf, 'PaperUnits', 'points', 'Units', 'points');
    set(gcf, 'Position', round([100 100 fac*3*d2 fac*d1]));

for t = 1:1:size(Y_ds,3)
    if (0)
        % plot filtered data
        subplot(131);imagesc(Y_ds(:,:,t),[nnY_ds,mmY_ds]); xlabel('Raw data (downsampled)','fontsize',14,'fontweight','bold'); axis equal; axis tight;
        colormap('bone');
        set(gca,'XTick',[],'YTick',[]);
        subplot(132);imagesc(M1_ds(:,:,t),[nnY_ds,mmY_ds]); xlabel('rigid corrected','fontsize',14,'fontweight','bold'); axis equal; axis tight;
        title(sprintf('Frame %i out of %i',t,size(Y_ds,3)),'fontweight','bold','fontsize',14);
        colormap('bone')
        set(gca,'XTick',[],'YTick',[]);
        subplot(133);imagesc(M2_ds(:,:,t),[nnY_ds,mmY_ds]); xlabel('non-rigid corrected','fontsize',14,'fontweight','bold'); axis equal; axis tight;
        colormap('bone')
        set(gca,'XTick',[],'YTick',[]);
    else
        % plot full data
        subplot(131);imagesc(Yf_ds(:,:,t),[nnYf_ds,mmYf_ds]); xlabel('Raw data (downsampled)','fontsize',14,'fontweight','bold'); axis equal; axis tight;
        colormap('bone');
        set(gca,'XTick',[],'YTick',[]);
        subplot(132);imagesc(M1f_ds(:,:,t),[nnYf_ds,mmYf_ds]); xlabel('rigid corrected','fontsize',14,'fontweight','bold'); axis equal; axis tight;
        title(sprintf('Frame %i out of %i',t,size(Y_ds,3)),'fontweight','bold','fontsize',14);
        colormap('bone')
        set(gca,'XTick',[],'YTick',[]);
        subplot(133);imagesc(M2f_ds(:,:,t),[nnYf_ds,mmYf_ds]); xlabel('non-rigid corrected','fontsize',14,'fontweight','bold'); axis equal; axis tight;
        colormap('bone')
        set(gca,'XTick',[],'YTick',[]);
    end
    drawnow;
    if make_avi
        currFrame = getframe(fig);
        writeVideo(vidObj,currFrame);
    end
end
if make_avi
    close(vidObj);
end

author: flatironinstitute

repo: NoRMCorre


text:
function [affine_matrix, transformed_matrix] = PTKSolveForRigidTranslation(image_to_transform, reference_image, reporting)
    % PTKSolveForRigidTranslation. Computes the transformation matrix to register one
    %     image segmentation to another based on an solving for a rigid
    %     translation.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2012.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    % Initial centroid registration
    [com_affine_matrix, com_affine_vector] = PTKRegisterCentroid(image_to_transform, reference_image, reporting);

    % Before computing a distance transform, we must resample the images so they are approximately isotropic
    % Find a voxel size to use for the registration image. We divide up thick
    % slices to give an approximately isotropic voxel size
    register_voxel_size = reference_image.VoxelSize;
    register_voxel_size = register_voxel_size./round(register_voxel_size/min(register_voxel_size));
    reference_image2 = reference_image.Copy;
    reference_image2.ResampleBinary(register_voxel_size);
    image_to_transform2 = image_to_transform.Copy;
    image_to_transform2.ResampleBinary(register_voxel_size);

    reference_image2.AddBorder(20);
    image_to_transform2.AddBorder(20);

    dt_float = PTKRunForEachComponentAndCombine(@MimImageUtilities.GetNormalisedDT, image_to_transform2, image_to_transform2, reporting);
    dt_ref = PTKRunForEachComponentAndCombine(@MimImageUtilities.GetNormalisedDT, reference_image2, reference_image2, reporting);

    dt_float.RescaleToMaxSize(128);

    dt_ref.RescaleToMaxSize(128);

    [affine_matrix, transformed_matrix] = Solve(dt_float, dt_ref, com_affine_vector(4:6), reporting);
end

function [affine_matrix, transformed_matrix] = Solve(image_to_transform, reference_image, starting_affine_vector, reporting)
    [i_o, j_o, k_o] = image_to_transform.GetCentredGlobalCoordinatesMm;
    [i_o, j_o, k_o] = ndgrid(i_o, j_o, k_o);

    [i_r, j_r, k_r] = reference_image.GetCentredGlobalCoordinatesMm;
    [i_r, j_r, k_r] = ndgrid(i_r, j_r, k_r);

    x_0 = starting_affine_vector';
    AnonFn = @(x) FnToMinimiseRigid(x, image_to_transform, reference_image, i_o, j_o, k_o, i_r, j_r, k_r, reporting);

    [x_vector, fval, exitflag, output] = fminsearch(AnonFn, x_0, optimset('TolX',0.0001, 'TolFun', 0.01, 'PlotFcns', @optimplotfval));

    affine_matrix = MimImageCoordinateUtilities.CreateAffineTranslationMatrix(x_vector);

    transformed_matrix = PTKRegisterImageAffineUsingCoordinates(image_to_transform, reference_image, affine_matrix, i_o, j_o, k_o, i_r, j_r, k_r, '*linear', reporting);
end

function closeness = FnToMinimiseRigid(x_vector, image_to_transform, reference_image, i_o, j_o, k_o, i_r, j_r, k_r, reporting)
    affine_matrix = MimImageCoordinateUtilities.CreateAffineTranslationMatrix(x_vector);
    transformed_image = PTKRegisterImageAffineUsingCoordinates(image_to_transform, reference_image, affine_matrix, i_o, j_o, k_o, i_r, j_r, k_r, '*linear', reporting);
    closeness = ComputeCloseness(transformed_image, reference_image);
end

function closeness = ComputeCloseness(image1, image2)
    diff = (image1.RawImage - image2.RawImage).^2;
    closeness = sqrt(sum(diff(:))/numel(image1.RawImage));
end

author: tomdoel

repo: pulmonarytoolkit


text:
function f_x = ParFor4(in1)
%PARFOR4
%    F_X = PARFOR4(IN1)

%    This function was generated by the Symbolic Math Toolbox version 8.2.
%    20-Sep-2019 09:35:35

u = in1(:,1);
ux = in1(:,4);
uxx = in1(:,5);
uxxx = in1(:,6);
f_x = ((u.*1.807785219463578e16+ux.*2.87844038948516e19+uxx.*1.111806867891814e17-uxxx.*7.848202243764388e17-u.*ux.*2.277629423136014e19-u.*uxxx.*3.623747593487516e18+4.461317998952448e15).*-1.0)./(u.*6.942494202789888e15-2.602608832707297e17);

author: dynamicslab

repo: SINDy-PI


text:
function [out1,out2,out3] = ginput2(arg1)
%GINPUT Graphical input from mouse.
%   [X,Y] = GINPUT(N) gets N points from the current axes and returns
%   the X- and Y-coordinates in length N vectors X and Y.  The cursor
%   can be positioned using a mouse (or by using the Arrow Keys on some
%   systems).  Data points are entered by pressing a mouse button
%   or any key on the keyboard except carriage return, which terminates
%   the input before N points are entered.
%
%   [X,Y] = GINPUT gathers an unlimited number of points until the
%   return key is pressed.
%
%   [X,Y,BUTTON] = GINPUT(N) returns a third result, BUTTON, that
%   contains a vector of integers specifying which mouse button was
%   used (1,2,3 from left) or ASCII numbers if a key on the keyboard
%   was used.

%   Copyright (c) 1984-96 by The MathWorks, Inc.
%   $Revision: 5.18 $  $Date: 1996/11/10 17:48:08 $

% Fixed version by Jean-Yves Bouguet to have a cross instead of 2 lines
% More visible for images

out1 = []; out2 = []; out3 = []; y = [];
c = computer;
if ~strcmp(c(1:2),'PC') & ~strcmp(c(1:2),'MA')
    tp = get(0,'TerminalProtocol');
else
    tp = 'micro';
end

if ~strcmp(tp,'none') & ~strcmp(tp,'x') & ~strcmp(tp,'micro'),
    if nargout == 1,
        if nargin == 1,
            eval('out1 = trmginput(arg1);');
        else
            eval('out1 = trmginput;');
        end
    elseif nargout == 2 | nargout == 0,
        if nargin == 1,
            eval('[out1,out2] = trmginput(arg1);');
        else
            eval('[out1,out2] = trmginput;');
        end
        if  nargout == 0
            out1 = [ out1 out2 ];
        end
    elseif nargout == 3,
        if nargin == 1,
            eval('[out1,out2,out3] = trmginput(arg1);');
        else
            eval('[out1,out2,out3] = trmginput;');
        end
    end
else

    fig = gcf;
    figure(gcf);

    if nargin == 0
        how_many = -1;
        b = [];
    else
        how_many = arg1;
        b = [];
        if  isstr(how_many) ...
            | size(how_many,1) ~= 1 | size(how_many,2) ~= 1 ...
            | ~(fix(how_many) == how_many) ...
            | how_many < 0
            error('Requires a positive integer.')
        end
        if how_many == 0
            ptr_fig = 0;
            while(ptr_fig ~= fig)
                ptr_fig = get(0,'PointerWindow');
            end
            scrn_pt = get(0,'PointerLocation');
            loc = get(fig,'Position');
            pt = [scrn_pt(1) - loc(1), scrn_pt(2) - loc(2)];
            out1 = pt(1); y = pt(2);
        elseif how_many < 0
            error('Argument must be a positive integer.')
        end
    end

    pointer = get(gcf,'pointer');
    set(gcf,'pointer','crosshair');
    fig_units = get(fig,'units');
    char = 0;

    while how_many ~= 0
        % Use no-side effect WAITFORBUTTONPRESS
        waserr = 0;
        eval('keydown = wfbp;', 'waserr = 1;');
        if(waserr == 1)
      if(ishandle(fig))
        set(fig,'pointer',pointer,'units',fig_units);
        error('Interrupted');
      else
        error('Interrupted by figure deletion');
      end
        end

        ptr_fig = get(0,'CurrentFigure');
        if(ptr_fig == fig)
            if keydown
                char = get(fig, 'CurrentCharacter');
                button = abs(get(fig, 'CurrentCharacter'));
                scrn_pt = get(0, 'PointerLocation');
                set(fig,'units','pixels')
                loc = get(fig, 'Position');
                pt = [scrn_pt(1) - loc(1), scrn_pt(2) - loc(2)];
                set(fig,'CurrentPoint',pt);
            else
                button = get(fig, 'SelectionType');
                if strcmp(button,'open')
                    button = 1; %b(length(b));
                elseif strcmp(button,'normal')
                    button = 1;
                elseif strcmp(button,'extend')
                    button = 2;
                elseif strcmp(button,'alt')
                    button = 3;
                else
                    error('Invalid mouse selection.')
                end
            end
            pt = get(gca, 'CurrentPoint');

            how_many = how_many - 1;

            if(char == 13) % & how_many ~= 0)
                % if the return key was pressed, char will == 13,
                % and that's our signal to break out of here whether
                % or not we have collected all the requested data
                % points.
                % If this was an early breakout, don't include
                % the <Return> key info in the return arrays.
                % We will no longer count it if it's the last input.
                break;
            end

            out1 = [out1;pt(1,1)];
            y = [y;pt(1,2)];
            b = [b;button];
        end
    end

    set(fig,'pointer',pointer,'units',fig_units);

    if nargout > 1
        out2 = y;
        if nargout > 2
            out3 = b;
        end
    else
        out1 = [out1 y];
    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function key = wfbp
%WFBP   Replacement for WAITFORBUTTONPRESS that has no side effects.

% Remove figure button functions
fprops = {'windowbuttonupfcn','buttondownfcn', ...
          'windowbuttondownfcn','windowbuttonmotionfcn'};
fig = gcf;
fvals = get(fig,fprops);
set(fig,fprops,{'','','',''})

% Remove all other buttondown functions
ax = findobj(fig,'type','axes');
if isempty(ax)
    ch = {};
else
    ch = get(ax,{'Children'});
end
for i=1:length(ch),
    ch{i} = ch{i}(:)';
end
h = [ax(:)',ch{:}];
vals = get(h,{'buttondownfcn'});
mt = repmat({''},size(vals));
set(h,{'buttondownfcn'},mt);

% Now wait for that buttonpress, and check for error conditions
waserr = 0;
eval(['if nargout==0,', ...
      '   waitforbuttonpress,', ...
      'else,', ...
      '   keydown = waitforbuttonpress;',...
      'end' ], 'waserr = 1;');

% Put everything back
if(ishandle(fig))
  set(fig,fprops,fvals)
  set(h,{'buttondownfcn'},vals)
end

if(waserr == 1)
  error('Interrupted');
end

if nargout>0, key = keydown; end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

repo: ekfmonoslam


text:
function [ind,indf] = get_index(d),
  ind= d.index;
  indf = d.findex;

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [ new_input_args ] = reorient_chipper_args( symmetry, datasize, varargin)
%REORIENT_CHIPPER_ARGS Applies a transform to the chipper indexing
%
% This function does nothing but transform the arguments into a chipper
% function based on the symmetry specified.
%
% A chipper function is a simplified function of a file reader.  It takes
% only three arguments: [minIndexInDimension1 maxIndexInDimension1],
% [minIndexInDimension2 maxIndexInDimension2], [subsampleInDimension1,
% subsampleInDimension2].  All other information such as the file handle
% and other auxilliary information must be contained within the chipper
% function.
%
% The SYMMETRY parameter is used to perform symmetry operations on the
% imported image data.  A combination of horizontal mirroring (hm), vertical
% mirroring (vm) and transpose(tr) is used to accomplish the symmetry
% operations.  The 'symmetry' input parameter is defined as [ hm vm tr ], where
% each of the elements can be either 0 or 1 to indicate whether to apply a
% transformation or not.
%
%   [ 0 0 0 ] - default setting; successive pixels on disk are interpreted
%               to fill the first dimension in a MATLAB array (vertical, if
%               viewed in imagesc/imshow without any reorienting)
%   [ 1 0 1 ] -  90 degrees CCW from [ 0 0 0 ]
%   [ 0 1 1 ] -  90 degrees  CW from [ 0 0 0 ]
%   [ 1 1 0 ] - 180 degrees     from [ 0 0 0 ]
%
%   [ 0 0 1 ] - transpose of [ 0 0 0 ];
%   [ 0 1 0 ] -  90 degrees CCW from [ 0 0 1 ]
%   [ 1 0 0 ] -  90 degrees  CW from [ 0 0 1 ]
%   [ 1 1 1 ] - 180 degrees     from [ 0 0 1 ]
%
% Written by: Wade Schwartzkopf, NGA/IDT
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

if symmetry(3)
    new_input_args={};
    if(nargin>=3), new_input_args{2}=varargin{1}; end;
    if(nargin>=4), new_input_args{1}=varargin{2}; end
    if(nargin>=5), new_input_args{3}=varargin{3}(end:-1:1); end;
else
    new_input_args=varargin;
end
if symmetry(1)
    if(nargin>=3)
        new_input_args{1}=double(datasize(1))-new_input_args{1}(end:-1:1)+1;
        % The double around datasize allows this to work on empty arguments.
    end
end
if symmetry(2)
    if(nargin>=4)
        new_input_args{2}=double(datasize(2))-new_input_args{2}(end:-1:1)+1;
    end
end

end

% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
function data = loadjson(fname,varargin)
%
% data=loadjson(fname,opt)
%    or
% data=loadjson(fname,'param1',value1,'param2',value2,...)
%
% parse a JSON (JavaScript Object Notation) file or string
%
% authors:Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
%            date: 2011/09/09
%         Nedialko Krouchev: http://www.mathworks.com/matlabcentral/fileexchange/25713
%            date: 2009/11/02
%         François Glineur: http://www.mathworks.com/matlabcentral/fileexchange/23393
%            date: 2009/03/22
%         Joel Feenstra:
%         http://www.mathworks.com/matlabcentral/fileexchange/20565
%            date: 2008/07/03
%
% $Id: loadjson.m 415 2013-10-07 16:38:31Z fangq $
%
% input:
%      fname: input file name, if fname contains "{}" or "[]", fname
%             will be interpreted as a JSON string
%      opt: a struct to store parsing options, opt can be replaced by
%           a list of ('param',value) pairs. The param string is equivallent
%           to a field in opt.
%
% output:
%      dat: a cell array, where {...} blocks are converted into cell arrays,
%           and [...] are converted to arrays
%
% license:
%     BSD license, see LICENSE_BSD.txt files for details
%
% -- this function is part of jsonlab toolbox (http://iso2mesh.sf.net/cgi-bin/index.cgi?jsonlab)
%

global pos inStr len  esc index_esc len_esc isoct arraytoken

if(regexp(fname,'[\{\}\]\[]','once'))
   string=fname;
elseif(exist(fname,'file'))
   fid = fopen(fname,'rt');
   string = fscanf(fid,'%c');
   fclose(fid);
else
   error('input file does not exist');
end

pos = 1; len = length(string); inStr = string;
isoct=exist('OCTAVE_VERSION');
arraytoken=find(inStr=='[' | inStr==']' | inStr=='"');
jstr=regexprep(inStr,'\\\\','  ');
escquote=regexp(jstr,'\\"');
arraytoken=sort([arraytoken escquote]);

% String delimiters and escape chars identified to improve speed:
esc = find(inStr=='"' | inStr=='\' ); % comparable to: regexp(inStr, '["\\]');
index_esc = 1; len_esc = length(esc);

opt=varargin2struct(varargin{:});
jsoncount=1;
while pos <= len
    switch(next_char)
        case '{'
            data{jsoncount} = parse_object(opt);
        case '['
            data{jsoncount} = parse_array(opt);
        otherwise
            error_pos('Outer level structure must be an object or an array');
    end
    jsoncount=jsoncount+1;
end % while

jsoncount=length(data);
if(jsoncount==1 && iscell(data))
    data=data{1};
end

if(~isempty(data))
      if(isstruct(data)) % data can be a struct array
          data=jstruct2array(data);
      elseif(iscell(data))
          data=jcell2array(data);
      end
end

%%
function newdata=parse_collection(id,data,obj)

if(jsoncount>0 && exist('data','var'))
    if(~iscell(data))
       newdata=cell(1);
       newdata{1}=data;
       data=newdata;
    end
end

%%
function newdata=jcell2array(data)
len=length(data);
newdata=data;
for i=1:len
      if(isstruct(data{i}))
          newdata{i}=jstruct2array(data{i});
      elseif(iscell(data{i}))
          newdata{i}=jcell2array(data{i});
      end
end

%%-------------------------------------------------------------------------
function newdata=jstruct2array(data)
fn=fieldnames(data);
newdata=data;
len=length(data);
for i=1:length(fn) % depth-first
    for j=1:len
        if(isstruct(getfield(data(j),fn{i})))
            newdata(j)=setfield(newdata(j),fn{i},jstruct2array(getfield(data(j),fn{i})));
        end
    end
end
if(~isempty(strmatch('x0x5F_ArrayType_',fn)) && ~isempty(strmatch('x0x5F_ArrayData_',fn)))
  newdata=cell(len,1);
  for j=1:len
    ndata=cast(data(j).x0x5F_ArrayData_,data(j).x0x5F_ArrayType_);
    iscpx=0;
    if(~isempty(strmatch('x0x5F_ArrayIsComplex_',fn)))
        if(data(j).x0x5F_ArrayIsComplex_)
           iscpx=1;
        end
    end
    if(~isempty(strmatch('x0x5F_ArrayIsSparse_',fn)))
        if(data(j).x0x5F_ArrayIsSparse_)
            if(~isempty(strmatch('x0x5F_ArraySize_',fn)))
                dim=data(j).x0x5F_ArraySize_;
                if(iscpx && size(ndata,2)==4-any(dim==1))
                    ndata(:,end-1)=complex(ndata(:,end-1),ndata(:,end));
                end
                if isempty(ndata)
                    % All-zeros sparse
                    ndata=sparse(dim(1),prod(dim(2:end)));
                elseif dim(1)==1
                    % Sparse row vector
                    ndata=sparse(1,ndata(:,1),ndata(:,2),dim(1),prod(dim(2:end)));
                elseif dim(2)==1
                    % Sparse column vector
                    ndata=sparse(ndata(:,1),1,ndata(:,2),dim(1),prod(dim(2:end)));
                else
                    % Generic sparse array.
                    ndata=sparse(ndata(:,1),ndata(:,2),ndata(:,3),dim(1),prod(dim(2:end)));
                end
            else
                if(iscpx && size(ndata,2)==4)
                    ndata(:,3)=complex(ndata(:,3),ndata(:,4));
                end
                ndata=sparse(ndata(:,1),ndata(:,2),ndata(:,3));
            end
        end
    elseif(~isempty(strmatch('x0x5F_ArraySize_',fn)))
        if(iscpx && size(ndata,2)==2)
             ndata=complex(ndata(:,1),ndata(:,2));
        end
        ndata=reshape(ndata(:),data(j).x0x5F_ArraySize_);
    end
    newdata{j}=ndata;
  end
  if(len==1)
      newdata=newdata{1};
  end
end

%%-------------------------------------------------------------------------
function object = parse_object(varargin)
    parse_char('{');
    object = [];
    if next_char ~= '}'
        while 1
            str = parseStr(varargin{:});
            if isempty(str)
                error_pos('Name of value at position %d cannot be empty');
            end
            parse_char(':');
            val = parse_value(varargin{:});
            eval( sprintf( 'object.%s  = val;', valid_field(str) ) );
            if next_char == '}'
                break;
            end
            parse_char(',');
        end
    end
    parse_char('}');

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

function object = parse_array(varargin) % JSON array is written in row-major order
global pos inStr isoct
    parse_char('[');
    object = cell(0, 1);
    dim2=[];
    if next_char ~= ']'
        [endpos e1l e1r maxlevel]=matching_bracket(inStr,pos);
        arraystr=['[' inStr(pos:endpos)];
        arraystr=regexprep(arraystr,'"_NaN_"','NaN');
        arraystr=regexprep(arraystr,'"([-+]*)_Inf_"','$1Inf');
        arraystr(find(arraystr==sprintf('\n')))=[];
        arraystr(find(arraystr==sprintf('\r')))=[];
        %arraystr=regexprep(arraystr,'\s*,',','); % this is slow,sometimes needed
        if(~isempty(e1l) && ~isempty(e1r)) % the array is in 2D or higher D
            astr=inStr((e1l+1):(e1r-1));
            astr=regexprep(astr,'"_NaN_"','NaN');
            astr=regexprep(astr,'"([-+]*)_Inf_"','$1Inf');
            astr(find(astr==sprintf('\n')))=[];
            astr(find(astr==sprintf('\r')))=[];
            astr(find(astr==' '))='';
            if(isempty(find(astr=='[', 1))) % array is 2D
                dim2=length(sscanf(astr,'%f,',[1 inf]));
            end
        else % array is 1D
            astr=arraystr(2:end-1);
            astr(find(astr==' '))='';
            [obj count errmsg nextidx]=sscanf(astr,'%f,',[1,inf]);
            if(nextidx>=length(astr)-1)
                object=obj;
                pos=endpos;
                parse_char(']');
                return;
            end
        end
        if(~isempty(dim2))
            astr=arraystr;
            astr(find(astr=='['))='';
            astr(find(astr==']'))='';
            astr(find(astr==' '))='';
            [obj count errmsg nextidx]=sscanf(astr,'%f,',inf);
            if(nextidx>=length(astr)-1)
                object=reshape(obj,dim2,numel(obj)/dim2)';
                pos=endpos;
                parse_char(']');
                return;
            end
        end
        arraystr=regexprep(arraystr,'\]\s*,','];');
        try
           if(isoct && regexp(arraystr,'"','once'))
                error('Octave eval can produce empty cells for JSON-like input');
           end
           object=eval(arraystr);
           pos=endpos;
        catch
         while 1
            val = parse_value(varargin{:});
            object{end+1} = val;
            if next_char == ']'
                break;
            end
            parse_char(',');
         end
        end
    end
    if(jsonopt('SimplifyCell',0,varargin{:})==1)
      try
        oldobj=object;
        object=cell2mat(object')';
        if(iscell(oldobj) && isstruct(object) && numel(object)>1 && jsonopt('SimplifyCellArray',1,varargin{:})==0)
            object=oldobj;
        elseif(size(object,1)>1 && ndims(object)==2)
            object=object';
        end
      catch
      end
    end
    parse_char(']');

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

function parse_char(c)
    global pos inStr len
    skip_whitespace;
    if pos > len || inStr(pos) ~= c
        error_pos(sprintf('Expected %c at position %%d', c));
    else
        pos = pos + 1;
        skip_whitespace;
    end

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

function c = next_char
    global pos inStr len
    skip_whitespace;
    if pos > len
        c = [];
    else
        c = inStr(pos);
    end

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

function skip_whitespace
    global pos inStr len
    while pos <= len && isspace(inStr(pos))
        pos = pos + 1;
    end

%%-------------------------------------------------------------------------
function str = parseStr(varargin)
    global pos inStr len  esc index_esc len_esc
 % len, ns = length(inStr), keyboard
    if inStr(pos) ~= '"'
        error_pos('String starting with " expected at position %d');
    else
        pos = pos + 1;
    end
    str = '';
    while pos <= len
        while index_esc <= len_esc && esc(index_esc) < pos
            index_esc = index_esc + 1;
        end
        if index_esc > len_esc
            str = [str inStr(pos:len)];
            pos = len + 1;
            break;
        else
            str = [str inStr(pos:esc(index_esc)-1)];
            pos = esc(index_esc);
        end
        nstr = length(str); switch inStr(pos)
            case '"'
                pos = pos + 1;
                if(~isempty(str))
                    if(strcmp(str,'_Inf_'))
                        str=Inf;
                    elseif(strcmp(str,'-_Inf_'))
                        str=-Inf;
                    elseif(strcmp(str,'_NaN_'))
                        str=NaN;
                    end
                end
                return;
            case '\'
                if pos+1 > len
                    error_pos('End of file reached right after escape character');
                end
                pos = pos + 1;
                switch inStr(pos)
                    case {'"' '\' '/'}
                        str(nstr+1) = inStr(pos);
                        pos = pos + 1;
                    case {'b' 'f' 'n' 'r' 't'}
                        str(nstr+1) = sprintf(['\' inStr(pos)]);
                        pos = pos + 1;
                    case 'u'
                        if pos+4 > len
                            error_pos('End of file reached in escaped unicode character');
                        end
                        str(nstr+(1:6)) = inStr(pos-1:pos+4);
                        pos = pos + 5;
                end
            otherwise % should never happen
                str(nstr+1) = inStr(pos), keyboard
                pos = pos + 1;
        end
    end
    error_pos('End of file while expecting end of inStr');

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

function num = parse_number(varargin)
    global pos inStr len isoct
    currstr=inStr(pos:end);
    numstr=0;
    if(isoct~=0)
        numstr=regexp(currstr,'^\s*-?(?:0|[1-9]\d*)(?:\.\d+)?(?:[eE][+\-]?\d+)?','end');
        [num, one] = sscanf(currstr, '%f', 1);
        delta=numstr+1;
    else
        [num, one, err, delta] = sscanf(currstr, '%f', 1);
        if ~isempty(err)
            error_pos('Error reading number at position %d');
        end
    end
    pos = pos + delta-1;

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

function val = parse_value(varargin)
    global pos inStr len
    % true = 1; false = 0;

    switch(inStr(pos))
        case '"'
            val = parseStr(varargin{:});
            return;
        case '['
            val = parse_array(varargin{:});
            return;
        case '{'
            val = parse_object(varargin{:});
            if isstruct(val)
                if(~isempty(strmatch('x0x5F_ArrayType_',fieldnames(val), 'exact')))
                    val=jstruct2array(val);
                end
            elseif isempty(val)
                val = struct;
            end
            return;
        case {'-','0','1','2','3','4','5','6','7','8','9'}
            val = parse_number(varargin{:});
            return;
        case 't'
            if pos+3 <= len && strcmpi(inStr(pos:pos+3), 'true')
                val = true;
                pos = pos + 4;
                return;
            end
        case 'f'
            if pos+4 <= len && strcmpi(inStr(pos:pos+4), 'false')
                val = false;
                pos = pos + 5;
                return;
            end
        case 'n'
            if pos+3 <= len && strcmpi(inStr(pos:pos+3), 'null')
                val = NaN;
                pos = pos + 4;
                return;
            end
    end
    error_pos('Value expected at position %d');
%%-------------------------------------------------------------------------

function error_pos(msg)
    global pos inStr len
    poShow = max(min([pos-15 pos-1 pos pos+20],len),1);
    if poShow(3) == poShow(2)
        poShow(3:4) = poShow(2)+[0 -1];  % display nothing after
    end
    msg = [sprintf(msg, pos) ': ' ...
    inStr(poShow(1):poShow(2)) '<error>' inStr(poShow(3):poShow(4)) ];
    error( ['JSONparser:invalidFormat: ' msg] );

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

function str = valid_field(str)
global isoct
% From MATLAB doc: field names must begin with a letter, which may be
% followed by any combination of letters, digits, and underscores.
% Invalid characters will be converted to underscores, and the prefix
% "x0x[Hex code]_" will be added if the first character is not a letter.
    pos=regexp(str,'^[^A-Za-z]','once');
    if(~isempty(pos))
        if(~isoct)
            str=regexprep(str,'^([^A-Za-z])','x0x${sprintf(''%X'',unicode2native($1))}_','once');
        else
            str=sprintf('x0x%X_%s',char(str(1)),str(2:end));
        end
    end
    if(isempty(regexp(str,'[^0-9A-Za-z_]', 'once' ))) return;  end
    if(~isoct)
        str=regexprep(str,'([^0-9A-Za-z_])','_0x${sprintf(''%X'',unicode2native($1))}_');
    else
        pos=regexp(str,'[^0-9A-Za-z_]');
        if(isempty(pos)) return; end
        str0=str;
        pos0=[0 pos(:)' length(str)];
        str='';
        for i=1:length(pos)
            str=[str str0(pos0(i)+1:pos(i)-1) sprintf('_0x%X_',str0(pos(i)))];
        end
        if(pos(end)~=length(str))
            str=[str str0(pos0(end-1)+1:pos0(end))];
        end
    end
    %str(~isletter(str) & ~('0' <= str & str <= '9')) = '_';

%%-------------------------------------------------------------------------
function endpos = matching_quote(str,pos)
len=length(str);
while(pos<len)
    if(str(pos)=='"')
        if(~(pos>1 && str(pos-1)=='\'))
            endpos=pos;
            return;
        end
    end
    pos=pos+1;
end
error('unmatched quotation mark');
%%-------------------------------------------------------------------------
function [endpos e1l e1r maxlevel] = matching_bracket(str,pos)
global arraytoken
level=1;
maxlevel=level;
endpos=0;
bpos=arraytoken(arraytoken>=pos);
tokens=str(bpos);
len=length(tokens);
pos=1;
e1l=[];
e1r=[];
while(pos<=len)
    c=tokens(pos);
    if(c==']')
        level=level-1;
        if(isempty(e1r)) e1r=bpos(pos); end
        if(level==0)
            endpos=bpos(pos);
            return
        end
    end
    if(c=='[')
        if(isempty(e1l)) e1l=bpos(pos); end
        level=level+1;
        maxlevel=max(maxlevel,level);
    end
    if(c=='"')
        pos=matching_quote(tokens,pos+1);
    end
    pos=pos+1;
end
if(endpos==0)
    error('unmatched "]"');
end

author: plotly


text:
function k = pascal_to_i4 ( i, j )

%*****************************************************************************80
%
%% PASCAL_TO_I4 converts Pacal triangle coordinates to a linear index.
%
%  Discussion:
%
%    We describe the grid points in a Pascal triangle in two ways:
%
%    As a linear index K:
%
%                     1
%                   2   3
%                 4   5   6
%               7   8   9   10
%
%    As elements (I,J) of Pascal's triangle:
%
%                     0,0
%                  1,0   0,1
%               2,0   1,1    0,2
%            3,0   2,1   1,2    0,3
%
%  Example:
%
%     K  I  J
%
%     1  0  0
%     2  1  0
%     3  0  1
%     4  2  0
%     5  1  1
%     6  0  2
%     7  3  0
%     8  2  1
%     9  1  2
%    10  0  3
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, J, the row and column indices.  I and J must
%    be nonnegative.
%
%    Output, integer K, the linear index of the (I,J) element.
%
  if ( i < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PASCAL_TO_I4 - Fatal error!\n' );
    fprintf ( 1, '  I < 0.\n' );
    fprintf ( 1, '  I = %d\n', i );
    error ( 'PASCAL_TO_I4 - Fatal error!' );
  elseif ( j < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PASCAL_TO_I4 - Fatal error!\n' );
    fprintf ( 1, '  J < 0.\n' );
    fprintf ( 1, '  J = %d\n', j );
    error ( 'PASCAL_TO_I4 - Fatal error!' );
  end

  d = i + j;

  k = ( d * ( d + 1 ) ) / 2 + j + 1;

  return
end


text:
function [matrix]=MakeCircle(Diameter,Size,edgewidth);
% MAKECIRCLE - make a circle with hard or soft edgesin a matrix of minimum size
% [matrix]=MakeCircle(Diameter,Size,edgewidth);
% Sizes are in pixels, edgewith > 1 = with of raised cosine function

if nargin < 1,
	help(mfilename);
	return;
end;
if nargin < 2,
	Size=Diameter;
end;
if nargin < 3,
	edgewidth = 0;
end;

% check
if Diameter+edgewidth>Size,
	disp(sprintf('Diameter + edgewidth are larger than Size'));
end;

% Make matrix coordinate system with 0 in center:
coords=-Size/2+.5:Size/2-.5;
[xx,yy]=meshgrid(coords);
dist_from_center = sqrt((xx.*xx)+(yy.*yy));clear xx yy;
fix_circ_calc = find(dist_from_center <= Diameter/2);
if edgewidth,
   edge_calc     = find(dist_from_center >  Diameter/2 & ...
                        dist_from_center <= Diameter/2 + edgewidth);
end;

%Make actual matrix
image_matrix=zeros(Size,Size);
image_matrix(fix_circ_calc)=1;
if edgewidth,
	dist_from_center = cos((dist_from_center-Diameter/2)./(edgewidth).*pi)./2+.5;
	image_matrix(edge_calc)=dist_from_center(edge_calc);
end;

matrix=image_matrix;
%ImageShow(matrix);


text:
classdef prtUiRocSelector < prtUiManagerPanel

    properties

        prtDs

        pf = [];
        pd = [];
        thresholds = [];
        selectedIndex = [];

        handleStruct
    end

    properties (Hidden, SetAccess='protected', GetAccess='protected')
        retainObsUpdateCallbackDepHelper = [];
        retainObsDepHelper = [];
    end
    properties (Dependent)
        retainObs
        retainObsUpdateCallback
    end

    methods
        function self = prtUiRocSelector(varargin)
            if nargin == 1
                self.prtDs = varargin{1};
            else
                self = prtUtilAssignStringValuePairs(self,varargin{:});
            end

            if nargin~=0 && ~self.hgIsValid
               self.create()
            end

            init(self);
        end

        function init(self)
            self.handleStruct.axes = axes('parent',self.managedHandle,'units','normalized','position',[0.1 0.1 0.85 0.8]);
            hold on
            grid on
            self.handleStruct.rocLine = plot(self.handleStruct.axes,nan,nan,'k');
            self.handleStruct.selectX = plot(self.handleStruct.axes,nan,nan,'kx','MarkerSize',12,'HitTest','off');
            hold off
            self.handleStruct.title = title(self.handleStruct.axes,'');
            set(self.handleStruct.rocLine,'HitTest','on','ButtonDownFcn',@(h,e)self.infoUpdate());

            self.updateRoc();
        end
        function updateRoc(self)

            self.selectedIndex = [];
            set(self.handleStruct.selectX,'XData',nan','YData',nan);
            set(self.handleStruct.title,'String','')

            if isempty(self.prtDs)
                set(self.handleStruct.rocLine,'XData',nan','YData',nan);
                self.selectedIndex = [];
                return
            end

            try
                if ~isempty(self.retainObs)
                    [self.pf ,self.pd, self.thresholds] = prtScoreRoc(self.prtDs.retainObservations(self.retainObs));
                else
                    [self.pf ,self.pd, self.thresholds] = prtScoreRoc(self.prtDs);
                end
            catch ME
                msgbox(ME.message,ME.identifier,'Error','Modal')
                set(self.handleStruct.rocLine,'XData',nan','YData',nan);
                return
            end
            set(self.handleStruct.rocLine,'XData',self.pf,'YData',self.pd);
            axis(self.handleStruct.axes,[0 1 0 1]);
        end
        function infoUpdate(self)
            cp = get(self.handleStruct.axes,'CurrentPoint');
            self.selectedIndex = mean(find(cp(1,1) > self.pf,1,'last'),find(cp(1,2) < self.pd,1,'first'));

            set(self.handleStruct.selectX,'XData',self.pf(self.selectedIndex),'YData',self.pd(self.selectedIndex));
            set(self.handleStruct.title,'String',sprintf('Pd = %0.3f, Pf = %0.3f, Theshold = %g',self.pd(self.selectedIndex), self.pf(self.selectedIndex),self.thresholds(self.selectedIndex)));

        end

        function val = get.retainObsUpdateCallback(self)
            val = self.retainObsUpdateCallbackDepHelper;
        end
        function set.retainObsUpdateCallback(self,val)
            assert(isempty(val) || (isa(val, 'function_handle') && nargin(val)==1),'retainObsUpdateCallback must be a function handle that accepts one input')

            self.retainObsUpdateCallbackDepHelper = val;
        end
        function val = get.retainObs(self)
            val = self.retainObsDepHelper;
        end
        function set.retainObs(self, val)
            if ~isempty(self.retainObsUpdateCallback)
                self.retainObsUpdateCallback(val)
            end
            self.retainObsDepHelper = val;
            self.updateRoc();
        end
        function updateRetainObs(self, val)
            self.retainObs = val;
        end

    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function varargout = flatmap(varargin)
% VL_FLATMAP Flatten a tree, assigning the label of the root to each node
%   [LABELS CLUSTERS] = VL_FLATMAP(MAP) labels each tree of the forest contained
%   in MAP. LABELS contains the linear index of the root node in MAP, CLUSTERS
%   instead contains a label between 1 and the number of clusters.
%
%   See also: VL_HELP().
[varargout{1:nargout}] = vl_flatmap(varargin{:});

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function [rev,file] = headRev(folders,strs,rec)
%HEADREV Retrieve the head revision for specified files
%
%   This function finds the most recently committed file in the specified
%   folders(s) and returns its revision and filename. The files must
%   contain the 'Revision' and 'Date' keywords.
%
%   REV = IOSR.SVN.HEADREV(FOLDERS) returns the head revision (i.e. the
%   highest revision) for all files in FOLDERS.
%
%   REV = IOSR.SVN.HEADREV(FOLDERS,STRS) allows additional filter strings
%   STRS to be specified. See GETCONTENTS for details of permitted filter
%   strings.
%
%   REV = IOSR.SVN.HEADREV(FOLDERS,STRS,REC), with REC = true, allows
%   folders to be searched recursively (default is false).
%
%   [REV,FILE] = IOSR.SVN.HEADREV(...) returns the filename to FILE of the
%   file with the highest revision.
%
%   FOLDERS and STRS may be a string specifying a single occurrence, or a
%   cell array of strings specifying multiple occurrences. This function
%   requires that the 'Revision' keyword is used in the searched functions.
%
%   Note that if the folders include files from an externally-defined
%   repository (which has been updated more recently than the native
%   repository), a misleading revision number may be presented.
%
%   See also IOSR.GENERAL.GETCONTENTS, IOSR.SVN.BUILDSVNPROFILE,
%            IOSR.SVN.READSVNKEYWORD.

%   Copyright 2016 University of Surrey.

    keyword = 'Date'; % use this data to sort files (doubtful that any others would work, without considerable parsing)

    if ischar(folders)
        folders = cellstr(folders);
    end

    if nargin > 1
        if ischar(strs)
            strs = cellstr(strs);
        end
        if isempty(strs)
            strs = {'files'};
        end
    else
        strs = {'files'};
    end

    % add subfolders if recursion is requested
    if nargin > 2
        if rec
            for f = folders
                folders = [folders; iosr.general.getContents(char(f),'filter','folders','rec',true,'path','full')]; %#ok<AGROW>
            end
        end
    end

    % get revision info for files
    svn_profile = iosr.svn.buildSvnProfile(folders,keyword,strs);

    % remove the keyword and convert revision strings to numbers
    svn_profile(:,2) = cellfun(@(x) (x(length(keyword)+3:end)),svn_profile(:,2),'uni',false);

    % remove NaN (files that don't have the revision keyword)
    IX1 = cellfun(@isempty,svn_profile(:,2));
    svn_profile = svn_profile(~IX1,:);

    % sort by revision number
    [~,IX2] = sort(svn_profile(:,2));
    svn_profile = svn_profile(IX2,:);

    % get highest revision number
    rev = iosr.svn.readSvnKeyword(svn_profile{end,1},'Revision');

    % return corresponding filename, if requested
    if nargout>1
        file = svn_profile{end,1};
    end

end

author: IoSR-Surrey

repo: MatlabToolbox


text:
function [mu, sigma, map] = correction_step(mu, sigma, z, map);
% Updates the belief, i.e., mu and sigma after observing landmarks,
% and augments the map with newly observed landmarks.
% The employed sensor model measures the range and bearing of a landmark
% mu: state vector containing robot pose and poses of landmarks obeserved so far.
% Current robot pose = mu(1:3)
% Note that the landmark poses in mu are stacked in the order by which they were observed
% sigma: the covariance matrix of the system.
% z: struct array containing the landmark observations.
% Each observation z(i) has an id z(i).id, a range z(i).range, and a bearing z(i).bearing
% The vector 'map' contains the ids of all landmarks observed so far by the robot in the order
% by which they were observed, NOT in ascending id order.

% For computing sigma
global scale;

% Number of measurements in this time step
m = size(z, 2);

% Measurement noise
Q = 0.01*eye(2);
n = length(mu);
zm = zeros(2,1);
S = zeros(2);
xbar = 0;
ybar = 0;

for i = 1:m

	% If the landmark is observed for the first time:
	if (isempty(find(map == z(i).id)))
	    % Add new landmark to the map
            [mu, sigma, map] = add_landmark_to_map(mu, sigma, z(i), map, Q);
	    % The measurement has been incorporated so we quit the correction step
	    continue;
	endif
	n=length(mu);
	sigma_x_z = zeros(n,2);
	% Compute sigma points from the predicted mean and covariance
        % This corresponds to line 6 on slide 32
	sigma_points = compute_sigma_points(mu, sigma);
        % Normalize!
	sigma_points(3,:) = normalize_angle(sigma_points(3,:));

	% Compute lambda
	n = length(mu);
	num_sig = size(sigma_points,2);
	lambda = scale - n;

        % extract the current location of the landmark for each sigma point
        % Use this for computing an expected measurement, i.e., applying the h function
	landmarkIndex = find(map==(z(i).id));
	landmarkXs = sigma_points(2*landmarkIndex + 2, :);
	landmarkYs = sigma_points(2*landmarkIndex + 3, :);

	% TODO: Compute z_points (2x2n+1), which consists of predicted measurements from all sigma points
        % This corresponds to line 7 on slide 32
	z_points(:,:) = [landmarkXs;landmarkYs] + repmat([z(i).range*cos(normalize_angle(z(i).bearing + mu(3)));z(i).range*sin(normalize_angle(z(i).bearing + mu(3)))],1,2*n+1);

        % setup the weight vector for mean and covariance
        wm = [lambda/scale, repmat(1/(2*scale), 1, 2*n)];
        wc = wm;

	% TODO: Compute zm, line 8 on slide 32
	% zm is the recovered expected measurement mean from z_points.
	% It will be a 2x1 vector [expected_range; expected_bearing].
        % For computing the expected_bearing compute a weighted average by
        % summing the sines/cosines of the angle
	for j=1:2*n+1
	    zm(1,1) = zm(1,1) + wm(j)*z_points(1,j);
	    xbar = xbar + wm(j)*cos(z_points(2,j));
	    ybar = ybar + wm(j)*sin(z_points(2,j));
	endfor
	    zm(2,1) = atan2(ybar,xbar);

	% TODO: Compute the innovation covariance matrix S (2x2), line 9 on slide 32
        % Remember to normalize the bearing after computing the difference
	for j=1:2*n+1
	    S = S + wc(j)*(z_points(:,j)-zm)*(z_points(:,j)-zm)' + Q;
	endfor

	% TODO: Compute Sigma_x_z, line 10 on slide 32
        % (which is equivalent to sigma times the Jacobian H transposed in EKF).
	% sigma_x_z is an nx2 matrix, where n is the current dimensionality of mu
        % Remember to normalize the bearing after computing the difference
	n = length(mu);
	for j=1:2*n+1
	    sigma_x_z = sigma_x_z + wc(j)*(sigma_points(:,j) - mu)*(z_points(:,j)-zm)';
	endfor

	% TODO: Compute the Kalman gain, line 11 on slide 32
	K = sigma_x_z/S;

	% Get the actual measurement as a vector (for computing the difference to the observation)
	z_actual = [z(i).range; z(i).bearing];

	% TODO: Update mu and sigma, line 12 + 13 on slide 32
        % normalize the relative bearing
	zdiff = z_actual - zm;
	zdiff(2) = normalize_angle(zdiff(2));
	mu = mu + K*zdiff;
	mu(3) = normalize_angle(mu(3));

	% TODO: Normalize the robot heading mu(3)

endfor

end

author: kiran-mohan

repo: SLAM-Algorithms-Octave


text:
function D = driving_function_mono_localwfs_sbl_pw(x0,nk,f,conf)
%DRIVING_FUNCTION_MONO_LOCALWFS_SBL_PW driving signal for a plane wave using
%local WFS with spatial bandwidth limitation
%
%   Usage: D = driving_function_mono_localwfs_sbl_pw(x0,nk,f,conf)
%
%   Input parameters:
%       x0          - position and direction of the secondary source / m [nx7]
%       nk          - propagation direction of plane wave / m [1x3]
%       f           - frequency of the monochromatic source / Hz
%       conf        - configuration struct (see SFS_config)
%
%   Output parameters:
%       D           - driving function [nx1]
%
%   See also: sound_field_mono_localwfs_sbl, driving_function_mono_localwfs_sbl
%
%   References:
%       Winter, Hahn, Spors (2017), "Time-Domain Realisation of Model-Based
%       Rendering for 2.5D Local Wave Field Synthesis Using Spatial
%       Bandwidth-Limitation", 25th European Signal Processing Conference
%       (EUSIPCO), pp. 718-722, https://bit.ly/2yMjlOw

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 4;
nargmax = 4;
narginchk(nargmin,nargmax);

%% ===== Configuration ========================================================
N0 = size(x0,1);
xref = conf.xref;
% Maximum order of circular basis expansion of sound field
if isempty(conf.localwfs_sbl.order)
    Nce = nfchoa_order(N0,conf);
else
    Nce = conf.localwfs_sbl.order;
end
% Resolution of plane wave decomposition
if isempty(conf.localwfs_sbl.Npw)
    Npw = 2*ceil(2*pi*0.9*f/conf.c*conf.secondary_sources.size/2);
else
    Npw = conf.localwfs_sbl.Npw;
end

%% ===== Computation ==========================================================
% Circular expansion coefficients, Winter et al. (2017), eq. (12)
Pm = circexp_mono_pw(nk,Nce,f,xref,conf);
% Modal window
wm = modal_weighting(Nce,conf);
Pm = bsxfun(@times,[wm(end:-1:2),wm],Pm);
% Plane wave decomposition, inverse FT of Winter et al. (2017), eq. (10)
Ppwd = pwd_mono_circexp(Pm,Npw);
% Driving signal, Winter et al. (2017), eq. (8)
D = driving_function_mono_wfs_pwd(x0,Ppwd,f,xref,conf);

author: sfstoolbox

repo: sfs-matlab


text:
function z = plus( x, y )
    %PLUS Addition of two TT/MPS tensors.
    %   Z = PLUS(X,Y) adds two TT/MPS tensors. The rank of the resulting
    %   tensor is 2*R.
    %
    %   See also MINUS, UMINUS.

    %   TTeMPS Toolbox.
    %   Michael Steinlechner, 2013-2016
    %   Questions and contact: michael.steinlechner@epfl.ch
    %   BSD 2-clause license, see LICENSE.txt

    % add sanity check...
    rx = x.rank;
    ry = y.rank;
    nx = x.size;

    z = TTeMPS( cell(1, x.order) );

    % first core:
    p = size(x.U{1},4);
    tmp = zeros( 1, nx(1), rx(2)+ry(2), p );
    tmp( 1, :, 1:rx(2), : ) = x.U{1};
    tmp( 1, :, rx(2)+1:end, : ) = y.U{1};
    z.U{1} = tmp;

    % central cores:
    for i = 2:x.order-1
        % possibility of block format:
        p = size(x.U{i},4);
        tmp = zeros( rx(i)+ry(i), nx(i), rx(i+1)+ry(i+1), p);
        tmp( 1:rx(i), :, 1:rx(i+1), :) = x.U{i};
        tmp( rx(i)+1:end, :, rx(i+1)+1:end, :) = y.U{i};
        z.U{i} = tmp;
    end

    % last core:
    p = size(x.U{end},4);
    tmp = zeros( rx(end-1)+ry(end-1), nx(end), 1, p );
    tmp( 1:rx(end-1), :, 1, : ) = x.U{end};
    tmp( rx(end-1)+1:end, :, 1, : ) = y.U{end};
    z.U{end} = tmp;
end

author: NicolasBoumal

repo: manopt


text:
function hermite_cubic_test12 ( )

%*****************************************************************************80
%
%% HERMITE_CUBIC_TEST12 tests HERMITE_CUBIC_LAGRANGE_INTEGRAL.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2011
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'HERMITE_CUBIC_TEST12:\n' );
  fprintf ( 1, '  HERMITE_CUBIC_LAGRANGE_INTEGRAL returns the integrals\n' );
  fprintf ( 1, '  of the four Lagrange basis functions associated \n' );
  fprintf ( 1, '  with F1, D1, F2 and D2 such that\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  P(X) = F1 * LF1(X) + D1 * LD1(X)\n' );
  fprintf ( 1, '       + F2 * LF2(X) + D2 * LD2(X).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The Lagrange basis function integrals:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '        X1          X2          LF1         LD1         LF2         LD2\n' );
  fprintf ( 1, '\n' );

  x2 = 1.0;
  for x1 = -6 : +2
    q = hermite_cubic_lagrange_integral ( x1, x2 );
    fprintf ( 1, '  %10.4f  %10.4f  %10.4f  %10.4f  %10.4f  %10.4f\n', ...
      x1, x2, q(1:4) )
  end

  return
end


text:
function [NewTessFile, iSurface] = tess_concatenate( TessFiles, NewComment, fileType, Labels)
% TESS_CONCATENATE: Concatenate various surface files into one new file.
%
% USAGE:  [NewTessFile, iSurface] = tess_concatenate(TessFiles, NewComment='New surface', fileType='Other', Labels)
%          [NewTessMat, iSurface] = tess_concatenate(TessMats,  NewComment='New surface', fileType='Other', Labels)
%
% INPUT:
%    - TessFiles   : Cell-array of paths to surfaces files to concatenate
%    - TessMats    : Array of loaded surface structures
%    - NewComment  : Name of the output surface
%    - fileType    : File type for the new file {'Cortex', 'InnerSkull', 'OuterSkull', 'Scalp', 'Other'}
% OUTPUT:
%    - NewTessFile : Filename of the newly created file
%    - iSurface    : Index of the new surface file

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

% Parse inputs
if (nargin < 4) || isempty(Labels)
    Labels = [];
end
if (nargin < 3) || isempty(fileType)
    fileType = [];
end
if (nargin < 2) || isempty(NewComment)
    NewComment = [];
end
% Save current scouts modifications
panel_scout('SaveModifications');
% Initialize new structure
NewTess = db_template('surfacemat');
isLeft = 0;
isRight = 0;
isWhite = 0;
isCortex = 0;
isVolAtlas = 0;
isSave = 1;

% Progress bar
bst_progress('start', 'Merge surfaces', 'Processing... ');
% Process all the files to concatenate
for iFile = 1:length(TessFiles)
    % Load tesselation
    if iscell(TessFiles)
        [oldTess, TessFiles{iFile}] = in_tess_bst(TessFiles{iFile});
        if isempty(oldTess)
            continue
        end
    % Files already loaded in calling function
    else
        oldTess = TessFiles(iFile);
        isSave = 0;
    end
    % Detect if right/left hemisphere
    if ~isempty(strfind(oldTess.Comment, 'lh.')) || ~isempty(strfind(oldTess.Comment, 'Lhemi')) || ~isempty(strfind(oldTess.Comment, 'Lwhite')) || ~isempty(strfind(oldTess.Comment, 'left'))
        isLeft = 1;
        scoutTag = ' L';
        scoutHemi = 'L';
        scoutComment = 'Cortex';
    elseif ~isempty(strfind(oldTess.Comment, 'rh.')) || ~isempty(strfind(oldTess.Comment, 'Rhemi')) || ~isempty(strfind(oldTess.Comment, 'Rwhite')) || ~isempty(strfind(oldTess.Comment, 'right'))
        isRight = 1;
        scoutTag = ' R';
        scoutHemi = 'R';
        scoutComment = 'Cortex';
    % Detect based on comment (tag ' L' or ' R' already present)
    elseif (length(oldTess.Comment) > 2) && strcmpi(oldTess.Comment(end-1:end), ' L')
        scoutTag = ' L';
        scoutHemi = 'L';
        scoutComment = oldTess.Comment;
    elseif (length(oldTess.Comment) > 2) && strcmpi(oldTess.Comment(end-1:end), ' R')
        scoutTag = ' R';
        scoutHemi = 'R';
        scoutComment = oldTess.Comment;
    % Guess based on the coordinates
    else
        if (nnz(oldTess.Vertices(:,2) > 0) > 5 * nnz(oldTess.Vertices(:,2) < 0))
            scoutTag = ' L';
            scoutHemi = 'L';
        elseif (nnz(oldTess.Vertices(:,2) < 0) > 5 * nnz(oldTess.Vertices(:,2) > 0))
            scoutTag = ' R';
            scoutHemi = 'R';
        else
            scoutTag = '';
            scoutHemi = 'U';
        end
        scoutComment = oldTess.Comment;
    end
    % Detect some specific types of surfaces
    if ~isempty(strfind(oldTess.Comment, 'white'))
        isWhite = 1;
    end
    if ~isempty(strfind(oldTess.Comment, 'cortex_'))
        isCortex = 1;
    end
    if ~isempty(strfind(oldTess.Comment, 'aseg')) || ~isempty(strfind(oldTess.Comment, 'subcortical'))
        isVolAtlas = 1;
    end
    % Concatenate current sub-tess to final tesselation structure
    offsetVertices   = size(NewTess.Vertices,1);
    NewTess.Faces    = [NewTess.Faces; oldTess.Faces + offsetVertices];
    NewTess.Vertices = [NewTess.Vertices; oldTess.Vertices];

    % History: Merged surface #i
    NewTess = bst_history('add', NewTess, 'merge', sprintf('Merge surface #%d: %s', iFile, oldTess.Comment));
    % History: Copy history of surface #i
    if isfield(oldTess, 'History') && ~isempty(oldTess.History)
        NewTess = bst_history('add', NewTess, oldTess.History, '  => ');
    end

    % Add an atlas "Structures" to reference the origins of each structure (if it does not exist)
    iAtlasStruct = find(strcmpi({oldTess.Atlas.Name}, 'Structures'));
    if isempty(iAtlasStruct)
        % Create new atlas
        iAtlasStruct = length(oldTess.Atlas) + 1;
        oldTess.Atlas(iAtlasStruct) = db_template('atlas');
        oldTess.Atlas(iAtlasStruct).Name = 'Structures';
        % Create one scout that describes all the structure
        oldTess.Atlas(iAtlasStruct).Scouts(1).Vertices = 1:length(oldTess.Vertices);
        oldTess.Atlas(iAtlasStruct).Scouts(1).Seed     = 1;
        oldTess.Atlas(iAtlasStruct).Scouts(1).Label    = scoutComment;
        oldTess.Atlas(iAtlasStruct).Scouts(1).Function = 'Mean';
        oldTess.Atlas(iAtlasStruct).Scouts(1).Region   = [scoutHemi 'U'];
        % If there are labels in inputs, try to get a color for the scout
        if ~isempty(Labels) && (size(Labels,2) == 3)
            iLabel = find(strcmpi(Labels(:,2), scoutComment));
            if ~isempty(iLabel)
                oldTess.Atlas(iAtlasStruct).Scouts(1).Color = Labels{iLabel,3} ./ 255;
            end
        end
        % Set scout color
        if isempty(oldTess.Atlas(iAtlasStruct).Scouts(1).Color)
            oldTess.Atlas(iAtlasStruct).Scouts(1) = panel_scout('SetColorAuto', oldTess.Atlas(iAtlasStruct).Scouts(1));
        end
    end

    % Concatenate atlases/scouts
    for iAtlasOld = 1:length(oldTess.Atlas)
        % Look for the same altas in the new surface
        iAtlasNew = find(strcmpi({NewTess.Atlas.Name}, oldTess.Atlas(iAtlasOld).Name));
        if isempty(iAtlasNew)
            iAtlasNew = length(NewTess.Atlas) + 1;
            NewTess.Atlas(iAtlasNew) = db_template('atlas');
            NewTess.Atlas(iAtlasNew).Name = oldTess.Atlas(iAtlasOld).Name;
        end
        % Loop over all scouts in the old surface to fix the indices of the vertices
        for iScout = 1:length(oldTess.Atlas(iAtlasOld).Scouts)
            % Adjust vertex indices
            oldTess.Atlas(iAtlasOld).Scouts(iScout).Vertices = oldTess.Atlas(iAtlasOld).Scouts(iScout).Vertices + offsetVertices;
            oldTess.Atlas(iAtlasOld).Scouts(iScout).Seed     = oldTess.Atlas(iAtlasOld).Scouts(iScout).Seed     + offsetVertices;
            % Add the first letter of the surface comment to the scout name
            if ~ismember(oldTess.Atlas(iAtlasOld).Scouts(iScout).Label(end), {'L', 'R'})
                oldTess.Atlas(iAtlasOld).Scouts(iScout).Label = [oldTess.Atlas(iAtlasOld).Scouts(iScout).Label, scoutTag];
            end
            % Set region name
            oldRegion = oldTess.Atlas(iAtlasOld).Scouts(iScout).Region;
            if isempty(oldRegion) || (length(oldRegion) < 2)
                oldTess.Atlas(iAtlasOld).Scouts(iScout).Region = [scoutHemi, 'U'];
            elseif (oldTess.Atlas(iAtlasOld).Scouts(iScout).Region(1) == 'U')
                oldTess.Atlas(iAtlasOld).Scouts(iScout).Region(1) = scoutHemi;
            end
        end
        % Add to surface
        if isempty(NewTess.Atlas(iAtlasNew).Scouts) && isempty(oldTess.Atlas(iAtlasOld).Scouts)
            NewTess.Atlas(iAtlasNew).Scouts = repmat(db_template('scout'), 0);
        elseif isempty(NewTess.Atlas(iAtlasNew).Scouts)
            NewTess.Atlas(iAtlasNew).Scouts = oldTess.Atlas(iAtlasOld).Scouts;
        else
            NewTess.Atlas(iAtlasNew).Scouts = [...
                struct_fix(db_template('scout'), NewTess.Atlas(iAtlasNew).Scouts), ...
                struct_fix(db_template('scout'), oldTess.Atlas(iAtlasOld).Scouts)];
        end
    end

    % Concatenate FreeSurfer registration spheres
    if isfield(oldTess, 'Reg') && isfield(oldTess.Reg, 'Sphere') && isfield(oldTess.Reg.Sphere, 'Vertices') && ~isempty(oldTess.Reg.Sphere.Vertices)
        if ~isfield(NewTess, 'Reg') || ~isfield(NewTess.Reg, 'Sphere') || ~isfield(NewTess.Reg.Sphere, 'Vertices')
            NewTess.Reg.Sphere.Vertices = oldTess.Reg.Sphere.Vertices;
        else
            NewTess.Reg.Sphere.Vertices = [NewTess.Reg.Sphere.Vertices; oldTess.Reg.Sphere.Vertices];
        end
    % Remove all spheres if they are not available in all the input files (disabled for now because cortex_cereb works well for projections)
    % elseif isfield(NewTess, 'Reg') &&  isfield(NewTess.Reg, 'Sphere') && isfield(NewTess.Reg.Sphere, 'Vertices') && ~isempty(NewTess.Reg.Sphere.Vertices)
    %     NewTess.Reg.Sphere.Vertices = [];
    end

    if isfield(oldTess, 'Reg') && isfield(oldTess.Reg, 'SphereLR') && isfield(oldTess.Reg.SphereLR, 'Vertices') && ~isempty(oldTess.Reg.SphereLR.Vertices)
        if ~isfield(NewTess, 'Reg') || ~isfield(NewTess.Reg, 'SphereLR') || ~isfield(NewTess.Reg.SphereLR, 'Vertices')
            NewTess.Reg.SphereLR.Vertices = oldTess.Reg.SphereLR.Vertices;
        else
            NewTess.Reg.SphereLR.Vertices = [NewTess.Reg.SphereLR.Vertices; oldTess.Reg.SphereLR.Vertices];
        end
    % Remove all spheres if they are not available in all the input files (disabled for now because cortex_cereb works well for projections)
    % elseif isfield(NewTess, 'Reg') &&  isfield(NewTess.Reg, 'SphereLR') && isfield(NewTess.Reg.SphereLR, 'Vertices') && ~isempty(NewTess.Reg.SphereLR.Vertices)
    %         NewTess.Reg.SphereLR.Vertices = [];
    end

    % Concatenate BrainSuite registration squares
    if isfield(oldTess, 'Reg') && isfield(oldTess.Reg, 'Square') && isfield(oldTess.Reg.Square, 'Vertices') && ~isempty(oldTess.Reg.Square.Vertices)
        if ~isfield(NewTess, 'Reg') || ~isfield(NewTess.Reg, 'Square') || ~isfield(NewTess.Reg.Square, 'Vertices')
            NewTess.Reg.Square.Vertices      = oldTess.Reg.Square.Vertices;
            NewTess.Reg.AtlasSquare.Vertices = oldTess.Reg.AtlasSquare.Vertices;
        else
            NewTess.Reg.Square.Vertices      = [NewTess.Reg.Square.Vertices;      oldTess.Reg.Square.Vertices];
            NewTess.Reg.AtlasSquare.Vertices = [NewTess.Reg.AtlasSquare.Vertices; oldTess.Reg.AtlasSquare.Vertices];
        end
    % Remove all spheres if they are not available in all the input files (disabled for now because cortex_cereb works well for projections)
    % elseif isfield(NewTess, 'Reg') && isfield(NewTess.Reg, 'Square') && isfield(NewTess.Reg.Square, 'Vertices') && ~isempty(NewTess.Reg.Square.Vertices)
    %     NewTess.Reg.Square.Vertices = [];
    %     NewTess.Reg.AtlasSquare.Vertices = [];
    end
end

% Sort scouts by name
for iAtlas = 1:length(NewTess.Atlas)
    [tmp, iSort] = sort({NewTess.Atlas(iAtlas).Scouts.Label});
    NewTess.Atlas(iAtlas).Scouts = NewTess.Atlas(iAtlas).Scouts(iSort);
end

% Detect surfaces types
if isLeft && isRight
    % File type: Cortex
    if isempty(fileType)
        fileType = 'Cortex';
    end
    % White matter
    if isWhite
        fileTag = 'cortex_white';
        if isempty(NewComment)
            NewComment = sprintf('white_%dV', length(NewTess.Vertices));
        end
    % Pial/cortex external envelope
    else
        fileTag = 'cortex_pial';
        if isempty(NewComment)
            NewComment = sprintf('cortex_%dV', length(NewTess.Vertices));
        end
    end
elseif isCortex && isVolAtlas
    fileTag = 'cortex_mixed';
    if isempty(fileType)
        fileType = 'Cortex';
    end
    if isempty(NewComment)
        NewComment = sprintf('cortex_mixed_%dV', length(NewTess.Vertices));
    end
else
    fileTag = 'concat';
    if isempty(fileType)
        fileType = 'Other';
    end
    if isempty(NewComment)
        NewComment = 'New surface';
    end
end
% Surface comments
NewTess.Comment = NewComment;
% History: Merge completed
NewTess = bst_history('add', NewTess, 'merge', 'Merge completed');

% ===== SAVE IN DATABASE =====
if isSave
    % Create new filename
    NewTessFile = bst_fullfile(bst_fileparts(TessFiles{1}), ['tess_' fileTag '.mat']);
    NewTessFile = file_unique(NewTessFile);
    % Save file
    bst_save(NewTessFile, NewTess, 'v7');
    % Make output filename relative
    NewTessFile = file_short(NewTessFile);
    % Get subject
    [sSubject, iSubject] = bst_get('SurfaceFile', TessFiles{1});
    % Register this file in Brainstorm database
    iSurface = db_add_surface(iSubject, NewTessFile, NewComment, fileType);
else
    NewTessFile = NewTess;
    iSurface = [];
end
% Close progress bar
bst_progress('stop');

author: brainstorm-tools


text:
function disp(L)
%DISP   Pretty-print a linop summary.
%
% See also DISPLAY.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Obtain size of the linop
[m, n] = size(L);

% Determine whether format is compact or loose:
loose = ~isequal(get(0, 'FormatSpacing'), 'compact');

if ( loose )
    fprintf('\n');
end

fprintf('   %ix%i linear operator', m,n)

if ( loose )
    fprintf('\n');
end

nc = size(L.constraint.functional, 1);
if ( nc == 1 )
    fprintf('\n   with 1 constraint/boundary condition.')
elseif ( nc > 0 )
    fprintf('\n   with %i constraints/boundary conditions.',nc)
end

if ( loose )
    fprintf('\n');
end

fprintf('\n');

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%CALIBRATECAMERACHARUCO  Calibrate a camera using ChArUco corners
%
%     [cameraMatrix, distCoeffs, reprojErr] = cv.calibrateCameraCharuco(charucoCorners, charucoIds, board, imageSize)
%     [cameraMatrix, distCoeffs, reprojErr, rvecs, tvecs, stdDevsIntrinsics, stdDevsExtrinsics, perViewErrors] = cv.calibrateCameraCharuco(...)
%     [...] = cv.calibrateCameraCharuco(..., 'OptionName',optionValue, ...)
%
% ## Input
% * __charucoCorners__ cell-array of detected charuco corners per frame
%   `{{[x,y],..}, ..}`.
% * __charucoIds__ cell-array of identifiers for each corner in
%   `charucoCorners` per frame (0-based) `{[id,..], ..}`
% * __board__ Marker board layout.
% * __imageSize__ input image size `[w,h]`.
%
% ## Output
% * __cameraMatrix__ Output 3x3 floating-point camera matrix
%   `A = [fx 0 cx; 0 fy cy; 0 0 1]`. If any of `UseIntrinsicGuess`,
%   `FixAspectRatio`, or `FixFocalLength` are specified, some or all of `fx`,
%   `fy`, `cx`, `cy` must be initialized before calling the function.
% * __distCoeffs__ Output vector of distortion coefficients
%   `[k1,k2,p1,p2,k3,k4,k5,k6,s1,s2,s3,s4]` of 4, 5, 8 or 12 elements.
% * __reprojErr__ the overall RMS re-projection error.
% * __rvecs__ Output cell-array of rotation vectors estimated for each board
%   view. That is, each k-th rotation vector together with the corresponding
%   k-th translation vector (see the next output parameter description) brings
%   the board pattern from the model coordinate space (in which object points
%   are specified) to the world coordinate space, that is, a real position of
%   the board pattern in the k-th pattern view (`k=0..M-1`).
% * __tvecs__ Output cell-array of translation vectors estimated for each
%   pattern view.
% * __stdDevsIntrinsics__ Output vector of standard deviations estimated for
%   intrinsic parameters. Order of deviations values:
%   `(fx,fy,cx,cy,k1,k2,p1,p2,k3,k4,k5,k6,s1,s2,s3,s4,taux,tauy)`. If one of
%   parameters is not estimated, its deviation is equals to zero.
% * __stdDevsExtrinsics__ Output vector of standard deviations estimated for
%   extrinsic parameters. Order of deviations values: `(R1, T1, ..., RM, TM)`
%   where `M` is number of pattern views, `Ri, Ti` are concatenated 1x3
%   vectors.
% * __perViewErrors__ Output vector of average re-projection errors estimated
%   for each pattern view.
%
% ## Options
% * __CameraMatrix__, __DistCoeffs__, __UseIntrinsicGuess__,
%   __FixPrincipalPoint__, __FixFocalLength__, __FixAspectRatio__,
%   __ZeroTangentDist__, __FixTangentDist__, __FixK1__, ..., __FixK6__,
%   __RationalModel__, __ThinPrismModel__, __FixS1S2S3S4__, __TiltedModel__,
%   __FixTauXTauY__, __UseLU__, __UseQR__
%   Different flags for the calibration process. See cv.calibrateCamera for
%   details.
% * __Criteria__ Termination criteria for the iterative optimization algorithm.
%   default `struct('type','Count+EPS', 'maxCount',30, 'epsilon',eps)`
%
% This function calibrates a camera using a set of corners of a ChArUco Board.
% The function receives a list of detected corners and its identifiers from
% several views of the Board.
%
% See also: cv.interpolateCornersCharuco, cv.calibrateCameraAruco,
%  cv.calibrateCamera, cv.Rodrigues
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function clusterQuantity = clusterAverage(clu, spikeQuantity)
% function clusterQuantity = clusterAverage(clu, spikeQuantity)
%
% get the average of some quantity across spikes in each cluster, given the
% quantity for each spike
%
% e.g.
% > clusterDepths = clusterAverage(clu, spikeDepths);
%
% clu and spikeQuantity must be vector, same size

% using a super-tricky algorithm for this - when you make a sparse
% array, the values of any duplicate indices are added. So this is the
% fastest way I know to make the sum of the entries of spikeQuantity for each of
% the unique entries of clu
[~, spikeCounts] = countUnique(clu);

% convert clu to indices, i.e. just values between 1 and nClusters.
[~,~,cluInds] = unique(clu);

% summation
q = full(sparse(cluInds, ones(size(clu)), double(spikeQuantity)));

% had sums, so dividing by spike counts gives the mean depth of each cluster
clusterQuantity = q./spikeCounts;

author: MouseLand

repo: Kilosort

save_path: github-repos/MATLAB/MouseLand-Kilosort


text:
function G = gradient( f )
%GRADIENT   Numerical gradient of a DISKFUN.
%   G = GRADIENT(F) returns the numerical gradient of the
%   DISKFUN F as a DISKFUNV G.
%
% See also DISKFUNV/DIV, DISKFUNV/CURL, DISKFUN/CURL, DISKFUNV/VORTICITY

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Empty check.
if isempty( f )
    G = diskfunv;
    return
end

fx = diff(f, 1);   % diff in x-variable
fy = diff(f, 2);   % diff in y-variable

G = diskfunv(fx,fy);

end


text:
%
% Make a string by combining a sequence of
% both strings and numbers (rounded), where each
% number is converted to a string
%
function str = makestr(varargin)

    n = numel(varargin);
    str = [];
    for i=1:n
        A = varargin{i};
        if ischar(A)
            str = [str A];
        elseif isnumeric(A)
            str = [str num2str(A,3)];
        else
            error('invalid input type');
        end
    end

end

author: brendenlake

repo: BPL

save_path: github-repos/MATLAB/brendenlake-BPL


text:
classdef IN1KMOP8 < EvoXBenchProblem
% <multi> <real> <large/none> <expensive/none>

%------------------------------- Reference --------------------------------
% Z. Lu, R. Cheng, Y. Jin, K. C. Tan, and K. Deb, Neural architecture
% search as multiobjective optimization benchmarks: Problem formulation and
% performance assessment, IEEE Transactions on Evolutionary Computation,
% 2023.
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            config.name = 'mnv3';
            config.args.objs = 'err&params&flops';
            config.args.normalized_objectives = false;
            obj.Setting@EvoXBenchProblem(config);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [seq, results] = get_sequence_results(seq)

if strcmpi(seq.format, 'otb')
    results.type = 'rect';
    results.res = seq.rect_position;
elseif strcmpi(seq.format, 'vot')
    seq.handle.quit(seq.handle);
else
    error('Uknown sequence format');
end

if isfield(seq, 'time')
    results.fps = seq.num_frames / seq.time;
else
    results.fps = NaN;
end

author: martin-danelljan

repo: ECO


text:
function element_node = grid_q9_element ( nelemx, nelemy )

%*****************************************************************************80
%
%% GRID_Q9_ELEMENT produces a grid of 9 node quadrilaterals.
%
%  Example:
%
%    Input:
%
%      NELEMX = 3, NELEMY = 2
%
%    Output:
%
%      ELEMENT_NODE =
%         1,  3, 17, 15,  2, 10, 16,  8,  9;
%         3,  5, 19, 17,  4, 12, 18, 10, 11;
%         5,  7, 21, 19,  6, 14, 20, 12, 13;
%        15, 17, 31, 29, 16, 24, 30, 22, 23;
%        17, 19, 33, 31, 18, 26, 32, 24, 25;
%        19, 21, 35, 33, 20, 28, 34, 26, 27.
%
%  Grid:
%
%   29---30---31---32---33---34---35
%    |    .    |    .    |    .    |
%    |    .    |    .    |    .    |
%   22 . 23 . 24 . 25 . 26 . 27 . 28
%    |    .    |    .    |    .    |
%    | 4  .    | 5  .    | 6  .    |
%   15---16---17---18---19---20---21
%    |    .    |    .    |    .    |
%    |    .    |    .    |    .    |
%    8 .  9 . 10 . 11 . 12 . 13 . 14
%    |    .    |    .    |    .    |
%    | 1  .    | 2  .    | 3  .    |
%    1----2----3----4----5----6----7
%
%  Reference Element Q9:
%
%    |
%    1  4--7--3
%    |  |     |
%    |  |     |
%    S  8  9  6
%    |  |     |
%    |  |     |
%    0  1--5--2
%    |
%    +--0--R--1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NELEMX, NELEMY, the number of elements along the
%    X and Y directions.  The number of elements generated will be
%    NELEMX * NELEMY.
%
%    Output, integer ELEMENT_NODE(9,NELEMX*NELEMY), the nodes that form
%    each element.
%

%
%  Node labeling:
%
%    NW----N----NE
%     |          |
%     W    C     E
%     |          |
%    SW----S----SE
%
  element = 0;

  for j = 1 : nelemy
    for i = 1 : nelemx

      sw = 2 * ( j - 1 )  * ( 2 * nelemx + 1 ) + 2 * ( i - 1 ) + 1;
      w  = sw +               2 * nelemx + 1;
      nw = sw +         2 * ( 2 * nelemx + 1 );

      s  = sw + 1;
      c  = sw + 1 +               2 * nelemx + 1;
      n  = sw + 1 +         2 * ( 2 * nelemx + 1 );

      se = sw + 2;
      e  = sw + 2 +               2 * nelemx + 1;
      ne = sw + 2 +         2 * ( 2 * nelemx + 1 );

      element = element + 1;

      element_node(1,element) = sw;
      element_node(2,element) = se;
      element_node(3,element) = ne;
      element_node(4,element) = nw;
      element_node(5,element) = s;
      element_node(6,element) = e;
      element_node(7,element) = n;
      element_node(8,element) = w;
      element_node(9,element) = c;

    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%PLAYREC Multichannel non-blocking audio recording and playback
%
% Run:
%     playrec('help')
%
% Tho get the actual help.
%
% PLAYREC is a MEX file which means it has to be compiled in order to work.
% Please see INSTALL-Matlab or INSTALL-Octave for details.
%
% Originally created by Robert Humphrey (see license_playrec.txt) and
% currently hosted at https://github.com/PlayrecForMatlab/playrec
% For further details visit the webpage or run playrec('about') and
% playrec('overview')
%
% This is however a sightly modified version doing an on-the-fly
% sample rate conversion if the target sample rate is not supported by
% the device.
%

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%--- help for transform_parameters ---
%
%  transform the parameters under linear restrictions and or dirichlet priors
%
%  ::
%
%    [obj,x0,lb,ub,vcov]=transform_parameters(obj,x0,lb,ub)
%    [obj,x0,lb,ub,vcov]=transform_parameters(obj,x0,lb,ub,vcov)
%
%  Args:
%
%     obj (rise | dsge | rfvar | svar): model object
%
%     x0 (empty | n x 1 vector): initial conditions of estimation in a space that
%       the user understands
%
%     lb (n x 1 vector): lower bound of the search space
%
%     ub (n x 1 vector): upper bound of the search space
%
%     vcov (empty | n x n matrix): variance covariance of the parameters
%
%  Returns:
%     :
%
%     - **obj** [rise|dsge|rfvar|svar]: model object with all restrictions set
%       up
%
%     - **x0** [empty|m x 1 vector]: transform initial conditions with m<=n
%
%     - **lb** [m x 1 vector]: transformed lower bound of the search space with
%       m<=n
%
%     - **ub** [m x 1 vector]: transformed upper bound of the search space with
%       m<=n
%
%     - **vcov** [empty|m x m matrix]: transformed variance covariance of the
%       parameters with m<=n
%
%  Note:
%
%     - Some checks have to be made after transformation in order to insure
%       that no element in the transformed lower bound exceeds its upper bound
%       counterpart.
%
%  See also:
%     - untransform_parameters
%
%

author: jmaih


text:
function prob_test150 ( )

%*****************************************************************************80
%
%% TEST150 tests UNIFORM_MEAN, UNIFORM_SAMPLE, UNIFORM_VARIANCE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST150\n' );
  fprintf ( 1, '  For the Uniform PDF:\n' );
  fprintf ( 1, '  UNIFORM_MEAN computes mean;\n' );
  fprintf ( 1, '  UNIFORM_SAMPLE samples;\n' );
  fprintf ( 1, '  UNIFORM_VARIANCE computes variance.\n' );

  a = 1.0;
  b = 10.0;

  check = uniform_check ( a, b );

  if ( ~check );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'TEST150 - Fatal error!\n' );
    fprintf ( 1, '  The parameters are not legal.\n' );
    return
  end

  mean = uniform_mean ( a, b );
  variance = uniform_variance ( a, b );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  PDF parameter A =     %14f\n', a );
  fprintf ( 1, '  PDF parameter B =     %14f\n', b );
  fprintf ( 1, '  PDF mean =            %14f\n', mean );
  fprintf ( 1, '  PDF variance =        %14f\n', variance );

  for i = 1 : nsample
    [ x(i), seed ] = uniform_sample ( a, b, seed );
  end

  mean = r8vec_mean ( nsample, x );
  variance = r8vec_variance ( nsample, x );
  xmax = max ( x(1:nsample) );
  xmin = min ( x(1:nsample) );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Sample size =     %6d\n', nsample );
  fprintf ( 1, '  Sample mean =     %14f\n', mean );
  fprintf ( 1, '  Sample variance = %14f\n', variance );
  fprintf ( 1, '  Sample maximum =  %14f\n', xmax );
  fprintf ( 1, '  Sample minimum =  %14f\n', xmin );

  return
end