text:
%% rrt_star_fn
% *RRT*FN* is a memory efficient variant of RRT*. It inherents all
% the optimization features of RRT*, but in contrast using limited
% number of nodes i.e. memory is limited.
%
%% Syntax
%  problem = rrt_star(map, max_iter, max_nodes, is_benchmark, rand_seed, variant)
%  function returns the object of the respective class with the result
%
%% Input
%  map           -- struct with appropriate fields (developer of
%                the class provides more information on this topic)
%  max_iter      -- number of iteration to solve the problem
%  max_nodes     -- the maximum allowed number of nodes added to the tree
%  is_benchmark  -- if true saves snapshots of the tree in a special directory
%                boolean variable
%  rand_seed     -- a random seed
%  variant       -- what class to choose, class used defines the problem space
%% Output
% *problem* is an object of an appropriate problem space representing class
%

author: olzhas


text:
function b = isPointOnLine(point, line, varargin)
%ISPOINTONLINE Test if a point belongs to a line.
%
%   B = isPointOnLine(POINT, LINE)
%   with POINT being [xp yp], and LINE being [x0 y0 dx dy].
%   Returns 1 if point lies on the line, 0 otherwise.
%
%   If POINT is an N-by-2 array of points, B is a N-by-1 array of booleans.
%
%   If LINE is a N-by-4 array of line, B is a 1-by-N array of booleans.
%
%   B = isPointOnLine(POINT, LINE, TOL)
%   Specifies the tolerance used for testing location on 3D line. Default value is 1e-14.
%
%   See also
%   lines2d, points2d, isPointOnEdge, isPointOnRay, isLeftOriented
%

% ------
% Author: David Legland
% E-mail: david.legland@inrae.fr
% Created: 2003-10-31
% Copyright 2003-2022 INRA - TPV URPOI - BIA IMASTE

% extract computation tolerance
tol = 1e-14;
if ~isempty(varargin)
    tol = varargin{1};
end

% test if lines are colinear, using third coordinate of 3D cross-product
% same test as:
% b = abs((xp-x0).*dy-(yp-y0).*dx)./hypot(dx, dy).^2 < tol;
b = bsxfun(...
    @rdivide, abs(...
    bsxfun(@times, bsxfun(@minus, point(:,1), line(:,1)'), line(:,4)') - ...
    bsxfun(@times, bsxfun(@minus, point(:,2), line(:,2)'), line(:,3)')), ...
    (line(:,3).^2 + line(:,4).^2)') < tol;

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function tt=simular(robot,t,y,mode)

    global radio
    global precision

    drawrobot3d(robot,[0 0 0 0 0 0])
    axis([-1.5 1.50 -1.5 0.75 -0.25 1.75]);
    adjust_view()
    f=length(t);
    d=1;

    pause(3)
    %Recorrido desde posicion inicial hasta inicio pintado
     tin=inicial(y,t,robot);
     in=length(tin);
     for i=1:in;
        T=tin{i};
        q = inversekinematic_kuka_kr160_r1570_nanoC(robot, T) ;
        drawrobot3d(robot, q)

        axis([-1.5 1.50 -1.5 0.75 -0.25 1.75]);
        hold on

         treal=directkinematic(robot,q(:,1));
         pin(:,i)=treal(1:3,4);
         tt{d}=treal;
         d=d+1;

         Representacion_en_esfera(y);
         hold on
        axis([-1.5 1.50 -1.5 0.75 -0.25 1.75]);
         switch mode
             case 'trayectoria'
                 plot3(pin(1,:),pin(2,:),pin(3,:));

             case 'letras'

         end
       % pause(0.005);
    end

    %Pintado esfera
    let=1;
    frase=[0; 0; 0];
    for i=1:f;
        T=y*t{i};
        q = inversekinematic_kuka_kr160_r1570_nanoC(robot, T) ;
        drawrobot3d(robot, q)

        axis([-1.5 1.50 -1.5 0.75 -0.25 1.75]);
        hold on
        treal=directkinematic(robot,q(:,1));
        tt{d}=treal;
        d=d+1;

            if abs((treal(1,4)-y(1,4))^2+(treal(2,4)-y(2,4))^2+(treal(3,4)-y(3,4))^2-radio^2) < precision
                frase(:,let)=treal(1:3,4);
                let=let+1;
            end
            p(:,i)=treal(1:3,4);

        %end
         Representacion_en_esfera(y);
         hold on
        axis([-1.5 1.50 -1.5 0.75 -0.25 1.75]);
         switch mode
             case 'trayectoria'
                  plot3(p(1,:),p(2,:),p(3,:));
                  plot3(pin(1,:),pin(2,:),pin(3,:));
             case 'letras'
                  plot3(frase(1,:),frase(2,:),frase(3,:));
         end

       % pause(0.005);
    end

    %Recorrido desde posicion final del pintado hasta posicion inicial
     tfin=final(y,t,robot);
     fin=length(tfin);
     for i=1:fin;
        T=tfin{i};
        q = inversekinematic_kuka_kr160_r1570_nanoC(robot, T) ;
        drawrobot3d(robot, q)

        axis([-1.5 1.50 -1.5 0.75 -0.25 1.75]);
        hold on
        treal=directkinematic(robot,q(:,1));
        tt{d}=treal;
        d=d+1;

            pfin(:,i)=treal(1:3,4);
            %axis([-1 1 -1 1]);

       % end
        Representacion_en_esfera(y);
        hold on
        axis([-1.5 1.50 -1.5 0.75 -0.25 1.75]);
         switch mode
             case 'trayectoria'
                  plot3(p(1,:),p(2,:),p(3,:));
                  plot3(pin(1,:),pin(2,:),pin(3,:));
                  plot3(pfin(1,:),pfin(2,:),pfin(3,:));
             case 'letras'
                 plot3(frase(1,:),frase(2,:),frase(3,:),'r');

         end

       % pause(0.005);
    end

    drawrobot3d(robot, [0 0 0 0 0 0])
    hold on

    Representacion_en_esfera(y);

        axis([-1.5 1.50 -1.5 0.75 -0.25 1.75]);
         switch mode
             case 'trayectoria'
                  plot3(p(1,:),p(2,:),p(3,:));
                  plot3(pin(1,:),pin(2,:),pin(3,:));
                  plot3(pfin(1,:),pfin(2,:),pfin(3,:));
             case 'letras'
                 %figure
                  %Representacion_en_esfera(y);
                 %hold on
                  plot3(frase(1,:),frase(2,:),frase(3,:),'r');

         end
end

repo: ARTE


text:
function v = tfocs_size( x )

% SIZE   TFOCS-friendly size operator.

v = { @zeros, m, n };

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function wishart_test06 ( )

%*****************************************************************************80
%
%% WISHART_TEST06 compares the Wishart and Bartlett sample matrices.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 July 2013
%
%  Author:
%
%    John Burkardt
%

%
%  Access the PDFLIB and RNGLIB libraries.
%
  addpath ( '../pdflib' );
  addpath ( '../rnglib' );
%
%  Initialize the RNGLIB library.
%  Normally, we would do this just once, here at the beginning.
%  In this example, however, we really want to do it just before
%  we call each of the sampling routines, so that they both access
%  the same set of random numbers...
%
  initialize ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'WISHART_TEST06:\n' );
  fprintf ( 1, '  Verify that, if using the same set of random numbers,\n' );
  fprintf ( 1, '    W = T'' * T,\n' );
  fprintf ( 1, '  where\n' );
  fprintf ( 1, '    W = wishart_sample ( n, df, sigma );\n' );
  fprintf ( 1, '    T = bartlett_sample ( n, df, sigma );\n' );
%
%  Set the parameters.
%
  n = 3;
  df = 5;
  r = [ 5.0, 1.0, 3.0; ...
        0.0, 4.0, 2.0; ...
        0.0, 0.0, 6.0 ];
  sigma = r' * r;
  r8mat_print ( n, n, sigma, '  Covariance SIGMA:' );
%
%  Initialize the random number package and compute W.
%
  initialize ( );
  w = wishart_sample ( n, df, sigma );
%
%  Initialize the random number package again, and compute T.
%
  initialize ( );
  t = bartlett_sample ( n, df, sigma );
%
%  Compute T' * T and compare it to W.
%
  tt = t' * t;

  diff = r8mat_norm_fro_affine ( n, n, w, tt );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Frobenius norm of error is %g\n', diff );
%
%  Release the PDFLIB and RNGLIB libraries.
%
  rmpath ( '../pdflib' );
  rmpath ( '../rnglib' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
#!/usr/bin/env octave
%% Machine Learning Online Class - Exercise 2: Logistic Regression
%
%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the second part
%  of the exercise which covers regularization with logistic regression.
%
%  You will need to complete the following functions in this exericse:
%
%     sigmoid.m
%     costFunction.m
%     predict.m
%     costFunctionReg.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% Load Data
%  The first two columns contains the exam scores and the third column
%  contains the label.

data = load('ex2data2.txt');
X = data(:, [1, 2]); y = data(:, 3);

plotData(X, y);

% Put some labels
hold on;

% Labels and Legend
xlabel('Microchip Test 1')
ylabel('Microchip Test 2')

% Specified in plot order
legend('y = 1', 'y = 0')
hold off;

%% =========== Part 1: Regularized Logistic Regression ============
%  In this part, you are given a dataset with data points that are not
%  linearly separable. However, you would still like to use logistic
%  regression to classify the data points.
%
%  To do so, you introduce more features to use -- in particular, you add
%  polynomial features to our data matrix (similar to polynomial
%  regression).
%

% Add Polynomial Features

% Note that mapFeature also adds a column of ones for us, so the intercept
% term is handled
X = mapFeature(X(:,1), X(:,2));

% Initialize fitting parameters
initial_theta = zeros(size(X, 2), 1);

% Set regularization parameter lambda to 1
lambda = 1;

% Compute and display initial cost and gradient for regularized logistic
% regression
[cost, grad] = costFunctionReg(initial_theta, X, y, lambda);

fprintf('Cost at initial theta (zeros): %f\n', cost);

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ============= Part 2: Regularization and Accuracies =============
%  Optional Exercise:
%  In this part, you will get to try different values of lambda and
%  see how regularization affects the decision coundart
%
%  Try the following values of lambda (0, 1, 10, 100).
%
%  How does the decision boundary change when you vary lambda? How does
%  the training set accuracy vary?
%

% Initialize fitting parameters
initial_theta = zeros(size(X, 2), 1);

% Set regularization parameter lambda to 1 (you should vary this)
lambda = 1;

% Set Options
options = optimset('GradObj', 'on', 'MaxIter', 400);

% Optimize
[theta, J, exit_flag] = ...
	fminunc(@(t)(costFunctionReg(t, X, y, lambda)), initial_theta, options);

% Plot Boundary
plotDecisionBoundary(theta, X, y);
hold on;
title(sprintf('lambda = %g', lambda))

% Labels and Legend
xlabel('Microchip Test 1')
ylabel('Microchip Test 2')

legend('y = 1', 'y = 0', 'Decision boundary')
hold off;

% Compute accuracy on our training set
p = predict(theta, X);

fprintf('Train Accuracy: %f\n', mean(double(p == y)) * 100);

pause;

author: SaveTheRbtz

repo: ml-class


text:
function f = createConvexHullMethod(p,xL,xU)

if ~isa(p.properties.derivative,'function_handle')
    f = [];
    return
end

if isa(p.properties.convexity,'char') && ~(isequal(p.properties.convexity,'none'))
    vexity = p.properties.convexity;
elseif isa(p.properties.convexity,'function_handle')
    % User-supplied method to derive convexity in region
    vexity = p.properties.convexity(xL,xU);
elseif ~isempty(p.properties.inflection)
    % Derive convexity by information about inflection
    vexity = DeriveVexityFromInflection(p.properties,xL,xU);
else
    vexity = 'none';
end

if isequal(vexity,'convex')
    if strcmpi(p.fcn,'blackbox')
        f0 = @(x)real(blackbox(x,p.arg{2}));
    else
        f0 = @(x)real(eval([p.fcn '(x)']));
    end
    f = @(xL,xU)createConvexHullMethodConvex(xL,xU,f0,p.properties.derivative);
elseif isequal(vexity,'concave')
    f0 = @(x)real(eval([p.fcn '(x)']));
    f = @(xL,xU)createConvexHullMethodConcave(xL,xU,f0,p.properties.derivative);
else
    f = [];
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function L = probitNoiseLikelihood(noise, mu, varsigma, y)

% PROBITNOISELIKELIHOOD Likelihood of the data under the PROBIT noise model.
% FORMAT
% DESC returns the likelihood of a data set under the  probit based classification noise model.
% ARG noise : the noise structure for which the likelihood is required.
% ARG mu : input mean locations for the likelihood.
% ARG varSigma : input variance locations for the likelihood.
% ARG y : target locations for the likelihood.
%
% SEEALSO : probitNoiseParamInit, probitNoiseLogLikelihood, noiseLikelihood
%
% COPYRIGHT : Neil D. Lawrence, 2004, 2005

% NOISE

D = size(y, 2);
for i = 1:D
  mu(:, i) = mu(:, i) + noise.bias(i);
end
L = cumGaussian((y.*mu)./(sqrt(noise.sigma2+varsigma)));

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function weight = gradAscent
%%
%
%   Description : LogisticRegression using gradAsscent
%   Author : Liulongpo
%   Time：2015-4-18 10:57:25
%
%%
clc
close all
clear
%%

data = load('testSet.txt');
[row , col] = size(data);
dataMat = data(:,1:col-1);
dataMat = [ones(row,1) dataMat] ;
labelMat = data(:,col);
alpha = 0.001;
maxCycle = 500;
weight = ones(col,1);
for i = 1:maxCycle
    h = sigmoid((dataMat * weight)');
    error = (labelMat - h');
    weight = weight + alpha * dataMat' * error;
end

figure
scatter(dataMat(find(labelMat(:) == 0),2),dataMat(find(labelMat(:) == 0),3),3);
hold on
scatter(dataMat(find(labelMat(:) == 1),2),dataMat(find(labelMat(:) == 1),3),5);
hold on
x = -3:0.1:3;
y = (-weight(1)-weight(2)*x)/weight(3);
plot(x,y)
hold off

end

function returnVals = sigmoid(inX)
    % 注意这里的sigmoid函数要用点除
    returnVals = 1.0./(1.0+exp(-inX));
end

repo: MachineLearning


text:
% convert mel frequency to linear frequency

function [linear_freq] = mel2linear(mel_freq)

% slow implementation
% for i=1:length(mel_freq)
%     linear_freq(i) = 700*( 10^(mel_freq(i)/2595)-1 );
% end

% fast implementation
linear_freq = 700*( 10.^(mel_freq/2595)-1 );

author: singaxiong

repo: SignalGraph


text:
function [Btu] = kJ2Btu(kJ)
% Convert energy or work from kilojoules to British thermal units.
% Chad A. Greene 2012
Btu = kJ*0.94781707775;

author: Sable

repo: mcbench-benchmarks


text:
function r = showrateh(h,err,k,opt,str)
%% SHOWRATEH rate of an err sequence
%
%  r = SHOWRATEH(N,err) finds the number r such that err = N^r and plots the
%  err vs N in loglog scale.
%
%  r = SHOWRATEH(N,err,k) finds the number r such that err(k:end)=N(k:end)^r.
%
%  The function accepts standard plotting properting. For example, r =
%  showrate(N,err,[],'r') will plot the error curve in red.
%
% See also showrate2, showresult, showmesh, showsolution
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

N = 1./h;
if (nargin<=2)
    k = 1; opt = '-*';
end
r = showrate(N,err,k,opt);
h_legend = legend(str,['C_1h^{' num2str(-r) '}'],'LOCATION','Best');
set(h_legend,'FontSize',12);
xlabel('log(1/h)');

repo: Aero-matlab


text:
function [w_upd, P_upd, mean_debug, cov_debug] = learnRLSModel(v, y_meas, w_old, P_old, acc_control_learning)

%__________________________________________________________________________
%% Documentation
%
% Authors:      Alexander Wischnewski (alexander.wischnewski@tum.de)
%
%
% Description:
%   learns a linaer model for the features inserted
%
% Inputs:
%   v:                      feature vector
%   y_meas:                 measured value
%   w_old:                  previous iteration weights
%   P_old:                  previous iteration covariance matrix
%   acc_control_learning:   Parameter structure for the learning algorithm
%
%
% Outputs:

% parameters
Q = 1e-8*eye(length(w_old));
R = 2;

% prediction step
y_pred = v'*w_old;
res = y_meas - y_pred;
P_pred = P_old + Q;

% update step
S = v'*P_pred*v + R;
K = P_pred*v/S;
w_upd = w_old + K*res;
P_upd = (eye(length(w_old)) - K*v')*P_pred;

% sample model and covariance at basis function points for debugging purposes
mean_debug = zeros(length(acc_control_learning.bf_vx_mps), 1);
cov_debug = zeros(length(acc_control_learning.bf_vx_mps), 1);
% for i = 1:1:length(acc_control_learning.bf_vx_mps)
%     dist_ay_feat = ((acc_control_learning.bf_vx_mps' - acc_control_learning.bf_vx_mps(i))/acc_control_learning.vx_width_mps).^2 + ...
%         ((acc_control_learning.bf_ay_req_mps2' - acc_control_learning.bf_ay_req_mps2(i))/acc_control_learning.ay_width_mps2).^2;
%     ay_feat = [exp(-dist_ay_feat), acc_control_learning.bf_ay_req_mps2(i)];
%     mean_debug(i) = ay_feat*w_upd;
%     cov_debug(i) = sqrt(ay_feat*P_upd*ay_feat');
% end

end

author: TUMFTM


text:
function s = QuinticTimeScaling(Tf, t)
% *** CHAPTER 9: TRAJECTORY GENERATION ***
% Takes Tf: Total time of the motion in seconds from rest to rest,
%       t: The current time t satisfying 0 < t < Tf.
% Returns s: The path parameter s(t) corresponding to a fifth-order
%            polynomial motion that begins and ends at zero velocity and
%            zero acceleration.
% Example Input:
%
% clear; clc;
% Tf = 2;
% t = 0.6;
% s = QuinticTimeScaling(Tf,t)
%
% Output:
% s =
%    0.1631

s = 10 * (t / Tf) ^ 3 - 15 * (t / Tf) ^ 4 + 6 * (t / Tf) ^ 5;
end

author: ShuoYangRobotics

repo: QuadrupedSim


text:
function x = pwPoly4(tGrid,xGrid,dxGrid,t)
% x = pwPoly4(tGrid,xGrid,t)
%
% This function does piece-wise quadratic interpolation of a set of data.
%
% INPUTS:
%   tGrid = [1, 2*n-1] = time at each grid point
%   xGrid = [m, 2*n-1] = function value at each grid point
%   dxGrid = [m, 2*n-1] = function slope at each grid point
%   t = [1, k] = vector of query times (must be contained within tGrid)
%
% OUTPUTS:
%   x = [m, k] = function value at each query time
%
% NOTES:
%   If t is out of bounds, then all corresponding values for x are replaced
%   with NaN
%

nGrid = length(tGrid);
if mod(nGrid-1,2)~=0 || nGrid < 3
    error('The number of grid-points must be odd and at least 3');
end

% Figure out sizes
n = floor((length(tGrid)-1)/2);
m = size(xGrid,1);
k = length(t);
x = zeros(m, k);

% Figure out which segment each value of t should be on
edges = [-inf, tGrid(1:2:end), inf];
[~, bin] = histc(t,edges);

% Loop over each quartic segment
for i=1:n
    idx = bin==(i+1);
    if sum(idx) > 0
        gridIdx = 2*(i-1) + [1,2,3];
        x(:,idx) = quartInterp(...
            tGrid(gridIdx), ...
            xGrid(:,gridIdx), ...
            dxGrid(:,gridIdx), ...
            t(idx));
    end
end

% Replace any out-of-bounds queries with NaN
outOfBounds = bin==1 | bin==(n+2);
x(:,outOfBounds) = nan;

end

function x = quartInterp(tGrid,xGrid,dxGrid,t)
%
% This function computes the interpolant over a single interval
%
% INPUTS:
%   tGrid = [1, 3] = time at endpoints and midpoint
%   xGrid = [m, 3] = function at endpoints and midpoint
%   dxGrid = [m, 3] = derivative at endpoints and midpoint
%   t = [1, p] = query times, spanned by tGrid
%
% OUTPUTS:
%   x = [m, p] = function at query times
%

% Rescale the query points to be on the domain [-1,1]
t = 2*(t-tGrid(1))/(tGrid(3)-tGrid(1)) - 1;

% Unpack function and derivative:
xLow = xGrid(:,1);
xMid = xGrid(:,2);
xUpp = xGrid(:,3);
dxLow = dxGrid(:,1);
dxUpp = dxGrid(:,3);

% Compute the coefficients:
a = dxUpp/4 - dxLow/4 - xLow/2 + xMid - xUpp/2;
b = dxLow/4 + dxUpp/4 + xLow/4 - xUpp/4;
c = dxLow/4 - dxUpp/4 + xLow - 2*xMid + xUpp;
d = (3*xUpp)/4 - dxUpp/4 - (3*xLow)/4 - dxLow/4;
e = xMid;

% Evaluate the polynomial for each dimension of the function:
p = length(t);
m = size(xGrid,1);
x = zeros(m,p);
for i=1:m
    x(i,:) = e(i,:) + t.*(d(i,:) + t.*(c(i,:) + t.*(b(i,:) + t.*a(i,:))));
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function test_failed = test_libltfat_long2fir(varargin)
test_failed = 0;

fprintf(' ===============  %s ================ \n',upper(mfilename));

definput.flags.complexity={'double','single'};
[flags]=ltfatarghelper({},definput,varargin);
dataPtr = [flags.complexity, 'Ptr'];

Larr =    [1, 9,11,110,1, 9, 8, 8,   11, 11, 10];
Loutarr = [11,9,12,221,1,10, 10, 11, 15, 16, 1000];

for do_complex = 0:1
    complexstring = '';
    if do_complex, complexstring = 'complex'; end
    funname = makelibraryname('long2fir',flags.complexity,do_complex);

    for Lidx = 1:numel(Larr)
        Lout = Larr(Lidx);
        L = Loutarr(Lidx);

            if do_complex
                z = cast((1:L)'+1i*(L:-1:1)',flags.complexity);
                zi = complex2interleaved(z);
                zout = randn(2*Lout,1,flags.complexity);

                ziPtr = libpointer(dataPtr,zi);
                zoutPtr = libpointer(dataPtr,zout);
            else
                z = cast((1:L)',flags.complexity);
                zi = z;
                zout = randn(Lout,1,flags.complexity);

                ziPtr = libpointer(dataPtr,zi);
                zoutPtr = libpointer(dataPtr,zout);
            end

            trueres = long2fir(z,Lout);

            status = calllib('libltfat',funname,ziPtr,L,Lout,zoutPtr);

            if do_complex
                res = norm(trueres - interleaved2complex(zoutPtr.Value));
            else
                res = norm(trueres - zoutPtr.Value);
            end

            [test_failed,fail]=ltfatdiditfail(res+status,test_failed,0);
            fprintf(['FIR2LONG OP L:%3i, Lout:%3i, %s %s %s %s\n'],L,Lout,flags.complexity,complexstring,ltfatstatusstring(status),fail);

            status = calllib('libltfat',funname,ziPtr,L,Lout,ziPtr);

            if do_complex
                res = norm(trueres - postpad(interleaved2complex(ziPtr.Value),numel(trueres)));
            else
                res = norm(trueres - postpad(ziPtr.Value,numel(trueres)));
            end

            [test_failed,fail]=ltfatdiditfail(res+status,test_failed,0);
            fprintf(['FIR2LONG IP L:%3i, Lout:%3i, %s %s %s %s\n'],L,Lout,flags.complexity,complexstring,ltfatstatusstring(status),fail);
    end
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function [list, path]=sct_tools_ls(fname, keeppath, keepext, folders,arborescence,select)
% [list, path]=sct_tools_ls(fname, keeppath?, keepext?, folders?,recursive?)
% Example: sct_tools_ls('ep2d*')
% example 2: sct_tools_ls('*',[],[],1) --> folders only
% example 3: sct_tools_ls('*',[],[],2) --> files only

if nargin < 2, keeppath=0; end
if nargin < 3, keepext=1; end
if nargin < 4 || isempty(folders), folders=0; end
if nargin < 6 || isempty(select), select=0; end
if nargin < 5, arborescence=0; end

% [list, path]=sct_tools_ls('*T.txt);
list=dir(fname);
[path,name,ext]= fileparts(fname);
path=[path filesep]; name = [name ext];
if strcmp(path,filesep)
    path=['.' filesep];
end

if folders==1
    list=list(cat(1,list.isdir));
elseif folders==2
    list=list(~cat(1,list.isdir));
end

% sort by name
list=sort_nat({list.name})';

% remove files starting with .
list(cellfun(@(x) strcmp(x(1),'.'), list))=[];
if keeppath
    for iL=1:length(list)
        list{iL}=[path list{iL}];
    end
end
pathcur = path;
path = repmat({path},[length(list),1]);

if ~keepext
    list=cellfun(@(x) sct_tool_remove_extension(x,keeppath),list,'UniformOutput',false);
end

if arborescence
    listdir = sct_tools_ls(pathcur,1,1,1);
    for idir = 1:length(listdir)
        [listidir, pathidir]=sct_tools_ls([listdir{idir} filesep name], keeppath, keepext, folders,arborescence,0);
        list = [list; listidir];
        path = cat(1,path, pathidir{:});
    end
end

if select, list=list{select}; end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function S=capitalize(S,skipset)
%function S=capitalize(S)
% Capitalize the first character in every word (excluding short words) in input
% string S, or in every word in every string in input cell array S.  The first word
% is always capitalized.
%
% Default set of words to skip is {'a','an','and','or','in','the'}
%
% function S=capitalize(S,skipset)
% User-supplied cell array of words to skip

if nargin<2
    skipset={'a','an','and','or','in','the'};
end
skipset=lower(skipset);

if iscell(S)
    S=cellfun(@(x) capitalize(x, skipset),S,'uniformoutput',false);
else
    skip=0;
    Sw=regexp(S,'[ \t]','split');
    skip(2:length(Sw))=ismember(lower(Sw(2:end)),skipset); %skip(1) is always 0
    Sw(~skip)=cellfun(@(x) [upper(x(1)) x(2:end)],Sw(~skip),'uniformoutput',false);
    S=sprintf('%s ',Sw{:});
    S(end)='';
end


text:
function r = subsasgn(r,s,b)
%SUBSASGN     Implements subscripted assignments for Taylor
%
%  example  r(2,:) = b
%

% written  05/21/09     S.M. Rump
% modified 02/28/10     S.M. Rump  multiple index
% modified 08/26/12     S.M. Rump  global variables removed
% modified 09/26/12     S.M. Rump  index handling (thanks to Matthew Weinstein)
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  K1 = getappdata(0,'INTLAB_TAYLOR_ORDER') + 1;

  if length(s)>1
    error('multiple indexing for Taylor assignment not allowed')
  end

  if strcmp(s.type,'()')     % assignment r(i) = b

    rEmpty = isempty(r);    % assignment r(i) = b for empty r
    if isempty(b)
      % does not work in Matlab 5.3 for sparse r
      index = reshape(1:prod(r.size),r.size);
      index(s.subs{:}) = [];
      r.size = size(index);
      r.t = r.t(:,index);
      if rndold
        setround(rndold)
      end
      return
    end

    if ~isa(b,'taylor')
      b = taylor(b);
    end

    if ~rEmpty
      resultIsintval = isa(r.t,'intval');
      if ~resultIsintval & isa(b.t,'intval')
        r = intval(r);
        resultIsintval = 1;
      end
      sizeincreased = 0;
      rsizeold = r.size;
      if length(s.subs)==1               % single index
        if ~isequal(s.subs{1},':')       % not call r(:)=...
          if r.size(1)==1                % row vector
            M = s.subs{1};
            sizeincreased = ( M>r.size(2) );
            if sizeincreased
              r.size(2) = M;
            end
          else
            sizeincreased = ( s.subs{1} > prod(r.size) );
            if sizeincreased
              error('In an assignment A(I) = B a Taylor matrix cannot be resized')
            end
          end
          if sizeincreased
            srx = r.size;
            if length(srx)==2
              if all(prod(srx)~=srx)
                error('matrix cannot be resized by assignment a(I) = b')
              end
            else
              error('attempt to grow size of array along ambiguous dimension')
            end
          end
        end
      else                            % multiple index
        for i=1:length(s.subs)
          if ~isequal(s.subs{i},':')
            M = max(s.subs{i});
            if ( M > r.size(i) )
              sizeincreased = 1;
              r.size(i) = M;
            end
          end
        end
      end
      if sizeincreased                % size increased, adapt .t
        rt = r.t;
        value = ones(rsizeold);
        value( s.subs{:} ) = 0;
        r.t = zeros(K1,prod(size(value)));
        if resultIsintval
          r.t = intval(r.t);
        end
        r.t( : , value==1 ) = rt;
      end
    else                     % assignment r(i) = b for empty r
      resultIsintval = isa(b.t,'intval');
      for i=1:length(s.subs)
        r.size(i) = max(s.subs{i});
      end
      if length(r.size)==1
        r.size = [1 r.size];
      end
      N = prod(r.size);
      if N>10
        r.t = sparse([],[],[],K1,N);
      else
        r.t = zeros(K1,N);
      end
      index = reshape(1:prod(r.size),r.size);
      r.t(:,index(s.subs{:})) = 0;
      if resultIsintval
        r.t = intval(r.t);
      end
    end
    value = reshape(1:prod(r.size),r.size);
    index = value( s.subs{:} );
    if ( prod(b.size)==1 ) & ( length(index)~=1 )
      r.t(:,index) = repmat(b.t,1,length(index));
    else
      r.t(:,index) = b.t;
    end
    if rEmpty
      r = class(r,'taylor');
    end
  else
    error('invalid index reference for taylor')
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
function ft_write_json(filename, json, varargin)

% FT_WRITE_JSON writes a MATLAB structure to a JSON file. Compared to the builtin
% MATLAB function, this implementation deals a bit different with missing values,
% booleans, and NaNs, and results in a more human-readable file.
%
% Use as
%   ft_write_json(filename, struct)
%
% See also FT_READ_JSON, JSONDECODE, JSONENCODE

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

% get the options
sort = ft_getopt(varargin, 'sort', true); % can be true/false, or 'yes'/'no'

ft_info('writing ''%s''\n', filename);
json = remove_empty(json);
if istrue(sort)
 json = sort_fields(json);
end
json = ft_struct2char(json); % convert strings into char-arrays
ft_hastoolbox('jsonlab', 1);
% write nan as 'n/a'
% write boolean as True/False
str = savejson('', json, 'NaN', '"n/a"', 'ParseLogical', true);
fid = fopen_or_error(filename, 'w');
fwrite(fid, str);
fclose(fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION this is shared with DATA2BIDS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = remove_empty(s)
if isempty(s)
  return
elseif isstruct(s)
  fn = fieldnames(s);
  fn = fn(structfun(@isempty, s));
  s = removefields(s, fn);
elseif istable(s)
  remove = false(1,size(s,2));
  for i=1:size(s,2)
    % find columns that are non-numeric and where all elements are []
    remove(i) = ~isnumeric(s{:,i}) && all(cellfun(@isempty, s{:,i}));
  end
  s = s(:,~remove);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION this is shared with DATA2BIDS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = sort_fields(x)
fn = fieldnames(x);
fn = sort(fn);
y = struct();
for i=1:numel(fn)
  y.(fn{i}) = x.(fn{i});
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


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

    xc = 60; yc = 90;

    sensor = 200 ./ ((x-xc).^2 + (y-yc).^2 + 200);

author: petercorke

repo: robotics-toolbox-matlab


text:
function [ element_node ] = grid_t3_element ( nelemx, nelemy )

%*****************************************************************************80
%
%% GRID_T3_ELEMENT produces a grid of pairs of 3 node triangles.
%
%  Example:
%
%    Input:
%
%      NELEMX = 3, NELEMY = 2
%
%    Output:
%
%      ELEMENT_NODE =
%         1,  2,  5;
%         6,  5,  2;
%         2,  3,  6;
%         7,  6,  3;
%         3,  4,  7;
%         8,  7,  4;
%         5,  6,  9;
%        10,  9,  6;
%         6,  7, 10;
%        11, 10,  7;
%         7,  8, 11;
%        12, 11,  8.
%
%  Grid:
%
%    9---10---11---12
%    |\ 8 |\10 |\12 |
%    | \  | \  | \  |
%    |  \ |  \ |  \ |
%    |  7\|  9\| 11\|
%    5----6----7----8
%    |\ 2 |\ 4 |\ 6 |
%    | \  | \  | \  |
%    |  \ |  \ |  \ |
%    |  1\|  3\|  5\|
%    1----2----3----4
%
%  Reference Element T3:
%
%    |
%    1  3
%    |  |\
%    |  | \
%    S  |  \
%    |  |   \
%    |  |    \
%    0  1-----2
%    |
%    +--0--R--1-->
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NELEMX, NELEMY, the number of elements along the
%    X and Y directions.  The number of elements generated will be
%    2 * NELEMX * NELEMY.
%
%    Output, integer ELEMENT_NODE(3,2*NELEMX*NELEMY), the nodes that form
%    each element.
%

%
%  Node labeling:
%
%    NW--NE
%     |\ |
%     | \|
%    SW--SE
%
  element = 0;

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

      sw = i     + ( j - 1 ) * ( nelemx + 1 );
      se = i + 1 + ( j - 1 ) * ( nelemx + 1 );
      nw = i     +   j       * ( nelemx + 1 );
      ne = i + 1 +   j       * ( nelemx + 1 );

      element = element + 1;

      element_node(1,element) = sw;
      element_node(2,element) = se;
      element_node(3,element) = nw;

      element = element + 1;

      element_node(1,element) = ne;
      element_node(2,element) = nw;
      element_node(3,element) = se;

    end
  end

  return
end


text:
function y = vertcat(varargin)
%VERTCAT (overloaded)

prenargin = nargin;
% Fast exit
if prenargin<2
    y=varargin{1};
    return
end

% Get dimensions
n = zeros(prenargin,1);
m = zeros(prenargin,1);
for i = 1:prenargin
    if isa(varargin{i},'blkvar')
        varargin{i} = sdpvar(varargin{i});
    end
    [n(i),m(i)]=size(varargin{i});
end

% Keep only non-empty
keep_these = find((n.*m)~=0);
if length(keep_these)<length(n)
    varargin = {varargin{keep_these}};
    n = n(keep_these);
    m = m(keep_these);
end;

% All widths should be equal
if any(m~=m(1))
    error('All rows in the bracketed expression must have the same number of columns.');
end

nblocks = size(varargin,2);
isasdpvar = zeros(nblocks,1);
% Get dimensions
for i = 1:nblocks
    isasdpvar(i) = isa(varargin{i},'sdpvar');
end

% Find all free variables used
all_lmi_variables = [];
for i = 1:nblocks
    if isasdpvar(i)
        all_lmi_variables = [all_lmi_variables varargin{i}.lmi_variables];
    end
end
all_lmi_variables = uniquestripped(all_lmi_variables);

% Pick one of the sdpvar objects to build on...
y = varargin{min(find(isasdpvar))};

% Some indexation tricks...
basis_i = [];
basis_j = [];
basis_s = [];
m = m(1);
sumn=sum(n);
shift=[0;cumsum(n)];
for j = 1:nblocks
    if isasdpvar(j)
        if length(all_lmi_variables)==length(varargin{j}.lmi_variables) && all_lmi_variables(1)==varargin{j}.lmi_variables(1) &&  all_lmi_variables(end)==varargin{j}.lmi_variables(end)
            % Avoid call to ismember and find
            in_this = 1:length(all_lmi_variables);
        else
            members = ismembcYALMIP(all_lmi_variables,varargin{j}.lmi_variables);
            in_this = find(members);
        end
        dummy = [1 1+in_this];
        [i2,j2,s2] = find(varargin{j}.basis);
        j2 = dummy(j2);
    else
        [i2,j2,s2] = find(varargin{j}(:));
    end
    if ~isempty(j2)
    col = 1+floor((i2-1)/n(j));
    row = i2-(col-1)*n(j);
    basis_i = [basis_i;row(:)+shift(j)+sumn*(col(:)-1)];
    basis_j = [basis_j;j2(:)];
    basis_s = [basis_s;s2(:)];
    end
end
basis = sparse(basis_i,basis_j,basis_s,sum(n)*m,1+length(all_lmi_variables));

y.dim(1) = sum(n);
y.dim(2) = m;
y.basis = basis;
y.lmi_variables = all_lmi_variables;
% Reset info about conic terms
y.conicinfo = [0 0];
y.extra.opname='';
y.extra.createTime = definecreationtime;


text:
function [W] = hp2W(hp)
% Convert power from mechanical horsepower to watts.
% Chad A. Greene 2012
W = hp*745.699871582;


text:
function asa159_test ( )

%*****************************************************************************80
%
%% ASA159_TEST tests ASA159.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 March 2009
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA159_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Test the ASA159 routines.\n' );

  asa159_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ASA159_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function A=ip2mat(B,m)
A=flipud(reshape(B,m)');

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function bench_plot(times,res,solvers,prob)
%BENCH_PLOT  Plot the results of a benchmark test
%
% bench_plot(times,res,solvers,prob) not intended to be called by a user -
% use optiBench.

%   Copyright (C) 2011 Jonathan Currie (I2C2)

nS = length(solvers);   % Number of Solvers
nP = length(times{1});  % Number of Test problems

tol = 1e-4;

%Determine maximum time for normalizing
maxT = 0;
for i = 1:nS
    s = sum(times{i});
    if(s > maxT), maxT = s; end
end
% maxT = 1;

[~,srtInd] = sort(times{1});

%Generate Plot Data
X = zeros(2*nP+3,nS); Y = X;
str = cell(nS,1);
for i = 1:nS
    acc = res{i}.acc(srtInd); ind = isinf(acc);
    if(any(ind))
        acc(ind) = 1;
    end
    ind = acc < tol;
    fM = sum(ind)/nP;
    fC = [0;cumsum(ind)/nP;fM];
    fT = [0;cumsum(times{i}(srtInd));maxT];
    [X(:,i),Y(:,i)] = stairs(fT,fC);
    str{i} = [upper(solvers{i}) ': ' num2str(fM*100,'%2.2f') '%'];
end

% %Determine maximum time for normalizing
% maxT = 0;
% for i = 1:nS
%     s = sum(times{i});
%     if(s > maxT), maxT = s; end
% end
%
% %Generate Plot Data
% X = zeros(2*nP+3,nS); Y = X;
% str = cell(nS,1);
% for i = 1:nS
%     acc = res{i}.acc; ind = isinf(acc);
%     if(any(ind))
%         acc(ind) = 1;
%     end
%     ind = acc < tol;
%     fM = sum(ind)/nP;
%     fC = [0;cumsum(ind)/nP;fM];
%     fT = [0;cumsum(times{i})/maxT;1];
%     [X(:,i),Y(:,i)] = stairs(fT,fC);
%     str{i} = [upper(solvers{i}) ': ' num2str(fM*100,'%2.2f') '%'];
% end

%Plot
plot(X,Y);
xlabel('Time [s]');
ylabel('Fraction of Problems Solved');
title([upper(prob) ' Benchmark Results for ' num2str(nP) ' Problems']);
axis([0 maxT 0 1.1]);
legend(str,'location','SE');

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function varargout=cosmo_warning(message, varargin)
% show a warning message; by default just once for each message
%
% cosmo_warning(message, ...)
% cosmo_warning(state)
% state=cosmo_warning()
%
% Inputs:
%   message      warning message to be shown, or one of:
%                'on'   : show all warning messages
%                'off'  : show no warning messages
%                'once' : show each warning message once [default]
%                'reset':
%   ...          if a warning message is provided according with
%                placeholders as used in sprintf, then the subsequent
%                arguments should contain their values
%   state        if a struct, then this queries or sets the state of
%                cosmo_warning.
%
% Notes:
%   - this function works more or less like matlab's warning function,
%     except that by default each warning is just shown once.
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    if isempty(get_from_state('when'))
        set_default_state();
    end

    if nargin==0
        varargout={get_state()};
        return
    end

    if isstruct(message)
        set_state(message);
        return
    end

    lmessage=lower(message);

    switch lmessage
        case {'on','off','once'}
            update_state('when',lmessage);
            return;

        case 'reset'
            set_default_state();
            return;

        otherwise
            show_warning(message,varargin{:});
    end

function show_warning(message,varargin)
    [identifier,full_message]=get_identifier_and_message(...
                                    message,varargin{:});

    shown_warnings=get_from_state('shown_warnings');
    has_warning=cosmo_match({full_message},shown_warnings);
    if ~has_warning
        shown_warnings{end+1}=full_message;
        update_state('shown_warnings',shown_warnings);
    end

    when=get_from_state('when');
    switch when
        case 'once'
            do_show_warning=~has_warning;

            me=mfilename();
            postfix=sprintf(['\n\nThis warning is shown only once, '...
                           'but the underlying issue may occur '...
                           'multiple times. To show each warning:\n'...
                           ' - every time:   %s(''on'')\n'...
                           ' - once:         %s(''once'')\n'...
                           ' - never:        %s(''off'')\n'],me,me,me);
            full_message=[full_message postfix];

        case 'off'
            do_show_warning=false;
        case 'on'
            do_show_warning=true;
        otherwise
            assert(false);
    end

    if do_show_warning
        state=warning(); % store state
        state_resetter=onCleanup(@()warning(state));

        warning('on','all');
        % avoid extra entry on the stack
        has_identifier=~isempty(identifier);

        if has_identifier
            warning(identifier,'%s',full_message);
        else
            warning('%s',full_message);
        end
    end

function [identifier,full_message]=get_identifier_and_message(...
                                            message,varargin)
    args=varargin;
    has_identifier=numel(args)>0 && has_warning_identifier(message);
    if has_identifier
        identifier=message;
        message=args{1};
        args=args(2:end);
    else
        identifier='';
    end

    if numel(args)>0
        full_message=sprintf(message, args{:});
    else
        full_message=message;
    end

function tf=has_warning_identifier(s)
    alpha_num='([a-z_A-Z0-9]+)';
    pat=sprintf('^%s(:%s)?:%s$',alpha_num,alpha_num,alpha_num);
    tf=~isempty(regexp(s,pat,'once'));

function s=get_state()
    s=get_or_set_state();

function set_state(s)
    get_or_set_state(s);

function set_default_state()
    s=struct();
    s.when='once';
    s.shown_warnings=cell(0);
    set_state(s);

function value=get_from_state(key)
    s=get_state();
    value=s.(key);

function update_state(key, value)
    s=get_state();
    s.(key)=value;
    set_state(s);

function varargout=get_or_set_state(s)
    persistent state;
    switch nargin
        case 0
            % get state
            if isempty(state)
                set_default_state();
            end

            varargout={state};

        case 1
            % set state
            state=s;
            varargout={};

        otherwise
            assert(false);
    end

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function [tm,im] = mapvert(tr,pi)
%MAPVERT find the tree-to-vertex mappings.
%   [TM,IM] = MAPVERT(TR,PI) returns the tree-to-vertex and
%   vertex-to-tree mappings for a given aabb-tree TR and a
%   collection of query vertices PI.
%
%   The tree-to-item mapping TM is a structure representing
%   the intersection of the items PI with the tree TR. TM.II
%   is an M-by-1 array of tree indices and TM.LL is an
%   M-by-1 cell array of item lists. Specifically, items in
%   the list TM.LL{JJ} intersect with the node TM.II(JJ).
%
%   The item-to-tree mapping IM is a structure representing
%   the inverse mapping. IM.II is an N-by-1 array of item
%   indices and IM.LL is an N-by-1 cell array of node lists.
%   Specifically, nodes in the list IM.LL{JJ} intersect with
%   the item IM.II(JJ).
%
%   See also QUERYSET, MAPRECT, MAKETREE

%   Darren Engwirda : 2014 --
%   Email           : de2363@columbia.edu
%   Last updated    : 06/04/2017

%----------------------- call SCANTREE to do the actual work
    if (nargout == +1)
       [tm   ] = scantree(tr,pi,@partvert);
    else
       [tm,im] = scantree(tr,pi,@partvert);
    end

end

function [j1,j2] = partvert(pi,b1,b2)
%PARTVERT partition points between boxes B1,B2 for SCANTREE.

    j1 = true(size(pi,1),1);
    j2 = true(size(pi,1),1);

    nd = size(b1,2) / +2;

    for ax = +1 : nd
%--------------- remains TRUE if inside bounds along axis AX
    j1 = j1 & pi(:,ax)>=b1(ax+nd*0) ...
            & pi(:,ax)<=b1(ax+nd*1) ;

%--------------- remains TRUE if inside bounds along axis AX
    j2 = j2 & pi(:,ax)>=b2(ax+nd*0) ...
            & pi(:,ax)<=b2(ax+nd*1) ;
    end

end

author: dengwirda


text:
function M = productmanifold(elements)
% Returns a structure describing a product manifold M = M1 x M2 x ... x Mn.
%
% function M = productmanifold(elements)
%
% Input: an elements structure such that each field contains a manifold
% structure.
%
% Output: a manifold structure M representing the manifold obtained by
% taking the Cartesian product of the manifolds described in the elements
% structure, with the metric obtainded by element-wise extension. Points
% and vectors are stored as structures with the same fieldnames as in
% elements.
%
% Example:
% M = productmanifold(struct('X', spherefactory(3), 'Y', spherefactory(4)))
% disp(M.name());
% x = M.rand()
%
% Points of M = S^2 x S^3 are represented as structures with two fields, X
% and Y. The values associated to X are points of S^2, and likewise points
% of S^3 for the field Y. Tangent vectors are also represented as
% structures with two corresponding fields X and Y.
%
% See also: powermanifold

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, Dec. 30, 2012.
% Contributors:
% Change log:
%
%   July  4, 2013 (NB):
%       Added support for vec, mat, tangent.
%       Added support for egrad2rgrad and ehess2rhess.
%       Modified hash function to make hash strings shorter.
%
%   Dec. 17, 2018 (NB):
%       Added check all_elements_provide() to many functions, so that if,
%       for example, one of the elements does not provide exp(), then the
%       product manifold also won't provide exp(). This makes it easier for
%       tools such as, for example, checkgradient, to determine whether exp
%       is available or not.
%
%   Feb. 10, 2020 (NB):
%       Added warnings about calling egrad2rgrad and ehess2rhess without
%       storedb and key, even if some base manifolds allow them.
%
%   Jan. 4, 2021 (NB):
%       Changes for compatibility with Octave 6.1.0: by introducing a
%       "helper" function, we separate out the pre-computations. This way,
%       all pre-computed quantities are passed as input to the helper
%       function. This makes them available to nested subfunctions.
%       The extra step is not necessary in Matlab.

    elems = fieldnames(elements);
    nelems = numel(elems);

    assert(nelems >= 1, ...
           'elements must be a structure with at least one field.');

    % Below are some precomputations for the mat/vec pair.
    %
    % Gather the length of the column vector representations of tangent
    % vectors for each of the manifolds. Raise a flag if any of the base
    % manifolds has no vec function available.
    vec_available = true;
    vec_lens = zeros(nelems, 1);
    for ii = 1 : nelems
        Mi = elements.(elems{ii});
        if isfield(Mi, 'vec')
            rand_x = Mi.rand();
            zero_u = Mi.zerovec(rand_x);
            vec_lens(ii) = length(Mi.vec(rand_x, zero_u));
        else
            vec_available = false;
            break;
        end
    end
    vec_pos = cumsum([1 ; vec_lens]);
    %
    vecmatareisometries = vec_available;
    for ii = 1 : nelems
        if ~isfield(elements.(elems{ii}), 'vecmatareisometries') || ...
           ~elements.(elems{ii}).vecmatareisometries()
            vecmatareisometries = false;
            break;
        end
    end
    %
    % Above are some precomputations for the mat/vec pair.

    % The helper function is the actual factory.
    M = productmanifoldhelper(elements, elems, nelems, vec_available, ...
                              vec_pos, vecmatareisometries);

end

function M = productmanifoldhelper(elements, elems, nelems, ...
                                   vec_available, vec_pos, ...
                                   vecmatareisometries)

    % Handy function to check if all elements provide the necessary methods
    function answer = all_elements_provide(method_name)
        answer = false;
        for i = 1 : nelems
            if ~isfield(elements.(elems{i}), method_name)
                return;
            end
        end
        answer = true;
    end

    M.name = @name;
    function str = name()
        str = 'Product manifold: ';
        str = [str sprintf('[%s: %s]', ...
                           elems{1}, elements.(elems{1}).name())];
        for i = 2 : nelems
            str = [str sprintf(' x [%s: %s]', ...
                   elems{i}, elements.(elems{i}).name())]; %#ok<AGROW>
        end
    end

    M.dim = @dim;
    function d = dim()
        d = 0;
        for i = 1 : nelems
            d = d + elements.(elems{i}).dim();
        end
    end

    M.inner = @inner;
    function val = inner(x, u, v)
        val = 0;
        for i = 1 : nelems
            val = val + elements.(elems{i}).inner(x.(elems{i}), ...
                                               u.(elems{i}), v.(elems{i}));
        end
    end

    M.norm = @(x, d) sqrt(M.inner(x, d, d));

    if all_elements_provide('dist')
        M.dist = @dist;
    end
    function d = dist(x, y)
        sqd = 0;
        for i = 1 : nelems
            sqd = sqd + elements.(elems{i}).dist(x.(elems{i}), ...
                                                 y.(elems{i}))^2;
        end
        d = sqrt(sqd);
    end

    if all_elements_provide('typicaldist')
        M.typicaldist = @typicaldist;
    end
    function d = typicaldist
        sqd = 0;
        for i = 1 : nelems
            sqd = sqd + elements.(elems{i}).typicaldist()^2;
        end
        d = sqrt(sqd);
    end

    M.proj = @proj;
    function v = proj(x, u)
        for i = 1 : nelems
            v.(elems{i}) = elements.(elems{i}).proj(x.(elems{i}), ...
                                                    u.(elems{i}));
        end
    end

    M.tangent = @tangent;
    function v = tangent(x, u)
        for i = 1 : nelems
            v.(elems{i}) = elements.(elems{i}).tangent(x.(elems{i}), ...
                                                       u.(elems{i}));
        end
    end

    % True by default, false if any false encountered
    M.tangent2ambient_is_identity = true;
    for k = 1 : nelems
        if isfield(elements.(elems{k}), 'tangent2ambient_is_identity')
            if ~elements.(elems{k}).tangent2ambient_is_identity
                M.tangent2ambient_is_identity = false;
                break;
            end
        end
    end

    M.tangent2ambient = @tangent2ambient;
    function v = tangent2ambient(x, u)
        for i = 1 : nelems
            if isfield(elements.(elems{i}), 'tangent2ambient')
                v.(elems{i}) = ...
                    elements.(elems{i}).tangent2ambient( ...
                                               x.(elems{i}), u.(elems{i}));
            else
                v.(elems{i}) = u.(elems{i});
            end
        end
    end

    M.egrad2rgrad = @egrad2rgrad;
    function g = egrad2rgrad(x, g)
        for i = 1 : nelems
            g.(elems{i}) = elements.(elems{i}).egrad2rgrad(...
                                               x.(elems{i}), g.(elems{i}));
        end
    end
    for ii = 1 : nelems
        if nargin(elements.(elems{ii}).egrad2rgrad) > 2
            warning('manopt:productmanifold:egrad2rgrad', ...
                   ['Product manifolds call M.egrad2rgrad with only two ', ...
                    'inputs:\nstoredb and key won''t be available.']);
            break;
        end
    end

    M.ehess2rhess = @ehess2rhess;
    function h = ehess2rhess(x, eg, eh, h)
        for i = 1 : nelems
            h.(elems{i}) = elements.(elems{i}).ehess2rhess(...
                 x.(elems{i}), eg.(elems{i}), eh.(elems{i}), h.(elems{i}));
        end
    end
    for ii = 1 : nelems
        if nargin(elements.(elems{ii}).ehess2rhess) > 4
            warning('manopt:productmanifold:ehess2rhess', ...
                   ['Product manifolds call M.ehess2rhess with only two ', ...
                    'inputs:\nstoredb and key won''t be available.']);
            break;
        end
    end

    if all_elements_provide('exp')
        M.exp = @exp;
    end
    function y = exp(x, u, t)
        if nargin < 3
            t = 1.0;
        end
        for i = 1 : nelems
            y.(elems{i}) = elements.(elems{i}).exp(x.(elems{i}), ...
                                                   u.(elems{i}), t);
        end
    end

    M.retr = @retr;
    function y = retr(x, u, t)
        if nargin < 3
            t = 1.0;
        end
        for i = 1 : nelems
            y.(elems{i}) = elements.(elems{i}).retr(x.(elems{i}), ...
                                                    u.(elems{i}), t);
        end
    end

    if all_elements_provide('log')
        M.log = @log;
    end
    function u = log(x1, x2)
        for i = 1 : nelems
            u.(elems{i}) = elements.(elems{i}).log(x1.(elems{i}), ...
                                                   x2.(elems{i}));
        end
    end

    M.hash = @hash;
    function str = hash(x)
        str = '';
        for i = 1 : nelems
            str = [str elements.(elems{i}).hash(x.(elems{i}))]; %#ok<AGROW>
        end
        str = ['z' hashmd5(str)];
    end

    M.lincomb = @lincomb;
    function v = lincomb(x, a1, u1, a2, u2)
        if nargin == 3
            for i = 1 : nelems
                v.(elems{i}) = elements.(elems{i}).lincomb(x.(elems{i}), ...
                                                        a1, u1.(elems{i}));
            end
        elseif nargin == 5
            for i = 1 : nelems
                v.(elems{i}) = elements.(elems{i}).lincomb(x.(elems{i}), ...
                                     a1, u1.(elems{i}), a2, u2.(elems{i}));
            end
        else
            error('Bad usage of productmanifold.lincomb');
        end
    end

    M.rand = @rand;
    function x = rand()
        for i = 1 : nelems
            x.(elems{i}) = elements.(elems{i}).rand();
        end
    end

    M.randvec = @randvec;
    function u = randvec(x)
        for i = 1 : nelems
            u.(elems{i}) = elements.(elems{i}).randvec(x.(elems{i}));
        end
        u = M.lincomb(x, 1/sqrt(nelems), u);
    end

    M.zerovec = @zerovec;
    function u = zerovec(x)
        for i = 1 : nelems
            u.(elems{i}) = elements.(elems{i}).zerovec(x.(elems{i}));
        end
    end

    if all_elements_provide('transp')
        M.transp = @transp;
    end
    function v = transp(x1, x2, u)
        for i = 1 : nelems
            v.(elems{i}) = elements.(elems{i}).transp(x1.(elems{i}), ...
                                              x2.(elems{i}), u.(elems{i}));
        end
    end

    if all_elements_provide('pairmean')
        M.pairmean = @pairmean;
    end
    function y = pairmean(x1, x2)
        for i = 1 : nelems
            y.(elems{i}) = elements.(elems{i}).pairmean(x1.(elems{i}), ...
                                                        x2.(elems{i}));
        end
    end

    if vec_available
        M.vec = @vec;
        M.mat = @mat;
    end

    function u_vec = vec(x, u_mat)
        u_vec = zeros(vec_pos(end)-1, 1);
        for i = 1 : nelems
            range = vec_pos(i) : (vec_pos(i+1)-1);
            u_vec(range) = elements.(elems{i}).vec(x.(elems{i}), ...
                                                   u_mat.(elems{i}));
        end
    end

    function u_mat = mat(x, u_vec)
        u_mat = struct();
        for i = 1 : nelems
            range = vec_pos(i) : (vec_pos(i+1)-1);
            u_mat.(elems{i}) = elements.(elems{i}).mat(x.(elems{i}), ...
                                                       u_vec(range));
        end
    end

    M.vecmatareisometries = @() vecmatareisometries;

    if all_elements_provide('lie_identity')
        M.lie_identity = @lie_identity;
    end

    function I = lie_identity()
        I = struct();
        for i = 1 : nelems
            Mi = elements.(elems{i});
            Ii = Mi.lie_identity();
            I.(elems{i}) = Ii;
        end
    end
end

author: NicolasBoumal

repo: manopt


text:
function [scr_lin,scr_qual] = linear_and_quality_fusion_scores(key,scores_obj_array,qual_obj_array,Qtrans,...
                                                    prior,niters,obj_func,cstepHessian,quiet)
% Fuses scores (using linear fusion) and then fuses in quality
% measures.  The fusion is applied to the same scores as those used
% to train the fusion function.  The input scores and key must
% already be loaded.  Use 'linear_and_quality_fusion_scores_from_files' if you
% want the script to do loading and saving for you.
% Inputs:
%   key: The Key containing information about the trials.
%   scores_obj_array: An array of Scores objects containing
%     system scores for the trials.
%   qual_obj_array: An array of Quality objects containing
%     quality measures for the segments.
%   Qtrans: A matrix indicating how the quality measures should be
%     combined.  Use an identity matrix of
%     size(numsystems,numsystems) if you want the quality measures
%     to be used as is.
%   prior: The effective target prior.
%   niters: The maximum number of training iterations.
%   obj_func: The objective function for the fusion training.  If [],
%     cllr objective is used.
%   cstepHessian: Boolean.  If true, the training algorithm will
%     calculate the Hessian matrix using complex step
%     differentiation which should make training faster.
%   quiet: Boolean.  If true, the training algorithm outputs fewer
%     messages describing its progress.
% Outputs:
%   scr_lin: The scores resulting from the linear fusion.
%   scr_qual: The scores resulting from the quality fusion.

assert(nargin==9)
assert(isa(key,'Key'))
assert(~isempty(Qtrans))

logprint(Logger.Info,'training linear fuser\n')
[flin,wlin] = train_linear_fusion_function(key,scores_obj_array,prior,obj_func,niters,cstepHessian,quiet);

logprint(Logger.Info,'linear fusing scores\n')
scr_lin = apply_linear_fusion_function(key,scores_obj_array,flin);

logprint(Logger.Info,'training quality fuser\n')
fqual = train_quality_fusion_function(key,scores_obj_array,qual_obj_array,Qtrans,wlin,prior,obj_func,niters,cstepHessian,quiet);

logprint(Logger.Info,'quality fusing scores\n')
scr_qual = apply_quality_fusion_function(key,scores_obj_array,qual_obj_array,Qtrans,fqual,scr_lin);

end

author: nesl


text:
function fd1d_heat_implicit_test ( )

%*****************************************************************************80
%
%% FD1D_HEAT_IMPLICIT_TEST tests the FD1D_HEAT_IMPLICIT library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 January 2012
%
%  Author:
%
%    John Burkardt
%
  timestamp ( )
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD1D_HEAT_IMPLICIT_TEST:\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '  Test the FD1D_HEAT_IMPLICIT library.\n' );

  fd1d_heat_implicit_test01 ( );
  fd1d_heat_implicit_test02 ( );
  fd1d_heat_implicit_test03 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'FD1D_HEAT_IMPLICIT_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function pass = test_determineDiscretization(pref)
% Test DETERMINEDISCRETIZATION method:

if ( nargin == 0 )
    pref = cheboppref();
end

%% Test the VALUES/COEFFS syntax:

% Test when values is passed with periodic boundary conditions:
dom = [0 2*pi];
N = chebop(@(x,u) diff(u,2) + cos(u), dom);
u0 = chebfun('0',dom);
x = chebfun('x',dom);
L = linearize(N, u0, x);
N.bc = 'periodic';
options = cheboppref();
options.discretization = 'values';
out = determineDiscretization(N, length(L.domain), options);
pass(1) = isequal(out.discretization, @trigcolloc);

% Test when coeffs is passed with periodic boundary conditions:
options.discretization = 'coeffs';
out = determineDiscretization(N, length(L.domain), options);
pass(2) = isequal(out.discretization, @trigspec);

% Test when values is passed with periodic boundary conditions and breakpoints:
dom = [0 pi 2*pi];
N = chebop(@(x,u) diff(u,2) + cos(u), dom);
u0 = chebfun('0',dom);
x = chebfun('x',dom);
L = linearize(N, u0, x);
options.discretization = 'values';
out = determineDiscretization(N, length(L.domain), options);
pass(3) = isequal(out.discretization, @chebcolloc2);

% Test when coeffs is passed with periodic boundary conditions and breakpoints:
options.discretization = 'coeffs';
out = determineDiscretization(N, length(L.domain), options);
pass(4) = isequal(out.discretization, @ultraS);

% Test when values is passed with dirichlet boundary conditions:
dom = [-1 1];
N = chebop(@(x,u) diff(u) + exp(u), dom);
u0 = chebfun('0',dom);
x = chebfun('x',dom);
L = linearize(N, u0, x);
N.bc = 'dirichlet';
options.discretization = 'values';
out = determineDiscretization(N, length(L.domain), options);
pass(5) = isequal(out.discretization, @chebcolloc2);

% Test when coeffs is passed with dirichlet boundary conditions:
options.discretization = 'coeffs';
out = determineDiscretization(N, length(L.domain), options);
pass(6) = isequal(out.discretization, @ultraS);

%% Test default:

% Default with dirichlet:
dom = [-1 1];
N = chebop(@(x,u) diff(u) + sin(u), dom);
u0 = chebfun('0',dom);
x = chebfun('x',dom);
L = linearize(N, u0, x);
N.bc = 'dirichlet';
out = determineDiscretization(N, length(L.domain), pref); % use pref
pass(7) = isequal(out.discretization, @chebcolloc2);

% Default with periodic:
N.bc = 'periodic';
out = determineDiscretization(N, length(L.domain), pref); % use pref
pass(8) = isequal(out.discretization, @trigcolloc);

%% Test CHEBCOLLOC1/ULTRAS:

options.discretization = @chebcolloc1;
out = determineDiscretization(N, length(L.domain), options);
pass(9) = isequal(out.discretization, @chebcolloc1);

options.discretization = @ultraS;
out = determineDiscretization(N, length(L.domain), options);
pass(10) = isequal(out.discretization, @ultraS);

%% Test a discontinuous RHS:

% Zero boundary condition at the left:
dom = [-1 1];
L = chebop(dom);
L.op = @(x,u) diff(u) + 2*u;
L.lbc = @(u) u - 0;
rhs = chebfun(@(x) abs(x), 'splitting', 'on');
lengthDom = max(length(L.domain),length(rhs.domain));
out = determineDiscretization(L, lengthDom, pref);
pass(11) = isequal(out.discretization, @chebcolloc2);

% Periodic boundary condition. The rhs is discontinuous so it should use
% CHEBCOLLOC2:
L.bc = 'periodic';
out = determineDiscretization(L, lengthDom, pref);
pass(12) = isequal(out.discretization, @chebcolloc2);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [vert,conn,tria,tnum] = smooth2(varargin)
%SMOOTH2 "hill-climbing" mesh-smoothing for two-dimensional,
%2-simplex triangulations.
%   [VERT,EDGE,TRIA,TNUM] = SMOOTH2(VERT,EDGE,TRIA,TNUM) re-
%   turns a "smoothed" triangulation {VERT,TRIA}, incorpora-
%   ting "optimised" vertex coordinates and mesh topology.
%
%   VERT is a V-by-2 array of XY coordinates in the triangu-
%   lation, EDGE is an array of constrained edges, TRIA is a
%   T-by-3 array of triangles, and TNUM is a T-by-1 array of
%   part indices. Each row of TRIA and EDGE define an eleme-
%   nt. VERT(TRIA(II,1),:), VERT(TRIA(II,2),:) and VERT(TRIA
%   (II,3),:) are the coordinates of the II-TH triangle. The
%   edges in EDGE are defined in a similar manner. NUM is an
%   array of part indexing, such that TNUM(II) is the index
%   of the part in which the II-TH triangle resides.
%
%   [VERT,EDGE,TRIA,TNUM] = SMOOTH2(... ,OPTS) passes an ad-
%   ditional options structure OPTS, containing user-defined
%   parameters, including:
%
% - OPTS.VTOL = {+1.0E-02} -- relative vertex movement tole-
%   rance, smoothing is converged when (VNEW-VERT) <= VTOL *
%   VLEN, where VLEN is a local effective length-scale.
%
% - OPTS.ITER = {+32} -- max. number of smoothing iterations
%
% - OPTS.DISP = {+ 4} -- smoothing verbosity. Set to INF for
%   quiet execution.
%
%   See also REFINE2, TRICOST, TRIDEMO

%   This routine is loosely based on the DISTMESH algorithm,
%   employing a "spring-based" analogy to redistribute mesh
%   vertices. Such an approach is described in: P.O. Persson
%   and Gilbert Strang. "A simple mesh generator in MATLAB."
%   SIAM review 46(2) 2004, pp: 329--345. Details of the al-
%   gorithm used here are somewhat different, with an alter-
%   ative spring-based update employed, in addition to hill-
%   climbing element quality guarantees, and vertex density
%   controls.

%-----------------------------------------------------------
%   Darren Engwirda : 2017 --
%   Email           : de2363@columbia.edu
%   Last updated    : 13/02/2020
%-----------------------------------------------------------

    vert = []; conn = []; tria = [] ;
    tnum = [];
    opts = []; hfun = []; harg = {} ;

%---------------------------------------------- extract args
    if (nargin>=+1), vert = varargin{1}; end
    if (nargin>=+2), conn = varargin{2}; end
    if (nargin>=+3), tria = varargin{3}; end
    if (nargin>=+4), tnum = varargin{4}; end
    if (nargin>=+5), opts = varargin{5}; end

   [opts] = makeopt(opts) ;

%---------------------------------------------- default CONN
    if (isempty(conn))

       [edge] = tricon2(tria);

        ebnd = edge(:,4) < +1;              %-- use bnd edge
        conn = edge(ebnd,1:2);

    end

%---------------------------------------------- default TNUM
    if (isempty(tnum))
        tnum = ones(size(tria, 1), 1) ;
    end

%---------------------------------------------- basic checks
    if ( ~isnumeric(vert) || ...
         ~isnumeric(conn) || ...
         ~isnumeric(tria) || ...
         ~isnumeric(tnum) || ...
         ~isstruct (opts) )
        error('smooth2:incorrectInputClass' , ...
            'Incorrect input class.') ;
    end

%---------------------------------------------- basic checks
    if (ndims(vert) ~= +2 || ...
        ndims(conn) ~= +2 || ...
        ndims(tria) ~= +2 || ...
        ndims(tnum) ~= +2 )
        error('smooth2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

    if (size(vert,2)~= +2 || ...
        size(conn,2)~= +2 || ...
        size(tria,2)~= +3 || ...
        size(tnum,2)~= +1 || ...
        size(tria,1)~= size(tnum,1) )
        error('smooth2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

    nvrt = size(vert,1) ;

%---------------------------------------------- basic checks
    if (min(min(conn(:,1:2))) < +1 || ...
            max(max(conn(:,1:2))) > nvrt )
        error('smooth2:invalidInputs', ...
            'Invalid EDGE input array.') ;
    end

    if (min(min(tria(:,1:3))) < +1 || ...
            max(max(tria(:,1:3))) > nvrt )
        error('smooth2:invalidInputs', ...
            'Invalid TRIA input array.') ;
    end

%---------------------------------------------- output title
    if (~isinf(opts.disp))
        fprintf(1,'\n') ;
        fprintf(1,' Smooth triangulation...\n') ;
        fprintf(1,'\n') ;
        fprintf(1,[...
' -------------------------------------------------------\n', ...
'      |ITER.|          |MOVE(X)|          |DTRI(X)|     \n', ...
' -------------------------------------------------------\n', ...
             ] ) ;
    end

%---------------------------------------------- polygon bnds
    node = vert; PSLG = conn; part = {};

    pmax = max(tnum(:));
    for ppos = +1 : pmax

        tsel = tnum == ppos ;
        tcur = tria(tsel,:) ;

       [ecur,tcur] ...
           = tricon2 (tcur) ;

        ebnd = ecur(:,4)==0 ;

        same = setset2( ...
            PSLG,ecur(ebnd,1:2));

        part{ppos} = find(same) ;

    end

%---------------------------------------------- inflate bbox
    vmin = min(vert,[],1);
    vmax = max(vert,[],1);

    vdel = vmax - 1.*vmin;
    vmin = vmin - .5*vdel;
    vmax = vmax + .5*vdel;

    vbox = [
        vmin(1), vmin(2)
        vmax(1), vmin(2)
        vmax(1), vmax(2)
        vmin(1), vmax(2)
           ] ;
    vert = [vert ; vbox] ;

%---------------------------------------------- DO MESH ITER
    tnow =  tic ;

    tcpu = struct('full',0.,'dtri',0., ...
        'tcon',0.,'iter',0.,'undo',0., ...
        'keep',0.) ;

    for iter = +1 : opts.iter

    %------------------------------------------ inflate adj.
        ttic = tic ;

       [edge,tria] = tricon2(tria,conn) ;

        tcpu.tcon = ...
            tcpu.tcon + toc(ttic) ;

    %------------------------------------------ compute scr.
        oscr = triscr2(vert,tria) ;

    %------------------------------------------ vert. iter's
        ttic =  tic ;

        nvrt = size(vert,1);
        nedg = size(edge,1);

        IMAT = sparse( ...
           edge(:,1),(1:nedg)',+1,nvrt,nedg) ;
        JMAT = sparse( ...
           edge(:,2),(1:nedg)',+1,nvrt,nedg) ;

        EMAT = IMAT + JMAT ;

        vdeg = sum(EMAT,2) ;                %-- vertex |deg|
        free = (vdeg == 0) ;

        vold = vert ;
        for isub = +1 : max(+2,min(+8,iter))

        %-- compute HFUN at vert/midpoints
            hvrt = evalhfn(vert, ...
                edge,EMAT,hfun,harg) ;

            hmid = hvrt(edge(:,1),:) ...
                 + hvrt(edge(:,2),:) ;
            hmid = hmid * +.5 ;

        %-- calc. relative edge extensions
            evec = vert(edge(:,2),:) ...
                 - vert(edge(:,1),:) ;
            elen = ...
                sqrt(sum(evec.^2,2)) ;

            scal = +1.0-elen./hmid ;
            scal = min (+1.0, scal);
            scal = max (-1.0, scal);

        %-- projected points from each end
            ipos = vert(edge(:,1),:) ...
                 -.67*[scal,scal].*evec;
            jpos = vert(edge(:,2),:) ...
                 +.67*[scal,scal].*evec;

           %scal = ...                      %-- nlin. weight
           %   max(abs(scal).^.5,eps^.75);
            scal = ...
               max(abs(scal).^ 1,eps^.75);

        %-- sum contributions edge-to-vert
            vnew = ...
            IMAT*([scal,scal] .* ipos) ...
          + JMAT*([scal,scal] .* jpos) ;

            vsum = max(EMAT*scal,eps^.75);

            vnew = vnew ./ [vsum,vsum] ;

        %-- fixed points. edge projection?
            vnew(conn(:),1:2) = ...
            vert(conn(:),1:2) ;

            vnew(vdeg==0,1:2) = ...
            vert(vdeg==0,1:2) ;

        %-- reset for the next local iter.
            vert = vnew ;

        end

        tcpu.iter = ...
            tcpu.iter + toc(ttic) ;

    %------------------------------------------ hill-climber
        ttic = tic ;

    %-- unwind vert. upadte if score lower
        nscr = ones(size(tria,1),1);
        btri = true(size(tria,1),1);

        umax = + 8 ;

        for undo = +1 : umax

            nscr(btri) = triscr2( ...
                vert,tria(btri,:)) ;

        %-- TRUE if tria needs "unwinding"
            smin = +.70 ;
            smax = +.90 ;
            sdel = .025 ;

            stol = smin+iter*sdel;
            stol = min (smax,stol) ;

            btri = nscr <= stol ...
                 & nscr <  oscr ;

            if (~any(btri)), break; end

        %-- relax toward old vert. coord's
            ivrt = ...
               unique(tria(btri,1:3));

            bvrt = ...
               false(size(vert,1),1) ;
            bvrt(ivrt) = true;

            if (undo ~= umax)
                bnew =  +.75 ^ undo ;
                bold =  +1.0 - bnew ;
            else
                bnew =  +0.0 ;
                bold =  +1.0 - bnew ;
            end

            vert(bvrt,:) = ...
                bold * vold(bvrt,:) ...
              + bnew * vert(bvrt,:) ;

            btri = any( ...
               bvrt(tria(:,1:3)),2) ;

        end

        oscr = nscr ;

        tcpu.undo = ...
            tcpu.undo + toc(ttic) ;

    %------------------------------------- test convergence!
        ttic = tic ;

        vdel = ...
            sum((vert-vold).^2,2) ;

        evec = vert(edge(:,2),:) ...
             - vert(edge(:,1),:) ;
        elen = ...
            sqrt(sum(evec.^2,2)) ;

        hvrt = evalhfn(vert, ...
            edge,EMAT,hfun,harg) ;

        hmid = hvrt(edge(:,1),:) ...
             + hvrt(edge(:,2),:) ;
        hmid = hmid * 0.5 ;
        scal = elen./hmid ;

        emid = vert(edge(:,1),:) ...
             + vert(edge(:,2),:) ;
        emid = emid * 0.5 ;

    %------------------------------------- |deg|-based prune
        keep = false(size(vert,1),1);
        keep(vdeg>+4) = true ;

        keep(conn(:)) = true ;
        keep(free(:)) = true ;

    %------------------------------------- 'density' control
        lmax = +5. / +4.  ;
        lmin = +1. / lmax ;

        less = scal<=lmin ;
    	more = scal>=lmax ;

        vbnd = false(size(vert,1),1);
        vbnd(conn(:,1)) = true ;
        vbnd(conn(:,2)) = true ;

        ebad = vbnd(edge(:,1)) ...     %-- not at boundaries
             | vbnd(edge(:,2)) ;

        less(ebad(:)) = false;
        more(ebad(:)) = false;

    %------------------------------------- force as disjoint
        lidx = find (less) ;

        for lpos = 1 : length(lidx)
            epos = lidx(lpos,1);
            inod = edge(epos,1);
            jnod = edge(epos,2);
        %--------------------------------- if still disjoint
            if (keep(inod) && ...
                keep(jnod) )

            keep(inod) = false ;
            keep(jnod) = false ;

            else

            less(epos) = false ;

            end
        end

        ebad = ...
           keep(edge(less,1)) ...
         & keep(edge(less,2)) ;

        more(ebad(:)) = false ;

    %------------------------------------- reindex vert/tria
        redo = ...
            zeros(size(vert,1),1);
        itop = ...
            length(find(keep));
        iend = ...
            length(find(less));

        redo(keep) = (1:itop)';

        redo(edge(less,1)) = ...        %-- to new midpoints
            (itop+1 : itop+iend)';
        redo(edge(less,2)) = ...
            (itop+1 : itop+iend)';

        vnew =[vert(keep,:) ;
               emid(less,:) ;
              ] ;
        tnew = reshape( ...
            redo(tria(:,+1:3)),[],3) ;

        ttmp = sort(tnew,2) ;           %-- filter collapsed
        okay = all( ...
            diff(ttmp,1,2)~=0,2) ;
        okay = ...
            okay & ttmp(:,1) > 0 ;
        tnew = tnew(okay,:) ;

    %------------------------------------- quality preserver
        nscr = ...
            triscr2  (vnew,tnew) ;

        stol = +0.80 ;

        tbad = nscr < stol ...
             & nscr < oscr(okay) ;

        vbad = ...
            false(size(vnew,1),1);
        vbad(tnew(tbad,:)) = true;

    %------------------------------------- filter edge merge
        lidx = find (less) ;

        ebad = ...
        vbad(redo(edge(lidx,1))) | ...
        vbad(redo(edge(lidx,2))) ;

        less(lidx(ebad)) = false ;

        keep(edge(...
        lidx(ebad),1:2)) =  true ;

    %------------------------------------- reindex vert/conn
        redo = ...
            zeros(size(vert,1),1);
        itop = ...
            length(find(keep));
        iend = ...
            length(find(less));

        redo(keep) = (1:itop)';

        redo(edge(less,1)) = ...
            (itop+1 : itop+iend)';
        redo(edge(less,2)) = ...
            (itop+1 : itop+iend)';

        vert =[vert(keep,:);
               emid(less,:);
               emid(more,:);
              ] ;
        conn = reshape( ...
            redo(conn(:,+1:2)),[],2) ;

        tcpu.keep = ...
            tcpu.keep + toc(ttic) ;

    %------------------------------------- build current CDT
        ttic = tic ;

       [vert,conn,tria,tnum] = ...
            deltri2 (vert, ...
                conn,node,PSLG,part) ;

        tcpu.dtri = ...
            tcpu.dtri + toc(ttic) ;

    %------------------------------------- dump-out progess!
        vdel = vdel./(hvrt.*hvrt) ;
        move = vdel > opts.vtol^2 ;
        nmov = ...
            length(find(move));

        ntri = size(tria,1) ;

        if (mod(iter,opts.disp)==+0)
            fprintf(+1, ...
            '%11i %18i %18i\n', ...
            [iter,nmov,ntri]) ;
        end

    %------------------------------------- loop convergence!
        if (nmov == +0), break; end

    end

    tria = tria( :,1:3);

%----------------------------------------- prune unused vert
    keep = false(size(vert,1),1);
    keep(tria(:)) = true ;
    keep(conn(:)) = true ;

    redo = zeros(size(vert,1),1);
    redo(keep) = ...
        (+1:length(find(keep)))';

    conn = ...
        reshape(redo(conn),[],2);
    tria = ...
        reshape(redo(tria),[],3);

    vert = vert(keep,:);

    tcpu.full = ...
        tcpu.full + toc(tnow) ;

    if (opts.dbug)
%----------------------------------------- print debug timer
        fprintf(1,'\n') ;
        fprintf(1,' Mesh smoothing timer...\n');
        fprintf(1,'\n') ;
        fprintf(1, ...
        ' FULL: %f \n', tcpu.full);
        fprintf(1, ...
        ' DTRI: %f \n', tcpu.dtri);
        fprintf(1, ...
        ' TCON: %f \n', tcpu.tcon);
        fprintf(1, ...
        ' ITER: %f \n', tcpu.iter);
        fprintf(1, ...
        ' UNDO: %f \n', tcpu.undo);
        fprintf(1, ...
        ' KEEP: %f \n', tcpu.keep);
        fprintf(1,'\n') ;
    end

    if (~isinf(opts.disp)), fprintf(1,'\n'); end

end

function [hvrt] = evalhfn(vert,edge,EMAT,hfun,harg)
%EVALHFN eval. the spacing-fun. at mesh vertices.

    if (~isempty (hfun))
        if (isnumeric(hfun))
            hvrt = hfun * ...
              ones(size(vert,1),1) ;
        else
            hvrt = feval( ...
                hfun,vert,harg{:}) ;
        end
    else

%-- no HFUN - HVRT is mean edge-len. at vertices!
        evec = vert(edge(:,2),:) - ...
               vert(edge(:,1),:) ;
        elen = sqrt(sum(evec.^2,2)) ;

        hvrt = (EMAT*elen) ...
            ./ max(sum(EMAT,2),eps) ;

        free = true(size(vert,1),1) ;
        free(edge(:,1)) = false;
        free(edge(:,2)) = false;

        hvrt(free)      =  +inf;

    end

end

function [opts] = makeopt(opts)
%MAKEOPT setup the options structure for SMOOTH2.

    if (~isfield(opts,'iter'))
        opts.iter = +32;
    else
    if (~isnumeric(opts.iter))
        error('smooth2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.iter)~= +1)
        error('smooth2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.iter <= +0)
        error('smooth2:invalidOptionValues', ...
            'Invalid OPT.ITER selection.') ;
    end
    end

    if (~isfield(opts,'disp'))
        opts.disp = + 4;
    else
    if (~isnumeric(opts.disp))
        error('smooth2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.disp)~= +1)
        error('smooth2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.disp <= +0)
        error('smooth2:invalidOptionValues', ...
            'Invalid OPT.DISP selection.') ;
    end
    end

    if (~isfield(opts,'vtol'))
        opts.vtol = +1.0E-02;
    else
    if (~isnumeric(opts.vtol))
        error('smooth2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.vtol)~= +1)
        error('smooth2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    if (opts.vtol <= 0.)
        error('smooth2:invalidOptionValues', ...
            'Invalid OPT.VTOL selection.') ;
    end
    end

    if (~isfield(opts,'dbug'))
        opts.dbug = false;
    else
    if (~islogical(opts.dbug))
        error('refine2:incorrectInputClass', ...
            'Incorrect input class.');
    end
    if (numel(opts.dbug)~= +1)
        error('refine2:incorrectDimensions', ...
            'Incorrect input dimensions.') ;
    end
    end

end


text:
classdef BT8 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP with bias feature

%------------------------------- Reference --------------------------------
% H. Li, Q. Zhang, and J. Deng, Biased multiobjective optimization and
% decomposition algorithm, IEEE Transactions on Cybernetics, 2017, 47(1):
% 52-66.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            [N,D] = size(X);
            I1    = 2 : 2 : D;
            I2    = 3 : 2 : D;
            Y     = X - repmat(X(:,1),1,D).^(0.5+1.5*repmat(0:D-1,N,1)/(D-1));
            DY    = Y.^2 + (1-exp(-Y.^2/1e-3))/5;
            PopObj(:,1) = X(:,1)         + sum(4*DY(:,I1).^2-cos(8*pi*DY(:,I1))+1,2);
            PopObj(:,2) = 1-sqrt(X(:,1)) + sum(4*DY(:,I2).^2-cos(8*pi*DY(:,I2))+1,2);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - R(:,1).^0.5;
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function A = spm_fnirs_sensitivity(DCM)
% Calculate sensitivity matrix which corresponds to the effective
% pathlength of detected photons for the channel measurements in the
% hemodynamic source.
% FORMAT [A] = spm_fnirs_sensitivity(DCM)
%
% DCM - DCM structure or its filename
%
% A - sensitivity matrix
%
% Green's function (see \dcm_fnirs\mmclab\estimate_greens_mmclab.m)
%--------------------------------------------------------------------------
% G.s - estimated Green's function from sensor (light emitter) positions
% into source positions [# sensor x # voxels x # wavelengths]
% G.d - estimated Green's function from sensor (light detector) positions
% into source positions [# sensor x # voxels x # wavelengths]
% G.xyz - MNI locations [3 x # voxels]
% G.elem - tissue types of voxels [3 x # voxels]
% 1-scalp, 2-CSF, 3-gray matter, 4-white matter
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Will Penny & Sungho Tak
% $id$

%--------------------------------------------------------------------------
if ~nargin || isempty(DCM)
    [DCM, sts] = spm_select(1, '^DCM*.*\.mat$', 'Select DCM.mat');
    if ~sts, S = []; return; end
elseif ischar(DCM)
    load(DCM);
end

%--------------------------------------------------------------------------
load(DCM.Y.P.fname.pos); % load channel positions
load(DCM.Y.P.fname.g);   % load Green's function

xyzh = [DCM.xY.xyz];     % hemodynamic source positions
xyzd = R.d.xyzH;         % detector positions
ns = R.s.ns;             % number of sources (emitter)
nd = R.d.nd;             % number of detectors
nh = size(xyzh, 2);      % number of hemodynamic sources

nnode = size(G.xyz, 2);  % number of nodes
r = DCM.options.rs .* 3; % maximum distance between point and distributed sources
rois = DCM.Y.P.rois;

%--------------------------------------------------------------------------
%-Find indices of Green's function of detector and hemodynamic source
node_d = zeros(nd, 1);
for i = 1:nd
    diff = G.xyz - xyzd(:, i) * ones(1, nnode);
    diff = sum(diff.^2);
    node_d(i, 1) = find(diff == min(diff));
end

if isfield(G, 'elem')
    [indx] = find(G.elem(5,:) == 3); % grey matter
    elem = G.elem(1:4, indx);
    node_g = sort(unique(elem(:))); % indx_node_grey
else
    node_g = 1:size(G.xyz, 2);
end

ng = length(node_g);
node_h = [];
A = [];

for i = 1:nh
    diff = G.xyz(:, node_g) - xyzh(:,i) * ones(1, ng);
    diff = sqrt(sum(diff.^2));

    if r == 0,
        A.node_h(i, 1) = node_g(find(diff == min(diff)));
    else
        indx_r = find(diff < r);
        [dist_h, indx_h] = sort(diff(indx_r));
        node_h = node_g(indx_r(indx_h));

        A(i).dist_h = dist_h;
        A(i).node_h = node_h;
    end
end

% Calculate sensitivity matrix using Green's function
%--------------------------------------------------------------------------
nwav = size(G.s, 3); % number of wavelengths
nroi = length(rois);

for i = 1:nwav
    niter = size(A, 2);
    for j = 1:niter
        S = [];
        for k = 1:nroi,
            indx_s = find(R.s.label == R.ch_sd(rois(k), 2));
            indx_d = find(R.d.label == R.ch_sd(rois(k), 3));
            S(k,:) = (G.s(indx_s, A(j).node_h,i) .* G.d(indx_d, A(j).node_h,i))./(G.s(indx_s, node_d(indx_d),i));
        end
        A(j).(sprintf('S%d',i)) = S;
    end
end

author: spm


text:
function pred_value = xval_simple_ols_loo(X, Y, pthreshold)
% Check: a very simple leave-one-out cross-validated regression
% pred_value = xval_simple_ols_loo(X, Y, p-value selection for univariate feature selection)
%
% pred_value: cross-validated predictions of outcome data
% X: n x variables matrix of predictors
% Y: n x 1 vector of outcomes
%
% Tor Wager, June 2010
%
% Go to any LASSO output directory and run this:
%
% maskInfo = iimg_read_img(fullfile(pwd, 'mask.img'), 2);
% dat{1} =  iimg_get_data(maskInfo, anticimages);
% pred_value = xval_simple_ols_loo(X, Y)

[N, k] = size(X);
pred_value = zeros(N, 1);

create_figure('test', 1, 2);

fprintf('Fold: %03d', 0);
for i = 1:N

    fprintf('\b\b\b%3.0f ', i);

    % select training data
    Xi = X; Xi(i, :) = [];
    Yi = Y; Yi(i) = [];

    % select features based on univariate correlations
    [r, p, Tstat] = correlation_fast_series(Xi, Yi);
    wh_features = p <= pthreshold;
    nfeatures(i) = sum(wh_features);
    if nfeatures(i) == 0
        disp('Warning: no features pass threshold');
        wh_features = p <= prctile(p, 10);
    else
        Xi = Xi(:, wh_features);
    end

    Xi = [ones(N-1, 1) Xi]; % add intercept

    % Make prediction
    b = pinv(Xi) * Yi;
    pred_value(i, 1) = [1 X(i, wh_features)] * b;

    create_figure('test', 1, 2, 1); subplot(1, 2, 1); plot(Xi*b, Yi, 'ko');
    plot(pred_value(i, 1), Y(i), 'ro', 'MarkerFaceColor', 'r');
    subplot(1, 2, 2);
   title('Black circles = training set; Red = holdout obs');
    drawnow;

end

cm = colormap(jet(N));
figure; hold on;
for i = 1:N
    plot(pred_value(i), Y(i), 'ko', 'MarkerFaceColor', cm(i, :));
end
xlabel('Predicted outcome (xval)'); ylabel('Outcome');
title('Color = order in data series');

% figure; hold on;
% for i = 1:N
%     plot(intercept_vals(i), pred_value(i), 'ko', 'MarkerFaceColor', cm(i, :));
% end

%
% figure; hold on;
%     plot(pred_value(whorder(1:24)), Y(whorder(1:24)), 'bo', 'MarkerFaceColor', 'b');
% plot(pred_value(whorder(25:end)), Y(whorder(25:end)), 'ro', 'MarkerFaceColor', 'r');

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore