text:
function  Metric_table  = MFSK_DEMOD(Data_in,M,E_s,SNR_dB)
%*** MFSK DEMOD ***
%function  Metric_table  = MFSK_DEMOD(Data_in,M,E_s,SNR_dB)
% Create a (M x n) matrix containing the metrics of the M-FSK demodulator,
% simulating an AWGN channel
% Data_in is random values [0 .. M-1], representing the data that must be
% transmitted
% M is the number of frequencies used (M-FSK)
% E_s is the transmitted signal energy
% SNR_dB is the signal-to-noise ratio, measured in dBs

    n = length(Data_in);

    %Calculate the mu and sigma_2 values from the SNR

    SNR = Convert_SNR(SNR_dB);
    N_o = E_s / SNR;        % Noise Power on channel in dB (Proakis - MFSK)

    mu = 0;
    sigma_2 = (1/2)*N_o;
    sigma = sqrt(sigma_2);

    %******************************
    %*** Create Decision matrix ***
    %******************************

    Metric_table = [];
    T_Metrics =[];

    %Create 2 x M metric outputs for timeslot i (M equals number of frequencies in M-FSK)
    for num = 1:n
        Signal = [];
        for i = 0:M-1
            if (i == Data_in(num))
                %It is the symbol that was transmitted

                %Determine a value for phi
                phi = 2 * pi * rand;

                % Compute the COS Component:
                Signal(i+1,1) = sqrt(E_s) * cos(phi) + normrnd(mu,sigma);

                % Compute the SIN Component
                Signal(i+1,2) = sqrt(E_s) * sin(phi) + normrnd(mu,sigma);

            else
                Signal(i+1,1) = normrnd(mu,sigma);
                Signal(i+1,2) = normrnd(mu,sigma);
            end
        end

        % Compute the Envelope metrics for each symbol
        Metrics = [];  % Vector containing the computed envelope metrics

        for i = 1:M
            Metrics(i) = Decision_metric(Signal(i,1),Signal(i,2));
        end

        Metric_table = [Metric_table transpose(Metrics)];
    end


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function createfigure_returns(x,y1, y2, y3,tit,legend1,legend2,legend3)
%CREATEFIGURE(X1,YMATRIX1)
%  x:  vector of xvalue data
%  y1,y2,y3:  vectors of densities yvalue data

% Create figure
figure('PaperSize',[20.98 29.68],...
    'Name',tit,'Color',[1 1 1]);

% Create axes
box('on');
hold('all');

ylow = min(min(y2),min(y3)); ylow = min(ylow,min(y1));
yhigh = max(max(y2),max(y3)); yhigh = max(yhigh, max(y1));

%ylow = -0.1;
%yhigh = 0.1;

subplot(3,1,1); plot(x,y1, 'black','DisplayName',legend1); xlabel('Time'); ylabel('Return Value');ylim([ylow yhigh]);title(legend1);
subplot(3,1,2); plot(x,y2, 'black', 'DisplayName',legend2); xlabel('Time'); ylabel('Return Value');ylim([ylow yhigh]);title(legend2);
subplot(3,1,3); plot(x,y3, 'black', 'DisplayName',legend3); xlabel('Time'); ylabel('Return Value');ylim([ylow yhigh]);title(legend3);

author: Sable

repo: mcbench-benchmarks


text:
%% HELLO_SCRIPT is a script to call HELLO_FUN.
%
%  Discussion:
%
%    The BATCH command runs scripts, not functions.  So we have to write
%    this short script if we want to work with BATCH!
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 May 2012
%
%  Author:
%
%    John Burkardt
%
  n = 10;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HELLO_SCRIPT\n' );
  fprintf ( 1, '  Call HELLO_FUN to say hello %d times\n', n );

  tic
  hello_fun ( n );
  toc

author: johannesgerer

repo: jburkardt-m


text:
%R2T Convert rotation matrix to a homogeneous transform
%
% T = R2T(R) is a homogeneous transform equivalent to an orthonormal
% rotation matrix R with a zero translational component.
%
% Notes::
% - Works for T in either SE(2) or SE(3)
%  - if R is 2x2 then T is 3x3, or
%  - if R is 3x3 then T is 4x4.
% - Translational component is zero.
% - For a rotation matrix sequence returns a homogeneous transform
%   sequence.
%
% See also T2R.

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

function T = r2t(R)

    % check dimensions: R is SO(2) or SO(3)
    d = size(R);
    if d(1) ~= d(2)
        error('matrix must be square');
    end
    if ~any(d(1) == [2 3])
        error('argument is not a rotation matrix (sequence)');
    end

    Z = zeros(d(1),1);
    B = [Z' 1];

    if numel(d) == 2
        % single matrix case
        T = [R Z; B];
    else
        %  matrix sequence case
        T = zeros(4,4,d(3));
        for i=1:d(3)
            T(:,:,i) = [R(:,:,i) Z; B];
        end
    end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function gX = lmcKernDiagGradX(kern, X)

% LMCKERNDIAGGRADX Gradient of LMC kernel's diagonal with respect to X.
% FORMAT
% DESC computes the gradient of the diagonal of the LMC kernel matrix with
% respect to the elements of the design matrix given in X.
% RETURN gX : the gradients of the diagonal with respect to each element of
% X. The returned matrix has the same dimensions as X.
% ARG kern : the kernel structure for which gradients are being
% computed.
% ARG X : the input data in the form of a design matrix.
%
% SEEALSO : lmcKernParamInit, kernDiagGradX
%
% COPYRIGHT : Mauricio A. Alvarez, 2010

% KERN

gX = zeros([kern.nout*size(X, 1) size(X,2)]);

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [diff_data] = timeSeriesComparison(t1,data1,t2,data2,mode)
% compare two time series the first one with lower spacing between epochs the second one with finer spacing
if strcmp(mode,'aggregate')
    edges = zeros(1,length(t2)+2+100);
    j = 1;
    rate = round(nan_mean(diff(t1))*86400)/86400;
    empty_idx = [];
    for i = 1:(length(t1) -1)
        edges(j) = t1(i) - (rate/2);
        j = j + 1;
        if ((t1(i+1) - t1(i)) > 1.5*rate)
            edges(j) = t1(i) + (rate/2);
            j = j + 1;
            empty_idx = [empty_idx; j];
        end
    end
    edges(j) = t1(end) + (rate/2);
    edges(j+1:end) = [];

    %edges = (data1.time.first.getMatlabTime-eps()) : 1:(data1.time.last.getMatlabTime+1);
    Y = discretize(t2,edges,'IncludedEdge','right');
    idx = false(size(Y));
    for ei = empty_idx'
        idx = idx | Y == ei;
    end
    idx = idx | isnan(Y);
    Y(idx) = [];
    ztds = zero2nan(data2(~idx));
    avg_data = accumarray(Y,ztds,[],@nan_mean);
    avg_data_time = edges(1:end-1) + (edges(2:end) - edges(1:end-1))/2;
    if length(avg_data_time) > max(Y)
    avg_data_time(max(Y)+1:end) = [];
    end
    avg_data_time = avg_data_time(avg_data~=0);
    avg_data = avg_data(avg_data~=0);
    [LIA,LocB] = ismembertol(avg_data_time, t1,1e-8);
    diff_data = nan(size(data1));
    diff_data(LocB(LocB~=0)) = zero2nan(avg_data(LIA)) - zero2nan(data1(LocB(LocB~=0)));

elseif strcmp(mode,'interpolate')
    diff_data = data2 - interp1(t1,data1,t2,'linear');
elseif strcmp(mode,'spline')
     edges = zeros(1,length(t2)+2+100);
    j = 1;
    rate = round(nan_mean(diff(t1))*86400)/86400;
    empty_idx = [];
    for i = 1:(length(t1) -1)
        edges(j) = t1(i) - (rate/2);
        j = j + 1;
        if ((t1(i+1) - t1(i)) > 1.5*rate)
            edges(j) = t1(i) + (rate/2);
            j = j + 1;
            empty_idx = [empty_idx; j];
        end
    end
    edges(j) = t1(end) + (rate/2);
    edges(j+1:end) = [];

    %edges = (data1.time.first.getMatlabTime-eps()) : 1:(data1.time.last.getMatlabTime+1);
    Y = discretize(t2,edges,'IncludedEdge','right');
    idx = false(size(Y));
    for ei = empty_idx'
        idx = idx | Y == ei;
    end
    idx = idx | isnan(Y);
    Y(idx) = [];
    avg_data_time = edges(1:end-1) + (edges(2:end) - edges(1:end-1))/2;
    if length(avg_data_time) > max(Y)
    avg_data_time(max(Y)+1:end) = [];
    end
    avg_data_time = avg_data_time(unique(Y));
    [LIA,LocB] = ismembertol(avg_data_time, t1,1e-7);

    [~,~,~, splined2] = splinerMat(t2, data2,rate,1e-5, avg_data_time);
    diff_data = nan(size(data1));
    diff_data(LocB(LocB~=0)) = zero2nan(splined2(LIA)) - zero2nan(data1(LocB(LocB~=0)));

end

author: goGPS-Project


text:
function [f] = VBA_psi(z)
% psi(x) = d[log(gamma(x))]/dx

try
    f = psi(z+1)-1./z; % for numerical purposes
catch
    siz = size(z);
    z=z(:);
    zz=z;
    f = 0.*z; % reserve space in advance
    %reflection point
    p=find(real(z)<0.5);
    if ~isempty(p)
        z(p)=1-z(p);
    end
    %Lanczos approximation for the complex plane
    g=607/128; % best results when 4<=g<=5
    c = [  0.99999999999999709182;
        57.156235665862923517;
        -59.597960355475491248;
        14.136097974741747174;
        -0.49191381609762019978;
        .33994649984811888699e-4;
        .46523628927048575665e-4;
        -.98374475304879564677e-4;
        .15808870322491248884e-3;
        -.21026444172410488319e-3;
        .21743961811521264320e-3;
        -.16431810653676389022e-3;
        .84418223983852743293e-4;
        -.26190838401581408670e-4;
        .36899182659531622704e-5];
    n=0;
    d=0;
    for k=size(c,1):-1:2
        dz=1./(z+k-2);
        dd=c(k).*dz;
        d=d+dd;
        n=n-dd.*dz;
    end
    d=d+c(1);
    gg=z+g-0.5;
    %log is accurate to about 13 digits...
    f = log(gg) + (n./d - g./gg) ;
    if ~isempty(p)
        f(p) = f(p)-pi*cot(pi*zz(p));
    end
    p=find(round(zz)==zz & real(zz)<=0 & imag(zz)==0);
    if ~isempty(p)
        f(p) = Inf;
    end
    f=reshape(f,siz);
end
return

author: MBB-team

repo: VBA-toolbox


text:
function [L,U,P] = lu_rightpr (A)
%LU_RIGHTPR recursive right-looking LU, with partial pivoting.
%
% Example:
%   [L,U,P] = lu_rightpr (A)
% See also: cs_demo

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

n = size (A,1) ;
if (n == 1)
    P = 1 ;
    L = 1 ;
    U = A ;
else
    [x,i] = max (abs (A (1:n,1))) ;                           % partial pivoting
    P1 = eye (n) ;
    P1 ([1 i],:) = P1 ([i 1], :) ;
    A = P1*A ;
    u11 = A (1,1) ;                                           % (6.10)
    u12 = A (1,2:n) ;                                         % (6.11)
    l21 = A (2:n,1) / u11 ;                                   % (6.12)
    [L22,U22,P2] = lu_rightpr (A (2:n,2:n) - l21*u12) ;       % (6.9) or (6.13)
    o = zeros(1,n-1) ;
    L = [ 1 o ; P2*l21 L22 ] ;                                % (6.14)
    U = [ u11 u12 ; o' U22 ] ;
    P = [ 1 o ; o' P2] * P1 ;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [net, info] = cnn_imagenet(varargin)
%CNN_IMAGENET   Demonstrates training a CNN on ImageNet
%  This demo demonstrates training the AlexNet, VGG-F, VGG-S, VGG-M,
%  VGG-VD-16, and VGG-VD-19 architectures on ImageNet data.

run(fullfile(fileparts(mfilename('fullpath')), ...
  '..', '..', 'matlab', 'vl_setupnn.m')) ;

opts.dataDir = fullfile(vl_rootnn, 'data','ILSVRC2012') ;
opts.modelType = 'alexnet' ;
opts.network = [] ;
opts.networkType = 'simplenn' ;
opts.batchNormalization = true ;
opts.weightInitMethod = 'gaussian' ;
[opts, varargin] = vl_argparse(opts, varargin) ;

sfx = opts.modelType ;
if opts.batchNormalization, sfx = [sfx '-bnorm'] ; end
sfx = [sfx '-' opts.networkType] ;
opts.expDir = fullfile(vl_rootnn, 'data', ['imagenet12-' sfx]) ;
[opts, varargin] = vl_argparse(opts, varargin) ;

opts.numFetchThreads = 12 ;
opts.lite = false ;
opts.imdbPath = fullfile(opts.expDir, 'imdb.mat');
opts.train = struct() ;
opts = vl_argparse(opts, varargin) ;
if ~isfield(opts.train, 'gpus'), opts.train.gpus = []; end;

% -------------------------------------------------------------------------
%                                                              Prepare data
% -------------------------------------------------------------------------

if exist(opts.imdbPath)
  imdb = load(opts.imdbPath) ;
  imdb.imageDir = fullfile(opts.dataDir, 'images');
else
  imdb = cnn_imagenet_setup_data('dataDir', opts.dataDir, 'lite', opts.lite) ;
  mkdir(opts.expDir) ;
  save(opts.imdbPath, '-struct', 'imdb') ;
end

% Compute image statistics (mean, RGB covariances, etc.)
imageStatsPath = fullfile(opts.expDir, 'imageStats.mat') ;
if exist(imageStatsPath)
  load(imageStatsPath, 'averageImage', 'rgbMean', 'rgbCovariance') ;
else
  train = find(imdb.images.set == 1) ;
  images = fullfile(imdb.imageDir, imdb.images.name(train(1:100:end))) ;
  [averageImage, rgbMean, rgbCovariance] = getImageStats(images, ...
                                                    'imageSize', [256 256], ...
                                                    'numThreads', opts.numFetchThreads, ...
                                                    'gpus', opts.train.gpus) ;
  save(imageStatsPath, 'averageImage', 'rgbMean', 'rgbCovariance') ;
end
[v,d] = eig(rgbCovariance) ;
rgbDeviation = v*sqrt(d) ;
clear v d ;

% -------------------------------------------------------------------------
%                                                             Prepare model
% -------------------------------------------------------------------------

if isempty(opts.network)
  switch opts.modelType
    case 'resnet-50'
      net = cnn_imagenet_init_resnet('averageImage', rgbMean, ...
                                     'colorDeviation', rgbDeviation, ...
                                     'classNames', imdb.classes.name, ...
                                     'classDescriptions', imdb.classes.description) ;
      opts.networkType = 'dagnn' ;

    otherwise
      net = cnn_imagenet_init('model', opts.modelType, ...
                              'batchNormalization', opts.batchNormalization, ...
                              'weightInitMethod', opts.weightInitMethod, ...
                              'networkType', opts.networkType, ...
                              'averageImage', rgbMean, ...
                              'colorDeviation', rgbDeviation, ...
                              'classNames', imdb.classes.name, ...
                              'classDescriptions', imdb.classes.description) ;
  end
else
  net = opts.network ;
  opts.network = [] ;
end

% -------------------------------------------------------------------------
%                                                                     Learn
% -------------------------------------------------------------------------

switch opts.networkType
  case 'simplenn', trainFn = @cnn_train ;
  case 'dagnn', trainFn = @cnn_train_dag ;
end

[net, info] = trainFn(net, imdb, getBatchFn(opts, net.meta), ...
                      'expDir', opts.expDir, ...
                      net.meta.trainOpts, ...
                      opts.train) ;

% -------------------------------------------------------------------------
%                                                                    Deploy
% -------------------------------------------------------------------------

net = cnn_imagenet_deploy(net) ;
modelPath = fullfile(opts.expDir, 'net-deployed.mat')

switch opts.networkType
  case 'simplenn'
    save(modelPath, '-struct', 'net') ;
  case 'dagnn'
    net_ = net.saveobj() ;
    save(modelPath, '-struct', 'net_') ;
    clear net_ ;
end

% -------------------------------------------------------------------------
function fn = getBatchFn(opts, meta)
% -------------------------------------------------------------------------

if numel(meta.normalization.averageImage) == 3
  mu = double(meta.normalization.averageImage(:)) ;
else
  mu = imresize(single(meta.normalization.averageImage), ...
                meta.normalization.imageSize(1:2)) ;
end

useGpu = numel(opts.train.gpus) > 0 ;

bopts.test = struct(...
  'useGpu', useGpu, ...
  'numThreads', opts.numFetchThreads, ...
  'imageSize',  meta.normalization.imageSize(1:2), ...
  'cropSize', meta.normalization.cropSize, ...
  'subtractAverage', mu) ;

% Copy the parameters for data augmentation
bopts.train = bopts.test ;
for f = fieldnames(meta.augmentation)'
  f = char(f) ;
  bopts.train.(f) = meta.augmentation.(f) ;
end

fn = @(x,y) getBatch(bopts,useGpu,lower(opts.networkType),x,y) ;

% -------------------------------------------------------------------------
function varargout = getBatch(opts, useGpu, networkType, imdb, batch)
% -------------------------------------------------------------------------
images = strcat([imdb.imageDir filesep], imdb.images.name(batch)) ;
if ~isempty(batch) && imdb.images.set(batch(1)) == 1
  phase = 'train' ;
else
  phase = 'test' ;
end
data = getImageBatch(images, opts.(phase), 'prefetch', nargout == 0) ;
if nargout > 0
  labels = imdb.images.label(batch) ;
  switch networkType
    case 'simplenn'
      varargout = {data, labels} ;
    case 'dagnn'
      varargout{1} = {'input', data, 'label', labels} ;
  end
end

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function y = boolean(x)

y = logical(x);

repo: gspbox


text:
function i4mat_write ( output_filename, m, n, table )

%*****************************************************************************80
%
%% I4MAT_WRITE writes an I4MAT file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 August 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the output filename.
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of points.
%
%    Input, integer TABLE(M,N), the points.
%
%    Input, logical HEADER, is TRUE if the header is to be included.
%

%
%  Open the file.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'I4MAT_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'I4MAT_WRITE - Error!' );
    return;
  end
%
%  Write the data.
%
  for j = 1 : n
    for i = 1 : m
      fprintf ( output_unit, '  %12d', round ( table(i,j) ) );
    end
    fprintf ( output_unit, '\n' );
  end
%
%  Close the file.
%
  fclose ( output_unit );

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   PARAMETERS Returns a data structure containing the parameters of the
%   ABB IRB140.
%
%   Author: Arturo Gil. Universidad Miguel Hernandez de Elche.
%   email: arturo.gil@umh.es date:   09/01/2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

robot.name= 'KUKA_KR20_3';

robot.DH.theta= '[q(1) q(2)-pi/2 q(3) q(4) q(5) q(6)+pi]';
robot.DH.d='[0.675 0 0 0.67 0 0.158]';
robot.DH.a='[0.26 0.68 -0.035 0 0 0]';
robot.DH.alpha= '[-pi/2 0 -pi/2 pi/2 -pi/2 0]';
robot.J=[];

robot.inversekinematic_fn = 'inversekinematic_kuka_kr20_3(robot, T)';
robot.directkinematic_fn = 'directkinematic(robot, q)';

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

%rotational: 0, translational: 1
robot.kind=['R' 'R' 'R' 'R' 'R' 'R'];

%minimum and maximum rotation angle in rad
robot.maxangle =[deg2rad(-185) deg2rad(185); %Axis 1, minimum, maximum
                deg2rad(-155) deg2rad(35); %Axis 2, minimum, maximum
                deg2rad(-130) deg2rad(154); %Axis 3
                deg2rad(-350) deg2rad(350); %Axis 4: Unlimited (400� default)
                deg2rad(-130) deg2rad(130); %Axis 5
                deg2rad(-350) deg2rad(350)]; %Axis 6: Really Unlimited to (800� default)

%maximum absolute speed of each joint rad/s or m/s
robot.velmax = [deg2rad(156); %Axis 1, rad/s
                deg2rad(156); %Axis 2, rad/s
                deg2rad(156); %Axis 3, rad/s
                deg2rad(330); %Axis 4, rad/s
                deg2rad(332); %Axis 5, rad/s
                deg2rad(616)];%Axis 6, rad/s

robot.accelmax=robot.velmax/0.1; % 0.1 is here an acceleration time

% end effectors maximum velocity
robot.linear_velmax = 2.5; %m/s

%base reference system
robot.T0 = eye(4);

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

% GRAPHICS
robot.graphical.has_graphics=1;
robot.graphical.color = [255 102 51]./255;
%for transparency
robot.graphical.draw_transparent=0;
%draw DH systems
robot.graphical.draw_axes=0;
%DH system length and Font size, standard is 1/10. Select 2/20, 3/30 for
%bigger robots
robot.graphical.axes_scale=1;
%adjust for a default view of the robot

robot.axis=[-1.1 1.7 -1.75 1.75 -0.5 2.2];
%read graphics files
robot = read_graphics(robot);

%DYNAMICS
robot.has_dynamics=0;

repo: ARTE


text: %OR Datafile overload

author: marianux

repo: ecg-kit

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


text:
report_this_filefun(mfilename('fullpath'));

% make a rate per year plot:
ed = (floor(min(newt2.Date)):1:ceil(max(newt2.Date)));

[ny,hy ] = hist(newt2.Date,ed-0.5);

figure
bar(hy,ny)

author: CelsoReyes


text:
% VL_TEST_INIT  Intialize unit test file
%   This is a VLFeat internal function used to initialize an .m file
%   containing a unit test.
%
%   See also: VL_TEST(), VL_HELP().

% Author:: Andrea Vedaldi

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% Copyright (C) 2013 Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

% Obtain the path to the m-file that called VL_TEST_INIT()
stack = dbstack(1) ;
mfilePath = stack(1).file ;
mfileName = stack(1).name ;
fid = fopen(mfilePath, 'r') ;

% Parse the source code of the caller to get the names of its
% subfunctions. One could use which('-subfun', stack(1).name) here,
% but this feature is undocumented.

setupFunc = [] ;
teardownFunc = [] ;

tests = {} ;
l = 0 ;
t = 0 ;
while 1
  tline = fgetl(fid) ;
  l = l + 1 ;
  if ~ischar(tline), break ; end
  mo = regexp(tline, '^\s*function\s+(\w\s*=\s*)?test_(?<name>[\w_]+)\s*\(.*\).*$', ...
              'names') ;
  if length(mo) > 0
    t = t + 1 ;
    tests{t}.name = mo(1).name ;
    tests{t}.func = str2func(['test_' mo(1).name]) ;
    tests{t}.lineNumber = l ;
  end

  mo = regexp(tline, '^\s*function\s+(\w+\s*=\s*)?(?<name>setup|teardown)\s*\(.*\).*$', ...
              'names') ;
  if length(mo) > 0
    switch mo(1).name
      case 'setup'
        if isempty(setupFunc)
          setupFunc = @setup ;
        else
          warning('%s:%d: ignoring redundant definition of the setup function', ...
                  mfilePath) ;
        end
      case 'teardown'
        if isempty(teardownFunc)
          teardownFunc = @teardown ;
        else
          warning('%s:%d: ignoring redundant definition of the setup function', ...
                  mfilePath) ;
        end
      otherwise
        assert(0) ;
    end
  end
end

if ~isempty(setupFunc), state = setupFunc() ; end

if length(varargin) > 0
  testFunc = str2func(sprintf('test_%s', varargin{1})) ;
  if exist('state','var')
    testFunc(state) ;
  else
    testFunc() ;
  end
  if ~isempty(teardownFunc), teardownFunc(state) ; end
  return ;
end

% Run all tests collecting the results.
if length(tests) == 0
  error('No tests found in ''%s''', mfileName)  ;
end

for t = 1:length(tests)
  subTestName = sprintf('test_%s', tests{t}.name) ;
  results{t}.testName = sprintf('%s:%s', mfileName(9:end), tests{t}.name) ;
  results{t}.succeded = true ;
  results{t}.exception = [] ;
  try
    if exist('state', 'var')
      tests{t}.func(state) ;
    else
      tests{t}.func() ;
    end
  catch exception
    results{t}.succeded = false ;
    results{t}.exception = exception ;
  end
end

if ~isempty(teardownFunc), teardownFunc(state) ; end

results = cat(2, results{:}) ;

author: jianxiongxiao

repo: ProfXkit


text:
function polyOut = parsePolygon(poly, format, varargin)
%PARSEPOLYGON Conversion between different polygon formats.
%
%   POLYOUT = parsePolygon(POLY, FORMAT)
%   Converts POLY to the specified format while FORMAT must be one of the
%   following options: 'cell','nan','repetition','polyshape'. Mainly used
%   in other polygon functions to bring the polygon into the required
%   format.
%
%   Example
%     r = [2.5, 2, 1];
%     poly = flipud(circleToPolygon([0 0 r(1)], round(2*pi*r(1))));
%     midCircle = circleToPolygon([0 0 r(2)], round(2*pi*r(2)));
%     innerCircle = flipud(circleToPolygon([0 0 r(3)], round(2*pi*r(3))));
%     poly = {poly, midCircle, innerCircle};
%     polyOut = parsePolygon(poly, 'repetition');
%     polyOut = parsePolygon(polyOut, 'nan');
%     polyOut = parsePolygon(polyOut, 'polyshape');
%     polyOut = parsePolygon(polyOut, 'cell');
%     assert(isequal(poly, polyOut))
%
%   See also
%   polygonToPolyshape
%

% ------
% Author: oqilipo
% E-mail: N/A
% Created: 2023-01-01, using MATLAB 9.13.0.2080170 (R2022b) Update 1
% Copyright 2023

% Parsing
p = inputParser;
addRequired(p,'format',...
    @(x) any(validatestring(x,{'cell','nan','repetition','polyshape'})));
parse(p,format,varargin{:});

format = p.Results.format;

% If polygon is a polyshape object
if isa(poly,'polyshape')
    polyOut = parsePolygon(poly.Vertices, format);
    return
end

% If polygon is a cell
if iscell(poly)
    [r, c] = size(poly);
    if r > 1 && c > 1
        error(['Non-supported polygon cell format with ' num2str(r) ' rows and ' num2str(c) 'columns!'])
    elseif r > 1 && c == 1
        polyCell = poly';
    elseif r == 1 && c > 1
        polyCell = poly;
    elseif r == 1 && c == 1
        polyCell = poly;
    end
    if any(cellfun(@(x) any(isnan(x(:))), poly))
        error('Non-supported polygon cell format. Polygons in cell format should not contain NaN values!')
    end
end

if ~iscell(poly) && any(isnan(poly(:)))
    % If polygon contains nan
    polyCell = splitPolygons(poly)';
elseif ~iscell(poly)
    % If polygon contains vertex repetitions
    polyBackup = poly;
    polyCell = {};
    count = 0;
    while ~isempty(poly)
        count = count+1;
        repLIdx = ismember(poly, poly(1,:), 'rows');
        repIdx = find(repLIdx);
        if sum(repLIdx) == 1
            polyCell = {poly};
            if count == 1
                break
            else
                warning('Contains vertex repetitions but the last vertex is not a repetiton!')
                polyCell = {polyBackup};
            end
        elseif sum(repLIdx)>1
            polyCell = [polyCell, {poly(1:repIdx(2)-1,:)}]; %#ok<AGROW>
            poly(1:repIdx(2),:) = [];
        end
    end
end

% Convert to output format
switch format
    case 'cell'
        polyOut = polyCell;
    case 'nan'
        polyOut = cell2mat((cellfun(@(x) [x; nan(1, size(x,2))], polyCell, 'UniformOutput',false))');
        polyOut(end,:) = [];
    case 'repetition'
        polyOut = cell2mat((cellfun(@(x) [x; x(1,:)], polyCell, 'UniformOutput',false))');
    case 'polyshape'
        polyOut = polyshape(polyCell{1});
        NoP = length(polyCell);
        if NoP > 1
            for p=2:NoP
                polyOut = addboundary(polyOut, polyCell{p});
            end
        end
end

end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
im = iread('tomato_124.jpg', 'gamma', 'sRGB', 'double');
figure(1); idisp(im)
randinit
[cls, cxy] = colorkmeans(im, 4);
figure(2); idisp(cls)
cxy

cls2 = (cls == 2);
figure(1); idisp(cls2);
binary_im = iclose(cls2, kcircle(15));
figure(2); idisp(binary_im)

f = iblobs(binary_im, 'boundary', 'class', 1)

idisp(im)
f.plot_boundary('r.')
f.plot_box('g')

author: petercorke

repo: machinevision-toolbox-matlab


text:
function plotnsgtf(c,shift,sr,fmin,fmax,bins,cutout,dynrange)
%PLOTNSGTF  Plot nonstationary Gabor filterbank coefficients
%   Usage:  plotnsgtf(c,shift,sr,fmin,fmax,bins,cutout,dynrange)
%           plotnsgtf(c,shift,sr,fmin,fmax,bins,cutout)
%           plotnsgtf(c,shift,sr,fmin,fmax,bins)
%           plotnsgtf(c,shift,sr,cutout,dynrange)
%           plotnsgtf(c,shift,sr,cutout)
%           plotnsgtf(c,shift,sr)
%           plotnsgtf(c,shift)
%           plotnsgtf(c)
%
%   Input parameters:
%         c        : Array of coefficients.
%         shift    : Vector of frequency shifts
%         sr       : signal sample rate in Hz (default 1 Hz)
%         fmin     : Minimum frequency used in the transform
%         fmax     : Maximum frequency used in the transform
%         bins     : Bins per octave (in constant or vector form)
%         cutout   : Desired part of the spectrogram, e.g.
%                    choice of 2 shows frequencies up to Nyquist
%                    (X shows the number_of_bins/X lowest frequency
%                    bins)
%         dynrange : Colorscale dynamic range in dB (default 60 dB)
%
%   Given a coefficient array c and the frequency shift vector shift,
%   this function plots the dB-scaled nonstationary Gabor filterbank
%   spectrogram corresponding to c. The vector shift and sampling rate
%   sr are used to determine the correct time axis labels. The frequency
%   axis is by default labeled by bin number.
%
%   For constant-Q nonstationary Gabor filterbanks, labeling with the
%   actual center frequencies is supported, requiring the filterbank
%   parameters fmin, fmax and bins as additional input.
%
%   The shown frequency range can be controlled with the cutout parameter
%   (default: 2) and the dynamic range of the spectrogram can be adjusted
%   with dynrange.
%
%   See also:  nsgtf, plotnsgt, plotslicq
%
%   Url: http://nsg.sourceforge.net/doc/plotting/plotnsgtf.php

% Copyright (C) 2013 Nicki Holighaus.
% This file is part of NSGToolbox version 0.1.0
%
% This work is licensed under the Creative Commons
% Attribution-NonCommercial-ShareAlike 3.0 Unported
% License. To view a copy of this license, visit
% http://creativecommons.org/licenses/by-nc-sa/3.0/
% or send a letter to
% Creative Commons, 444 Castro Street, Suite 900,
% Mountain View, California, 94041, USA.

% Author:  Gino Velasco, Nicki Holighaus and Radu C. Frunza
% Original code by Florent Jaillet
% Date: 26.04.13

ticklabels = 0;

if nargin < 8
    % Default value for colorscale dynamic.
    dynrange = 60;
    if nargin < 7 % Default value for frequency cutout.
        cutout = 2;
        if nargin < 6
            ticklabels = 0;
            if nargin == 5
                cutout = fmin;
                dynrange = fmax;
            elseif nargin < 4
                cutout = 2;
                if nargin < 3
                    % Default sampling rate
                    sr = 1;
                    if nargin < 2
                        error('Not enough input arguments');
                    end
                end
            else
                cutout = fmin;
            end
        end
    end
end

N = length(shift)-2;
L=sum(shift);

% For the coefficient output of nsgtf_real, adjust the frequeny
% channels shown, if necessary

if N > size(c,1) && cutout < 2
    cutout = cutout/2;
end

cla %Clear previous axes

% Compute maximum of the representation for colorscale dynamic handling.
if iscell(c) == 1
    if size(c{1},2) > 1
        error(['Multichannel spectrograms are not supported. Please ',...
            'use ''cellfun(@(x) x(:,k),c,''UniformOutput'',0)'' to ',...
            'select the k-th channel.']);
    end
    temp=cell2mat(c);
else
    if size(c,3) > 1
        error(['Multichannel spectrograms are not supported. Please ',...
            'use ''c(:,:,k)'' to select the k-th channel.']);
    end
    temp=c;
end
ma=20*log10(max(abs(temp(:))));

% Plot the representation: as the sampling grid in the time frequency plane
% is irregular, the representation by done by plotting many images next to
% each other, with one image for each window
if iscell(c) == 1
    hold('on');
    for ii=1:floor((N-1)/cutout)+1
        temp = 20*log10(abs(c{ii})+eps);
        % +eps is here to avoid log of 0
        % Octave cannot plot images that are only one point wide, so we use
        % images that are to points wide
        imagesc([0,L-1]/sr,[ii,ii+1],[temp,temp].',...
            [ma-dynrange,ma]);
    end
    hold('off');
    axis('tight');
else
    imagesc([0, (L-1)/sr],[1 size(c,2)],...
        20*log10(abs(c)'+eps),[ma-dynrange, ma]);
    select=[0,(L-1)/sr,1,round(size(c,2)/cutout)+1];
    axis(select);
    set(gca,'YDir','normal');
end

% Compute YTickLabels for the correct frequencies
if ticklabels == 1

    vfq = [0, fmin*2.^(0:log2(fmax/fmin))];

    if length(bins) < length(vfq) - 1
        xbins = [bins, bins(end)*ones(1,length(vfq)-length(bins)-1)];
    else
        xbins = bins(1:length(vfq)-1);
    end;

    sbins = [0,2,1+cumsum(xbins)];

    yTick = [2;2+2*cumsum(xbins(1:end-1))';(1+N/2)];
    yTick = unique([yTick(yTick <= floor((length(shift)-1)/cutout)+1);...
        floor((N-1)/cutout)+1]);

    yTickLabel = zeros(length(yTick),1);

    for kk = 1:length(yTick)
        ind = find(yTick(kk) < sbins,1);
        yTickLabel(kk) = vfq(ind-1)*...
            2^(((yTick(kk)-1)-sbins(ind-1))/(sbins(ind)-sbins(ind-1)));
    end;

    yTickLabel = num2str(round(yTickLabel),5);
    set(gca,'YTick',yTick,'YTickLabel',yTickLabel);
end

author: nesl


text:
function mesh = mmgRunLs(mmg)
% Copyright (c) 2018-2019
% Matthieu Aussal, CMAP, Ecole Polytechnique
% Algiane Froehly, CARDAMOME, INRIA-SOFT
% LGPL Lesser General Public License v3.0.
% Remeshing using Mmg tools : https://www.mmgtools.org

% Current and mmg directory
%here  = pwd;
there = which('mmg.m');
there = there(1:end-6);

oldMeshFileName = mmg.oldMeshFileName;
newMeshFileName = mmg.newMeshFileName;
levelSetOldFileName = mmg.levelSetOldFileName;
levelSetNewFileName = mmg.levelSetNewFileName;

% Move to mmg directory
%cd(there)

% Check mesh (3D nodes, triangle or tetra, colours are int)
mesh = mmg.Mesh;
if (size(mesh.vtx,2) ~= 3) || (size(mesh.elt,2) < 3) || (size(mesh.elt,2) > 4)
    error('mmgRun.m : unavailable case, please use 3D vertices.')
end

% Check colours are int
if (norm(mesh.col-floor(mesh.col),'inf') > 1e-12) || (min(mesh.col)<0)
    error('mmgRun.m : only integer values for colours.')
end

% Delete .mesh and .sol files
deleteFiles(oldMeshFileName,newMeshFileName,levelSetOldFileName,levelSetNewFileName)

%figure(1)
s.connec  = mesh.elt();
s.coord   = mesh.vtx(:,1:2);
mN1 = Mesh().create(s);
%mN1.plot();
msRelator = MasterSlaveRelator(s.coord);
masterSlave = msRelator.getRelation();
cell = CellNodesDescriptor(s.coord);
n1 = masterSlave(:);
n2 = cell.cornerNodes;
nT = [n1;n2];
mN1.computeEdges();
nodeInBoundaryEdges = mN1.edges.computeBoundaryEdges(nT);
edges = nodeInBoundaryEdges;
edges(:,3) = 0;

nodes = find(edges(:,3) ~= 10);

% Write original mesh in .mesh format
mmgMeshWrite(oldMeshFileName,mesh,mmg.Req,nodes,edges);

% Write size map in .sol format
if ~isempty(mmg.Map)
    mmgSolWrite(levelSetOldFileName,mmg.Map,is2d(mesh));
end

% Operating system
if ismac
    os = './mac_';
elseif ispc
    os = 'win_';
elseif isunix
    os = '/./uni_';
else
    disp('mmg.m : unavailable case')
end
os = [there,os];

% Element form
if (size(mesh.elt,2) == 3)
    if is2d(mesh)
        bin = 'mmg2d_O3';
    else
        bin = 'mmgs_O3';
    end
elseif (size(mesh.elt,2) == 4)
    bin = 'mmg3d_O3';
end

% Add extension (for windows only)
if ispc
    bin = [bin,'.exe '];
else
    bin = [bin,' '];
end

% Define temporary mesh files
file = ['-in ',oldMeshFileName,' -out ',newMeshFileName,' '];

% Convert input to command
field = fieldnames(mmg);
opt   = '';
for i = 1:length(field)
    cmd = getfield(mmg,field{i});
    if ~isempty(cmd) && ~strcmp(field{i},'Mesh') && ~strcmp(field{i},'Req')...
       && ~strcmp(field{i},'Map') ...
       && ~strcmp(field{i},'oldMeshFileName') && ~strcmp(field{i},'newMeshFileName') ...
       && ~strcmp(field{i},'levelSetOldFileName') && ~strcmp(field{i},'levelSetNewFileName')
       %&& ~strcmp(field{i},'original') && ~strcmp(field{i},'refined') ...
        opt = [opt,cmd,' '];
    end
end

if ~isempty(mmg.Map)
    cmd = ['-sol ',levelSetOldFileName];
    opt = [opt,cmd,' '];
    cmd = ['-ls 0'];
    opt = [opt,cmd,' '];
end

% Execute mmg binaries
command = [os,bin,file,opt];
[~,~] = system(command);

% Read refined mesh
[remesh,req,label,edges] = mmgMeshRead(newMeshFileName);

% Replace required entites by original ones
if ~isempty(req)
    % Read original mesh for double/single troubles
    origin = mmgMeshRead(oldMeshFileName);

    % Vertex indices from required elements
    I = mesh.elt(mmg.Req,:);
    J = remesh.elt(req,:);

    % Unicity of vertices
    I = unique(I(:));
    J = unique(J(:));

    % Appairement using knnsearch with foat original mesh (to be improved)
    idx = knnsearch(origin.vtx(I,:),remesh.vtx(J,:),'K',1);

    % Replace remesh required vertices using original mesh vertices
    remesh.vtx(J,:) = mesh.vtx(I(idx),:);
end

% Save output
mesh = remesh;

% Delete .mesh and .sol files
deleteFiles(oldMeshFileName,newMeshFileName,levelSetOldFileName,levelSetNewFileName);

% Back to current directory
%cd(here)
end

function deleteFiles(oldMeshFileName,newMeshFileName,levelSetOldFileName,levelSetNewFileName)
if exist(oldMeshFileName,'file')
    delete(oldMeshFileName);
end
if exist(newMeshFileName,'file')
    delete(newMeshFileName);
end
if exist(levelSetOldFileName,'file')
    delete(levelSetOldFileName);
end
if exist(levelSetNewFileName,'file')
    delete(levelSetNewFileName);
end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [s] = readstruct(filename, defaults)
% readstruct Read a key-value file to a structure
%
% This function reads a key-value file to a structure. If possible, the
% values are converted to numbers, otherwise the value is kept as string.
%
% Input:
% - filename (string): Path to the file.
% - defaults (struct): A structure merged with read data, can be used to
%   set defaults.
%
% Output:
% - s (structure): Resulting structure.
%

fid = fopen(filename);
lines = textscan(fid, '%s', 'Delimiter', '\n');
fclose(fid);

lines = lines{1};

s = struct();

for i = 1:numel(lines)
	[key, value] = parse_line(lines{i});

	key = strrep(strrep(key, ' ', '_'), '-', '_');

	if ~isempty(key)

        tokens = strsplit(key, '.');

        s = update_struct(s, tokens, value);

	end

end;

if nargin > 1
	s = struct_merge(s, defaults);
end

end

function s = update_struct(s, keys, value)

    if numel(keys) == 1
        s.(keys{1}) = value;
    else
        if ~isfield(s, keys{1}) || ~isstruct(s.(keys{1}))
            s.(keys{1}) = struct();
        end
        s.(keys{1}) = update_struct(s.(keys{1}), keys(2:end), value);
    end

end

function [key, value] = parse_line(line)

    delimiter = find(line == '=', 1);

	if isnan(delimiter)
		key = []; value = [];
		return;
	end;

	key = line(1:delimiter-1);
	value = line(delimiter+1:end);

	[numvalue, status] = str2num(value); %#ok<ST2NM>

	if status
		value = numvalue;
	end

end

author: votchallenge

repo: toolkit-legacy


text:
function y = callback_doubleoperator(x,dir,options)

% callback_doubleoperator - apply two calbacks
%
%   y = callback_doubleoperator(x,dir,options);
%
%   Copyright (c) 2008 Gabriel Peyre

c1 = getoptions(options, 'callback1', 0, 1);
c2 = getoptions(options, 'callback2', 0, 1);

if dir==1
    y = feval(c1, x, dir, options);
    y = feval(c2, y, dir, options);
else
    y = feval(c2, x, dir, options);
    y = feval(c1, y, dir, options);
end

author: gpeyre

repo: matlab-toolboxes


text:
function [mri] = gui_avw_open(mri,command,parent)

% gui_avw_open - Load & Display Analyze 7.5 data
%
% Useage: [mri] = gui_avw_open(mri,[command],[parent])
%
% mri is a structure, generated by 'mri_toolbox_defaults'
% command is either 'init' or 'load'
% parent is a handle to the gui that calls this gui, useful
% for updating the UserData field of the parent from this gui.
% The mri structure may be returned to the parent when the parent
% handle is given.
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:56 $

% Licence:  GNU GPL, no express or implied warranties
% History:  08/2002, Darren.Weber@flinders.edu.au
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~exist('mri','var'),
    mri = mri_toolbox_defaults;
elseif isempty(mri),
    mri = mri_toolbox_defaults;
end

if ~exist('command','var'),
    command = 'init';
elseif isempty(command),
    command = 'init';
end

command = lower(command);

switch command,
case 'init',
    if exist('parent','var'),
        MRIOpen = INIT(mri,parent);
    else
        MRIOpen = INIT(mri,'');
    end
otherwise,
    MRIOpen = get(gcbf,'Userdata');
    set(MRIOpen.gui,'Pointer','watch');
end

switch command,

case 'plot',

    MRIOpen.mri.plot = 1;
    MRIOpen.mri = mri_open(MRIOpen.mri);

case 'return',

    MRIOpen.mri.plot = 0;
    MRIOpen.mri = mri_open(MRIOpen.mri);

case 'save',

    fprintf('\ngui_avw_open: Save As not implemented yet.\n');

otherwise,

end

% -- tidy up & return

switch command,
case 'init',
case 'cancel',
    GUI.parent = MRIOpen.parent;
    mri_updateparent(GUI);
    close gcbf;
otherwise,
    set(MRIOpen.gui,'Pointer','arrow');
    set(MRIOpen.gui,'Userdata',MRIOpen);

    mri = mri_updateparent(MRIOpen,0);

    if isequal(get(MRIOpen.handles.Bhold,'Value'),0),
        close gcbf;
        if isfield(MRIOpen,'parent'),
            parent = MRIOpen.parent.gui;
        else
            parent = [];
        end
    else
        parent = MRIOpen.gui;
    end
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [MRIOpen] = INIT(mri,parent)
    % GUI General Parameters

    GUIwidth  = 500;
    GUIheight = 120;

    version = '$Revision: 1.1 $';
    name = sprintf('AVW File Open [v %s]\n',version(11:15));

    GUI = figure('Name',name,'Tag','AVW_OPEN',...
                 'NumberTitle','off',...
                 'MenuBar','none','Position',[1 1 GUIwidth GUIheight]);
    movegui(GUI,'center');

    Font.FontName   = 'Helvetica';
    Font.FontUnits  = 'Pixels';
    Font.FontSize   = 12;
    Font.FontWeight = 'normal';
    Font.FontAngle  = 'normal';

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Voltage Data Selection and Parameters

    G.Title_data = uicontrol('Parent',GUI,'Style','text','Units','Normalized',Font, ...
        'Position',[.01 .75 .17 .2],...
        'String','Data Type:','HorizontalAlignment','left');

    mri.type = 'Analyze';
    G.PmriType = uicontrol('Tag','PmriType','Parent',GUI,'Style','edit',...
        'Units','Normalized',Font,  ...
        'Position',[.20 .75 .25 .2],...
        'String',mri.type);

    switch mri.orient
        case 'auto',                mriOrient = 1;
        case 'axial unflipped',     mriOrient = 2;
        case 'coronal unflipped',   mriOrient = 3;
        case 'sagittal unflipped',  mriOrient = 4;
        case 'axial flipped',       mriOrient = 5;
        case 'coronal flipped',     mriOrient = 6;
        case 'sagittal flipped',    mriOrient = 7;
        otherwise,                  mriOrient = 1;
    end

    G.PmriOrient = uicontrol('Tag','PmriOrient','Parent',GUI,'Style','popupmenu',...
        'Units','Normalized',Font,  ...
        'Position',[.50 .75 .25 .2],...
        'TooltipString','Force loading with this orientation, otherwise use auto detection',...
        'String',{'auto' 'axial unflipped' 'coronal unflipped' 'sagittal unflipped' ...
                  'axial flipped'  'coronal flipped'  'sagittal flipped'},'Value',mriOrient,...
        'Callback',strcat('MRIOpen = get(gcbf,''Userdata'');',...
                          'MRIOpen.mri.orient = popupstr(MRIOpen.handles.PmriOrient);',...
                          'set(gcbf,''Userdata'',MRIOpen); clear MRIOpen;'));

    G.Title_path = uicontrol('Parent',GUI,'Style','text','Units','Normalized',Font, ...
        'Position',[.01 .50 .17 .2],...
        'String','Path','HorizontalAlignment','left');
    G.EmriPath = uicontrol('Parent',GUI,'Style','edit','Units','Normalized',Font,  ...
    	'Position',[.20 .50 .58 .2], 'String',mri.path,...
        'Callback',strcat('MRIOpen = get(gcbf,''Userdata'');',...
                          'MRIOpen.mri.path = get(MRIOpen.handles.EmriPath,''String'');',...
                          'set(gcbf,''Userdata'',MRIOpen); clear MRIOpen;'));

    G.Title_file = uicontrol('Parent',GUI,'Style','text','Units','Normalized',Font, ...
        'Position',[.01 .25 .17 .2],...
        'String','File','HorizontalAlignment','left');
    G.EmriFile = uicontrol('Parent',GUI,'Style','edit','Units','Normalized',Font,  ...
    	'Position',[.20 .25 .58 .2], 'String',mri.file,...
        'Callback',strcat('MRIOpen = get(gcbf,''Userdata'');',...
                          'MRIOpen.mri.file = get(MRIOpen.handles.EmriFile,''String'');',...
                          'set(gcbf,''Userdata'',MRIOpen); clear MRIOpen;'));

    Font.FontWeight = 'bold';

    % BROWSE: Look for the data
    browsecommand = strcat('MRIOpen = get(gcbf,''Userdata'');',...
        'cd(MRIOpen.mri.path);',...
        '[file, path] = uigetfile(',...
        '{''*.hdr;*.img'', ''Analyze Files (*.hdr,*.img)'';', ...
        ' ''*.*'',   ''All Files (*.*)''},', ...
        '''Select MRI File'');',...
        'if ~isequal(path,0), MRIOpen.mri.path = path; end;',...
        'if ~isequal(file,0), MRIOpen.mri.file = file; end;',...
        'set(MRIOpen.handles.EmriPath,''String'',MRIOpen.mri.path);',...
        'set(MRIOpen.handles.EmriFile,''String'',MRIOpen.mri.file);',...
        'if ~isempty(findstr(file,''.hdr'')) | ~isempty(findstr(file,''.img'')), ',...
            'set(MRIOpen.handles.PmriType,''Value'',1); ',...
            'MRIOpen.mri.type = popupstr(MRIOpen.handles.PmriType);',...
        'else, ',...
            'set(MRIOpen.handles.PmriType,''Value'',3); ',...
            'MRIOpen.mri.type = popupstr(MRIOpen.handles.PmriType);',...
        'end;',...
        'set(gcbf,''Userdata'',MRIOpen); clear MRIOpen file path;');
    G.BmriFile = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized',Font, ...
        'Position',[.01 .01 .17 .2], 'String','BROWSE',...
        'BackgroundColor',[0.8 0.8 0.0],...
        'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
        'Callback', browsecommand );

    % PLOT: Load & plot the data!
    G.Bplot = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
        'Position',[.20 .01 .18 .2],...
        'String','PLOT','BusyAction','queue',...
        'TooltipString','Plot the MRI data and return p struct.',...
        'BackgroundColor',[0.0 0.5 0.0],...
        'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
        'Callback',strcat('MRIOpen = get(gcbf,''Userdata'');',...
                          'mri = gui_avw_open(MRIOpen.mri,''plot'');',...
                          'clear MRIOpen;'));

    % Save As
	G.Bsave = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
        'Position',[.40 .01 .18 .2],'HorizontalAlignment', 'center',...
        'String','SAVE AS','TooltipString','MRI File Conversion Tool (not implemented yet)',...
        'BusyAction','queue',...
        'Visible','on',...
        'BackgroundColor',[0.0 0.0 0.75],...
        'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
        'Callback',strcat('MRIOpen = get(gcbf,''Userdata'');',...
                          'mri = gui_avw_open(MRIOpen.mri,''save'');',...
                          'clear MRIOpen;'));

    % Quit, return file parameters
    G.Breturn = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
        'Position',[.60 .01 .18 .2],...
        'String','RETURN','BusyAction','queue',...
        'TooltipString','Return p struct to workspace and parent GUI.',...
        'BackgroundColor',[0.75 0.0 0.0],...
        'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
        'Callback',strcat('MRIOpen = get(gcbf,''Userdata'');',...
                          'mri = gui_avw_open(MRIOpen.mri,''return'');',...
                          'clear MRIOpen;'));

    % Cancel
    G.Bcancel = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
        'Position',[.80 .01 .18 .2],...
        'String','CANCEL','BusyAction','queue',...
        'TooltipString','Close, do not return parameters.',...
        'BackgroundColor',[0.75 0.0 0.0],...
        'ForegroundColor', [1 1 1], 'HorizontalAlignment', 'center',...
        'Callback',strcat('MRIOpen = get(gcbf,''Userdata'');',...
                          'mri = gui_avw_open(MRIOpen.mri,''cancel'');',...
                          'clear MRIOpen;'));

	% Help
	G.Bhelp = uicontrol('Parent',GUI,'Style','pushbutton','Units','Normalized', Font, ...
        'Position',[.80 .25 .18 .2],'String','Help','BusyAction','queue',...
        'BackgroundColor',[1 1 0],...
        'ForegroundColor',[0 0 0], 'HorizontalAlignment', 'center',...
        'Callback','doc mri_toolbox;');

    % Hold GUI Open checkbox
	G.Bhold = uicontrol('Parent',GUI,'Style','checkbox','Units','Normalized', Font, ...
        'Position',[.80 .50 .18 .2],'String','Hold GUI','BusyAction','queue',...
        'TooltipString','MRI File Load GUI remains open after ''Plot'' or ''Return'' commands.',...
        'Value',mri.hold,'HorizontalAlignment', 'center');

    % Store userdata
    if exist('parent','var'), MRIOpen.parent.gui = parent; end
    MRIOpen.gui = GUI;
    MRIOpen.handles = G;
    MRIOpen.mri = mri;
    set(GUI,'Userdata',MRIOpen);
    set(GUI,'HandleVisibility','callback');

return

author: PatternRecognition

repo: OpenBMI


text:
function [U,UT] = isometric_curve_flow(V,varargin)
  % ISOMETRIC_CURVE_FLOW Isometric flow for curves as described by "Robust
  % Fairing via Conformal Curvature Flow" by [Crane et al. 2013]
  %
  % U = isometric_curve_flow(V)
  % U = isometric_curve_flow(V,'ParameterName',ParameterValue,...)
  %
  % Inputs:
  %   V  #V by 2 list of vertex positions
  %   Optional:
  %     'Tao' followed by timestep tao
  % Outputs:
  %   U  #V by 2 list of vertex positions
  %

  % time step
  tao = 1e-2;

  ii = 1;
  while ii < numel(varargin)
    switch varargin{ii}
    case 'Tao'
      assert(ii+1<=numel(varargin));
      ii = ii+1;
      tao = varargin{ii};
    otherwise
      error('Unsupported parameter: %s',varargin{ii});
    end
    ii = ii+1;
  end

  % number of points
  n = size(V,1);

  % Evaluate curvature and edge vectors
  [kappa,alpha,ev,l] = curvature(V);

  % desired flow direction
  kappa_dot = - 2 * kappa;

  % defining edge length at each point
  egde_vector = zeros(n,1);
  for k=2:n-1
    edge_vector(k) = 0.5*(norm(V(k+1,:)-V(k,:))+norm(V(k-1,:)-V(k,:)));
  end
  edge_vector(1) = 0.5*(norm(V(2,:)-V(1,:))+norm(V(size(V,1),:)-V(1,:)));
  edge_vector(n) = 0.5*(norm(V(1,:)-V(n,:))+norm(V(n,:)-V(n-1,:)));
  % mass matrix
  M = diag(sparse(0.5*(l([end 1:end-1]) + l)));

  % re-defining orthonormal basis
  C(:,1) = ones(n,1)/sqrt(ones(n,1)'*M*ones(n,1));

  C(:,2) = V(:,1) - ((V(:,1)'*M*C(:,1)))*C(:,1);
  C(:,2) = C(:,2)/sqrt(C(:,2)'*M*C(:,2));

  C(:,3) = V(:,2) - ((V(:,2)'*M*C(:,1)))*C(:,1) - ((V(:,2)'*M*C(:,2)))*C(:,2) ;
  C(:,3) = C(:,3)/sqrt(C(:,3)'*M*C(:,3));

  % kappa_dot = kappa_dot - ...
  %  (kappa_dot'*M*C(:,1))*C(:,1) - ...
  %  (kappa_dot'*M*C(:,2))*C(:,2) - ...
  %  (kappa_dot'*M*C(:,3))*C(:,3) ;
  kappa_dot = kappa_dot - sum(bsxfun(@times,sum(bsxfun(@times,kappa_dot,M*C)),C),2);

  % Take explicit euler step
  kappa = kappa + tao * kappa_dot;

  % arbitrarily let theta_0 = 0
  theta_0 = 0;
  % % angle between last segment and -x-axis
  %theta_0 = atan2( ...
  %  ev(end,1).*0 - ev(end,2).*-1, ...
  %  ev(end,1).*-1 + ev(end,2).*0 );
  theta_0 = atan2( ...
    -1.*ev(end,2) - 0.*ev(end,1), ...
    -1.*ev(end,1) + 0.*ev(end,2));
  theta = theta_0 + cumsum(kappa .* (0.5 * (l([end 1:end-1]) + l)));

  % recover tangents: T(i,:) is along V(i,:) to V(i+1,:)
  T = bsxfun(@times,l,[cos(theta) sin(theta)]);
  UT = [V(1,:); bsxfun(@plus,V(1,:),cumsum(-T))];

  % Recover positions
  % build laplacian
  L = sparse(1:n,[2:n 1],-1./l,n,n);
  % symmetric
  L = L+L';
  % diagonal diagonals
  L = L - diag(sum(L,2));
  % Bulid rhs
  b = bsxfun(@rdivide,T,l);
  b = b-b([end 1:end-1],:);
  % % solve
  %U = L\b;
  % Fix U(1,:) and solve
  U = V;
  U(2:end,:) = L(2:end,2:end) \ (-L(2:end,1) * U(1,:) + b(2:end,:));

  % % fit a rigid transformation
  %[R,T] = fit_rigid(V,U);
  %U = bsxfun(@plus,U*R,T);

end

author: alecjacobson

repo: gptoolbox


text:
function E = r8mat_expm3 ( A )

%*****************************************************************************80
%
%% R8MAT_EXPM3 approximates the matrix exponential using an eigenvalue approach.
%
%  Discussion:
%
%    exp(A) = V * D * V
%
%    where V is the matrix of eigenvectors of A, and D is the diagonal matrix
%    whose i-th diagonal entry is exp(lambda(i)), for lambda(i) an eigenvalue
%    of A.
%
%    This function is accurate for matrices which are symmetric, orthogonal,
%    or normal.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 March 2011
%
%  Author:
%
%    Cleve Moler, Charles Van Loan
%
%  Reference:
%
%    Cleve Moler, Charles VanLoan,
%    Nineteen Dubious Ways to Compute the Exponential of a Matrix,
%    Twenty-Five Years Later,
%    SIAM Review,
%    Volume 45, Number 1, March 2003, pages 3-49.
%
%  Parameters:
%
%    Input, real A(N,N), the matrix.
%
%    Output, real E(N,N), the estimate for exp(A).
%
  [ V, D ] = eig ( A );
  E = V * diag ( exp ( diag ( D ) ) ) / V;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [engine, loglik] = enter_evidence(engine, evidence, varargin)
% ENTER_EVIDENCE Add the specified evidence to the network (jtree)
% [engine, loglik] = enter_evidence(engine, evidence, ...)
%
% evidence{i} = [] if X(i) is hidden, and otherwise contains its observed value (scalar or column vector).
%
% The following optional arguments can be specified in the form of name/value pairs:
% [default value in brackets]
%
% soft    - a cell array of soft/virtual evidence;
%           soft{i} is a prob. distrib. over i's values, or [] [ cell(1,N) ]
%
% e.g., engine = enter_evidence(engine, ev, 'soft', soft_ev)

bnet = bnet_from_engine(engine);
ns = bnet.node_sizes(:);
N = length(bnet.dag);

engine.evidence = evidence; % store this for marginal_nodes with add_ev option
engine.maximize = 0;

% set default params
exclude = [];
soft_evidence = cell(1,N);

% parse optional params
args = varargin;
nargs = length(args);
for i=1:2:nargs
  switch args{i},
   case 'soft',    soft_evidence = args{i+1};
   otherwise,
    error(['invalid argument name ' args{i}]);
  end
end

onodes = find(~isemptycell(evidence));
hnodes = find(isemptycell(evidence));
pot_type = determine_pot_type(bnet, onodes);
 if strcmp(pot_type, 'cg')
  check_for_cd_arcs(onodes, bnet.cnodes, bnet.dag);
end

% Evaluate CPDs with evidence, and convert to potentials
pot = cell(1, N);
for n=1:N
  fam = family(bnet.dag, n);
  e = bnet.equiv_class(n);
  if isempty(bnet.CPD{e})
    error(['must define CPD ' num2str(e)])
  else
    pot{n} = convert_to_pot(bnet.CPD{e}, pot_type, fam(:), evidence);
  end
end
clqs = engine.clq_ass_to_node(1:N);

% soft evidence
soft_nodes = find(~isemptycell(soft_evidence));
S = length(soft_nodes);
if S > 0
  assert(pot_type == 'd');
  assert(mysubset(soft_nodes, bnet.dnodes));
end
for i=1:S
  n = soft_nodes(i);
  pot{end+1} = dpot(n, ns(n), soft_evidence{n});
end
clqs = [clqs engine.clq_ass_to_node(soft_nodes)];

[clpot, seppot] = init_pot(engine, clqs, pot, pot_type, onodes);
[clpot, seppot] = collect_evidence(engine, clpot, seppot);
[clpot, seppot] = distribute_evidence(engine, clpot, seppot);

C = length(clpot);
ll = zeros(1, C);
for i=1:C
  [clpot{i}, ll(i)] = normalize_pot(clpot{i});
end
loglik = ll(1); % we can extract the likelihood from any clique

engine.clpot = clpot;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [u,eqn,info] = fracLapP1P2(node,elem,pde,option)
%% FRACLAPP1P2 solves fractional Laplacian equation using P1-P2 element
%
% [u,eqn,info] = fracLapP1P1(node,elem,pde,option) solves the fractional
% Laplacian equation
%
%  (-\Delta^^s u = f in \Omega with u = 0 on \partial \Omega.

%% Options
if ~exist('option','var'), option = []; end
if ~isfield(option,'plotflag'), option.plotflag = 0; end

%% Parameters
global s;
alpha = 1-2*s;
if s == 0.5
    Gamma = 1;
else
    Gamma = 5/(2*s)+0.01;  % a different grading factor
end

%% Mesh
N = size(node,1); NT = size(elem,1);
% My = round(sqrt(NT));
L = 1+log10(NT)/3; % to control the truncated error, L should increase as logN
My = round(L*2*NT^(1/4)); % number of elements in y direction
% My = round(pde.L*sqrt(sqrt(NT/2)));
y = gradmap(0,L,Gamma,My);
Ny = length(y); % number of vertices in y-direction
NTy = Ny - 1;   % number of elements in y-direction
NTtotal = NT*NTy;
Nxdof = N;      % dof in x direction is number of vertices
Nydof = (Ny + NTy); % dof in y direction is number of vertices plus elements
Ndof = Nxdof*Nydof;

tic; % start the timing for the assembling
%% Stiffness matrix and mass matrix in the extended direction
% quantities in y direction
hy = diff(y);
a = zeros(NTy,5);
for i = 1:5
    a(:,i) = diff(y.^(alpha+i)/(alpha+i));
end
y1 = y(1:end-1);
y2 = y(2:end);
y12 = y1 + y2;
y1y2 = y1.*y2;
ym = y12/2;
% stiffness matrix in y direction
Ay = zeros(NTy,3,3);
Ay(:,1,1) = a(:,1)./hy.^2;
Ay(:,1,2) = -Ay(:,1,1);
Ay(:,2,1) = -Ay(:,1,1);
Ay(:,2,2) = Ay(:,1,1);
Ay(:,1,3) = 8*(a(:,2) - ym.*a(:,1))./hy.^3;
Ay(:,3,1) = Ay(:,1,3);
Ay(:,2,3) = -Ay(:,1,3);
Ay(:,3,2) = Ay(:,2,3);
Ay(:,3,3) = 64*(a(:,3) - 2*ym.*a(:,2) + ym.^2.*a(:,1))./hy.^4;
% mass matrix in y direction
My = zeros(NTy,3,3);
My(:,1,1) = (a(:,3) - 2*y1.*a(:,2) + y1.^2.*a(:,1))./hy.^2;
My(:,1,2) = (-a(:,3) + (y1+y2).*a(:,2) - y1.*y2.*a(:,1))./hy.^2;
My(:,2,1) = My(:,1,2);
My(:,2,2) = (a(:,3) - 2*y2.*a(:,2) + y2.^2.*a(:,1))./hy.^2;
My(:,1,3) = 4*(a(:,4) - (y12+y2).*a(:,3) + (y2.*y12 + y1y2).*a(:,2) ...
             - y2.*y1y2.*a(:,1))./hy.^3;
My(:,3,1) = My(:,1,3);
My(:,2,3) = 4*(-a(:,4) + (y12+y1).*a(:,3) - (y1.*y12 + y1y2).*a(:,2) ...
             + y1.*y1y2.*a(:,1))./hy.^3;
My(:,3,2) = My(:,2,3);
My(:,3,3) = 16*(a(:,5) - 2*y12.*a(:,4) + (y12.^2 + 2*y1y2).*a(:,3)...
               -2*y12.*y1.*y2.*a(:,2) + y1y2.^2.*a(:,1))./hy.^4;

%% Stiffness matrix and mass matrix in the original direction
[Dphi,area] = gradbasis(node,elem);
% Compute a piecewise constant diffusion coefficient
if ~isfield(pde,'d'), pde.d = []; end
if ~isfield(option,'dquadorder'), option.dquadorder = 1; end
if ~isempty(pde.d) && isnumeric(pde.d)
   K = pde.d;                                 % d is an array
end
if ~isempty(pde.d) && ~isnumeric(pde.d)       % d is a function
    [lambda,weight] = quadpts(option.dquadorder);
    nQuad = size(lambda,1);
    K = zeros(NT,1);
    for p = 1:nQuad
		pxy = lambda(p,1)*node(elem(:,1),:) ...
			+ lambda(p,2)*node(elem(:,2),:) ...
			+ lambda(p,3)*node(elem(:,3),:);
        K = K + weight(p)*pde.d(pxy);
   end
end
if ~isempty(pde.d) % build the coefficients into the scaled area
    areaK = K.*area;
else
    areaK = area;
end
At = zeros(NT,3,3);
Mt = zeros(NT,3,3);
for i = 1:3
    for j = 1:3
        At(:,i,j) = (Dphi(:,1,i).*Dphi(:,1,j) + Dphi(:,2,i).*Dphi(:,2,j)).*areaK;
        Mt(:,i,j) = areaK*((i==j)+1)/12;
    end
end

%% Assemble stiffness matrix
dofMap = repmat(1:Nxdof,Nydof,1) + repmat((0:Nydof-1)'*Nxdof,1,Nxdof);
ii = zeros(9*9*NTtotal,1);
jj = zeros(9*9*NTtotal,1);
sA = zeros(9*9*NTtotal,1);
index = 0;
for m = 1:3
    for n = 1:3
        for i = 1:3
            for j = 1:3
                if m < 3
                    ii(index+1:index+NTtotal) = dofMap(m:Ny+m-2,elem(:,i));
                else % m = 3: middle points of elements in y direction
                    ii(index+1:index+NTtotal) = dofMap(Ny+1:Nydof,elem(:,i));
                end
                if n < 3
                    jj(index+1:index+NTtotal) = dofMap(n:Ny+n-2,elem(:,j));
                else % n = 3: middle points of elements in y direction
                    jj(index+1:index+NTtotal) = dofMap(Ny+1:Nydof,elem(:,j));
                end
                sA(index+1:index+NTtotal) = kron(At(:,i,j),My(:,m,n)) + ...
                                            kron(Mt(:,i,j),Ay(:,m,n));
                index = index + NTtotal;
            end
        end
    end
end
A = sparse(ii,jj,sA,Ndof,Ndof);
clear ii jj sA At Mt Ay My

%% Assemble the right hand side
b = zeros(Ndof,1);
u = zeros(Ndof,1);

%% Set up boundary conditions
% find Dirichlet boundary nodes
bdNode = findboundary(elem);  % boundary vertices of mesh in x-direction
lateralbd = dofMap(:,bdNode);
topbd = dofMap(Ny,:);

% Modify the matrix to include the Dirichlet boundary condition
bdidx = zeros(Ndof,1);
bdidx(lateralbd) = 1;
bdidx(topbd) = 1;
Tbd = spdiags(bdidx,0,Ndof,Ndof);
T = spdiags(1-bdidx,0,Ndof,Ndof);
A = T*A*T + Tbd;

% Compute boundary integral over the original domain
if ~isfield(option,'fquadorder')
    option.fquadorder = 4;
end
ft = zeros(NT,3);
ds = 2^(1-2*s)*gamma(1-s)/gamma(s);
if isfield(pde,'f') && ~ischar(pde.f)
    [lambdaf,weightf] = quadpts(option.fquadorder);
    phif = lambdaf;                 % linear bases
    nQuadf = size(lambdaf,1);
    for p = 1:nQuadf
        % quadrature points in the x-y coordinate
        pxy = lambdaf(p,1)*node(elem(:,1),:) ...
            + lambdaf(p,2)*node(elem(:,2),:) ...
            + lambdaf(p,3)*node(elem(:,3),:);
        fp = pde.f(pxy);
        for i = 1:3
            ft(:,i) = ft(:,i) + weightf(p)*phif(p,i)*fp;
        end
    end
elseif isfield(pde,'f') && ischar(pde.f)
    if strcmp(pde.f,'intx') % f is a distribution
        % compute (intx, d_x phi)
        [lambdaf,weightf] = quadpts(option.fquadorder);
        nQuad = size(lambdaf,1);
        bt = zeros(NT,1);
        for p = 1:nQuad
            pxy = lambdaf(p,1)*node(elem(:,1),:) ...
                + lambdaf(p,2)*node(elem(:,2),:) ...
                + lambdaf(p,3)*node(elem(:,3),:);
            bt = bt + weightf(p)*pde.intx(pxy);
        end
%         [Dphi,area] = gradbasis(node,elem);
        for i = 1:3
            ft(:,i) = Dphi(:,1,i).*bt; % (intx, d_xphi)
        end
    end
end
ft = ds*ft.*repmat(area,1,3);
b = b + accumarray(elem(:),ft(:),[Ndof,1]);
% Neumann edfts are considered as open set. So the corner points should be
% set as Dirichlet boundary condition!
b(bdNode) = 0;
clear ft

%% Record assembling time
assembleTime = toc;
if ~isfield(option,'printlevel'), option.printlevel = 1; end
if option.printlevel >= 2
    fprintf('Time to assemble matrix equation %4.2g s\n',assembleTime);
end

%% Solve the system of linear equations
freeDof = find(bdidx==0);
if isempty(freeDof), return; end
% Set up solver type
if isempty(option) || ~isfield(option,'solver')    % no option.solver
    if Ndof <= 2e3  % Direct solver for small size systems
        option.solver = 'direct';
    else            % MGCG  solver for large size systems
        option.solver = 'mg';
    end
end
solver = option.solver;
% solve
switch solver
    case 'direct'
        tic;
        u(freeDof) = A(freeDof,freeDof)\b(freeDof);
        residual = norm(b - A*u);
        info = struct('solverTime',toc,'itStep',0,'err',residual,'flag',2,'stopErr',residual);
    case 'none'
        info = struct('solverTime',[],'itStep',0,'err',[],'flag',3,'stopErr',[]);
    case 'mg'
        option.x0 = u;
        option.solver = 'VCYCLE';
        option.freeDof = freeDof;
        [u,info] = mgfracLapP1P2(A,b,elem,option);
    case 'amg'
        option.solver = 'CG';
        [u(freeDof),info] = amg(A(freeDof,freeDof),b(freeDof),option);
end

%% Compute error using boundary integral
if isfield(pde,'exactu')
    err = zeros(NT,1);
    [lambda,weight] = quadpts(8);
    phi = lambda;                 % linear bases
    nQuadf = size(lambda,1);
    for p = 1:nQuadf
        % quadrature points in the x-y coordinate
        pxy = lambda(p,1)*node(elem(:,1),:) ...
            + lambda(p,2)*node(elem(:,2),:) ...
            + lambda(p,3)*node(elem(:,3),:);
        fp = pde.f(pxy);
        uhp = u(elem(:,1))*phi(p,1) + u(elem(:,2))*phi(p,2) + u(elem(:,3))*phi(p,3);
        up = pde.exactu(pxy);
        err = err + weight(p)*fp.*(up - uhp);
%         err = err + weight(p)*fp.*(up - 2*uhp);
    end
    err = ds*sum(err.*area);
%     err = ds*sum(err.*area) + u'*(A*u);
    err = sqrt(abs(err));
else
    err = 0;
end
info.errH1 = err;

%% Output information
eqn = struct('A',A,'b',b,'freeDof',freeDof);
info.assembleTime = assembleTime;

repo: ifem


text:
% STATCOND  - compare two or more data conditions statistically using
%               standard parametric or nonparametric permutation-based ANOVA
%               (1-way or 2-way) or t-test methods. Parametric testing uses
%               FCDF from the Matlab Statistical Toolbox.
% Usage:
%          >> [stats, df, pvals, surrog] = statcond( data, 'key','val'... );
%
% Inputs:
%   data       = one-or two-dimensional cell array of data matrices.
%                   For nonparametric, permutation-based testing, the
%                last dimension of the data arrays (which may be of up to
%                4 dimensions) is permuted across conditions, either in
%                a 'paired' fashion (not changing the, e.g., subject or
%                trial order in the last dimension) or in an umpaired
%                fashion (not respecting this order). If the number of
%                elements in the last dimension is not the same across
%                conditions, the 'paired' option is turned 'off'.  Note:
%                All other dimensions MUST be constant across conditions.
%                   For example, consider a (1,3) cell array of matrices
%                of size (100,20,x) each holding a (100,20) time/frequency
%                transform from each of x subjects. Only the last dimension
%                (here x, the number of subjects) may differ across the
%                three conditions.
%                   The test used depends on the size of the data array input.
%                When the data cell array has 2 columns and the data are
%                paired, a paired t-test is performed; when the data are
%                unpaired, an unpaired t-test is performed. If 'data'
%                has only one row (paired or unpaired) and more than 2
%                columns, a one-way ANOVA is performed. If the data cell
%                array contains several rows and columns, and the data is
%                paired, a two-way repeated measure ANOVA is performed.
%                NOTE THAT IF THE DATA is unpaired, EEGLAB will use a
%                balanced 1 or 2 way ANOVA and parametric results might not
%                be meaningful (bootstrap and permstatcondutation should be fine).
%
% Optional inputs:
%   'paired'   = ['on'|'off'] pair the data array {default: 'on' unless
%                the last dimension of data array is of different lengths}.
%   'method'   = ['perm'|'bootstrap'|'param'] method for computing the p-values:
%                 'param' or 'parametric' = parametric testing (standard ANOVA
%                                           or t-test);
%                 'perm' or 'permutation' = non-parametric testing using
%                                           surrogate data
%                 'bootstrap' = non-parametric bootstrap
%                  made by permuting the input data {default: 'param'}
%   'naccu'    = [integer] Number of surrogate data copies to use in 'perm'
%                 or 'bootstrap' method estimation (see above) {default: 200}.
%   'verbose'  = ['on'|'off'] print info on the command line {default: 'on'}.
%   'variance' = ['homegenous'|'inhomogenous'] this option is exclusively
%                for parametric statistics using unpaired t-test. It allows
%                to compute a more accurate value for the degree of freedom
%                using the formula for inhomogenous variance (see
%                ttest2_cell function). Default is 'inhomegenous'.
%   'surrog'   = surrogate data array (see output).
%   'stats'    = F- or T-value array (see output).
%   'tail'     = ['one'|'two'] run one-tailed (F-test) or two tailed
%                (T-test). This option is only relevant when using the
%                'surrog' input. Otherwise it is ignored.
%   'forceanova' = ['on'|'off'] force the use of ANOVA calculation even
%                for 2x1 designs. Default is 'off'.
%   'alpha'    = [float] p-value threshold value. Allow returning
%                confidence intervals and mask (requires structoutput below).
%   'structoutput' = ['on'|'off'] return an output structure instead of
%                the regular output. Allow to output mask and confidence
%                intervals.
%   'cluster'  = ['on'|'off'] cluster correction for multiple comparison.
%                Only functional when alpha is set (if alpha is NaN, it
%                sets it to 0.05) and for 1-way Anova or t-test.
%
% Legacy parameters:
%   'threshold' - now 'alpha'
%   'mode'      - now 'method'
%
% Outputs:
%   stats      = F- or T-value array of the same size as input data without
%                the last dimension. A T value is returned only when the data
%                includes exactly two conditions.
%   df         = degrees of freedom, a (2,1) vector, when F-values are returned
%   pvals      = array of p-values. Same size as input data without the last
%                data dimension. All returned p-values are two-tailed.
%   surrog     = surrogate statistic values (same size as stats output with the last
%                dimension filled with a number ('naccu') of surrogate data sets.
%
% Important note: When a two-way ANOVA is performed, outputs are cell arrays
%                 with three elements: output(1) = row effects;
%                 output(2) = column effects; output(3) = interactions
%                 between rows and columns.
%
% Examples:
%      >> a = { rand(1,10) rand(1,10)+0.5 }; % pseudo 'paired' data vectors
%         [t df pvals] = statcond(a);        % perform paired t-test
%           pvals =
%              5.2807e-04 % standard t-test probability value
%         % Note: for different RAND outputs, results will differ.
%
%         [t df pvals surog] = statcond(a, 'method', 'perm', 'naccu', 2000);
%           pvals =
%              0.0065 % nonparametric t-test using 2000 permuted data sets
%
%         a = { rand(2,11) rand(2,10) rand(2,12)+0.5 }; % pseudo 'unpaired'
%         [F df pvals] = statcond(a); % perform an unpaired ANOVA
%           pvals =
%              0.00025 % p-values for difference between columns
%              0.00002 % for each data row
%
%         a = { rand(3,4,10) rand(3,4,10) rand(3,4,10); ...
%               rand(3,4,10) rand(3,4,10) rand(3,4,10)+0.5 };
%         % pseudo (2,3)-condition data array, each entry containing
%         %                                    ten (3,4) data matrices
%         [F df pvals] = statcond(a);  % perform a paired 2-way ANOVA
%         % Output:
%           pvals{1} % a (3,4) matrix of p-values; effects across rows
%           pvals{2} % a (3,4) matrix of p-values; effects across columns
%           pvals{3} % a (3,4) matrix of p-values; interaction effects
%                                      % across rows and columns
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2005-
%         With thanks to Robert Oostenveld for fruitful discussions
%         and advices on this function.
%
% See also: ANOVA1_CELL, ANOVA2_CELL, ANOVA2RM_CELL, FCDF

% perform a paired t-test
% -----------------------
% a = { rand(2,10) rand(2,10) };
% [t df pval] = statcond(a); pval
% [h p t stat] = ttest( a{1}(1,:), a{2}(1,:)); p
% [h p t stat] = ttest( a{1}(2,:), a{2}(2,:)); p
%
% compare significance levels
% --------------------------
% a = { rand(1,10) rand(1,10) };
% [F df pval] = statcond(a, 'method', 'perm', 'naccu', 200); pval
% [h p t stat] = ttest( a{1}(1,:), a{2}(1,:)); p

% Copyright (C) Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. 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 HOLDER 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.

function [ ori_vals, df, pvals, surrogval ] = statcond( data, varargin );

    if nargin < 1
        help statcond;
        return;
    end
    try, warning('off', 'MATLAB:divideByZero'); catch, end;

    if exist('finputcheck')
        g = finputcheck( varargin, { 'naccu'      'integer'   [1 Inf]             200;
                                     'method'     'string'    { 'param','parametric','perm','permutation','bootstrap' }  'param';
                                     'mode'       'string'    { }                 '';
                                     'paired'     'string'    { 'on','off','auto' }      'auto';
                                     'surrog'     { 'real','cell' }      []       [];
                                     'stats'      { 'real','cell' }      []       [];
                                     'structoutput' 'string'  { 'on','off' }      'off';
                                     'forceanova'   'string'  { 'on','off' }      'off';
                                     'arraycomp'  'string'    { 'on','off' }      'on';
                                     'cluster'    'string'    { 'on','off' }      'off';
                                     'alpha'      'real'      []                  NaN;
                                     'tail'       'string'    { 'one','both','upper','lower'}    'both';
                                     'variance'   'string'    { 'homogenous','inhomogenous' }    'inhomogenous';
                                     'returnresamplingarray' 'string'    { 'on','off' }      'off';
                                     'verbose'    'string'    { 'on','off' }      'on' }, 'statcond');
        if ischar(g), error(g); end
    else
        g = struct(varargin{:});
        if ~isfield(g, 'naccu'),     g.naccu = 200; end
        if ~isfield(g, 'method'),    g.method  = 'param'; end
        if ~isfield(g, 'paired'),    g.paired = 'auto'; end
        if ~isfield(g, 'surrog'),    g.surrog = []; end
        if ~isfield(g, 'orivals'),   g.orivals = []; end
        if ~isfield(g, 'arraycomp'), g.arraycomp = 'on'; end
        if ~isfield(g, 'verbose'),   g.verbose = 'on'; end
        if ~isfield(g, 'tail'),      g.tail = 'both'; end
        if ~isfield(g, 'variance'),  g.variance = 'homogenous'; end
        if ~isfield(g, 'structoutput'), g.structoutput = 'on'; end
        if ~isfield(g, 'returnresamplingarray'),   g.returnresamplingarray = 'off'; end
    end
    if ~isempty(g.mode), g.method = g.mode; end
    if size(data,2) == 1, data  = transpose(data); end % cell array transpose

    % other settings
    % --------------
    if strcmpi(g.method, 'parametric'), g.method = 'param'; end
    if strcmpi(g.method, 'permutation'), g.method = 'perm'; end
    if strcmpi(g.verbose, 'on'), verb = 1; else verb = 0; end
    if strcmp(g.method, 'param' ) && exist('fcdf') ~= 2
      myfprintf('on',['statcond(): parametric testing requires fcdf() \n' ...
               '            from the Matlab StatsticaL Toolbox.\n' ...
               '            Running nonparametric permutation tests\n.']);
      g.method = 'perm';
    end
    g.naccu = round(g.naccu);

    % pairing
    % -------
    paired{1} = 'on';
    if size(data,1) == 1 && length(unique(cellfun('size', data, ndims(data{1}) ))) > 1
        paired{1} = 'off';
    else
        paired{1} = 'on';
        for iCol = 1:size(data,2)
            if length(unique(cellfun('size', data(:,iCol), ndims(data{1}) ))) > 1
                paired{1} = 'off';
            end
        end
        paired{2} = 'on';
        for iRow = 1:size(data,1)
            if length(unique(cellfun('size', data(iRow,:), ndims(data{1}) ))) > 1
                paired{2} = 'off';
            end
        end
    end
    if length(paired) > 1
        if (strcmpi(paired{1}, 'off') && strcmpi(paired{2}, 'on')) || ...
                (strcmpi(paired{1}, 'on') && strcmpi(paired{2}, 'off'))
            myfprintf(g.verbose, 'Possible mixed paired and unpaired independent variables, using balanced Anova (which assumes all unpaired)\n');
            paired{1} = 'off';
        end
    end
    if strcmpi(g.paired, 'auto')
        g.paired = paired{1};
    else
        if strcmpi(g.paired, 'on') && strcmpi(paired{1}, 'off')
            myfprintf(g.verbose, 'You set to use paired statistics but the number of cases differs\n');
            g.paired = 'off';
        end
    end

    % reshape matrices
    % ----------------
    nd = size(data{1});
    nd = nd(1:end-1);
    for index = 1:prod(size(data))
        data{index} = reshape(data{index}, [prod(nd) size(data{index},myndims(data{index}))]);
    end

    if ~strcmpi(g.method, 'param') && isempty(g.surrog)
         tmpsize   = size(data{1});
         surrogval = zeros([ tmpsize(1:end-1) g.naccu ], 'single');
    else surrogval = [];
    end

    % check for NaNs or Inf
    % ---------------------
    for iDat = 1:length(data(:))
        if any(isnan(reshape(data{iDat}, prod(size(data{iDat})),1))) || ...
                any(isinf(reshape(data{iDat}, prod(size(data{iDat})),1)))
            error('Statcond: One of the input array contains NaNs or Infinite values');
        end
    end

    % bootstrap flag
    % --------------
    if strcmpi(g.method, 'bootstrap'), bootflag = 1;
    else                               bootflag = 0;
    end

    if isempty(g.surrog)

        % return resampling array
        % -----------------------
        if strcmpi(g.returnresamplingarray, 'on')
            if strcmpi(g.arraycomp, 'on')
                ori_vals = surrogdistrib( data, 'method', g.method, 'pairing', g.paired, 'naccu', g.naccu);
            else
                ori_vals = surrogdistrib( data, 'method', g.method, 'pairing', g.paired);
            end
            return;
        end

        % text output
        % -----------
        myfprintf(verb,'%d x %d, ', size(data,1), size(data,2));
        if strcmpi(g.paired, 'on')
             myfprintf(verb,'paired data, ');
        else myfprintf(verb,'unpaired data, ');
        end
        if size(data,1) == 1 && size(data,2) == 2
             myfprintf(verb,'computing T values\n');
        else myfprintf(verb,'computing F values\n');
        end
        if size(data,1) > 1
            if strcmpi(g.paired, 'on')
                 myfprintf(verb,'Using 2-way repeated measure ANOVA\n');
            else myfprintf(verb,'Using balanced 2-way ANOVA (not suitable for parametric testing, only bootstrap)\n');
            end
        elseif size(data,2) > 2
            if strcmpi(g.paired, 'on')
                 myfprintf(verb,'Using 1-way repeated measure ANOVA\n');
            else myfprintf(verb,'Using balanced 1-way ANOVA (equivalent to Matlab anova1)\n');
            end
        else
            if strcmpi(g.paired, 'on')
                 myfprintf(verb,'Using paired t-test\n');
            else myfprintf(verb,'Using unpaired t-test\n');
            end
        end
        if ~strcmpi(g.method, 'param')
            if bootflag, myfprintf(verb,'Bootstraps (of %d):', g.naccu);
            else         myfprintf(verb,'Permutations (of %d):', g.naccu);
            end
        end
    end

    tail = g.tail;
    if isempty(g.surrog)
        if size(data,1) == 1 % only one row

            if size(data,2) == 2 && strcmpi(g.forceanova, 'off')

                % paired t-test (very fast)
                % -------------
                [ori_vals, df] = ttest_cell_select(data, g.paired, g.variance);

                if strcmpi(g.method, 'param')

                    % Check if exist tcd.m file from the Statistics Toolbox (Bug 1352 )
                    if exist('tcdf','file') == 2  && license('test', 'Statistics_Toolbox')
                        pvals = 2*tcdf(-abs(ori_vals), df);
                    else
                        pvals = 2*mytcdf(-abs(ori_vals), df);
                    end

                    pvals = reshape(pvals, size(pvals));
                else
                    if strcmpi(g.arraycomp, 'on')
                        try
                            myfprintf(verb,'...');
                            res = surrogdistrib( data, 'method', g.method, 'pairing', g.paired, 'naccu', g.naccu);
                            surrogval = ttest_cell_select( res, g.paired, g.variance);
                        catch
                           lasterr
                           myfprintf(verb,'\nFast computation failed because of memory limitation, reverting to standard computation');
                           g.arraycomp = 'off';
                        end
                    end
                    if strcmpi(g.arraycomp, 'off')
                        [res, precomp] = surrogdistrib( data, 'method', g.method, 'pairing', g.paired);
                        for index = 1:g.naccu
                            res = surrogdistrib( {}, 'precomp', precomp);
                            if mod(index, 10) == 0, myfprintf(verb,'%d ', index); end
                            if mod(index, 100) == 0, myfprintf(verb,'\n'); end
                            if myndims(res{1}) == 1
                                 surrogval(index)     = ttest_cell_select(res, g.paired, g.variance);
                            else surrogval(:,index)   = ttest_cell_select(res, g.paired, g.variance);
                            end
                        end
                    end
                end
            else
                % one-way ANOVA (paired) this is equivalent to unpaired t-test
                % -------------
                tail = 'one';
                [ori_vals, df] = anova1_cell_select( data, g.paired );
                if strcmpi(g.method, 'param')
                    pvals = 1-fcdf(ori_vals, df(1), df(2));
                else
                    if strcmpi(g.arraycomp, 'on')
                        try
                            myfprintf(verb,'...');
                            res = surrogdistrib( data, 'method', g.method, 'pairing', g.paired, 'naccu', g.naccu);
                            surrogval = anova1_cell_select( res, g.paired );
                        catch,
                            myfprintf(verb,'\nFast computation failed because of memory limitation, reverting to standard computing');
                            g.arraycomp = 'off';
                        end
                    end
                    if strcmpi(g.arraycomp, 'off')
                        [res, precomp] = surrogdistrib( data, 'method', g.method, 'pairing', g.paired);
                        for index = 1:g.naccu
                            if mod(index, 10) == 0, myfprintf(verb,'%d ', index); end
                            if mod(index, 100) == 0, myfprintf(verb,'\n'); end

                            res = surrogdistrib( {}, 'precomp', precomp);
                            if myndims(data{1}) == 1
                            	 surrogval(index)     = anova1_cell_select( res, g.paired );
                            else surrogval(:,index)   = anova1_cell_select( res, g.paired );
                            end
                        end
                    end
                end
            end
        else
            % two-way ANOVA (paired or unpaired)
            % ----------------------------------
            tail = 'one';
            [ ori_vals{1}, ori_vals{2}, ori_vals{3}, df{1}, df{2}, df{3} ] = anova2_cell_select( data, g.paired );
            if strcmpi(g.method, 'param')
                pvals{1} = 1-fcdf(ori_vals{1}, df{1}(1), df{1}(2));
                pvals{2} = 1-fcdf(ori_vals{2}, df{2}(1), df{2}(2));
                pvals{3} = 1-fcdf(ori_vals{3}, df{3}(1), df{3}(2));
            else
                surrogval = { surrogval surrogval surrogval };
                dataori   = data;
                if strcmpi(g.arraycomp, 'on')
                    try
                        myfprintf(verb,'...');
                        res = surrogdistrib( data, 'method', g.method, 'pairing', g.paired, 'naccu', g.naccu);
                        [ surrogval{1}, surrogval{2}, surrogval{3} ] = anova2_cell_select( res, g.paired );
                    catch
                        myfprintf(verb,'\nFast computation failed because of memory limitation, reverting to standard computing');
                        g.arraycomp = 'off';
                    end
                end
                if strcmpi(g.arraycomp, 'off')
                    [res, precomp] = surrogdistrib( data, 'method', g.method, 'pairing', g.paired);
                    for index = 1:g.naccu
                        if mod(index, 10) == 0, myfprintf(verb,'%d ', index); end
                        if mod(index, 100) == 0, myfprintf(verb,'\n'); end

                        res = surrogdistrib( {}, 'precomp', precomp);
                        if myndims(data{1}) == 1
                         	 [ surrogval{1}(index),     surrogval{2}(index),     surrogval{3}(index)     ] = anova2_cell_select( res, g.paired );
                        else [ surrogval{1}(:,index),   surrogval{2}(:,index),   surrogval{3}(:,index)   ] = anova2_cell_select( res, g.paired );
                        end
                    end
                end
            end
        end
        myfprintf(verb,'\n');
    else
        surrogval = g.surrog;
        ori_vals  = g.stats;
        df        = [];
    end

    % compute p-values
    % ----------------
    if ~strcmpi(g.method, 'param')
        if iscell( surrogval )
            pvals{1} = stat_surrogate_pvals(surrogval{1}, ori_vals{1}, tail);
            pvals{2} = stat_surrogate_pvals(surrogval{2}, ori_vals{2}, tail);
            pvals{3} = stat_surrogate_pvals(surrogval{3}, ori_vals{3}, tail);
        else
            if strcmpi(g.cluster, 'on')
                if isnan(g.alpha)
                    g.alpha = 0.05;
                    disp('Alpha value set to 0.05 automatically');
                end
                if size(surrogval,3) > 1 || size(surrogval,4) > 1 || isnan(g.alpha)
                    error('Cluster method not implemented for alpha NaN or for more than 2 dims');
                else
                    tmpPvals = 2*tcdf(-abs(surrogval), size(surrogval,1)) < g.alpha;
                    largestCluster = zeros(1, size(surrogval,2));
                    signifPos = zeros(1, size(surrogval,1));
                    for iSurog = 1:size(surrogval,2)
                        signifPos = tmpPvals(:,iSurog);
                        while any(signifPos)
                            tmpInd = find(signifPos);
                            currentInd = tmpInd(1);
                            currentCount = 0;
                            while currentInd+currentCount <= length(signifPos) && signifPos(currentInd+currentCount) == 1
                                signifPos(currentInd+currentCount) = 0;
                                currentCount = currentCount+1;
                            end
                            if currentCount > largestCluster(iSurog), largestCluster(iSurog) = currentCount; end
                        end
                    end
                    largestCluster = sort(largestCluster);
                    thresholdVal = largestCluster(round(length(largestCluster)*(1-g.alpha)));

                    % assign p-val status
                    signifPos = 2*tcdf(-abs(ori_vals), length(ori_vals)) < g.alpha;
                    pvals     = ones(size(signifPos));
                    while any(signifPos)
                        tmpInd = find(signifPos);
                        currentInd = tmpInd(1);
                        currentCount = 0;
                        while currentInd+currentCount <= length(signifPos) && signifPos(currentInd+currentCount) == 1
                            signifPos(currentInd+currentCount) = 0;
                            currentCount = currentCount+1;
                        end
                        if currentCount >= thresholdVal, pvals(currentInd:currentInd+currentCount-1) = 0; end
                    end

                end
            else
                pvals = stat_surrogate_pvals(surrogval, ori_vals, tail);
            end
        end
        try, warning('on', 'MATLAB:divideByZero'); catch, end
    end

    [ ori_vals, pvals ] = reshape_results( nd, ori_vals, pvals);
    [ surrogval ]       = reshape_results( [nd g.naccu], surrogval);

    % confidence intervals
    % --------------------
    if ~isnan(g.alpha)
        outputstruct.ci = stat_surrogate_ci(surrogval, g.alpha, tail);
        if strcmpi(g.structoutput, 'off')
            disp('Warning: returning confidence interval requires an output structure');
        end
        if iscell(pvals)
            for ind = 1:length(pvals)
                outputstruct.mask{ind} = pvals{ind} < g.alpha;
            end
        else
            outputstruct.mask = pvals < g.alpha;
        end
    end

    % create a structure for outputting values
    % ---------------------------------------
    if strcmpi(g.structoutput, 'on')
        outputstruct.method = g.method;
        outputstruct.pval   = pvals;
        outputstruct.df     = df;
        outputstruct.surrog = surrogval;
        if length(data(:)) == 2
             outputstruct.t = ori_vals;
        else outputstruct.f = ori_vals;
        end
        outputstruct.stat   = ori_vals;
        ori_vals = outputstruct;
    end

% compute ANOVA 2-way
% -------------------
function [f1, f2, f3, df1, df2, df3] = anova2_cell_select( res, paired)
    if strcmpi(paired,'on')
        [f1, f2, f3, df1, df2, df3] = anova2rm_cell( res );
    else
        [f1, f2, f3, df1, df2, df3] = anova2_cell( res );
    end

% compute ANOVA 1-way
% -------------------
function [f, df] = anova1_cell_select( res, paired)
    if strcmpi(paired,'on')
        [f, df] = anova1rm_cell( res );
    else
        [f, df] = anova1_cell( res );
    end

% compute t-test
% -------------------
function [t, df] = ttest_cell_select( res, paired, homogenous)
    if strcmpi(paired,'on')
        [t, df] = ttest_cell( res{1}, res{2});
    else
        [t, df] = ttest2_cell( res{1}, res{2}, homogenous);
    end

% function to compute the number of dimensions
% --------------------------------------------
function val = myndims(a)
    if ndims(a) > 2
        val = ndims(a);
    else
        if size(a,1) == 1
            val = 2;
        elseif size(a,2) == 1
            val = 1;
        else
            val = 2;
        end
    end

% function for verbose messages
% -----------------------------
function myfprintf(verb, varargin)
    if verb
        fprintf(varargin{:});
    end

% function to replace tcdf
% ------------------------
function p = mytcdf(x,v)

if length(v) == 1
    v = repmat(v, size(x));
end

x2 = x.^2;
inds1 = (v < x2);
inds2 = (v >= x2);
if any(inds1(:)), p(inds1) = betainc(v(inds1) ./ (v(inds1) + x2(inds1)), v(inds1)/2, 0.5, 'lower') / 2; end
if any(inds2(:)), p(inds2) = betainc(x2(inds2) ./ (v(inds2) + x2(inds2)), 0.5, v(inds2)/2, 'upper') / 2; end
inds = (x > 0);
if any(inds)
    p(inds) = 1 - p(inds);
end

inds = (v > 1e7);
if any(inds(:)), p(inds) = normcum(x(inds)); end

p(x == 0) = 0.5;
if isempty(p)
    p = ones(size(x));
else
    p = reshape(p, size(x));
end
function [p] = normcum(z)
p = 0.5 * erfc(-z ./ sqrt(2));

% reshape results
% ---------------
function varargout = reshape_results(nd, varargin)
    if length(varargin) > 1
        for index = 1:length(varargin)
            varargout{index} = reshape_results(nd, varargin{index});
        end
    elseif iscell(varargin{1})
        for index = 1:length(varargin{1})
            varargout{1}{index} = reshape_results(nd, varargin{1}{index});
        end
    else
        if ~isempty(varargin{1})
            if length(nd) == 1, nd = [ nd 1 ]; end
            varargout{1} = reshape(varargin{1}, nd);
        else varargout{1} = [];
        end
    end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function r8btree_interp_test ( )

%*****************************************************************************80
%
%% R8BTREE_INTERP_TEST demonstrates R8BTREE_INTERP.
%
%  Discussion:
%
%    We sample the function P77_FUN over the interval [0,2].
%
%    Then we compare interpolated to exact values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 November 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8BTREE_INTERP_TEST\n' );
  fprintf ( 1, '  Build a BTREE from sample values of a function.\n' );
  fprintf ( 1, '  Then compare exact to interpolated values.\n' );
%
%  Initialize the BTREE.
%
  node_num = 0;
  tree = [];
  data_num = 2;
  tree_data = [];
%
%  Generate 33 nodes in [0,2];
%
  for log_bottom = 1 : 5

    bottom = 2^log_bottom;

    for top = 1 : 2 : bottom - 1

      x = 2 * top / bottom;
      fx = p77_fun ( x );
      node_data = [ x, fx ];

      [ node_num, tree, tree_data ] = r8btree_node_add ( node_num, tree, ...
        data_num, tree_data, node_data );

    end

    if ( log_bottom == 1 )

      x = 0.0;
      fx = p77_fun ( x );
      node_data = [ x, fx ];

      [ node_num, tree, tree_data ] = r8btree_node_add ( node_num, tree, ...
        data_num, tree_data, node_data );

      x = 2.0;
      fx = p77_fun ( x );
      node_data = [ x, fx ];

      [ node_num, tree, tree_data ] = r8btree_node_add ( node_num, tree, ...
        data_num, tree_data, node_data );

    end

  end
%
%  Print the BTREE.
%
  r8btree_print_ordered ( node_num, tree, data_num, tree_data );
%
%  Do interpolation.
%
  int_num = 21;
  x = linspace ( 0.0, 2.0, int_num );
  fx = r8btree_interp ( node_num, tree, data_num, tree_data, int_num, x );

  fprintf ( 1, '\n' );
  fprintf ( 1, '      X          Exact F    Interp F\n' );
  fprintf ( 1, '\n' );
  for i = 1 : int_num
    fprintf ( 1, '  %10f  %10f  %10f\n', x(i), p77_fun ( x(i) ), fx(i) );
  end

  return
end


text:
function [bw] = kernelbw(xi)
%% Silverman bandwidth

N      = length(xi);
prop   = 1.0;
sig    = std(xi);
iqrSig = 0.7413 * iqr(xi);

if (max(iqrSig) == 0)
   iqrSig = sig;
end

bw = prop * min(sig, iqrSig) * N^(-1 / (4 + 1));

end

author: Sable

repo: mcbench-benchmarks


text:
function C = tolUnion(A, B, tol)
%TOLUNION   UNION with a tolerance for checking floating-point equality.
%   C = TOLUNION(A, B, TOL) for real vectors A and B is the same as UNION(A, B)
%   except that two values are deemed "equal" if they are within an absolute
%   difference TOL of one another.  If this happens, the mean of the two values
%   is placed in C.
%
%   C = TOLUNION(A, B) uses a default tolerance of 100*EPS*MAX(NORM(A, Inf),
%   NORM(B, Inf)).

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

if ( nargin < 3 )
    tol = 100*eps*max(norm(A, Inf), norm(B, Inf));
end

% TODO:  This might not do the right thing if there are multiple values that
% are all close to one another, say, 2, 2 + 2*eps, 2 + 4*eps, and 2 + 6*eps,
% and tol = 3*eps.  This code will merge all of those into one value:  2 + eps
% (the average of 2 and 2 + 2*eps), and one can argue that behavior is wrong.
% We can resolve this if this ever becomes an issue in practice.
C = union(A, B);
ind = find(diff(C) < tol);
C(ind) = (C(ind) + C(ind+1))/2;
C(ind+1) = [];

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun