text:
function gmst=utc2gmst(t,ut1_utc)

ep2000=[2000 1 1 12 0 0];

tut=timeadd(t,ut1_utc);
[ut,tut0]=time2sec(tut);
t1=timediff(tut0,epoch2time(ep2000))/86400.0/36525.0;
t2=t1^2; t3=t1^2;
gmst0=24110.54841+8640184.812866*t1+0.093104*t2-6.2E-6*t3;
gmst=gmst0+1.002737909350795*ut;

if gmst<0
    gmst=-mod(abs(gmst),86400.0)*pi/43200.0;
else
    gmst=mod(abs(gmst),86400.0)*pi/43200.0;
end

return

repo: GINav


text:
classdef testTopOptLevelSetPrinting < testTopOptPrinting

    properties (Access = protected)
        fileOutputName = 'testTopOptLevelSetPrinting';
        printMode = 'DesignVariable';
    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function outface=outersurf(node,face)
%
% outface=outersurf(node,face)
%
% extract the out-most shell of a complex surface mesh
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%    node:  node coordinates
%    face:  surface triangle list
%
% output:
%    outface: the out-most shell of the surface mesh
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

face=face(:,1:3);

ed=surfedge(face);
if(~isempty(ed))
   error('open surface is detected, you have to close it first, consider meshcheckrepair() with meshfix option');
end

[no,el]=fillsurf(node,face);

outface=volface(el);

[no,outface]=removeisolatednode(no,outface);
maxfacenode=max(outface(:));

[I,J]=ismember(round(no(1:maxfacenode,:)*1e10),round(node*1e10),'rows');

% if(sum(I)~=maxfacenode)
%     error('mesh tessellation failed');
% end
outface=J(outface);
[ii,jj]=find(outface==0);
outface(ii,:)=[];

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function kern = sqexpKernParamInit(kern)

% SQEXPKERNPARAMINIT SQEXP kernel parameter initialisation.
% This kernel is a 'pre-packaged' compound kernel of the form
% {'rbf', 'lin', 'bias', 'white'}. Using this kernel removes
% the overhead of mutliple calls through the 'cmpnd' kernel.
%
% SEEALSO sqexpKernParamInit
%
% FORMAT
% DESC initialises the pre-built compound squared exponential
%  kernel structure with some default parameters.
% ARG kern : the kernel structure which requires initialisation.
% RETURN kern : the kernel structure with the default parameters placed in.
%
% SEEALSO : kernCreate, kernParamInit
%
% COPYRIGHT : Neil D. Lawrence, 2004

% KERN

kern.inverseWidth = 1;
kern.rbfVariance = 1;
kern.whiteVariance = 1;
kern.biasVariance = 1;
kern.nParams = 4;

kern.transforms(1).index = [1 2 3 4];
kern.transforms(1).type = optimiDefaultConstraint('positive');

kern.isStationary = false;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function dicomrt_showabout(item)
% dicomrt_showabout
%
% Display information, copyright and discalimer.
%
% Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org)

if strcmp(item,'showabout_dicomrt_explore')
    what2say = {'This is dicomrt_explore.','','',...
            'Dicomrt_explore is a wrapper for a custom version of sliceomatic.',...
            'This function is part of the DICOM-RT Toolbox.',...
            'Copyright (C) 2002 Emiliano Spezi (emiliano.spezi@physics.org).',...
            'The DICOM-RT Toolbox is distributed under the General Public License (www.gnu.org/copyleft/gpl.html).'};
    helpdlg(what2say,'About dicomrt_explore')
elseif strcmp(item,'showabout_sliceomatic')
    what2say = {'The DICOM-RT Toolbox incorporates a custom version of sliceomatic v. 2.0.','','',...
            'Sliceomatic is a slice and isosurface volume exploration GUI.',...
            'Sliceomatic was written by Eric Ludlam <eludlam@mathworks.com>.',...
            'Copyright 2000, 2001, 2002 The MathWorks Inc.',...
            'Sliceomatic can be downloaded from Matlab-Central (www.mathworks.com/matlabcentral).','',...
            'Original code modified by Emiliano Spezi <emiliano.spezi@physics.org> ',...
            'to include the following capabilities:',...
            '1) axes limits control',...
            '2) support for all DICOM patient orientations',...
            '3) contour level selection',...
            '4) superposition of DICOM RTSTRUCT objects'};
    helpdlg(what2say,'About sliceomatic')
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [nu, U, llh, Ezz, Ezy] = kalmanSmoother(model, X)
% Kalman smoother (forward-backward algorithm for linear dynamic system)
% NOTE: This is the exact implementation of the Kalman smoother algorithm in PRML.
% However, this algorithm is not practical. It is numerical unstable.
% Input:
%   X: d x n data matrix
%   model: model structure
% Output:
%   nu: q x n matrix of latent mean mu_t=E[z_t] w.r.t p(z_t|x_{1:T})
%   U: q x q x n latent covariance U_t=cov[z_t] w.r.t p(z_t|x_{1:T})
%   Ezz: q x q matrix E[z_tz_t^T]
%   Ezy: q x q matrix E[z_tz_{t-1}^T]
%   llh: loglikelihood
% Written by Mo Chen (sth4nth@gmail.com).
A = model.A; % transition matrix
G = model.G; % transition covariance
C = model.C; % emission matrix
S = model.S;  % emision covariance
mu0 = model.mu0; % prior mean
P0 = model.P0;  % prior covairance

n = size(X,2);
q = size(mu0,1);
mu = zeros(q,n);
V = zeros(q,q,n);
P = zeros(q,q,n); % C_{t+1|t}
Amu = zeros(q,n); % u_{t+1|t}
llh = zeros(1,n);

% forward
PC = P0*C';
R = C*PC+S;
K = PC/R;
mu(:,1) = mu0+K*(X(:,1)-C*mu0);
V(:,:,1) = (eye(q)-K*C)*P0;
P(:,:,1) = P0;  % useless, just make a point
Amu(:,1) = mu0; % useless, just make a point
llh(1) = logGauss(X(:,1),C*mu0,R);
for i = 2:n
    [mu(:,i), V(:,:,i), Amu(:,i), P(:,:,i), llh(i)] = ...
        forwardUpdate(X(:,i), mu(:,i-1), V(:,:,i-1), A, G, C, S);
end
llh = sum(llh);
% backward
nu = zeros(q,n);
U = zeros(q,q,n);
Ezz = zeros(q,q,n);
Ezy = zeros(q,q,n-1);

nu(:,n) = mu(:,n);
U(:,:,n) = V(:,:,n);
Ezz(:,:,n) = U(:,:,n)+nu(:,n)*nu(:,n)';
for i = n-1:-1:1
    [nu(:,i), U(:,:,i), Ezz(:,:,i), Ezy(:,:,i)] = ...
        backwardUpdate(nu(:,i+1), U(:,:,i+1), mu(:,i), V(:,:,i), Amu(:,i+1), P(:,:,i+1), A);
end

function [mu1, V1, Amu, P, llh] = forwardUpdate(x, mu0, V0, A, G, C, S)
k = numel(mu0);
P = A*V0*A'+G;                                               % 13.88
PC = P*C';
R = C*PC+S;
K = PC/R;                                                    % 13.92
Amu = A*mu0;
CAmu = C*Amu;
mu1 = Amu+K*(x-CAmu);                                        % 13.89
V1 = (eye(k)-K*C)*P;                                         % 13.90
llh = logGauss(x,CAmu,R);                                    % 13.91

function [nu0, U0, E00, E10] = backwardUpdate(nu1, U1, mu, V, Amu, P, A)
J = V*A'/P;                                                  % 13.102
nu0 = mu+J*(nu1-Amu);                                        % 13.100
U0 = V+J*(U1-P)*J';                                          % 13.101
E00 = U0+nu0*nu0';                                           % 13.107
E10 = U1*J'+nu1*nu0';                                        % 13.106

author: PRML

repo: PRMLT

save_path: github-repos/MATLAB/PRML-PRMLT


text:
%firws() - Designs windowed sinc type I linear phase FIR filter
%
% Usage:
%   >> b = firws(m, f);
%   >> b = firws(m, f, w);
%   >> b = firws(m, f, t);
%   >> b = firws(m, f, t, w);
%
% Inputs:
%   m - filter order (mandatory even)
%   f - vector or scalar of cutoff frequency/ies (-6 dB;
%       pi rad / sample)
%
% Optional inputs:
%   w - vector of length m + 1 defining window {default blackman}
%   t - 'high' for highpass, 'stop' for bandstop filter {default low-/
%       bandpass}
%
% Output:
%   b - filter coefficients
%
% Example:
%   fs = 500; cutoff = 0.5; df = 1;
%   m  = firwsord('hamming', fs, df);
%   b  = firws(m, cutoff / (fs / 2), 'high', windows('hamming', m + 1));
%
% References:
%   Smith, S. W. (1999). The scientist and engineer's guide to digital
%   signal processing (2nd ed.). San Diego, CA: California Technical
%   Publishing.
%
% Author: Andreas Widmann, University of Leipzig, 2005
%
% See also:
%   firwsord, invfirwsord, kaiserbeta, windows

%123456789012345678901234567890123456789012345678901234567890123456789012

% Copyright (C) 2005 Andreas Widmann, University of Leipzig, widmann@uni-leipzig.de
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
%
% $Id$

function [b, a] = firws(m, f, t, w)

    a = 1;

    if nargin < 2
        ft_error('Not enough input arguments');
    end
    if length(m) > 1 || ~isnumeric(m) || ~isreal(m) || mod(m, 2) ~= 0 || m < 2
        ft_error('Filter order must be a real, even, positive integer.');
    end
    f = f / 2;
    if any(f <= 0) || any(f >= 0.5)
        ft_error('Frequencies must fall in range between 0 and 1.');
    end
    if nargin < 3 || isempty(t)
        t = '';
    end
    if nargin < 4 || isempty(w)
        if ~isempty(t) && ~ischar(t)
            w = t;
            t = '';
        else
            w = windows('blackman', (m + 1));
        end
    end
    w = w(:)'; % Make window row vector

    b = fkernel(m, f(1), w);

    if length(f) == 1 && strcmpi(t, 'high')
        b = fspecinv(b);
    end

    if length(f) == 2
        b = b + fspecinv(fkernel(m, f(2), w));
        if isempty(t) || ~strcmpi(t, 'stop')
            b = fspecinv(b);
        end
    end

% Compute filter kernel
function b = fkernel(m, f, w)
    m = -m / 2 : m / 2;
    b(m == 0) = 2 * pi * f; % No division by zero
    b(m ~= 0) = sin(2 * pi * f * m(m ~= 0)) ./ m(m ~= 0); % Sinc
    b = b .* w; % Window
    b = b / sum(b); % Normalization to unity gain at DC

% Spectral inversion
function b = fspecinv(b)
    b = -b;
    b(1, (length(b) - 1) / 2 + 1) = b(1, (length(b) - 1) / 2 + 1) + 1;


text:
function tests = PlanTest
  tests = functiontests(localfunctions);
end

function setupOnce(tc)
    clc
    load map1          % load map
    tc.TestData.map = map;
    tc.TestData.goal = [50 30];
    tc.TestData.start = [20 10];
end

function teardownOnce(tc)
    close all
end

function bug2_test(tc)

    nav = Bug2(tc.TestData.map);
    tc.verifyInstanceOf(nav, 'Bug2');
    s = nav.char();
    verifyTrue(tc, ischar(s) );

    nav.plot();

    p = nav.query(tc.TestData.start, tc.TestData.goal);

    tc.verifyTrue(size(p,2) == 2, 'plan must have 2 columns');
    tc.verifyEqual(p(1,:), tc.TestData.start, 'start point must be on path');
    tc.verifyEqual(p(end,:), tc.TestData.goal, 'goal point must be on path');
    tc.verifyFalse( any( nav.isoccupied(p') ), 'path must have no occupied cells');

    nav.plot()
    nav.plot(p);

    tc.assumeTrue(ispc || ismac);  % FILTER video generation for Travis
    fname = fullfile(tempdir, 'bug.mp4');
    p = nav.query(tc.TestData.start, tc.TestData.goal, 'animate', 'movie', fname);
    tc.verifyTrue(exist(fname, 'file') == 2);
    delete(fname);

    tc.verifyError( @() nav.plan(), 'RTB:Bug2:badcall');

    map = zeros(10,10);
    map(3:7,3:7) = 1;
    map(4:6,4:6) = 0;
    nav = Bug2(map);
    tc.verifyError( @() nav.query([5 5], [2 2]), 'RTB:bug2:noplan');

end

function dxform_test(tc)

    nav = DXform(tc.TestData.map);
    tc.verifyInstanceOf(nav, 'DXform');

    s = nav.char();
    verifyTrue(tc, ischar(s) );

    nav.plot();
    nav.plan(tc.TestData.goal);
    nav.plot();
    nav.plan(tc.TestData.goal, 'animate');

    p = nav.query(tc.TestData.start);

    tc.verifyTrue(size(p,2) == 2, 'plan must have 2 columns');
    tc.verifyEqual(p(1,:), tc.TestData.start, 'start point must be on path');
    tc.verifyEqual(p(end,:), tc.TestData.goal, 'goal point must be on path');
    tc.verifyFalse( any( nav.isoccupied(p') ), 'path must have no occupied cells');

    nav.query(tc.TestData.start, 'animate');

    nav.plot();
    nav.plot(p);

    nav.plot3d();
    nav.plot3d(p);
end

function distancexform_test(tc)
    map = zeros(10,10);
    map(4:6,4:6) =1;
    %args= {'verbose'}
    args = {};

    dx1 = distancexform(map, [5 8], 'fast', 'noipt', args{:});
    tc.verifyClass(dx1, 'double');
    tc.verifyEqual(dx1(8,5), 0);
    tc.verifyTrue(all(all(isnan(dx1(map==1)))));
    i=sub2ind(size(dx1), 8, 5);
    tc.verifyTrue(all(dx1(i+[-11 -10 -9 -1 1 11 10 9])) > 0);

    tc.verifySize(dx1, size(map));
    dx2 = distancexform(map, [5 8], 'nofast', 'ipt', args{:});
    tc.verifyClass(dx1, 'double');
    tc.verifySize(dx1, size(map));
    dx3 = distancexform(map, [5 8], 'nofast', 'noipt', args{:});
    tc.verifyClass(dx1, 'double');
    tc.verifySize(dx1, size(map));

    tc.verifyEqual(dx1, dx2, 'absTol', 1e-6, 'MEX ~= bwdist');
    tc.verifyEqual(dx1, dx3, 'MEX ~= MATLAB');

    dx1 = distancexform(map, [5 8], 'cityblock', 'fast', 'noipt', args{:});
    tc.verifyClass(dx1, 'double');
    tc.verifyTrue(all(all(isnan(dx1(map==1)))));
    i=sub2ind(size(dx1), 8, 5);
    tc.verifyTrue(all(dx1(i+[-11 -10 -9 -1 1 11 10 9])) > 0);
    tc.verifySize(dx1, size(map));
    dx2 = distancexform(map, [5 8], 'cityblock', 'nofast', 'ipt', args{:});
    tc.verifyClass(dx1, 'double');
    tc.verifySize(dx1, size(map));
    dx3 = distancexform(map, [5 8], 'cityblock', 'nofast', 'noipt', args{:});
    tc.verifyClass(dx1, 'double');
    tc.verifySize(dx1, size(map));

    tc.verifyEqual(dx1, dx2, 'absTol', 1e-6, 'MEX ~= bwdist');
    tc.verifyEqual(dx1, dx3, 'MEX ~= MATLAB');

    dx1 = distancexform(map, [5 8], 'cityblock', 'animate');
    dx1 = distancexform(map, [5 8], 'cityblock', 'animate', 'delay', 0.1);

    tc.assumeTrue(ispc || ismac);  % FILTER video generation for Travis
    fname = fullfile(tempdir, 'bug.mp4');
    dx1 = distancexform(map, [5 8], 'cityblock', 'animate', 'movie', fname);
    tc.verifyTrue(exist(fname, 'file') == 2);
    delete(fname)
end

function dstar_test(tc)

    % create a planner
    nav = Dstar(tc.TestData.map, 'quiet');
    tc.verifyInstanceOf(nav, 'Dstar');

    s = nav.char();
    verifyTrue(tc, ischar(s) );

    nav.plot();

    % plan path to goal
    nav.plan(tc.TestData.goal);
    nav.plan(tc.TestData.goal, 'animate');

    % execute it
    p = nav.query(tc.TestData.start);

        tc.verifyTrue(size(p,2) == 2, 'plan must have 2 columns');
    tc.verifyEqual(p(1,:), tc.TestData.start, 'start point must be on path');
    tc.verifyEqual(p(end,:), tc.TestData.goal, 'goal point must be on path');
    tc.verifyFalse( any( nav.isoccupied(p') ), 'path must have no occupied cells');

    nav.query(tc.TestData.start, 'animate');

    nav.plot();
    nav.plot(p);

    % add a swamp
    for r=78:85
        for c=12:45
            nav.modify_cost([c;r], 2);
        end
    end

    % replan
    nav.plan();

    % show new path
    nav.query(tc.TestData.start);

    p = nav.query(tc.TestData.start);
    tc.verifyTrue(size(p,2) == 2, 'plan must have 2 columns');
    tc.verifyEqual(p(1,:), tc.TestData.start, 'start point must be on path');
    tc.verifyEqual(p(end,:), tc.TestData.goal, 'goal point must be on path');
    tc.verifyFalse( any( nav.isoccupied(p') ), 'path must have no occupied cells');

    nav.plot(p);

    nav.modify_cost([12 45; 78 85], 2);
    nav.modify_cost([12 13 14; 78 79 80], [2 3 4]);
end

function prm_test(tc)

    randinit
    nav = PRM(tc.TestData.map);
        tc.verifyInstanceOf(nav, 'PRM');

    s = nav.char();
    verifyTrue(tc, ischar(s) );

    nav.plot();

    nav.plan();
    nav.plot();

    p = nav.query(tc.TestData.start, tc.TestData.goal);

    tc.verifyTrue(size(p,2) == 2, 'plan must have 2 columns');
    tc.verifyEqual(p(1,:), tc.TestData.start, 'start point must be on path');
    tc.verifyEqual(p(end,:), tc.TestData.goal, 'goal point must be on path');
    tc.verifyFalse( any( nav.isoccupied(p') ), 'path must have no occupied cells');

    nav.query(tc.TestData.start, tc.TestData.goal);
    nav.plot(p);
end

function lattice_test(tc)

    nav = Lattice();
    tc.verifyInstanceOf(nav, 'Lattice');

    nav.plan('iterations', 8);
    nav.plot()

    start = [1 2 pi/2]; goal = [2 -2 0];
    p = nav.query( start, goal );
    verifyEqual(tc, size(p,1), 7);
    tc.verifyTrue(size(p,2) == 3, 'plan must have 3 columns');
    tc.verifyEqual(p(1,:), start, 'AbsTol', 1e-10, 'start point must be on path');
    tc.verifyEqual(p(end,:), goal, 'AbsTol', 1e-10, 'goal point must be on path');

    nav.plot

    nav.plan('cost', [1 10 10])
    p = nav.query( start, goal );
    verifyEqual(tc, size(p,1), 9);
    tc.verifyTrue(size(p,2) == 3, 'plan must have 3 columns');
    tc.verifyEqual(p(1,:), start, 'AbsTol', 1e-10, 'start point must be on path');
    tc.verifyEqual(p(end,:), goal, 'AbsTol', 1e-10, 'goal point must be on path');

    load road
    nav = Lattice(road, 'grid', 5, 'root', [50 50 0])
    nav.plan();
    start = [30 45 0]; goal = [50 20 0];
    p = nav.query(start, goal);
    tc.verifyTrue(size(p,2) == 3, 'plan must have 3 columns');
    tc.verifyEqual(p(1,:), start, 'AbsTol', 1e-10, 'start point must be on path');
    tc.verifyEqual(p(end,:), goal, 'AbsTol', 1e-10, 'goal point must be on path');
    tc.verifyFalse( any( nav.isoccupied(p(:,1:2)') ), 'path must have no occupied cells');

end

function rrt_test(tc)
    randinit

    car = Bicycle('steermax', 0.5);
    %nav = RRT(car, 'goal', goal, 'range', 5);
    nav = RRT(car, 'npoints', 400);
    tc.verifyInstanceOf(nav, 'RRT');

    s = nav.char();

    nav.plan();
    nav.plot();

    start = [0 0 0]; goal = [0 2 0];
    p = nav.query(start, goal);

    tc.verifyTrue(size(p,2) == 3, 'plan must have 3 columns');
    tc.verifyFalse( any( nav.isoccupied(p(:,1:2)') ), 'path must have no occupied cells');

    % bigger example with obstacle
    load road

    randinit
    nav = RRT(car, road, 'npoints', 1000, 'root', [50 22 0], 'simtime', 4);
    nav.plan();
    p = nav.query([40 45 0], [50 22 0]);
    tc.verifyTrue(size(p,2) == 3, 'plan must have 3 columns');
    tc.verifyFalse( any( nav.isoccupied(p(:,1:2)') ), 'path must have no occupied cells');

end

author: petercorke

repo: robotics-toolbox-matlab


text: cd ~/project/data/qdots

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%% OPTI Toolbox Differentiation Demo
%
% This file demonstrates each of the different differentiation algorithms
% supplied with the OPTI Toolbox. Note you should be familiar with the
% operation of OPTI Toolbox by reading the accompanying examples, as well
% as you should have completed the NLP demo.
%
%   Copyright (C) 2011 Jonathan Currie (I2C2)

%% Gradient vs Jacobian
% The terms Gradient and Jacobian are associated with differentiation
% functions, and for the purposes of this toolbox they are basically
% interchangeable:
%
% Gradient - The first derivative of the objective function (1 x n)
% Jacobian - The first derivative of the constraint function (m x n)
%
% Where you will note the main difference is the gradient is a vector, and
% the Jacobian is a matrix. However this is just a rule of thumb and will
% not hold in all instances (constraint functions with one row for
% example!). All functions are named as returning a Jacobian, which could
% also be a Gradient Vector.

%% Problem 1
% For the following examples we are going to use this nonlinear objective
% function:

fun = @(x) sin(pi*x(1)/12) * cos(pi*x(2)/16);

x = [0.75 0.25]';

%% Example 1 - Automatic Differentiation (AD)
% Automatic differentiation is one of the most powerful differentiation
% strategies which can provide error free gradients of a function by
% applying the chain rule to each operation.
%
% The AD routines implemented in OPTI are supplied by adiff, a Matlab
% project by William McIlhagga. While most Matlab functions are overloaded,
% you cannot use external code (i.e. via MEX) or any toolbox or class
% functions.

dx = autoJac(fun,x)

%% Example 2 - Numerical Differentiation (ND)
% Numerical differentiation using finite differences is a computationally
% expensive procedure which can result in an inaccurate gradient if the
% internal perturbations are not chosen correctly.
%
% The ND routine implemented in OPTI is the Intel MKL djacobi function
% which approximates the derivative using central differences. This is
% implemented via a MEX function which repeatedly calls the function in
% order to close in on the gradient.

dx = mklJac(fun,x)

%% Example 3 - Symbolic Differentiation (SD)
% Symbolic differentiation analytically differentiates the function as a
% symbolic expression, resulting in a single expression for the gradient.
% Complications occur if the function cannot be analytically differentiated
% or the symbolic routine cannot find a derivative.
%
% The SD routine implemented in OPTI uses the Matlab Symbolic Toolbox as
% well as two wrapper functions in order to generate the gradient. The
% wrapper functions convert the function handle to a symbolic expression
% and vice-versa.

grad = symJac(fun)

if(~isempty(grad)) %don't run if Symbolic Toolbox not installed
    dx = grad(x)
end

%% Problem 2
% For the next few examples we will be solving this NLP:

obj = @(x) log(1+x(1)^2) - x(2);

lb = [-2 -2]';
ub = [2 2]';

%% Example 4 - Applying AD to NLP Solving
% To use AD to generate the objective gradient for IPOPT create a function
% handle which calls autoJac, and then pass this to the optiprob function:

grad = @(x) autoJac(obj,x);

opts = optiset('solver','ipopt');
Opt = opti('obj',obj,'grad',grad,'bounds',lb,ub,'options',opts)

[x,fval,exitflag,info] = solve(Opt,[0;0]);
fval
info

%% Example 5 - Applying ND to NLP Solving
% To use ND to generate the objective gradient for IPOPT just replace the
% above with mklJac:

grad = @(x) mklJac(obj,x);

Opt = opti('obj',obj,'grad',grad,'bounds',lb,ub,'options',opts)

[x,fval,exitflag,info] = solve(Opt,[0;0]);
fval
info

%% Example 6 - Applying SD to NLP Solving
% To use SD we can generate the analytical gradient and use this as our
% gradient function:

grad = symJac(obj)

if(~isempty(grad))
    Opt = opti('obj',obj,'grad',grad,'bounds',lb,ub,'options',opts)

    [x,fval,exitflag,info] = solve(Opt,[0;0]);
    fval
    info
else
    display('Cannot run example without Symbolic Toolbox!');
end

%% Conclusion
% By default OPTI Toolbox will always use Numerical Differentiation
% (mklJac) for gradients / Jacobians if one has not been provided. This
% provides the best combination of flexibility with complex objective /
% constraint functions and performance. However if you can use SD to
% generate a gradient this would be a preferred method.
%
% If you find the optimizer is failing with ND and AD is a suitable
% candidate you can try it to ensure your gradients are accurate.

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function cot_values_test ( )

%*****************************************************************************80
%
%% COT_VALUES_TEST demonstrates the use of COT_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 January 2015
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'COT_VALUES_TEST:\n' );
  fprintf ( 1, '  COT_VALUES stores values of \n' );
  fprintf ( 1, '  the cotangent function.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = cot_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% function to show limo results - deprecated

% Copyright (C) 2016 Ramon Martinez
%
% 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.

classdef pop_limoresults < handle

    %class properties - access is private so nothing else can access these
    properties (Access = public)
        gui_h;
        study;
        limofiles_path;
        limofiles_level1;
        limofiles_level2;
        regnames;
        catvarnames;
        datorica_indx;
        reg_indx;
        cat_indx;
        string_mplot;
        limostruct_indx
    end

    methods
        % =================================================================
        function obj = pop_limoresults(STUDY,analysis)

            if ~ismember(analysis,{'dat','ica'}), return; end

            % Check if STUDY.design.limo
            %--------------------------------------------------------------
            datorica_list = {'dat','ica'};
            obj.datorica_indx = find(strcmp(analysis,datorica_list));

            [obj.gui_h,STUDY] = guibuilder(STUDY,analysis);
            obj.string_mplot  = get(obj.gui_h.popupmenu_measure2plot,'String');
            obj.gui_h         = guihandles(obj.gui_h.fig);
            obj.study         = STUDY;

            level_tmp       = get(obj.gui_h.popupmenu_level         ,'Value');
            m2plot_tmp      = get(obj.gui_h.popupmenu_measure2plot  ,'Value');
            var2plot_indx   = get(obj.gui_h.popupmenu_modelvar2plot ,'Value');
            [var2plot,filespath] = getmeasures2plot(obj.study,level_tmp-1,m2plot_tmp,obj.datorica_indx);

            % Detecting Measure and updating limofiles_path and
            %--------------------------------------------------------------
            obj.limofiles_path            = filespath;
            obj.limofiles_level1.index    = var2plot_indx;
            obj.limofiles_level1.guiname  = var2plot.guiname;
            obj.limofiles_level1.pathname = var2plot.pathname;
            % Updating Electrode list
            %--------------------------------------------------------------
            try
                load(fullfile(filespath,'LIMO.mat'),'LIMO');
                electoplot_list = ['All Channels';{LIMO.data.chanlocs.labels}'];
                electoplot_indx = 1;
                set(obj.gui_h.listbox_elect2plot,'String',electoplot_list);
                set(obj.gui_h.listbox_elect2plot,'Value',electoplot_indx);
            catch
            end

            % Creating Name of Regressors
            %--------------------------------------------------------------
            try
                % --
                % Getting the Names of the variables
                tmpregnames  = {};
                vartype_tmp  = {};
                reg_indx     = {};
                for i = 1: length(STUDY.design(STUDY.currentdesign).variable)
                    if strcmp(STUDY.design(STUDY.currentdesign).variable(i).vartype,'categorical')
                        if isempty(tmpregnames)
                            for j = 1: length(STUDY.design(STUDY.currentdesign).variable(i).value)
                                if  ~iscell(STUDY.design(STUDY.currentdesign).variable(i).value{1})
                                    tmpregnames{j} = [STUDY.design(STUDY.currentdesign).variable(i).label '_' STUDY.design(STUDY.currentdesign).variable(i).value{j}];
                                else
                                    for jvars = 1:length(STUDY.design(STUDY.currentdesign).variable(i).value{j})
                                        if jvars == 1
                                            jointtmp = STUDY.design(STUDY.currentdesign).variable(i).value{j}{jvars};
                                        else
                                            jointtmp = strcat (jointtmp, '&',  STUDY.design(STUDY.currentdesign).variable(i).value{j}{jvars});
                                        end
                                    end
                                    tmpregnames{j} = [STUDY.design(STUDY.currentdesign).variable(i).label '_' jointtmp{:}];
                                end
                            end
                            [vartype_tmp{1:length(tmpregnames)}] = deal(STUDY.design(STUDY.currentdesign).variable(i).vartype);
                        else
                            for j = 1:length(STUDY.design(STUDY.currentdesign).variable(i).value)
                                if  ~iscell(STUDY.design(STUDY.currentdesign).variable(i).value{1})
                                    tmpregnames{end+1} = [STUDY.design(STUDY.currentdesign).variable(i).label '_' STUDY.design(STUDY.currentdesign).variable(i).value{j}];
                                else
                                    for jvars = 1:length(STUDY.design(STUDY.currentdesign).variable(i).value{j})
                                        if jvars == 1
                                            jointtmp = STUDY.design(STUDY.currentdesign).variable(i).value{j}{jvars};
                                        else
                                            jointtmp = strcat (jointtmp, '&',  STUDY.design(STUDY.currentdesign).variable(i).value{j}{jvars});
                                        end
                                    end
                                    tmpregnames{end+1} = [STUDY.design(STUDY.currentdesign).variable(i).label '_' jointtmp];
                                end

                                %
                                vartype_tmp{end+1} = STUDY.design(STUDY.currentdesign).variable(i).vartype;
                            end
                        end
                    else
                        if isempty(tmpregnames)
                            tmpregnames{1} = STUDY.design(STUDY.currentdesign).variable(i).label;
                            vartype_tmp{1} = STUDY.design(STUDY.currentdesign).variable(i).vartype;

                        else
                            tmpregnames{end+1} = STUDY.design(STUDY.currentdesign).variable(i).label;
                            vartype_tmp{end+1} = STUDY.design(STUDY.currentdesign).variable(i).vartype;
                        end
                    end
                end

                % Set regressors order as LIMO assumes it (first all categorical, then all continuous)
                %----------------------------------------------------------------------------------
                cont_indx   = find(strcmp('continuous',vartype_tmp));
                cat_indx    = find(strcmp('categorical',vartype_tmp));
                tmpregnames = tmpregnames([cat_indx cont_indx]);

                % Adding Baseline
                %----------------------------------------------------------
                tmpregnames{end+1} = 'Baseline';
                % --
                limo = load(fullfile(obj.limofiles_path,'LIMO.mat' ));
                Nreg = size(limo.LIMO.design.X,2);

                 % Rgressors names in case of split continuous ones
                 %---------------------------------------------------------
                if ~isempty(cont_indx) && (Nreg ~= length(tmpregnames))
                    display('Split regressor detected. Giving generic names to continuous regressors ...' );
                    if ~isempty(cat_indx), inintcont = length(cat_indx)+1; else inintcont = 1 ; end

                    for i = inintcont:(Nreg-1)
                        tmpregnames{i} = ['Cont_' num2str(i - (inintcont - 1))];
                        display(['Naming Regressor: ' tmpregnames{i}]);
                    end
                end

                %---------------------------------------------------------
                cat_comb = [];
                if ~isempty(cat_indx) && size(cat_indx,2) > 1
                    cat_comb = nchoosek(1:length(cat_indx),2);
                end

                % Generating cat index
                %---------------------------------------------------------
                if ~isempty(cat_indx)
                    allcat_indx  = [1: length(cat_indx),(length(tmpregnames) + 1):(length(tmpregnames) + size(cat_comb,1))];
                end

                % Adding Regressor to names
                %---------------------------------------------------------
                for i = 1:length(tmpregnames)
                    %tmpregnames{i} = ['Regressor_' tmpregnames{i}];
                    reg_indx{i}    = ['[' num2str(i) ']'];
                end

                % Adding Combinations of cat regressors
                %----------------------------------------------------------
                if ~isempty(cat_comb)
                    for j =1:size(cat_comb,1)
                        strtmp = tmpregnames(cat_comb(j,:));
                        tmpregnames{end+1} = [strtmp{1} sprintf(' - %s', strtmp{2:end})];
                        clear strtmp;
                        reg_indx{end+1}    = ['[' num2str(cat_comb(j,:) ) ']'];
                    end
                end

                % Retrieving list of independent variables
                if ~isempty(cat_indx)
                    tmpregnames2 = tmpregnames(allcat_indx);
                    cat_indx     =  reg_indx(allcat_indx);
                else
                    tmpregnames2 = [];
                end

                % Adding 'Regressor' to names
                %---------------------------------------------------------
                for i = 1:length(tmpregnames)
                    tmpregnames{i} = ['Regressor_' tmpregnames{i}];
                    plusindx = strfind(tmpregnames{i},' - ');
                    if ~isempty(plusindx)
                        tmpregnames{i} = [tmpregnames{i}(1:plusindx+2) 'Regressor_' tmpregnames{i}(plusindx+3:end)];
                    end
                end

                % Updating fields
                %----------------------------------------------------------
                obj.regnames    = tmpregnames;
                obj.catvarnames = tmpregnames2;
                obj.reg_indx    = reg_indx;
                obj.cat_indx    = cat_indx;

            catch
                % Case where model is not computed
                %----------------------------------------------------------
                display(['File LIMO.mat not founded in :'  fullfile(obj.limofiles_path,'LIMO.mat')]);
                set(obj.gui_h.popupmenu_modelvar2plot,'Enable','off');
                set(obj.gui_h.popupmenu_plottype     ,'Enable','off');
                set(obj.gui_h.pushbutton_plot        ,'Enable','off');
                eeglab_warning('Make sure to compute the model for this measure');
            end

            % CALLBACKS
            %--------------------------------------------------------------
            %set the callback function for button_plot
            set(obj.gui_h.pushbutton_plot, 'Callback', @obj.callback_plot);

            %set the callback function for popupmenu_level
            set(obj.gui_h.popupmenu_level,'Callback', @obj.callback_popupmenu_level);

            %set the callback function for popupmenu_measure2plot
            set(obj.gui_h.popupmenu_measure2plot,'Callback', @obj.callback_popupmenu_level);

            %set the callback function for popupmenu_plottype
            set(obj.gui_h.popupmenu_plottype,'Callback', @obj.callback_popupmenu_plottype);

            %set the callback function for popupmenu_dataorresult
            set(obj.gui_h.popupmenu_dataorresult,'Callback', @obj.callback_popupmenu_dataorresult);

            %set the callback function for popupmenu_modelvar2plot
            set(obj.gui_h.popupmenu_modelvar2plot,'Callback', @obj.callback_popupmenu_modelvar2plot);

            %set the callback function for checkbox_stats
            set(obj.gui_h.checkbox_stats,'Callback', @obj.callback_checkbox_stats);

            %set the callback function for checkbox_stats
            set(obj.gui_h.listbox_elect2plot,'Callback', @obj.callback_listbox_elect2plot);

            %set the callback function for buttom_designadd
            set(obj.gui_h.pushbutton_designadd, 'Callback', @obj.callback_pushbutton_designadd);

        end
        % =================================================================
        function obj = callback_popupmenu_level(obj,~,~)

            val_level = get(obj.gui_h.popupmenu_level,'Value');        % Getting value (subject from level)
            val_mplot = get(obj.gui_h.popupmenu_measure2plot,'Value'); % Getting value from popupmenu_measure2plot
            val_dor = get(obj.gui_h.popupmenu_dataorresult,'Value');   % Getting value from popupmenu_dataorresult

            % Getting values popupmenu_modelvar2plot
            string_vplot  = get(obj.gui_h.popupmenu_modelvar2plot ,'String');
            val_vplot     = get(obj.gui_h.popupmenu_modelvar2plot ,'Value' );
            val_typeplot  = get(obj.gui_h.popupmenu_plottype      ,'Value' );
            string_elec   = get(obj.gui_h.listbox_elect2plot      ,'String');
            val_elec      = get(obj.gui_h.listbox_elect2plot      ,'Value' );

            % Single subject level
            if val_level ~= 1

                [var2plot,filespath,tmp,obj.study] = getmeasures2plot(obj.study,val_level-1,val_mplot,obj.datorica_indx,1);

                if strcmp(get(obj.gui_h.popupmenu_measure2plot,'String'),{' '})
                    set(obj.gui_h.popupmenu_measure2plot  ,'String', obj.string_mplot);
                    set(obj.gui_h.popupmenu_measure2plot  ,'Value',1);
                end

                % Updating Electrode list
                % =========================================================

                if ~strcmp('No Variables Computed',var2plot.guiname{1})

                    load(fullfile(filespath,'LIMO.mat'),'LIMO');
                    electoplot_list = ['All Channels';{LIMO.data.chanlocs.labels}'];
                    set(obj.gui_h.listbox_elect2plot,'String',electoplot_list);

                    % Getting index for var2plot_list and setting electrode
                    % index
                    if val_typeplot == 1 || val_typeplot == 2

                        % Updating electrode list val
                        %----------------------------------------------
                        set(obj.gui_h.listbox_elect2plot,'Value',1);

                         % Just Checking if the selected variable exist in
                         % the sub (case of Results only)
                         %-------------------------------------------------
                         if val_dor == 2
                            if ~ismember(string_vplot(val_vplot),var2plot.guiname)
                                % Getting index of  'Condition_effect_1.mat' for default value
                                % 1 st pass
                                var2plot_indx = find(strcmp(var2plot.guiname,'Condition_effect_1'), 1);
                                % 2nd pass
                                if isempty(var2plot_indx)
                                    var2plot_indx = find(strcmp(var2plot.guiname,'Covariate_effect_1'),1);
                                end
                                % Setting to 1 if all fails
                                if isempty(var2plot_indx)
                                    var2plot_indx = 1;
                                end
                                set(obj.gui_h.popupmenu_modelvar2plot,'String',var2plot.guiname);
                                set(obj.gui_h.popupmenu_modelvar2plot,'Value',var2plot_indx);
                            else
                                var2plot_indx = find(strcmp(string_vplot(val_vplot),var2plot.guiname));
                                set(obj.gui_h.popupmenu_modelvar2plot,'String',var2plot.guiname);
                                set(obj.gui_h.popupmenu_modelvar2plot,'Value',var2plot_indx);
                            end
                         end
                        %--------------------------------------------------
                    elseif val_typeplot == 3
                        %var2plot_list = string_vplot;

                        % Updating electrode list val
                        %--------------------------------------------------
                        newval_elec = find(strcmp(string_elec(val_elec),electoplot_list));
                        if isempty(newval_elec)
                            newval_elec = 2;
                        end
                        set(obj.gui_h.listbox_elect2plot,'Value',newval_elec);
                    end

                    % Enabling on GUI features
                    %------------------------------------------------------
                    set(obj.gui_h.popupmenu_measure2plot  ,'Enable','on');
                    set(obj.gui_h.popupmenu_dataorresult  ,'Enable','on');
                    set(obj.gui_h.popupmenu_modelvar2plot ,'Enable','on');
                    set(obj.gui_h.popupmenu_plottype      ,'Enable','on');
                    set(obj.gui_h.pushbutton_plot         ,'Enable','on');
                    set(obj.gui_h.checkbox_stats	      ,'Enable','on');
                    set(obj.gui_h.edit_pval               ,'Enable','on');
                else
                    % Enabling off GUI features
                    %------------------------------------------------------
                    set(obj.gui_h.popupmenu_dataorresult  ,'Enable','off');
                    set(obj.gui_h.popupmenu_modelvar2plot ,'Enable','off');
                    set(obj.gui_h.popupmenu_plottype      ,'Enable','off');
                    set(obj.gui_h.pushbutton_plot         ,'Enable','off');

                    eeglab_warning('Make sure to compute the model for this measure');
                end

                % Updating limofiles_level1 structure
                %----------------------------------------------------------
                if strcmp(get(obj.gui_h.popupmenu_modelvar2plot,'Enable'),'off')
                    obj.limofiles_level1.guiname = [];
                else
                    obj.limofiles_path          = filespath;
                    obj.limofiles_level1.guiname  = var2plot.guiname;
                    obj.limofiles_level1.pathname = var2plot.pathname;
                    obj.limofiles_level1.index    = get(obj.gui_h.popupmenu_modelvar2plot,'Value');
                end

            % Group level
            else
                measure_list  = {'erp','spec'};
                datorica_list = {'dat', 'ica'};
                datatype = [datorica_list{obj.datorica_indx} measure_list{get(obj.gui_h.popupmenu_measure2plot,'Value')}];

                obj.limostruct_indx = find(~cellfun(@isempty,strfind({obj.study.design(obj.study.currentdesign).limo.datatype},datatype)));
                if ~isempty(obj.limostruct_indx) && isfield(obj.study.design(obj.study.currentdesign).limo, 'groupmodel') && ~isempty(obj.study.design(obj.study.currentdesign).limo(obj.limostruct_indx).groupmodel)
                    obj.limofiles_level2.guiname  = {obj.study.design(obj.study.currentdesign).limo(obj.limostruct_indx).groupmodel.guiname}';
                    obj.limofiles_level2.pathname = {obj.study.design(obj.study.currentdesign).limo(obj.limostruct_indx).groupmodel.filename}';
                    obj.limofiles_level2.index = 1;
                    set(obj.gui_h.popupmenu_modelvar2plot,'String',obj.limofiles_level2.guiname);
                    set(obj.gui_h.popupmenu_modelvar2plot,'Value',1);

                    % Loading expected electrodes
                    %----------------------------------------------------------
                    try
                        load(obj.study.design(obj.study.currentdesign).limo(obj.limostruct_indx).chanloc);
                        electoplot_list = ['All Channels';{expected_chanlocs.labels}'];
                        set(obj.gui_h.listbox_elect2plot,'String',electoplot_list);
                    catch
                        % implement this
                    end

                    if val_typeplot == 1 || val_typeplot == 2

                        % Updating electrode list val
                        %----------------------------------------------
                        set(obj.gui_h.listbox_elect2plot,'Value',1);

                    elseif val_typeplot == 3

                        % Updating electrode list val
                        %--------------------------------------------------
                        newval_elec = find(strcmp(string_elec(val_elec),electoplot_list));
                        if isempty(newval_elec)
                            newval_elec = 2;
                        end
                        set(obj.gui_h.listbox_elect2plot,'Value',newval_elec);
                    end

                    % Enabling on/off GUI features
                    %------------------------------------------------------
                    set(obj.gui_h.popupmenu_dataorresult  ,'Value', 2);
                    set(obj.gui_h.popupmenu_dataorresult  ,'Enable','off');
                    set(obj.gui_h.popupmenu_modelvar2plot ,'Enable','on');
                    set(obj.gui_h.popupmenu_plottype      ,'Enable','on');
                    set(obj.gui_h.pushbutton_plot         ,'Enable','on');
                    set(obj.gui_h.popupmenu_mcc	          ,'Enable','off');
                    set(obj.gui_h.popupmenu_compmethod    ,'Enable','off');
                    set(obj.gui_h.checkbox_stats	      ,'Enable','off');
                    set(obj.gui_h.edit_pval               ,'Enable','on');
                    set(obj.gui_h.pushbutton_designadd    ,'Enable','off');% Check this
                else
                    % Disabling GUI features
                    %------------------------------------------------------
                    set(obj.gui_h.popupmenu_dataorresult  ,'Enable','off');
                    set(obj.gui_h.popupmenu_modelvar2plot ,'Enable','off');
                    set(obj.gui_h.popupmenu_plottype      ,'Enable','off');
                    set(obj.gui_h.pushbutton_plot         ,'Enable','off');
                    set(obj.gui_h.pushbutton_designadd    ,'Enable','off');
                    set(obj.gui_h.edit_pval               ,'Enable','off');
                    eeglab_warning('No group level results founded');
                end
            end
        end
        % =================================================================
        function obj = callback_popupmenu_modelvar2plot(obj,~,~)

            var2plot_indx = get(obj.gui_h.popupmenu_modelvar2plot ,'Value');
            if get(obj.gui_h.popupmenu_level,'Value') ~= 1
                obj.limofiles_level1.index = var2plot_indx;
            else
                obj.limofiles_level2.index = var2plot_indx;
            end
        end
        % =================================================================
        function obj = callback_popupmenu_plottype(obj,~,~)

            val_typeplot = get( obj.gui_h.popupmenu_plottype    ,'Value');
            val_dor      = get(obj.gui_h.popupmenu_dataorresult ,'Value');
            val_level    = get( obj.gui_h.popupmenu_level       ,'Value');

            if (val_typeplot == 1 || val_typeplot == 2)
                if val_level ~= 1
                    % Call callback_popupmenu_level
                    obj = callback_popupmenu_level(obj);
                end

                % Updating listbox_elect2plot and disabling it
                set(obj.gui_h.listbox_elect2plot   ,'Value',1);
                set(obj.gui_h.listbox_elect2plot   ,'Enable','off');

                % Disabling 'Add design' button
                if  val_dor == 1
                    set(obj.gui_h.pushbutton_designadd ,'Enable','off');
                else
                    set(obj.gui_h.pushbutton_designadd ,'Enable','on');
                end

            elseif val_typeplot == 3

                if val_level ~= 1
                    if val_dor == 2
                        listtmp = obj.regnames;
                    elseif val_dor == 1
                        listtmp = obj.catvarnames;
                    end
                    set(obj.gui_h.popupmenu_modelvar2plot ,'String',listtmp);
                    set(obj.gui_h.popupmenu_modelvar2plot ,'Value',1);
                end
                set(obj.gui_h.listbox_elect2plot   ,'Enable','on');
                set(obj.gui_h.listbox_elect2plot   ,'Value' ,2);
                set(obj.gui_h.pushbutton_designadd ,'Enable','off');
            end

        end
        % =================================================================
        function obj = callback_popupmenu_dataorresult(obj,~,~)
            val_ptype = get( obj.gui_h.popupmenu_plottype   ,'Value');
            val_dor   = get(obj.gui_h.popupmenu_dataorresult,'Value');

            % Case for data
            if val_dor == 1

                if val_ptype == 1 || val_ptype == 2
                    set(obj.gui_h.popupmenu_modelvar2plot,'String',obj.catvarnames);
                    set(obj.gui_h.popupmenu_modelvar2plot,'Value',1);
                elseif val_ptype == 3
                    set(obj.gui_h.popupmenu_modelvar2plot,'String',obj.catvarnames); % THIS MUST INCLUDE CONT!!!!!!!!!!!!!?????
                    set(obj.gui_h.popupmenu_modelvar2plot,'Value',1);
                end

                % 'Add design' disable
                set(obj.gui_h.pushbutton_designadd ,'Enable','off');

            % Case for results
            elseif val_dor == 2
                if val_ptype == 1 || val_ptype == 2
                    % Call callback_popupmenu_level
                    obj = callback_popupmenu_level(obj);

                % 'Add design' enable
                set(obj.gui_h.pushbutton_designadd ,'Enable','on');

                elseif val_ptype == 3
                    set(obj.gui_h.popupmenu_modelvar2plot,'String',obj.regnames);
                    set(obj.gui_h.popupmenu_modelvar2plot,'Value',1);

                    % Call callback_popupmenu_level
                    obj = callback_popupmenu_level(obj);

                    % 'Add design' disable
                    set(obj.gui_h.pushbutton_designadd ,'Enable','off');

                end
            end
        end
        % =================================================================
        function obj = callback_checkbox_stats(obj,~,~)

            if get(obj.gui_h.checkbox_stats        ,'Value');
                set(obj.gui_h.popupmenu_mcc        ,'Enable','on');
                set(obj.gui_h.popupmenu_compmethod ,'Enable','on');
            else
                set(obj.gui_h.popupmenu_mcc        ,'Enable','off');
                set(obj.gui_h.popupmenu_mcc        ,'Value',1);
                set(obj.gui_h.popupmenu_compmethod ,'Enable','off');
                set(obj.gui_h.popupmenu_compmethod ,'Value',1);
            end
        end
        % =================================================================
        function obj = callback_listbox_elect2plot(obj,~,~)

            valtmp_plottype = get( obj.gui_h.popupmenu_plottype,'Value');
            valtmp_chan     = get( obj.gui_h.listbox_elect2plot,'Value');

            if (valtmp_plottype == 3) && (valtmp_chan == 1)
                eeglab_warning('Invalid selection for requested plot. Selecting 1st electrode instead ');
                set(obj.gui_h.listbox_elect2plot,'Value',2);
            end
        end
        % =================================================================
        function obj = callback_plot(obj,~, ~)

             % Prepare inputs for std_limoresults
            level       = get(obj.gui_h.popupmenu_level,'Value');
            if level == 1
                levelval = 2;
                testindxval = obj.limofiles_level2.index;
            else
                subjindxval = level - 1;
                levelval    = 1;
                testindxval = obj.limofiles_level1.index;
            end
            measure_list  = {'erp','spec'};
            datorica_list = {'dat', 'ica'};
            measure_val = [datorica_list{obj.datorica_indx} measure_list{get(obj.gui_h.popupmenu_measure2plot,'Value')}];

            % DEFS
            checkbox_stat_status = get(obj.gui_h.checkbox_stats,'Value');
            if checkbox_stat_status
                handles.p = str2double(get(obj.gui_h.edit_pval,'string'));
                indxtmp      = get(obj.gui_h.popupmenu_mcc,'Value');
                handles.MCC  = indxtmp;
                handles.tfce = 0;
                handles.dir  = pwd;
                indxtmp           = get(obj.gui_h.popupmenu_compmethod,'Value');
                if indxtmp == 1
                    handles.bootstrap = 0;
                else
                    handles.bootstrap = 1;
                end
            else
                handles.p         = 0.05;
                handles.MCC       = 1;
                handles.dir       = pwd;
                handles.bootstrap = 0;
                handles.tfce      = 0;
            end

            % Determine if 'Original Data' or 'Results'
            if get(obj.gui_h.popupmenu_dataorresult,'Value') == 1
                flagdataval = 1;
            else
                flagdataval = 0 ;
            end

            % Getting the plot type
            plottypeval   = get(obj.gui_h.popupmenu_plottype,'Value');
            if plottypeval == 3
                chanindxval = get( obj.gui_h.listbox_elect2plot,'Value')-1;
                if chanindxval == 0,
                    eeglab_warning('Must select a valid electrode');
                    return;
                end
            end
            % ---------------------
            % --- Data ---
            if flagdataval == 1
                % Level 1
                if levelval == 1
                    if plottypeval ~= 3
                        std_limoresults(obj.study,'plottype',plottypeval,'flagdata',flagdataval,'measure',measure_val,'level',1,'subjindx',subjindxval,'regressor',obj.cat_indx(get(obj.gui_h.popupmenu_modelvar2plot,'value')));
                    else
                        std_limoresults(obj.study,'plottype',3,'flagdata',flagdataval,'measure',measure_val,'level',1,'subjindx',subjindxval,'regressor',obj.cat_indx(get(obj.gui_h.popupmenu_modelvar2plot,'value')),'chanindx',chanindxval);
                    end
                    % Level 2
                else
                    fprintf('Under construction');

                end
            % --- Results ---
            else
                % Level 1
                if levelval == 1
                    if plottypeval ~= 3
                    std_limoresults(obj.study,'plottype',plottypeval,'flagdata',0,'measure',measure_val,'level',1,'subjindx',subjindxval,'testindx',testindxval);
                    else
                        std_limoresults(obj.study,'plottype',3,'flagdata',0,'measure',measure_val,'level',1,'subjindx',subjindxval,'regressor',obj.reg_indx(get(obj.gui_h.popupmenu_modelvar2plot,'value')),'chanindx',chanindxval);
                    end

                % Level 2
                else
                    if plottypeval ~= 3
                    std_limoresults(obj.study,'plottype',plottypeval,'flagdata',0,'measure',measure_val,'level',2,'testindx',testindxval);
                    else
                        std_limoresults(obj.study,'plottype',3,'flagdata',0,'measure',measure_val,'level',2,'testindx',testindxval,'chanindx',chanindxval);
                    end
                end
            end
        end
        % =================================================================
        function obj = callback_pushbutton_designadd(obj,~,~)
            clear -global LIMO;

            val_level  = get( obj.gui_h.popupmenu_level        ,'Value');
            val_mplot  = get(obj.gui_h.popupmenu_measure2plot  ,'Value');
            stringtmp  = get(obj.gui_h.popupmenu_modelvar2plot ,'String');
            pwdtmp = pwd;
            limo_contrast_manager(fullfile(obj.limofiles_path,'LIMO.mat'));
            cd(pwdtmp);
            htmp = findall(0,'Type','Figure','Tag','figure_limo_contrast_manager');
            if ~isempty(htmp)
                waitfor(htmp); clear htmp;
            end
            % Get new generated file and add it to the list as current
            % ---------------------------------------------------------
            [currentfiles,tmp,tmp,STUDY] = getmeasures2plot(obj.study,val_level-1,val_mplot,obj.datorica_indx,1);
            obj.study = STUDY;
            newfile = setdiff(currentfiles.guiname,stringtmp);
            if ~isempty(newfile)
                %stringtmp(end+1) = newfile;
                set(obj.gui_h.popupmenu_modelvar2plot,'String',currentfiles.guiname);
                set(obj.gui_h.popupmenu_modelvar2plot,'Value',length(currentfiles.guiname));
                %obj.limofiles_filename = newfile{1};
            else
                set(obj.gui_h.popupmenu_modelvar2plot,'Value',1);
            end
        end
    end
end
% =================================================================
%% =================================================================
function [handles,STUDY] = guibuilder(STUDY,analysis)

% Create list "Level of Analysis" (Using First Subject)
%--------------------------------------------------------------------------
subjnames = STUDY.design(STUDY.currentdesign).cases.value;

listlevel = {'Group (Level 2)'};
for i=1:length(subjnames)
    listlevel{i+1,1} = [subjnames{i} ' (Level 1)'];
end

% Create list for "Measures to Plot" (Using First Subject and first measure available)
%--------------------------------------------------------------------------
measure2plot_list     = {'ERP','Power Spectrum'};
datorica_list         = {'dat', 'ica'};
measure2plot_internal = {'erp','spec'};
datorica              =  find(strcmp(analysis,datorica_list));

for i = 1:length(measure2plot_internal)
    measure_flag{i} = any(cellfun(@(x) ~isempty(x),strfind({STUDY.design(STUDY.currentdesign).limo.datatype},measure2plot_internal{i})));
end
measure_flag = cell2mat(measure_flag);

% Create list "Variable to Plot" (Using First Subject)
%--------------------------------------------------------------------------
if ~logical(sum(measure_flag))
    var2plot_list = {'No Variable has been found'};
    var2plot_indx = 1;
else
    tmp = find(measure_flag);
    measure2plot_indx = tmp(1); % With this we are forcing to pick the first measure available
    [var2plot,tmp,tmp,STUDY] = getmeasures2plot(STUDY,1,measure2plot_indx,datorica,1);
    if ~iscell(var2plot)
        var2plot_list = var2plot.guiname;
    else
        var2plot_list = var2plot;
    end
    % Getting index of  'Condition_effect_1.mat' for default value
    var2plot_indx = find(strcmp(var2plot.guiname,'Condition_effect_1'), 1);
    if isempty(var2plot_indx)
        var2plot_indx = find(strcmp(var2plot.guiname,'Covariate_effect_1'),1);
    end
end

% Create list elect2plot_list
%--------------------------------------------------------------------------
tmpchan         =  struct2cell(STUDY.changrp(1,:));
electoplot_list = ['All Channels';{tmpchan{1,1,:}}'];
electoplot_indx = 1;

% TOOLTIPS texts
%--------------------------------------------------------------------------

ttiptext{1} = ['Group or individual analysis. Individual analysis' 10 ...          % ttiptext_level
    'is the same a single subject analysis.' 10 ...
    'Group analysis is the analysis of a group of subjects.'];
ttiptext{2} = ['EEG measure to plot.' 10 '(ERP/Log Spectrum)'];                    % ttiptext_m2plot
ttiptext{3} = 'Select to plot Original Data or Results from the Linear Model';     % ttiptext_dorresult
ttiptext{4} = ['Selection of experimental conditions, subject groups, other' 10 ...% ttiptext_mvar2plot
    'regressors and contrast. Some options might be shaded out' 10 ...
    'based on the ?Plot Type?selected below.'];
ttiptext{5} = 'Select type of plot';                                                % ttiptext_plottype
ttiptext{6} = 'Select p-value to compute statistic';                                % ttiptext_pval
ttiptext{7} = 'Select the channels to be used in the plot';                         % ttiptext_chan

ttiptext{8} = 'Select Multiple Comparisons Method';
ttiptext{9} = 'Select Computational Approach to Perform Multiple Comparisons';
% Stuff for texts
%--------------------------------------------------------------------------
textnames = {'text_level','text_measure2plot','text_plottype','text_dataorresult','text_modelvar2plot','text_pval','text_elect2plot','text_mcc','text_compmethod'};
textpos   = {[0.0349 0.941 0.333 0.0428],...
    [0.0349 0.832 0.306 0.0592],...
    [0.0349 0.737 0.186 0.0625],...
    [0.0349 0.655 0.422 0.0526],...
    [0.0349 0.562 0.422 0.0526],...
    [0.0349 0.48 0.399 0.0428 ],...
    [0.0349 0.378 0.345 0.0428],...
    [0.0426 0.591 0.705 0.242 ],...
    [0.0426 0.227 0.481 0.197]};
textstring = {'Level of Analysis','Data Measure','Plot Type','Data or Results','Variable to Plot','Threshold (p-value )','Electrodes to Plot','Method','Computational Approach'};
% Stuff for popupmenus
%--------------------------------------------------------------------------
plottype_list    = {'Combined Plot','Scalp Maps Series','Time Course'};
datorresult_list = {'Data';'Linear Model Results'};
mcc_list         = {'None';'Clustering';'TFCE';'Max'};
compmethod_list  = {'None','Default Method'};
xalign1          = 0.450; %0.481;
popupnames = {'popupmenu_level','popupmenu_measure2plot','popupmenu_plottype','popupmenu_dataorresult','popupmenu_modelvar2plot','popupmenu_mcc','popupmenu_compmethod'};
popuppos   = {[xalign1 0.921 0.550 0.0625],...           % level
              [xalign1 0.829 0.550 0.0625],...           % measure2plot
              [xalign1 0.737 0.550 0.0625],...           % plottype
              [xalign1 0.645 0.550 0.0625],...           % dataorresult
              [xalign1+0.08 0.553 0.550-0.08 0.0625],... % modelvar2plot
              [0.535 0.561 0.465 0.273 ],...             % mcc
              [0.535 0.152 0.465 0.273]};                % compmethod
popupstring = {listlevel,measure2plot_list,plottype_list,datorresult_list,var2plot_list,mcc_list,compmethod_list};
popupindex  = {2,measure2plot_indx,1,2,var2plot_indx,1,1};
%-------------------------------------------------------------------------
compmethod_list = {''};                                                    % UPDATE THIS
pval_default    = num2str(0.5);
icadefs;
color           = BACKEEGLABCOLOR;

% Figure
%--------------------------------------------------------------------------
handles.fig = figure('MenuBar','none',...
    'Name','pop_limoresults',...
    'NumberTitle','off',...
    'Units', 'Points',...
    'Color', color,...
    'Position',[549.072 185.686 283.538 466.212],...
    'Resize', 'off');

%Panel 1 (Plot Settings)
%--------------------------------------------------------------
% Panel
handles.panel_1 = uipanel('parent',handles.fig,...
    'units','normalized',...
    'position',[0.0182 0.306 0.956 0.683],...
    'title','Plot Settings',...
    'FontSize',12,...
    'backgroundcolor',color,...
    'FontWeight','bold',...
    'Tag', 'panel_1');

%Panel 2 (Compute Statistics)
%--------------------------------------------------------------------------
handles.panel_2 = uipanel('parent',handles.fig,...
    'units','normalized',...
    'position',[0.0182 0.12 0.956 0.173],...
    'title','Correct for Multiple Comparisons',... % 'Compute Statistics (optional)'
    'backgroundcolor',color,...
    'FontSize',12,...
    'FontWeight','bold',...
    'Tag','panel_2');

% Texts
%--------------------------------------------------------------------------
for i = 1:length(textnames)
    if i > 7 parent = handles.panel_2 ; else parent = handles.panel_1; end

    handles.(textnames{i}) = uicontrol('parent'              ,parent ,...
        'style'               ,'text',...
        'units'               ,'normalized',...
        'position'            ,textpos{i},...
        'string'              ,textstring{i},...
        'FontSize'            ,11,...
        'HorizontalAlignment' ,'left',...
        'backgroundcolor'     ,color,...
        'tooltipString'       ,ttiptext{i},...
        'Tag'                 ,textnames{i});
end

% Popups Menus
%--------------------------------------------------------------------------
ttip_popup = {ttiptext{[1:5,8,9]}};
for i = 1:length(popupnames)
    if i > 5 parent = handles.panel_2; else parent = handles.panel_1; end

    handles.(popupnames{i}) = uicontrol('parent'        ,parent,...
        'style'         ,'popupmenu',...
        'units'         ,'normalized',...
        'position'      ,popuppos{i},...
        'string'        ,popupstring{i},...
        'value'         ,popupindex{i},...
        'Fontsize'      ,10,...
        'tooltipString' ,ttip_popup{i},...
        'Tag'           ,popupnames{i});
end
set(handles.popupmenu_mcc,       'Enable','off');
set(handles.popupmenu_compmethod,'Enable','off');

% Edit
handles.edit_pval = uicontrol('parent',handles.panel_1,...
    'style','edit',...
    'units','normalized',...
    'position',[xalign1+0.02 0.451 0.271 0.0724],...
    'string',pval_default,...
    'backgroundcolor',[0.929 0.929 0.929],...
    'tooltipString',ttiptext{6},...
    'Tag','edit_pval');

% List
handles.listbox_elect2plot = uicontrol('parent',handles.panel_1,...
    'style','listbox',...
    'units','normalized',...
    'position',[xalign1+0.02 0.0263 0.488 0.395],...
    'string',electoplot_list,...
    'Value',electoplot_indx,...
    'backgroundcolor',[0.929 0.929 0.929],...
    'FontSize',11,...
    'enable','off',...
    'tooltipString',ttiptext{7},...
    'Tag','listbox_elect2plot');

%Panel 2 (Compute Statistics)
%--------------------------------------------------------------------------
% Checkbox
handles.checkbox_stats = uicontrol('parent',handles.fig,...
    'style','checkbox',...
    'units','normalized',...
    'position',[0.73 0.253 0.106 0.0493],...
    'string','',...
    'backgroundcolor',color,...
    'Tag','checkbox_stats');

% Button
handles.pushbutton_plot = uicontrol('parent',handles.fig,...
    'style','pushbutton',...
    'units','normalized',...
    'position',[0.113 0.0321 0.766 0.0664],...
    'string','PLOT',...
    'backgroundcolor',[0.929 0.929 0.929],...
    'Tag', 'pushbutton_plot');

% Button add design
handles.pushbutton_designadd = uicontrol('parent',handles.panel_1,...
    'style','pushbutton',...
    'units','normalized',...
    'position',[xalign1+0.02 0.553+0.005 0.07 0.0625 - 0.005],...
    'string','+',...
    'FontWeight','bold',...
    'FontSize', 12,...
    'tooltipString','Add contrast',...
    'backgroundcolor',[0.929 0.929 0.929],...
    'Tag', 'pushbutton_designadd');
end

%% =================================================================
function [var2plot,filespath,limoindx,STUDY] = getmeasures2plot(STUDY,subjN,measureindx,datoricaindx,checknewfile_flag)

% measureindx  is the index to {'erp','spec'}
% datoricaindx is the index to {'dat', 'ica'}

% Init
filespath = [];
limoindx = [];
var2plot.guiname{1}  = 'No Variables Computed';
measure_list         = {'erp','spec'};
datorica_list        = {'dat', 'ica'};
if ~exist('checknewfile_flag','var'), checknewfile_flag = 0; end

requested_datatype = [datorica_list{datoricaindx} measure_list{measureindx}];
measures_computed  = {STUDY.design(STUDY.currentdesign).limo.datatype};

for i = 1:length(measure_list)
    measure_flag{i} = any(cellfun(@(x) ~isempty(x),strfind(measures_computed,measure_list{i})));
end
measure_flag = cell2mat(measure_flag);

% Warning for gui
if ~logical(sum(measure_flag))
    eeglab_warning('No results have been found for this subject');
end
if ismember(requested_datatype,measures_computed)
    limoindx  = find(strcmp(requested_datatype,measures_computed));
    % Generating list
    %----------------
    filespath         = STUDY.design(STUDY.currentdesign).limo(limoindx).basefolder_level1{subjN};

    if checknewfile_flag
        tmp = dir(filespath);
        if ~isempty({tmp.name})
            var2plot_list        = {tmp.name}';

            % Cleaning out the list
            %----------------------
            list2cleanout = {'Betas.mat','LIMO.mat','Yr.mat','Yhat.mat','Res.mat','.','..'};
            for i = 1:length(list2cleanout)
                ind2delete{i} = find(strcmp(var2plot_list,list2cleanout{i}));
            end
            nonemptyvals              = find(cellfun(@(x) ~isempty(x), ind2delete));
            ind2delete                = [ind2delete{nonemptyvals}];
            var2plot_list(ind2delete) = [];

            %Update STUDY

            [newfiles, newfilesindx] = setdiff(cellfun(@(x) fullfile(filespath,x),var2plot_list,'UniformOutput', false),STUDY.design(STUDY.currentdesign).limo(limoindx).model(subjN).filename);
            if ~isempty(newfiles)
                STUDY.design(STUDY.currentdesign).limo(limoindx).model(subjN).guiname(newfilesindx) = cellfun(@(x) x(1:end-4),var2plot_list(newfilesindx),'UniformOutput', false);
                STUDY.design(STUDY.currentdesign).limo(limoindx).model(subjN).filename(newfilesindx) = cellfun(@(x) fullfile(filespath,x),var2plot_list(newfilesindx),'UniformOutput', false);
                STUDY = pop_savestudy( STUDY, [],'filepath', STUDY.filepath,'savemode','resave');
                display('pop_limoresults: New contrast files added to STUDY ');
            end
        end
    end

    var2plot.guiname  = (STUDY.design(STUDY.currentdesign).limo(limoindx).model(subjN).guiname);
    var2plot.pathname = (STUDY.design(STUDY.currentdesign).limo(limoindx).model(subjN).filename);
end
end

% =================================================================
function eeglab_warning(text)

icadefs;
eeglabcolor = BACKEEGLABCOLOR;
h           = warndlg(text,'EEGLAB Warning');
set(h,'color', eeglabcolor);

uiwait(h);
end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function [ generator_new, change_l2, seed ] = cvt_iteration ( m, n, ...
  generator, sample_num_cvt, sample_function_cvt, seed )

%*****************************************************************************80
%
%%  CVT_ITERATION takes one step of the CVT iteration.
%
%  Discussion:
%
%    The routine is given a set of points, called "generators", which
%    define a tessellation of the region into Voronoi cells.  Each point
%    defines a cell.  Each cell, in turn, has a centroid, but it is
%    unlikely that the centroid and the generator coincide.
%
%    Each time this CVT iteration is carried out, an attempt is made
%    to modify the generators in such a way that they are closer and
%    closer to being the centroids of the Voronoi cells they generate.
%
%    A large number of sample points are generated, and the nearest generator
%    is determined.  A count is kept of how many points were nearest to each
%    generator.  Once the sampling is completed, the location of all the
%    generators is adjusted.  This step should decrease the discrepancy
%    between the generators and the centroids.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 February 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of Voronoi cells.
%
%    Input, real GENERATOR(M,N), the Voronoi cell generators.
%
%    Input, integer SAMPLE_NUM_CVT, the number of sample points.
%
%    Input, integer SAMPLE_FUNCTION_CVT, region sampling function:
%    -1, sampling function is RAND (MATLAB intrinsic);
%    0, sampling function is UNIFORM;
%    1, sampling function is HALTON;
%    2, sampling function is GRID;
%
%    Input, integer SEED, the random number seed.
%
%    Output, real GENERATOR_NEW(M,N), the new Voronoi cell generators.
%
%    Output, integer SEED, the new random number seed.
%
%    Output, real CHANGE_L2, the L2 norm of the difference between
%    the input and output data.
%
  generator_new(1:m,1:n) = 0.0;
  tally(1:n) = 0;
  reset = 1;

  for j = 1 : sample_num_cvt
%
%  Generate a sampling point X.
%
    [ x(1:m), seed ] = region_sampler ( m, 1, sample_num_cvt, ...
       sample_function_cvt, reset, seed );
%
%  Ensure that X is a column vector.
%
    x = x(:)

    reset = 0;
%
%  Find the nearest cell generator.
%
    nearest = find_closest ( m, n, x, generator );
%
%  Add X to the averaging data for GENERATOR(*,NEAREST).
%
    for i = 1 : m
      generator_new(i,nearest) = generator_new(i,nearest) + x(i);
    end
    tally(nearest) = tally(nearest) + 1;

  end
%
%  Compute the new generators.
%
  for j = 1 : n
    if ( tally(j) ~= 0 )
      generator_new(1:m,j) = generator_new(1:m,j) / tally(j);
    end
  end
%
%  Determine the change.
%
  change_l2 = 0.0;
  for j = 1 : n
    for i = 1 : m
      change_l2 = change_l2 + ( generator_new(i,j) - generator(i,j) )^2;
    end
  end
  change_l2 = sqrt ( change_l2 );

  return
end


text:
function [S2G, W, M2] = quadratureS2Grid(bandwidth, varargin)
%
% Syntax
%   [S2G, W, M2] = quadratureS2Grid(M) quadrature grid of type chebyshev
%   [S2G, W, M2] = quadratureS2Grid(M, 'gauss') quadrature grid of type gauss
%

persistent S2G_p;
persistent W_p;
persistent M2_p;

if check_option(varargin, 'gauss')
  load(fullfile(mtexDataPath,'vector3d','quadratureS2Grid_gauss.mat'),'gridIndex');
else
  load(fullfile(mtexDataPath,'vector3d','quadratureS2Grid_chebyshev.mat'),'gridIndex');
end
index = find(gridIndex.bandwidth >= bandwidth, 1);
if isempty(index)
  index = size(gridIndex,1);
  warning('M is too large, instead we are giving you the largest quadrature grid we got.');
end

M2 = gridIndex.bandwidth(index);

if ~isempty(M2_p) && M2_p == M2
  S2G = S2G_p;
  W = W_p;
  M2 = M2_p;

else
  name = cell2mat(gridIndex.name(index));
  if check_option(varargin, 'gauss')
    data = load(fullfile(mtexDataPath,'vector3d','quadratureS2Grid_gauss.mat'),name);
  else
    data = load(fullfile(mtexDataPath,'vector3d','quadratureS2Grid_chebyshev.mat'),name);
  end

  data = data.(name);
  S2G = vector3d.byPolar(data(:, 1), data(:, 2));

  if check_option(varargin, 'gauss')
    W = data(:,3);
  else
    W = 4*pi/size(data, 1) .* ones(size(S2G));
  end

  % store the data
  S2G_p = S2G;
  W_p = W;
  M2_p = M2;
end

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
% Signal Processing Toolbox check
% axs Nov 2020
function SigProcT_here = SignalProcessingToolboxCheck

% Get Toolbox list from ver
tbxs = ver;
[a, tbxs_n] = size(tbxs);
tbx_cell = cell(1,tbxs_n);
tbx_cell = {tbxs.Name};

SigProcT_here = any(ismember(tbx_cell,'Signal Processing Toolbox'));

if SigProcT_here == 0
    % if the Signal Processing Toolbox is not found
    warning_text = 'Matlab Signal Processing Toolbox has not been found.';
    warning_text2= 'This is required for many ERPLAB functions, especially filters';
    warning(warning_text)
    warning(warning_text2)
    disp('See Signal Processing Toolbox help here')

end

author: ucdavis

repo: erplab

save_path: github-repos/MATLAB/ucdavis-erplab


text:
function t_miqps_master(quiet)
%T_MIQPS_MASTER  Tests of MILP/MIQP solvers via MIQPS_MASTER().

%   MP-Opt-Model
%   Copyright (c) 2010-2020, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MP-Opt-Model.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://github.com/MATPOWER/mp-opt-model for more info.

if nargin < 1
    quiet = 0;
end

algs = {'DEFAULT', 'CPLEX', 'MOSEK', 'GUROBI', 'GLPK', 'OT'};
names = {'DEFAULT', 'CPLEX', 'MOSEK', 'Gurobi', 'glpk', 'intlin/lin/quadprog'};
check = {@have_miqp_solver, 'cplex', 'mosek', 'gurobi', 'glpk', 'intlinprog'};
does_qp = [0 1 1 1 0 0];
if have_feature('gurobi') || have_feature('cplex') || have_feature('mosek')
    does_qp(1) = 1;
end

n = 48;
nqp = 28;
nmiqp = 6;
t_begin(n*length(algs), quiet);

diff_alg_warn_id = 'optim:linprog:WillRunDiffAlg';
if have_feature('quadprog') && have_feature('quadprog', 'vnum') == 7.005
    s1 = warning('query', diff_alg_warn_id);
    warning('off', diff_alg_warn_id);
end

for k = 1:length(algs)
    if ~isempty(check{k}) && ...
            (ischar(check{k}) && ~have_feature(check{k}) || ...
             isa(check{k}, 'function_handle') && ~check{k}())
        t_skip(n, sprintf('%s not installed', names{k}));
    else
        opt = struct('verbose', 0, 'alg', algs{k});
        mpopt = struct( ...
            'verbose', 0, ...
            'opf', struct( ...
                'violation', 5e-6 ), ...
            'cplex', struct( ...
                'lpmethod', 0, ...
                'qpmethod', 0, ...
                'opt', 0 ), ...
            'mosek', struct( ...
                'lp_alg', 0, ...
                'max_it', 0, ...
                'gap_tol', 0, ...
                'max_time', 0, ...
                'num_threads', 0, ...
                'opt', 0 ) ...
        );
        switch names{k}
            case 'CPLEX'
%               alg = 0;        %% default uses barrier method with NaN bug in lower lim multipliers
                alg = 2;        %% use dual simplex
                mpopt.cplex.lpmethod = alg;
                mpopt.cplex.qpmethod = min([4 alg]);
                opt.cplex_opt = cplex_options([], mpopt);
            case 'MOSEK'
%                 sc = mosek_symbcon;
%                 alg = sc.MSK_OPTIMIZER_DUAL_SIMPLEX;    %% use dual simplex
%                 alg = sc.MSK_OPTIMIZER_INTPNT;          %% use interior point
%                 mpopt.mosek.lp_alg = alg;
                mpopt.mosek.gap_tol = 1e-10;
%                 mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_PFEAS = 1e-10;
%                 mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_DFEAS = 1e-10;
%                 mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_INFEAS = 1e-10;
%                 mpopt.mosek.opts.MSK_DPAR_INTPNT_TOL_REL_GAP = 1e-10;
                vnum = have_feature('mosek', 'vnum');
                if vnum >= 8
%                     mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_PFEAS = 1e-10;
%                     mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_DFEAS = 1e-10;
%                     mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_INFEAS = 1e-10;
%                     mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_MU_RED = 1e-10;
                    mpopt.mosek.opts.MSK_DPAR_INTPNT_QO_TOL_REL_GAP = 1e-10;
                end
%                 opt.verbose = 3;
                opt.mosek_opt = mosek_options([], mpopt);
        end

        t = sprintf('%s - 3-d LP : ', names{k});
        %% based on example from 'doc linprog'
        c = [-5; -4; -6];
        A = [1 -1  1;
             -3  -2  -4;
             3  2  0];
        l = [-Inf; -42; -Inf];
        u = [20; Inf; 30];
        xmin = [0; 0; 0];
        x0 = [];
        [x, f, s, out, lam] = miqps_master([], c, A, l, u, xmin, [], [], [], opt);
        t_is(s, 1, 12, [t 'success']);
        t_is(x, [0; 15; 3], 6, [t 'x']);
        t_is(f, -78, 6, [t 'f']);
        t_is(lam.mu_l, [0;1.5;0], 9, [t 'lam.mu_l']);
        t_is(lam.mu_u, [0;0;0.5], 9, [t 'lam.mu_u']);
        t_is(lam.lower, [1;0;0], 9, [t 'lam.lower']);
        t_is(lam.upper, zeros(size(x)), 9, [t 'lam.upper']);

        if does_qp(k)
            t = sprintf('%s - unconstrained 3-d quadratic : ', names{k});
            %% from http://www.akiti.ca/QuadProgEx0Constr.html
            H = [5 -2 -1; -2 4 3; -1 3 5];
            c = [2; -35; -47];
            x0 = [0; 0; 0];
            [x, f, s, out, lam] = miqps_master(H, c, [], [], [], [], [], [], [], opt);
            t_is(s, 1, 12, [t 'success']);
            t_is(x, [3; 5; 7], 7, [t 'x']);
            t_is(f, -249, 11, [t 'f']);
            t_ok(isempty(lam.mu_l), [t 'lam.mu_l']);
            t_ok(isempty(lam.mu_u), [t 'lam.mu_u']);
            t_is(lam.lower, zeros(size(x)), 13, [t 'lam.lower']);
            t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);

            t = sprintf('%s - constrained 2-d QP : ', names{k});
            %% example from 'doc quadprog'
            H = [   1   -1;
                    -1  2   ];
            c = [-2; -6];
            A = [   1   1;
                    -1  2;
                    2   1   ];
            l = [];
            u = [2; 2; 3];
            xmin = [0; 0];
            x0 = [];
            [x, f, s, out, lam] = miqps_master(H, c, A, l, u, xmin, [], x0, [], opt);
            t_is(s, 1, 12, [t 'success']);
            t_is(x, [2; 4]/3, 7, [t 'x']);
            t_is(f, -74/9, 6, [t 'f']);
            t_is(lam.mu_l, [0;0;0], 13, [t 'lam.mu_l']);
            t_is(lam.mu_u, [28;4;0]/9, 4, [t 'lam.mu_u']);
            t_is(lam.lower, zeros(size(x)), 7, [t 'lam.lower']);
            t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);

            t = sprintf('%s - constrained 4-d QP : ', names{k});
            %% from https://v8doc.sas.com/sashtml/iml/chap8/sect12.htm
            H = [   1003.1  4.3     6.3     5.9;
                    4.3     2.2     2.1     3.9;
                    6.3     2.1     3.5     4.8;
                    5.9     3.9     4.8     10  ];
            c = zeros(4,1);
            A = [   1       1       1       1;
                    0.17    0.11    0.10    0.18    ];
            l = [1; 0.10];
            u = [1; Inf];
            xmin = zeros(4,1);
            x0 = [1; 0; 0; 1];
            [x, f, s, out, lam] = miqps_master(H, c, A, l, u, xmin, [], x0, [], opt);
            t_is(s, 1, 12, [t 'success']);
            t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
            t_is(f, 3.29/3, 6, [t 'f']);
            t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
            t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
            t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
            t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);

            t = sprintf('%s - (struct) constrained 4-d QP : ', names{k});
            p = struct('H', H, 'A', A, 'l', l, 'u', u, 'xmin', xmin, 'x0', x0, 'opt', opt);
            [x, f, s, out, lam] = miqps_master(p);
            t_is(s, 1, 12, [t 'success']);
            t_is(x, [0; 2.8; 0.2; 0]/3, 5, [t 'x']);
            t_is(f, 3.29/3, 6, [t 'f']);
            t_is(lam.mu_l, [6.58;0]/3, 6, [t 'lam.mu_l']);
            t_is(lam.mu_u, [0;0], 13, [t 'lam.mu_u']);
            t_is(lam.lower, [2.24;0;0;1.7667], 4, [t 'lam.lower']);
            t_is(lam.upper, zeros(size(x)), 13, [t 'lam.upper']);
        else
            t_skip(nqp, sprintf('%s does not handle QP problems', names{k}));
        end

        t = sprintf('%s - infeasible LP : ', names{k});
        p = struct('A', sparse([1 1]), 'c', [1;1], 'u', -1, 'xmin', [0;0], 'opt', opt);
        [x, f, s, out, lam] = miqps_master(p);
        t_ok(s <= 0, [t 'no success']);

% opt.verbose = 3;
        t = sprintf('%s - 2-d MILP : ', names{k});
        %% from MOSEK 6.0 Guided Tour, section  7.13.1, https://docs.mosek.com/6.0/toolbox/node009.html
        c = [-2; -3];
        A = sparse([195 273; 4 40]);
        u = [1365; 140];
        xmax = [4; Inf];
        vtype = 'I';
        p = struct('c', c, 'A', A, 'u', u, 'xmax', xmax, 'vtype', vtype, 'opt', opt);
        [x, f, s, out, lam] = miqps_master(p);
        t_is(s, 1, 12, [t 'success']);
        t_is(x, [4; 2], 12, [t 'x']);
        t_is(lam.mu_l, [0; 0], 12, [t 'lam.mu_l']);
        t_is(lam.mu_u, [0; 0], 12, [t 'lam.mu_u']);
        t_is(lam.lower, [0; 0], 12, [t 'lam.lower']);
        t_is(lam.upper, [2; 3], 12, [t 'lam.upper']);

        if does_qp(k)
            t = sprintf('%s - 4-d MIQP : ', names{k});
            %% from cplexmiqpex.m, CPLEX_Studio_Academic124/cplex/examples/src/matlab/cplexmiqpex.m
            H = sparse([ 33   6    0    0;
                          6  22   11.5  0;
                          0  11.5 11    0;
                          0   0    0    0]);
            c = [-1 -2 -3 -1]';
            Aineq = [-1  1  1 10;
               1 -3  1  0];
            bineq = [20  30]';
            Aeq   = [0  1  0 -3.5];
            beq   =  0;
            xmin    = [ 0;   0;   0; 2];
            xmax    = [40; Inf; Inf; 3];
            A = sparse([Aeq; Aineq]);
            l = [beq; -Inf; -Inf];
            u = [beq; bineq];
            vtype = 'CCCI';
            p = struct('H', H, 'c', c, 'A', A, 'l', l, 'u', u, ...
                'xmin', xmin, 'xmax', xmax, 'vtype', vtype, 'opt', opt);
            [x, f, s, out, lam] = miqps_master(p);
            t_is(s, 1, 12, [t 'success']);
            t_is(x, [7; 7; 0; 2], 7, [t 'x']);
            t_is(lam.mu_l, [466; 0; 0], 6, [t 'lam.mu_l']);
            t_is(lam.mu_u, [0; 272; 0], 6, [t 'lam.mu_u']);
            t_is(lam.lower, [0; 0; 349.5; 4350], 5, [t 'lam.lower']);
            t_is(lam.upper, [0; 0; 0; 0], 7, [t 'lam.upper']);
        else
            t_skip(nmiqp, sprintf('%s does not handle MIQP problems', names{k}));
        end
% opt.verbose = 0;
    end
end

if have_feature('quadprog') && have_feature('quadprog', 'vnum') == 7.005
    warning(s1.state, diff_alg_warn_id);
end

t_end;

function TorF = have_miqp_solver()
TorF = have_feature('cplex') || have_feature('glpk') || ...
    have_feature('gurobi') || have_feature('intlinprog') || ...
    have_feature('mosek');

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function msm_to_mm_coordinate ( output_filename, a, type, symmetry )

%*****************************************************************************80
%
%% MSM_TO_MM_COORDINATE writes a MATLAB Sparse Matrix to a Matrix Market file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    31 October 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the name of the file to which the information
%    should be written.
%
%    Input, sparse matrix A, the NROW by NCOL matrix, stored in MATLAB sparse
%    matrix format, which is to be written to the file.
%
%    Input, string TYPE, describes the arithmetic data type, as
%    'REAL', 'INTEGER', 'COMPLEX', or 'PATTERN'.
%
%    Input, string SYMMETRY, describes the symmetry type, as 'GENERAL',
%    'SYMMETRIC', 'SKEW-SYMMETRIC' or 'HERMITIAN'.
%
  if ( s_eqi ( type, 'REAL' ) )
    msm_to_mm_coordinate_real ( output_filename, a, symmetry );
  elseif ( s_eqi ( type, 'INTEGER' ) )
    msm_to_mm_coordinate_integer ( output_filename, a, symmetry );
  elseif ( s_eqi ( type, 'COMPLEX' ) )
    msm_to_mm_coordinate_complex ( output_filename, a, symmetry );
  elseif ( s_eqi ( type, 'PATTERN' ) )
    msm_to_mm_coordinate_pattern ( output_filename, a, symmetry );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MSM_TO_MM_COORDINATE - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of TYPE.\n' );
    error ( 'MSM_TO_MM_COORDINATE - Fatal error!' );
  end

  return
end


text:
%  thetaj : switching parameter
%  thetaj_bar: steady state value
%  pj : 1 if the parameter is "partitioned"
%  iota1 : 0 if FRWZ perturbation
%  sig : perturbation parameter
%

author: jmaih


text:
classdef MaxPropTimeReachedValidator < AbstractLaunchVehicleDataValidator
    %MaxPropTimeReachedValidator Summary of this class goes here
    %   Detailed explanation goes here

    properties
        lvdData LvdData
    end

    methods
        function obj = MaxPropTimeReachedValidator(lvdData)
            obj.lvdData = lvdData;
        end

        function [errors, warnings] = validate(obj)
            errors = LaunchVehicleDataValidationError.empty(0,1);
            warnings = LaunchVehicleDataValidationWarning.empty(0,1);

            maxPropTime = obj.lvdData.settings.maxScriptPropTime;
            lastRunExecTime = obj.lvdData.script.lastRunExecTime;

            if(lastRunExecTime >= maxPropTime)
                str = sprintf('Maximum script propagation time of %.3f sec reached or exceeded.  Propagation terminated.  Consider increasing the maximum script propagation time.', maxPropTime);
                warnings(end+1) = LaunchVehicleDataValidationWarning(str);
            end
        end
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [rsmooth,gsmooth,bsmooth,stdim] = flyproc(im)

% This function contains the same algorithm as the demo script flyexdemo.m.
% All graphical output steps have been removed, and the script has been
% converted into a function which takes as input the filename of a
% Drosophila image, and outputs the rotated and cropped image, along with
% curves fit to the red, green and blue channels. It is used as a utility
% function by localflyexdemo.m and dctflyexdemo.m.
%
% This function uses the Image Processing ans Curve Fitting Toolboxes.
%
% Sam Roberts

%   Copyright 2006 The MathWorks, Inc.

im = imread(im);

RGBmax = max(im,[],3);

mask = im2bw(RGBmax,20/255);

mask = medfilt2(mask,[3,3]);

mask = imclose(mask,strel('disk',5,0));

L = bwlabel(mask);

stats = regionprops(L,...
    {'Orientation','MajorAxisLength','MinorAxisLength'});

stdim = imrotate(im,-stats.Orientation);
warning off MATLAB:colon:nonIntegerIndex
stdim = stdim(size(stdim,1)/2-stats.MinorAxisLength/2:size(stdim,1)/2+stats.MinorAxisLength/2,...
    size(stdim,2)/2-stats.MajorAxisLength/2:size(stdim,2)/2+stats.MajorAxisLength/2,:);
warning on MATLAB:colon:nonIntegerIndex

RGBmax = max(stdim,[],3);

RGBequal = adapthisteq(RGBmax,'NumTiles',...
    [ceil(size(RGBmax,1)/13),ceil(size(RGBmax,2)/13)]);

RGBequal = medfilt2(RGBequal);

bw=im2bw(RGBequal,105/255);
[L,num] = bwlabel(bw,4);

stats = regionprops(L,'Centroid','PixelIdxList');

r = squeeze(stdim(:,:,1));
g = squeeze(stdim(:,:,2));
b = squeeze(stdim(:,:,3));
[rows,columns,tmp]=size(stdim);
for i=1:num
    data(i,1:2) = stats(i).Centroid./[columns,rows]*100;
    data(i,3)=mean(r(stats(i).PixelIdxList));
    data(i,4)=mean(g(stats(i).PixelIdxList));
    data(i,5)=mean(b(stats(i).PixelIdxList));
end

middlestrip = (data(:,2)>40 & data(:,2)<60);
middledata = data(middlestrip,:);

opts = fitoptions('smoothingspline','SmoothingParam',0.01);
fitresultr = fit(middledata(:,1),middledata(:,3),'smoothingspline',opts);
fitresultg = fit(middledata(:,1),middledata(:,4),'smoothingspline',opts);
fitresultb = fit(middledata(:,1),middledata(:,5),'smoothingspline',opts);

rsmooth = feval(fitresultr,0:100);
gsmooth = feval(fitresultg,0:100);
bsmooth = feval(fitresultb,0:100);

author: Sable

repo: mcbench-benchmarks


text:
clear
clc
close all

addpath('data')
addpath('src')
dataset = {'1_mECS', '2_Kolod', '3_Pollen', '4_Usoskin'}

for i = 1:4

    % perform the analysis for the current dataset
    load(['Test_' dataset{i}]);
    C = max(true_labs); %%% number of clusters
    rng(i,'twister'); %%% for reproducibility
    [y, S, F, ydata,alpha] = SIMLR(in_X,C,10);

    % report NMI values
    NMI_i = Cal_NMI(y,true_labs);
    fprintf(['The NMI value for dataset ' dataset{i} ' is %f\n'], NMI_i);

    % visualization
    figure;
    gscatter(ydata(:,1),ydata(:,2),true_labs);

end

author: BatzoglouLabSU

repo: SIMLR


text:
function [n1,val]=zerotest(y,m,w)
[m,n1]=min(abs(y-m));
if sum(w)==0
val=1;
else
val=w(n1);
end


text:
function h = hinton(w);
%HINTON	Plot Hinton diagram for a weight matrix.
%
%	Description
%
%	HINTON(W) takes a matrix W and plots the Hinton diagram.
%
%	H = HINTON(NET) also returns the figure handle H which can be used,
%	for instance, to delete the  figure when it is no longer needed.
%
%	To print the figure correctly in black and white, you should call
%	SET(H, 'INVERTHARDCOPY', 'OFF') before printing.
%
%	See also
%	DEMHINT, HINTMAT, MLPHINT
%

%	Copyright (c) Ian T Nabney (1996-2001)

% Set scale to be up to 0.9 of maximum absolute weight value, where scale
% defined so that area of box proportional to weight value.

% Use no more than 640x480 pixels
xmax = 640; ymax = 480;

% Offset bottom left hand corner
x01 = 40; y01 = 40;
x02 = 80; y02 = 80;

% Need to allow 5 pixels border for window frame: but 30 at top
border = 5;
top_border = 30;

ymax = ymax - top_border;
xmax = xmax - border;

% First layer

[xvals, yvals, color] = hintmat(w);
% Try to preserve aspect ratio approximately
if (8*size(w, 1) < 6*size(w, 2))
  delx = xmax; dely = xmax*size(w, 1)/(size(w, 2));
else
  delx = ymax*size(w, 2)/size(w, 1); dely = ymax;
end

h = figure('Color', [0.5 0.5 0.5], ...
  'Name', 'Hinton diagram', ...
  'NumberTitle', 'off', ...
  'Colormap', [0 0 0; 1 1 1], ...
  'Units', 'pixels', ...
  'Position', [x01 y01 delx dely]);
set(gca, 'Visible', 'off', 'Position', [0 0 1 1]);
hold on
patch(xvals', yvals', color', 'Edgecolor', 'none');
axis equal;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function chebyshev_polynomial_test01 ( )

%*****************************************************************************80
%
%% CHEBYSHEV_POLYNOMIAL_TEST01 tests T_PROJECT_COEFFICIENTS_DATA.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2012
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CHEBYSHEV_POLYNOMIAL_TEST01:\n' );
  fprintf ( 1, '  T_PROJECT_COEFFICIENTS_DATA estimates the Chebyshev polynomial\n' );
  fprintf ( 1, '  coefficients for a function given as data (x,fx).\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here, we use fx = f(x) = x^2 for the data.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Since T(0,x) = 1 and T(2,x) = 2*x^2 - 1, the correct expansion is\n' );
  fprintf ( 1, '  f(x) = 1/2 T(0,x) + 0 T(1,x) + 1/2 T(2,x) + 0 * all other polys.\n' );
%
%  Data in [0,1];
%
  a = 0.0;
  b = 1.0;
  m = 20;
  seed = 123456789;
  [ x, seed ] = r8vec_uniform_01 ( m, seed );
  d = x.^2;

  r8vec2_print ( m, x, d, '  Data ( X, D ):' );

  n = 4;
  c = t_project_coefficients_data ( a, b, m, n, x, d );

  r8vec_print ( n, c, '  Coefficients of Chebyshev expansion of degree 4.' );
%
%  Compare Chebyshev expansion and original function.
%
  d2 = t_project_value ( m, n, x, c );

  fprintf ( 1, '\n' );
  fprintf ( 1, '   I      X(I)     Data(I)      Chebyshev(X(I))\n' );
  fprintf ( 1, '\n' );
  for i = 1 : m
    fprintf ( 1, '  %2d  %12g  %12g  %12g\n', i, x(i), d(i), d2(i) );
  end

  return
end