text:
function c = tapas_ehgf_ar1_binary_config
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Contains the configuration for the enhanced Hierarchical Gaussian Filter (eHGF)
% for binary inputs in the absence of perceptual uncertainty.
%
% The HGF is the model introduced in
%
% Mathys C, Daunizeau J, Friston, KJ, and Stephan KE. (2011). A Bayesian foundation
% for individual learning under uncertainty. Frontiers in Human Neuroscience, 5:39.
%
% The binary HGF model has since been augmented with a positive factor kappa1 which
% scales the second level with respect to the first, i.e., the relation between the
% first and second level is
%
% p(x1=1|x2) = s(kappa1*x2), where s(.) is the logistic sigmoid.
%
% By default, kappa1 is fixed to 1, leading exactly to the model introduced in
% Mathys et al. (2011).
%
% This file refers to BINARY inputs (Eqs 1-3 in Mathys et al., (2011));
% for continuous inputs, refer to tapas_ehgf_config.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The HGF configuration consists of the priors of parameters and initial values. All priors are
% Gaussian in the space where the quantity they refer to is estimated. They are specified by their
% sufficient statistics: mean and variance (NOT standard deviation).
%
% Quantities are estimated in their native space if they are unbounded (e.g., the omegas). They are
% estimated in log-space if they have a natural lower bound at zero (e.g., the sigmas).
%
% Parameters can be fixed (i.e., set to a fixed value) by setting the variance of their prior to
% zero. Aside from being useful for model comparison, the need for this arises whenever the scale
% and origin at the j-th level are arbitrary. This is the case if the observation model does not
% contain the representations mu_j and sigma_j. A choice of scale and origin is then implied by
% fixing the initial value mu_j_0 of mu_j and either kappa_j-1 or omega_j-1.
%
% Fitted trajectories can be plotted by using the command
%
% >> tapas_ehgf_binary_plotTraj(est)
%
% where est is the stucture returned by tapas_fitModel. This structure contains the estimated
% perceptual parameters in est.p_prc and the estimated trajectories of the agent's
% representations (cf. Mathys et al., 2011). Their meanings are:
%
%         est.p_prc.mu_0       row vector of initial values of mu (in ascending order of levels)
%         est.p_prc.sa_0       row vector of initial values of sigma (in ascending order of levels)
%         est.p_prc.rho        row vector of rhos (representing drift; in ascending order of levels)
%         est.p_prc.ka         row vector of kappas (in ascending order of levels)
%         est.p_prc.om         row vector of omegas (in ascending order of levels)
%
% Note that the first entry in all of the row vectors will be NaN because, at the first level,
% these parameters are either determined by the second level (mu_0 and sa_0) or undefined (rho,
% kappa, and omega).
%
%         est.traj.mu          mu (rows: trials, columns: levels)
%         est.traj.sa          sigma (rows: trials, columns: levels)
%         est.traj.muhat       prediction of mu (rows: trials, columns: levels)
%         est.traj.sahat       precisions of predictions (rows: trials, columns: levels)
%         est.traj.v           inferred variance of random walk (rows: trials, columns: levels)
%         est.traj.w           weighting factors (rows: trials, columns: levels)
%         est.traj.da          volatility prediction errors  (rows: trials, columns: levels)
%         est.traj.ud          updates with respect to prediction  (rows: trials, columns: levels)
%         est.traj.psi         precision weights on prediction errors  (rows: trials, columns: levels)
%         est.traj.epsi        precision-weighted prediction errors  (rows: trials, columns: levels)
%         est.traj.wt          full weights on prediction errors (at the first level,
%                                  this is the learning rate) (rows: trials, columns: levels)
%
% Note that in the absence of sensory uncertainty (which is the assumption here), the first
% column of mu, corresponding to the first level, will be equal to the inputs. Likewise, the
% first column of sa will be 0 always.
%
% Tips:
% - When analyzing a new dataset, take your inputs u and use
%
%   >> est = tapas_fitModel([], u, 'tapas_ehgf_binary_config', 'tapas_bayes_optimal_binary_config');
%
%   to determine the Bayes optimal perceptual parameters (given your current priors as defined in
%   this file here, so choose them wide and loose to let the inputs influence the result). You can
%   then use the optimal parameters as your new prior means for the perceptual parameters.
%
% - If you get an error saying that the prior means are in a region where model assumptions are
%   violated, lower the prior means of the omegas, starting with the highest level and proceeding
%   downwards.
%
% - Alternatives are lowering the prior means of the kappas, if they are not fixed, or adjusting
%   the values of the kappas or omegas, if any of them are fixed.
%
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2020 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Config structure
c = struct;

% Model name
c.model = 'ehgf_ar1_binary';

% Number of levels (minimum: 3)
c.n_levels = 3;

% Input intervals
% If input intervals are irregular, the last column of the input
% matrix u has to contain the interval between inputs k-1 and k
% in the k-th row, and this flag has to be set to true
c.irregular_intervals = false;

% Sufficient statistics of Gaussian parameter priors

% Initial mus and sigmas
% Format: row vectors of length n_levels
% For all but the first two levels, this is usually best
% kept fixed to 1 (determines origin on x_i-scale). The
% first level is NaN because it is determined by the second,
% and the second implies neutrality between outcomes when it
% is centered at 0.
c.mu_0mu = [NaN, 0, 1];
c.mu_0sa = [NaN, 0, 0];

c.logsa_0mu = [NaN,   log(0.1), log(1)];
c.logsa_0sa = [NaN,          0,      0];

% Phis
% Format: row vector of length n_levels.
% Undefined (therefore NaN) at the first level.
% Fix this to zero (-Inf in logit space) to set to zero.
c.logitphimu = [NaN, -Inf, tapas_logit(0.5,1)]; %tapas_logit(0.1,1)
c.logitphisa = [NaN,    0,                  2];

% ms
% Format: row vector of length n_levels.
% This should be fixed for all levels where the omega of
% the next lowest level is not fixed because that offers
% an alternative parametrization of the same model.
c.mmu = [NaN, c.mu_0mu(2), c.mu_0mu(3)];
c.msa = [NaN,           0,           1];

% Rhos
% Format: row vector of length n_levels.
% Undefined (therefore NaN) at the first level.
% Fix this to zero to turn off drift.
c.rhomu = [NaN, 0, 0];
c.rhosa = [NaN, 0, 0];

% Kappas
% Format: row vector of length n_levels-1.
% Fixing log(kappa1) to log(1) leads to the original HGF model.
% Higher log(kappas) should be fixed (preferably to log(1)) if the
% observation model does not use mu_i+1 (kappa then determines the
% scaling of x_i+1).
c.logkamu = [log(1), log(1)];
c.logkasa = [     0,      0];

% Omegas
% Format: row vector of length n_levels.
% Undefined (therefore NaN) at the first level.
c.ommu = [NaN,  -3,   2];
c.omsa = [NaN,   4,   4];

% Gather prior settings in vectors
c.priormus = [
    c.mu_0mu,...
    c.logsa_0mu,...
    c.logitphimu,...
    c.mmu,...
    c.rhomu,...
    c.logkamu,...
    c.ommu,...
         ];

c.priorsas = [
    c.mu_0sa,...
    c.logsa_0sa,...
    c.logitphisa,...
    c.msa,...
    c.rhosa,...
    c.logkasa,...
    c.omsa,...
         ];

% Check whether we have the right number of priors
expectedLength = 5*c.n_levels+2*(c.n_levels-1)+1;
if length([c.priormus, c.priorsas]) ~= 2*expectedLength
    error('tapas:hgf:PriorDefNotMatchingLevels', 'Prior definition does not match number of levels.')
end

% Model function handle
c.prc_fun = @tapas_ehgf_ar1_binary;

% Handle to function that transforms perceptual parameters to their native space
% from the space they are estimated in
c.transp_prc_fun = @tapas_ehgf_ar1_binary_transp;

end

author: translationalneuromodeling

repo: tapas


text:
function sftpack_test05 ( )

%*****************************************************************************80
%
%% TEST05 tests R8VEC_SQSTB and R8VEC_SQSTF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2010
%
%  Author:
%
%    John Burkardt
%
  n = 256;
  alo = 0.0;
  ahi = 5.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST05\n' );
  fprintf ( 1, '  For real slow quarter wave sine transforms,\n' );
  fprintf ( 1, '  R8VEC_SQSTF does a forward transform;\n' );
  fprintf ( 1, '  R8VEC_SQSTB does a backward transform.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The number of data items is N = %d\n', n );
%
%  Set the data values.
%
  seed = 123456789;

  [ x, seed ] = r8vec_uniform ( n, alo, ahi, seed );

  r8vec_print_part ( n, x, 10, '  The original data:' );
%
%  Compute the coefficients.
%
  y = r8vec_sqstf ( n, x );

  r8vec_print_part ( n, y, 10, '  The sine coefficients:' );
%
%  Now compute inverse transform of coefficients.  Should get back the
%  original data.

  x = r8vec_sqstb ( n, y );

  r8vec_print_part ( n, x, 10, '  The retrieved data:' );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function m = uminus(k)
%  kronMatrix/uminus
%
%     -k negates the elements of k.
%
%     -kron(k) == kron(-k)

m = k;
if (length(m.a) > length(m.b))
   m.b = -(m.b);
else
   m.a = -(m.a);
end

repo: IRtools


text:
function segment_length = p09_boundary_segment_length ( segment_index, h )

%*****************************************************************************80
%
%% P09_BOUNDARY_SEGMENT_LENGTH returns boundary segment lengths in problem 09.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson and Gilbert Strang,
%    A Simple Mesh Generator in MATLAB,
%    SIAM Review,
%    Volume 46, Number 2, June 2004, pages 329-345.
%
%  Parameters:
%
%    Input, integer SEGMENT_INDEX, the index of one of the boundary segments.
%
%    Input, real H, the suggested spacing between points.
%
%    Output, integer SEGMENT_LENGTH, the number of points in the segment.
%
  if ( h <= 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P09_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );
    fprintf ( 1, '  Nonpositive H = %f\n', h );
    error ( 'P09_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );
  end

  if ( segment_index == 1 )

    n = round ( 4.0 / h );
    n = max ( n, 5 );
    segment_length = n + mod ( 4 - mod ( n - 1, 4 ), 4 );

  elseif ( segment_index == 2 )

    n = round ( 0.6 / h );
    n = max ( n, 7 );
    segment_length = n + mod ( 6 - mod ( n - 1, 6 ), 6 );

  elseif ( segment_index == 3 )

    n = round ( 0.6 / h );
    n = max ( n, 7 );
    segment_length = n + mod ( 6 - mod ( n - 1, 6 ), 6 );

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'P09_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );
    fprintf ( 1, '  Illegal SEGMENT_INDEX = %d\n', segment_index );
    error ( 'P09_BOUNDARY_SEGMENT_LENGTH - Fatal error!' );

  end

  return
end


text:
% Test file for chebtech/minandmax.m

function pass = test_minandmax(pref)

% Get preferences.
if ( nargin < 1 )
    pref = chebtech.techPref();
end

for n = 1:2
    if ( n == 1 )
        testclass = chebtech1();
    else
        testclass = chebtech2();
    end

    %%
    % Spot-check the extrema for a few functions.

    pass(n, 1) = test_spotcheck_minmax(testclass, @(x) ...
        ((x-0.2).^3 -(x-0.2) + 1).*sec(x-0.2), ...
        0.710869767377087, 1.884217141925336, pref);
    pass(n, 2) = test_spotcheck_minmax(testclass, @(x) sin(10*x), -1, 1, pref);
    pass(n, 3) = test_spotcheck_minmax(testclass, @airy, airy(1), airy(-1), ...
        pref);
    pass(n, 4) = test_spotcheck_minmax(testclass, @(x) -1./(1 + x.^2), -1, ...
        -0.5, pref);
    pass(n, 5) = test_spotcheck_minmax(testclass, ...
        @(x) (x - 0.25).^3.*cosh(x), ...
        (-1.25)^3*cosh(-1), 0.75^3*cosh(1), pref);

    %%
    % Check operation for array-valued inputs.

    fun_op = @(x) [sin(10*x) airy(x) (x - 0.25).^3.*cosh(x)];
    f = testclass.make(fun_op, [], pref);
    [y, x] = minandmax(f);
    y_exact = [-1 airy(1)  (-1.25)^3*cosh(-1);
                1 airy(-1) 0.75^3*cosh(1)];

    pass(n, 6) = all(abs(y(:) - y_exact(:)) < 10*max(vscale(f)*eps));

    % Check that the points x are indeed extreme points of the function
    % operator.
    for k = 1:1:size(f.coeffs, 2)
        fx = fun_op(x(:, k));
        if ( max(abs(fx(:, k) - y_exact(:, k))) > 1e1*max(vscale(f)*eps) )
            pass(n, 6) = 0;
            break;
        end
    end

    % Test complex-array-valued CHEBTECH objects.
    f = testclass.make(@(x) [exp(sin(2*x)), 1i*cos(20*x)]);
    [vals, pos] = minandmax(f);
    f1 = testclass.make(@(x) exp(sin(2*x)));
    [vals1, pos1] = minandmax(f1);
    f2 = testclass.make(@(x) 1i*cos(20*x));
    [vals2, pos2] = minandmax(f2);
    pass(n, 7) = norm(abs(vals) - abs([vals1 vals2]), inf) < ...
        1e2*max(vscale(f)*eps);
    % Note, we don't expect pos(:,2) = pos2 as the min and max are not unique.
    pass(n, 8) = norm(pos(:,1) - pos1, inf) < 500*max(vscale(f)*eps);

end

end

% Spot-check the results for a given function.
function result = test_spotcheck_minmax(testclass, fun_op, exact_min, ...
    exact_max, pref)

f = testclass.make(fun_op, [], pref);
[y, x] = minandmax(f);
y_exact = [exact_min ; exact_max];
fx = fun_op(x);
result = ((max(abs(y - y_exact)) < 10*vscale(f)*eps) && ...
          (max(abs(fx - y_exact)) < 10*vscale(f)*eps));

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function [f]=ref_idwilt_1(coef,g,a,M)
%REF_IDWILT_1  Reference IDWILT by IDGT
%

%   Author : Peter L. Søndergaard

L=size(g,1);
N=L/a;
W=size(coef,2);

coef=reshape(coef,M*2,N/2,W);

coef2=zeros(2*M,N,W);

if 1

  % ----- Loop version ---

  for n=0:N/2-1

    % m=0
    coef2(1,2*n+1,:) = coef(1,n+1,:);

    for m=0:M-1

      % m odd
      for m=1:2:M-1
	coef2(m+1,2*n+1,:)     = -i/sqrt(2)*coef(m+1,n+1,:);
	coef2(m+1,2*n+2,:)     =  1/sqrt(2)*coef(M+m+1,n+1,:);

	coef2(2*M-m+1,2*n+1,:) =  i/sqrt(2)*coef(m+1,n+1,:);
	coef2(2*M-m+1,2*n+2,:) =  1/sqrt(2)*coef(M+m+1,n+1,:);
      end;

      % m even
      for m=2:2:M-1
	coef2(m+1,2*n+1,:)     =  1/sqrt(2)*coef(m+1,n+1,:);
	coef2(m+1,2*n+2,:)     = -i/sqrt(2)*coef(M+m+1,n+1,:);

	coef2(2*M-m+1,2*n+1,:) =  1/sqrt(2)*coef(m+1,n+1,:);
	coef2(2*M-m+1,2*n+2,:) =  i/sqrt(2)*coef(M+m+1,n+1,:);
      end;

    end;

    % m=nyquest
    if mod(M,2)==0
      coef2(M+1,2*n+1,:) = coef(M+1,n+1,:);
    else
      coef2(M+1,2*n+2,:) = coef(M+1,n+1,:);
    end;

  end;

else
  % ----- Vector version ---

  % First and middle modulation are transferred unchanged.
  coef2(1,1:2:N,:) = coef(1,:,:);
  if mod(M,2)==0
    coef2(M+1,1:2:N,:) = coef(M+1,:,:);
  else
    coef2(M+1,2:2:N,:) = coef(M+1,:,:);
  end;

  if M>2
    coef2(3:2:M,1:2:N,:)        = 1/sqrt(2)*coef(3:2:M,:,:);
    coef2(3:2:M,2:2:N,:)        = -1/sqrt(2)*i*coef(M+3:2:2*M,:,:);

    coef2(2*M-1:-2:M+2,1:2:N,:) = 1/sqrt(2)*coef(3:2:M,:,:);
    coef2(2*M-1:-2:M+2,2:2:N,:) =  1/sqrt(2)*i*coef(M+3:2:2*M,:,:);
  end;

  % sine, first column.
  coef2(2:2:M,1:2:N,:)        = -1/sqrt(2)*i*coef(2:2:M,:,:);
  coef2(2:2:M,2:2:N,:)        = 1/sqrt(2)*coef(M+2:2:2*M,:,:);

  coef2(2*M:-2:M+2,1:2:N,:)   =  1/sqrt(2)*i*coef(2:2:M,:,:);
  coef2(2*M:-2:M+2,2:2:N,:)   = 1/sqrt(2)*coef(M+2:2:2*M,:,:);

end;

f=idgt(reshape(coef2,2*M,N,W),g,a);

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function obj = coef_Jacob3_qt_hand_eye_func_new(in1)
coef_J3 = in1(:,3);
coef_J6 = in1(:,6);
coef_J8 = in1(:,8);
coef_J9 = in1(:,9);
coef_J16 = in1(:,16);
coef_J18 = in1(:,18);
coef_J19 = in1(:,19);
coef_J25 = in1(:,25);
coef_J26 = in1(:,26);
coef_J27 = in1(:,27);
coef_J28 = in1(:,28);
coef_J29 = in1(:,29);
coef_J30 = in1(:,30);
coef_J31 = in1(:,31);
coef_J38 = in1(:,38);
coef_J40 = in1(:,40);
coef_J41 = in1(:,41);
coef_J47 = in1(:,47);
coef_J48 = in1(:,48);
coef_J49 = in1(:,49);
coef_J50 = in1(:,50);
coef_J51 = in1(:,51);
coef_J52 = in1(:,52);
coef_J53 = in1(:,53);
coef_J59 = in1(:,59);
coef_J60 = in1(:,60);
coef_J61 = in1(:,61);
coef_J62 = in1(:,62);
coef_J63 = in1(:,63);
coef_J64 = in1(:,64);
coef_J65 = in1(:,65);
coef_J66 = in1(:,66);
coef_J67 = in1(:,67);
coef_J68 = in1(:,68);
coef_J69 = in1(:,69);
coef_J70 = in1(:,70);
obj = [coef_J3,coef_J6,coef_J8.*2.0,coef_J9,coef_J16,coef_J18.*2.0,coef_J19,coef_J25.*3.0,coef_J26.*2.0,coef_J27,coef_J28,coef_J29,coef_J30,coef_J31,coef_J38,coef_J40.*2.0,coef_J41,coef_J47.*3.0,coef_J48.*2.0,coef_J49,coef_J50,coef_J51,coef_J52,coef_J53,coef_J59.*4.0,coef_J60.*3.0,coef_J61.*2.0,coef_J62.*2.0,coef_J63.*2.0,coef_J64.*2.0,coef_J65.*2.0,coef_J66,coef_J67,coef_J68,coef_J69,coef_J70];

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function I = eye(d)
%EYE       Identity operator.
%   This function is deprecated. Use OPERATORBLOCK.EYE.

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

I = linop( operatorBlock.eye(double(d)) );

end


text: %INVSIG The inverse sigmoid of a dataset

author: marianux

repo: ecg-kit

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


text:
classdef TestPerspectiveTransform
    %TestPerspectiveTransform

    methods (Static)
        function test_2d_numeric
            M = [cos(pi/4) sin(pi/4) 0 ; ...
                -sin(pi/4) cos(pi/4) 0 ; ...
                        0         0  1];
            % Nx2x1
            src = randn(10,2);
            dst = cv.perspectiveTransform(src, M);
            validateattributes(dst, {class(src)}, {'size',size(src)});
            % 1xNx2
            src = shiftdim(randn(10,2), -1);
            dst = cv.perspectiveTransform(src, M);
            validateattributes(dst, {class(src)}, {'size',size(src)});
            % Nx1x2
            src = permute(randn(10,2), [1 3 2]);
            dst = cv.perspectiveTransform(src, M);
            validateattributes(dst, {class(src)}, {'size',size(src)});
        end

        function test_3d_numeric
            M = eye(4);
            % Nx3
            src = randn(10,3);
            dst = cv.perspectiveTransform(src, M);
            validateattributes(dst, {class(src)}, {'size',size(src)});
            % 1xNx3
            src = permute(randn(10,3), [3 1 2]);
            dst = cv.perspectiveTransform(src, M);
            validateattributes(dst, {class(src)}, {'size',size(src)});
            % Nx1x3
            src = permute(randn(10,3), [1 3 2]);
            dst = cv.perspectiveTransform(src, M);
            validateattributes(dst, {class(src)}, {'size',size(src)});
        end

        function test_2d_cell
            M = eye(3);
            pts = num2cell(randn(10,2), 2);
            dst = cv.perspectiveTransform(pts, M);
            validateattributes(dst, {'cell'}, {'vector', 'numel',numel(pts)});
            cellfun(@(pt) validateattributes(pt, {'numeric'}, ...
                {'vector', 'numel',2}), dst);
        end

        function test_3d_cell
            M = eye(4);
            pts = num2cell(randn(10,3), 2);
            dst = cv.perspectiveTransform(pts, M);
            cellfun(@(pt) validateattributes(pt, {'numeric'}, ...
                {'vector', 'numel',3}), dst);
        end

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

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function prob_test155 ( )

%*****************************************************************************80
%
%% TEST155 tests VON_MISES_MEAN, VON_MISES_SAMPLE, VON_MISES_CIRCULAR_VARIANCE;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  nsample = 1000;
  seed = 123456789;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST155\n' );
  fprintf ( 1, '  For the Von Mises PDF:\n' );
  fprintf ( 1, '  VON_MISES_MEAN computes the mean;\n' );
  fprintf ( 1, '  VON_MISES_SAMPLE samples.\n' );
  fprintf ( 1, ...
    '  VON_MISES_CIRCULAR_VARIANCE computes the circular variance;\n' );

  a = 1.0;
  b = 2.0;

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

  check = von_mises_check ( a, b );

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

  mean = von_mises_mean ( a, b );
  variance = von_mises_circular_variance ( a, b );

  fprintf ( 1, '  PDF mean =              %14f\n', mean );
  fprintf ( 1, '  PDF circular variance = %14f\n', variance );

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

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

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

  return
end


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% function B = getElasticMatrixStg(omega,m,mu,lambda)
%
% Builds the elasticity matrix B for a domain defined by omega
% and a staggered grid discretization defined by m:
%
%     | a\nabla       0        0 |
% B = |       0 a\nabla        0 |
%     |       0       0  a\nabla |
%     |           b \div         |,
%
% where \nabla and \div are weighted by a=\sqrt(mu) and b=\sqrt(\lambda+\mu),
% respectively; defaults: \mu = 1, \lambda = 0.
% Note: the matrices are simple, it is size that matters.
%
% =============================================================================

function B = getElasticMatrixStg(omega,m,mu,lambda)
if nargin == 0,
  help(mfilename);
  runMinimalExample;
  B = 'endOfMinimalExample';
  return;
end;

                                        % set defautls
if ~exist('mu','var'),     mu     = [];  end;
if ~exist('lambda','var'), lambda = [];  end;
if isempty(mu),
  warning('mu has not yet been set! set mu=1');
  mu = 1;
end;
if isempty(lambda),
  warning('lambda has not yet been set! lambda=0');
  lambda = 0;
end;

dim = length(omega)/2;
h   = (omega(2:2:end)-omega(1:2:end))./m; % voxel size for integration
a   = sqrt(mu);
b   = sqrt(mu+lambda);

% for the dimension ms(i) of the i-th staggered grid, the i-th component
% of of m=[m(1),m(2),m(3)] has to be increased by one; the i-th staggered
% grid contains ns(i) points

ms  = @(i)   m + ((1:dim)==i);
ns  = @(i)   prod(ms(i));
% Example: m=(2,3), ms(1) = (3,3); ms(2)=(2,4), ns(1) = 9, ns(2) = 8.

% the dimension of \partial^{h,i}_k is p(i,k)-by-ns(i),
% where q(i)=prod(ms(i)) is the length of the i-th staggered grid
% and p(i,k) is the product of a modified m:
%   m(i) has to be increase by one and
%   m(k) has to be decrase  by one
p = @(i,k) prod(m + ((1:dim)==i) - ((1:dim==k)));

% Example:
% p(1,1) = 2*3         = 6,   p(1,2) = (2+1)*(3-1) = 6, ns(1) = 9
% p(2,1) = (2-1)*(3+1) = 4,   p(2,2) = 2*3         = 6, ns(2) = 8

% E is a q-by-q identity matrix, where q=m(k)+1 for i==k and q=m(k) (i~=k)
E = @(i,k) speye(m(k)+(i==k));

% Example:
% E(1,1) = speye(m(1)+1) = speye(3)
% E(1,2) = speye(m(2))   = speye(3)
% E(2,1) = speye(m(1))   = speye(2)
% E(2,2) = speye(m(2)+1) = speye(4)

% D is the 1D derivative matrix of size p-by-(p+1)
%        |-1 1     |
% D = |  .  .   |/h(k), where p=m(k) (i==k) and p=m(k)-1 (i~=k)
%        |     -1 1|

D = @(i,k) spdiags(ones(m(k)-(i~=k),1)*[-1,1],...
  [0,1],m(k)-(i~=k),m(k)+(i==k))/h(k);

% Example:
% D(1,1) is m(1)-by-(m(1)+1) = 2-by-3
% D(1,2) is (m(2)-1)-by-m(2) = 2-by-3
% D(2,1) is (m(1)-1)-by-m(1) = 1-by-2
% D(2,2) is m(2)-by-(m(2)+1) = 3-by-4

switch dim
  case 2
    % build the 2D elasticity operator
    %
    %      | a\nabla 0      |   |a D11  0     |
    %  B = |                | = |a D12  0     |
    %      | 0      a\nabla |   |0      a D21 |
    %      |                |   |0      a D22 |
    %      | b Div1 b Div2  |   |b D11  b D22 |

    D11 = kron(E(1,2),D(1,1));  D12 = kron(D(1,2),E(1,1));
    D21 = kron(E(2,2),D(2,1));  D22 = kron(D(2,2),E(2,1));

    % Example:
    % D11 = 3-by-3 \otimes 2-by-3 = 6-by-9
    % D12 = 2-by-3 \otimes 3-by-3 = 6-by-9
    % D21 = 4-by-4 \otimes 1-by-2 = 4-by-8
    % D22 = 3-by-4 \otimes 2-by-2 = 6-by-8

    B = [  a*D11,sparse(p(1,1),ns(2));
           a*D12,sparse(p(1,2),ns(2));
           sparse(p(2,1),ns(1)),a*D21;
           sparse(p(2,2),ns(1)),a*D22;
           b*D11,b*D22];
    % Example: B = 28-by-17
  case 3
    % build the 3D elasticity operator
    %
    %      |a D11  0      0    |
    %  B = |a D12  0      0    |
    %      |a D13  0      0    |
    %      |0      a D21  0    |
    %      |0      a D22  0    |
    %      |0      a D23  0    |
    %      |0      0      a D31|
    %      |0      0      a D32|
    %      |0      0      a D33|
    %      |b D11  b D22  b D33|

    % Example: m=[2,3,4] ns = (36,32,30)
    D11 = kron(kron(E(1,3),E(1,2)),D(1,1));
    D12 = kron(kron(E(1,3),D(1,2)),E(1,1));
    D13 = kron(kron(D(1,3),E(1,2)),E(1,1));
    D21 = kron(kron(E(2,3),E(2,2)),D(2,1));
    D22 = kron(kron(E(2,3),D(2,2)),E(2,1));
    D23 = kron(kron(D(2,3),E(2,2)),E(2,1));
    D31 = kron(kron(E(3,3),E(3,2)),D(3,1));
    D32 = kron(kron(E(3,3),D(3,2)),E(3,1));
    D33 = kron(kron(D(3,3),E(3,2)),E(3,1));

    % Example:
    % D11 = kron(kron(4-by-4,3-by-3),2-by-3) = 24-by-36
    % D12 = kron(kron(4-by-4,2-by-3),3-by-3) = 24-by-36
    % D13 = kron(kron(3-by-4,3-by-3),3-by-3) = 27-by-36
    % D21 = kron(kron(4-by-4,4-by-4),1-by-2) = 16-by-32
    % D22 = kron(kron(4-by-4,3-by-4),2-by-2) = 24-by-32
    % D23 = kron(kron(3-by-4,4-by-4),2-by-2) = 24-by-32
    % D31 = kron(kron(5-by-5,3-by-3),1-by-2) = 15-by-30
    % D32 = kron(kron(5-by-5,2-by-3),2-by-2) = 20-by-30
    % D32 = kron(kron(4-by-5,3-by-3),2-by-2) = 24-by-30

    B = [
      a * D11, sparse(p(1,1),ns(2)), sparse(p(1,1),ns(3)); ...
      a * D12, sparse(p(1,2),ns(2)), sparse(p(1,2),ns(3)); ...
      a * D13, sparse(p(1,3),ns(2)), sparse(p(1,3),ns(3)); ...
      sparse(p(2,1),ns(1)), a * D21, sparse(p(2,1),ns(3)); ...
      sparse(p(2,2),ns(1)), a * D22, sparse(p(2,2),ns(3)); ...
      sparse(p(2,3),ns(1)), a * D23, sparse(p(2,3),ns(3)); ...
      sparse(p(3,1),ns(1)), sparse(p(3,1),ns(2)), a * D31; ...
      sparse(p(3,2),ns(1)), sparse(p(3,2),ns(2)), a * D32; ...
      sparse(p(3,3),ns(1)), sparse(p(3,3),ns(2)), a * D33; ...
      b * D11,              b * D22,              b * D33
      ];
    % Example:
    % B = [75-by-36    0      0
    %         0     64-by-32  0
    %         0        0     59-by-30
    %      24-by-36 24-by-32 24-by-30 ]
    %   = 228-by-98

  otherwise,  error('nyi');
end;
%------------------------------------------------------------------------------
function runMinimalExample

omega = [0,1,0,2,0,3]; m = [4,5,6];
B2 = getElasticMatrixStg(omega(1:4),m(1:2),1,0);
B3 = getElasticMatrixStg(omega,m,1,0);
figure(1); clf;
subplot(1,2,1); spy(B2); title(sprintf('%s: B2',mfilename));
subplot(1,2,2); spy(B3); title(sprintf('%s: B3',mfilename));
%==============================================================================


text:
function [ state ] = nonBlocking_movePTPJointSpace( t , jPos, relVel)
%% This function is used for performing non-blocking point to point motion in joint space,
% for the KUKA iiwa 7 R 800.

%% Syntax:
% [ state ] = nonBlocking_movePTPJointSpace( t , jPos, relVel)

%% About:
% This function is used to move the robot point to point in joint space.

%% Arreguments:
% t: is the TCP/IP connection
% jPos: is the destanation position in joint space, it is 1x7 cell array,
% joint angles are in radians.
% relVel: is a double, from zero to one, specifies the over-ride relative
% velocity.

% Copy right, Mohammad SAFEEA, 9th of May 2017

% check also the function: nonBlocking_isGoalReached()

    global paramName;
    global paramVal;
    paramName='jpos';
    paramVal=jPos;

    theCommand=['jRelVel_',num2str(relVel),'_']; % set over ride.
    fprintf(t, theCommand);
    message=fgets(t);

    sendJointsPositions( t ,jPos); % send destination joint positions.

    theCommand='doPTPinJS';
    fprintf(t, theCommand); % start the point to point motion.
    message=fgets(t);
    state=checkAcknowledgment(message);

end

repo: KST-Kuka-Sunrise-Toolbox


text:
function pstruct = tapas_ehgf_namep(pvec)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2020 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

pstruct = struct;

l = length(pvec)/5;

if l ~= floor(l)
    error('tapas:hgf:UndetNumLevels', 'Cannot determine number of levels');
end

pstruct.mu_0      = pvec(1:l);
pstruct.sa_0      = pvec(l+1:2*l);
pstruct.rho       = pvec(2*l+1:3*l);
pstruct.ka        = pvec(3*l+1:4*l-1);
pstruct.om        = pvec(4*l:5*l-1);
pstruct.pi_u      = pvec(5*l);

end


text:
% SERVER.GETREPLY gets a reply structure stored by the slave with given index.
function result = getReply(server, slaveIdx)

  if (slaveIdx > server.slaveCount)
    error('Slave index is out of bounds.');
  end

  if (server.options.emulate)
    % emulate the following:
    % 1.  deserialization of data by slave
    % 2.  processing of request by slave
    % 3.  serialization of reply by slave.
    % 4.  de-serialization of reply by server.

    % Deserialize request
    reply = emulateSlaveProcessRequest(server, slaveIdx);
  else
    if (slaveIdx > length(server.slave_sock))
      error('No socket for slave');
    end
    sock=server.slave_sock{slaveIdx};

    % this blocks until slave creates output stream
    is=java.io.ObjectInputStream(sock.getInputStream());
    % block until we get a new message;  convert to MATLAB
    reply = java2matlab(is.readObject());
  end

  result=reply.result;
  if (isfield(reply, 'error'))
    rethrow(reply.error);
  end

  % wrap any slaveRef objects in a serverRef...
  result = serverRef.wrapSlaveRefs(result, server, slaveIdx);

author: huashiyiqike

repo: LSTM-MATLAB


text:
function boolVal=pointIsInUBeam(points,uBounds,boundType)
%%POINTISINUBEAM A region in a single direction cosine in 2D space is
%      defined either in terms of a central direction and a half beamwidth
%      (the beamwidth being the distance between the bounds of the region)
%      or in terms of a starting direction cosine and an ending direction
%      cosine, determine whether provided points are within that span.
%      Because beams and points that go past the -1 to 1 region alias back
%      into that region, one cannot simply compare the values in the points
%      to those of the bounds.
%
%INPUTS: points A vector or matrix of direction cosines values to test
%               whether or not they are in the beam. nonaliased values are
%               between -1 and 1.
%       uBounds A 2X1 vector defining the bounds of the beam. If
%               boundType=0, then uBounds(1) is the direction of the
%               center of the beam and uBounds(2) is the half
%               beamwidth of the beam. If boundType=1, then uBounds(1) is
%               the lower bound and azBounds(2) is the upper bound --where
%               the values are taken in increasing (aliased) order. If
%               bounds go outside of the beam span, then they are assumed
%               to alias back into the valid span (a beam split between two
%               directions 180 degrees apart).
%     boundType This is either 0 or 1 and affects the format of uBounds,
%               as described above. The default if omitted or an empty
%               matrix is passed is 0.
%
%OUTPUTS: boolVal This has the same size as the input and the entries are
%                 true if the corresponding entry in the input is in the
%                 beam and false if it is not in the beam.
%
%EXAMPLE:
%This plots a beam that is aliased across the -1,1 boundary. Then, points
%uniformly spaced across the -1 to 1 u span are tested for being in the
%beam. Those in the beam are plotted in red and all the other points are
%plotted in blue. One can see that the function correctly handles aliasing
%around the -1,1 discontinuity.
% range=[0;1];
% beamHalfwidth=0.1;
% closeEnds=true;
% figure(1)
% clf
% hold on
% uCenter=0.95;
% drawUBeam(uCenter,beamHalfwidth,range,[],[],closeEnds,'-k','linewidth',2)
%
% numPoints=1000;
% uPoints=linspace(-1,1,numPoints);
% boundType=0;
% uBounds=[uCenter;beamHalfwidth];
% boolVal=pointIsInUBeam(uPoints,uBounds,boundType);
%
% %Plot all the points at unit range in blue.
% points=[ones(1,numPoints);uPoints];
% useHalfRange=true;
% xyPts=ru2Cart2D(points,useHalfRange);
% scatter(xyPts(1,:),xyPts(2,:),100,'.b')
% %Plot the points that gate in red.
% points=points(:,boolVal);
% xyPts=ru2Cart2D(points,useHalfRange);
% scatter(xyPts(1,:),xyPts(2,:),200,'.r')
% axis([-1,1,-0.5,1])
%
%January 2023 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<3||isempty(boundType))
    boundType=0;
end

if(boundType==0)
    %If defined in terms of a beam center and a beamwidth....
    beamCenter=uBounds(1);
    beamHalfwidth=uBounds(2);

    if(beamHalfwidth>=1)
        boolVal=true(size(points));
        return;
    end

    u1=beamCenter-beamHalfwidth;
    u2=beamCenter+beamHalfwidth;
else
    %If the bounds are directly given.
    u1=uBounds(1);
    u2=uBounds(2);
end
%Alias everything so that u1=0. Everything is considered in terms of
%increasing u from 0.
u2=wrapRange(u2-u1,0,2);
points=wrapRange(points-u1,0,2);

%Test whether the direction cosines are in or on the edge of the beam.
boolVal=(points<=u2);

end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
%% (Internal) Get names of annotations from annotation structure
%
%   [AnnNames, all_annotations] = getAnnNames(aux_struct)
%
% Arguments:
%
%      + aux_struct:
%
%      + retries: times to check the existence
%
% Output:
%
%      + AnnNames:
%
%      + all_annotations:
%
% Example:
%
%
% See also ECGtask_QRS_detection
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 21/4/2015
% Copyright 2008-2015
%
function [AnnNames, all_annotations] = getAnnNames(aux_struct)

    AnnNames = [];

    for fname = rowvec(fieldnames(aux_struct))
        if( isfield(aux_struct.(fname{1}), 'time') )
            AnnNames = [AnnNames; cellstr(fname{1}) cellstr('time')];
        end
        if( isfield(aux_struct.(fname{1}), 'qrs') )
            AnnNames = [AnnNames; cellstr(fname{1}) cellstr('qrs')];
        end
    end

    cant_anns = size(AnnNames,1);

    all_annotations = cell(cant_anns,1);
    for ii = 1:cant_anns
        all_annotations{ii} = aux_struct.(AnnNames{ii,1}).(AnnNames{ii,2});
    end


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% Tutorial for FAIR: shows how to visulaize 3D data with FAIR.
%
% - load data             (here 3D MRI)
% - setup viewer          (viewImage2D)
% - view  data
%==============================================================================

clear, %close all,
help(mfilename),

% setup data
fprintf('%s\n','generate data (here: MRI)');
load mri;
omega = [0,20,0,20,0,10];    % specify physical domain
T     = double(squeeze(D));
m     = size(T);             % get size of data

%%
fprintf('%s\n','use image montage for visualization of data')
FAIRfigure(1); clf;
imgmontage(T,omega,m,'colormap',gray(100),'numbering','on');
title('visualization of 3D data - imgmontage','fontsize',30)

%%
fprintf('%s\n','use a volumetric view')
FAIRfigure(2); clf;
volView(T,omega,m,'facecolor',[240,140,100]/256,'facealpha',0.75); hold on;
colormap(gray(100))
%%
fprintf('%s\n','use a volumetric view + slices')
FAIRfigure(3); clf;
volView(T,omega,m,'facecolor',[240,140,100]/256,'facealpha',0.75); hold on;
colormap(gray(100))
vh = viewSlices(T,omega,m);

%%
FAIRfigure(4); clf; colordef(gcf,'black'); colormap(gray(100));
for j=1:m(3);
  vh = viewSlices(T,omega,m,'s1',64,'s2',64,'s3',[j]);
  FAIRpause(1/10);
end;
echo off
%==============================================================================


text:
function a = covar2 ( m, n, b, a )

%*****************************************************************************80
%
%% COVAR2 returns the covariance matrix for a rectangular matrix.
%
%  Discussion:
%
%    M should be at least N, and B has maximal rank N.
%    The covariance matrix is defined as
%
%      A = inverse ( B' * B ).
%
%  Properties:
%
%    A is symmetric: A' = A.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the order of B.
%    N <= M is required.
%
%    Input, real B(M,N), the matrix whose covariance matrix
%    is desired.
%
%    Output, real A(N,N), the covariance matrix of B.
%
%    Output, integer IERROR, error flag.
%    IERROR = 0 for no error,
%    IERROR is nonzero to mean that A is singular.
%
  if ( m < n )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'COVAR2 - Fatal error!\n' );
    fprintf ( 1, '  M < N\n' );o
    error ( 'COVAR2 - Fatal error!' );
  end
%
%  Compute B' * B.
%
  btb(1:n,1:n) = b(1:m,1:n)' * b(1:m,1:n);
%
%  Factor the matrix.
%
  [ btb, pivot, error ] = r8mat_gefa ( btb, n );

  if ( ierror ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'COVAR2 - Fatal error!\n' );
    fprintf ( 1, '  The matrix B'' * B does not have full rank.\n' );
    fprintf ( 1, '  This means that B does not have rank N.\n' );
    error ( 'COVAR2 - Fatal error!' );
  end
%
%  Solve for A by solving the equations ( B' * B ) * A = I
%  one column at a time.
%
  a = identity ( n, n );

  job = 0;

  for j = 1 : n
    a(1:n,j) = r8mat_gesl ( btb, n, pivot, a(1:n,j), job );
  end

  return
end


text:
function fiff_write_string(fid,kind,data)
%
% fiff_write_string(fid,kind,data)
%
% Writes a string tag
%
%     fid           An open fif file descriptor
%     kind          The tag kind
%     data          The string data to write
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%   Revision 1.2  2006/04/23 15:29:40  msh
%   Added MGH to the copyright
%
%   Revision 1.1  2006/04/10 23:26:54  msh
%   Added fiff reading routines
%
%

me='MNE:fiff_write_string';

if nargin ~= 3
        error(me,'Incorrect number of arguments');
end

FIFFT_STRING=10;
FIFFV_NEXT_SEQ=0;
datasize=size(data,2);
count = fwrite(fid,int32(kind),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(FIFFT_STRING),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(datasize),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,int32(FIFFV_NEXT_SEQ),'int32');
if count ~= 1
    error(me,'write failed');
end
count = fwrite(fid,data,'uchar');
if count ~= datasize
    error(me,'write failed');
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Spectral factorization using Kolmogorov 1939 approach.
% (code follows pp. 232-233, Signal Analysis, by A. Papoulis)
%
% Computes the minimum-phase impulse response which satisfies
% given auto-correlation.
%
% Input:
%   r: top-half of the auto-correlation coefficients
%      starts from 0th element to end of the auto-corelation
%      should be passed in as a column vector
% Output
%   h: impulse response that gives the desired auto-correlation

function h = spectral_fact(r)

% length of the impulse response sequence
n = length(r);

% over-sampling factor
mult_factor = 100;        % should have mult_factor*(n) >> n
m = mult_factor*n;

% computation method:
% H(exp(jTw)) = alpha(w) + j*phi(w)
% where alpha(w) = 1/2*ln(R(w)) and phi(w) = Hilbert_trans(alpha(w))

% compute 1/2*ln(R(w))
w = 2*pi*((0:m-1)/m);
R = [ ones(m,1) 2*cos(kron(w',1:n-1)) ]*r;
alpha = 1/2*log(R);

% find the Hilbert transform
alphatmp = fft(alpha);
alphatmp(floor(m/2)+1:m) = -alphatmp(floor(m/2)+1:m);
alphatmp(1) = 0;
alphatmp(floor(m/2)+1) = 0;
phi = real(ifft(1j*alphatmp));

% now retrieve the original sampling
index  = find(rem(0:m-1,mult_factor)==0);
alpha1 = alpha(index);
phi1   = phi(index);

% compute the impulse response (inverse Fourier transform)
h = real(ifft(exp(alpha1+1j*phi1),n));

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function mv = mv_exportMap(mv, param, saveFlag, varargin);
%
% mv = mv_exportMap(mv, [param], [saveFlag], [options]);
%
% For Multi-Voxel UI,  export a set of voxel parameter values
% as a map. [more]
%
%
% ras 08/05.
if notDefined('mv')
    mv = get(gcf, 'UserData');
end

if notDefined('param')
    % put up a dialog
    paramList = {'D-prime' 'D-prime-refcond' 'Preferred Condition' ...
                 'Selectivity Index' 'Voxel Reliability'  ...
                 'Mutual Information' 'Omnibus Reliability' ...
                 'Mean Voxel Amplitudes' 'Amplitudes for Each Condition' ...
                 'Contrast'};

    dlg.fieldName = 'param';
    dlg.style = 'popup';
    dlg.string = 'Export which parameter?';
    dlg.list = paramList;
    dlg.value = 1;

    dlg(2).fieldName = 'saveFlag';
    dlg(2).style = 'checkbox';
    dlg(2).string = 'Save Parameter To File';
    dlg(2).value = 1;

    resp = generalDialog(dlg, 'Export to Map...');
    if isempty(resp),  disp('mv_exportMap: User Aborted'); return; end
    param = resp.param;
    saveFlag = resp.saveFlag;
end

% Get the map name,  and values for the map
switch lower(param)
    case 'd-prime',
        mapName = 'd-prime';
        vals = mv_dprime(mv);
    case 'd-prime-refcond',
        mapName = 'd-prime-refcond';
        [vals mapName] = mv_dprime_cond(mv);
    case 'preferred condition',
        mapName = 'Preferred Condition';
        [dprime vals] = mv_dprime(mv);

    case 'voxel reliability',
        if ~isfield(mv, 'voxRSorting'),  mv = mv_sortByVoxR(mv);  end
        mapName = 'Voxel Reliability';
        vals = mv.voxRSorting.metric;

    case 'mutual information',
        if ~isfield(mv, 'wta'),   mv = mv_reliability(mv, 'plotFlag', 0);  end
        mv = mv_mutualInformation(mv, [], 'auto', 1);
        mapName = 'Mutual Information';
        vals = mv.mutualInf.Im;

    case 'omnibus reliability',
        if ~isfield(mv, 'MISorting'),  mv=mv_sortByOmniR(mv); end
        mapName = 'Omnibus Reliability';
        vals = mv.MISorting.metric;

    case 'selectivity index',
        % special function call for high-res / lo-res comparisons
        mv = mv_exportSelectivity(mv, saveFlag);
        return

    case 'mean voxel amplitudes',
        % ask user what sort of amplitudes to export
        ampNames = {'Peak-Baseline Amplitude' 'Beta Weights' ...
                    'Dot-product Amplitudes'};
        n = cellfind({'difference' 'betas' 'relamps'}, mv.params.ampType);
        mapName = ampNames{n};
        vals = mean(mv_amps(mv), 2)';

    case 'amplitudes for each condition',
        % special function call: multiple maps
        mv = mv_exportAmplitudeMap(mv);
        return

    case 'contrast',
        % contrast map values: conditions specified in varargin
        if length(varargin) < 2
            % get from dialog
           prompt = {'Active Conditions' 'Control Conditions' 'Map Name'};
           name = 'Input for Peaks function';
           resp = inputdlg(prompt, mfilename, 1, {'1' '0' ''});
           active = str2num(resp{1});
           control = str2num(resp{2});
		   mapName = resp{3};

        else
            active = varargin{1};
            control = varargin{2};
			if length(varargin) > 2
				mapName = varargin{3};
			end

        end
        if ~isfield(mv, 'glm'), mv = mv_applyGlm(mv); end
        vals = glm_contrast(mv.glm, active, control);

        activeName = [mv.trials.condNames{active+1}];
        ctrlName = [mv.trials.condNames{control+1}];
	   if isempty(mapName)
			mapName = sprintf('%sV%s', activeName, ctrlName);
	   end

	otherwise,
		myErrorDlg(sprintf('Parameter %s is not defined.', param));

end

% initialize a view to get parameters about the
% map volume
fn = sprintf('getSelected%s', mv.roi.viewType);
view = eval(fn);
if isempty(view)
    % no selected view of the proper type -- make a hidden one
    mrGlobals; loadSession; saveFlag = 1;
    fn = sprintf('initHidden%s', roi.viewType); view = eval(fn);
end
mapdims = viewGet(view, 'dataSize');
nScans = viewGet(view, 'numScans');
scan = mv.params.scans(1);

% plug in the values to the map volume:
mapvol = zeros(mapdims);
ind = roiIndices(view,mv.coords);
if ~exist('vals','var')
	fprintf('error: %s map does not exist\n',lower(param));
else
	mapvol(ind) = vals;
end

% plug the map volume into the view
fprintf('Assigning %s data to param map for scan %i\n', mapName, scan);
map = cell(1, nScans);
map{scan} = mapvol;
if ~isequal(view.name, 'hidden')
    view = setParameterMap(view, map, mapName);
    refreshScreen(view);
end

% evaluate this in the workspace,  so the view
% itself is updated
assignin('base', 'map', map);
evalin('base', sprintf('%s=setParameterMap(%s, map, ''%s'');', ...
    view.name, view.name, mapName));

% save if selected
if saveFlag==1
    mapName = sprintf('%s_%s', mapName, mv.roi.name);
    mapPath = fullfile(dataDir(view), mapName);
    if exist(mapPath, 'dir')
        load(mapPath, 'map', 'mapName');
        map{scan} = mapvol;
    end
    save(mapPath, 'map', 'mapName');
    fprintf('Saved map %s.\n', mapPath);
end

% this is a check for mrVista 2: if we're running a mrVista2 session GUI,
% update it to show the new map
global GUI
if ~isempty(GUI)
    sessionGUI_selectDataType;
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function dpts = move_pts(pts, By, Bx)
dpts = pts;
for j = 1 : size(pts, 2)
    r = ceil(pts(1, j));
    c = ceil(pts(2, j));
    dpts(1, j) = dpts(1, j) + By(r, c);
    dpts(2, j) = dpts(2, j) + Bx(r, c);
end

repo: pTVreg


text:
% p23 例5 埃特金逐步插值
clear;
format long;
x = 0.3 :0.1 : 0.7;
y = [0.29850, 0.39646, 0.49311, 0.58813, 0.68122];
xx = 0.462; % 插值点
fij(:,1) = y;
for i = 1 : 3
    for j = i+1 : 5
       fij(j,i+1) = fij(i,i)*(xx-x(j))/(x(i)-x(j)) + fij(j,i)*(xx-x(i))/(x(j)-x(i))  % 显示中间结果，体现迭代
    end
end
fij

repo: NumericalAnalysis


text: %UMINUS Mapping overload


text:
classdef MimMarkerCyan < MimGuiPlugin
    % MimMarkerCyan. Gui Plugin for switching transparency of zero
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimMarkerCyan is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    properties
        ButtonText = 'Cyan'
        SelectedText = 'Cyan'
        ToolTip = 'Place cyan markers'
        Category = 'Marker colour'
        Visibility = 'Dataset'
        Mode = 'Markers'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 5
        ButtonHeight = 1
        Icon = 'single_marker.png'
        IconColour = GemMarkerPoint.DefaultColours{4}
        Location = 4
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.NewMarkerColour = 4;
        end

        function enabled = IsEnabled(gui_app)
            enabled = gui_app.IsDatasetLoaded && isequal(gui_app.GetCurrentModeName, MimModes.MarkerMode);
        end

        function is_selected = IsSelected(gui_app)
            is_selected = gui_app.ImagePanel.NewMarkerColour == 4;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function [f, g] = fgplvmSequenceObjectiveGradient(xvec, model, Y)

% FGPLVMSEQUENCEOBJECTIVEGRADIENT Wrapper function for objective
% and gradient of a single sequence in latent space and the corresponding output sequence.
% FORMAT
% DESC provides a wrapper function for the negative log probability
% of a given data sequence under the posterior distribution of the
% Gaussian process induced by the training data..
% ARG vecx : time ordered locations in input space for the sequence
% placed as a vector using the matlab X(:)' notation.
% ARG model : the model structure for which the negative log
% probability of the given data under the posterior is to be computed.
% ARG Y : time ordered locations in data spaces for the sequence.
% RETURN f : the negative of the log probability of the given data
% sequence under the posterior distribution induced by the training data.
% RETURN g : the gradient of the negative of the returned log probability
% with respect to the latent sequence.
%
% SEEALSO : fgplvmCreate, fgplvmSequenceLogLikelihood, fgplvmOptimiseSequence
%
% COPYRIGHT : Neil D. Lawrence, 2006

% FGPLVM

% Check how the optimiser has given the parameters
if size(xvec, 1) > size(xvec, 2)
  % As a column vector ... transpose everything.
  transpose = true;
  X = reshape(xvec', size(Y, 1), model.q);
else
  transpose = false;
  X = reshape(xvec, size(Y, 1), model.q);
end
f = - fgplvmSequenceLogLikelihood(model, X, Y);

if nargout > 1
  g = - fgplvmSequenceLogLikeGradient(model, X, Y);
  g = g(:)';
end
if transpose
  g = g';
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function [ c ] = get_bpsk_cap( snr_db )
% assumes e[x^2] = 1;

x_vec = [-1, 1];
n_0 = 1/2*10^(-snr_db/10);
delta_y = sqrt(n_0) * 0.001;
max_value = min(10000, max(x_vec) + 3 + 3*sqrt(n_0));
y_vec = -max_value: delta_y : max_value;
p_y = zeros(length(y_vec), 1);
for y_index = 1 : length(y_vec)
    y = y_vec(y_index);
    for x_index = 1 : length(x_vec)
        x = x_vec(x_index);
        p_y(y_index) = p_y(y_index) + exp(-(y - x).^2/2/n_0)/sqrt(2*pi*n_0) * 0.5;
    end
end

p_y = p_y / (sum(p_y)*delta_y);

h = 0;
for y_index = 1 : length(y_vec)
    if p_y(y_index) > 0
        h = h + (-log2(p_y(y_index))) * p_y(y_index) * delta_y;
    end
end
c = h - 0.5 * (1 + log(2*pi*n_0))/log(2) ;

author: tavildar

repo: Polar

save_path: github-repos/MATLAB/tavildar-Polar


text:
function [out1,out2,out3] = myginput(arg1)
% almost the same as matlab's ginput, but does not change the mouse pointer

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

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

   fig = gcf;
   figure(gcf);

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

   % Remove figure button functions
   state = uisuspend(fig);
   pointer = get(gcf,'pointer');
   %set(gcf,'pointer','fullcrosshair');
   fig_units = get(fig,'units');
   char = 0;

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

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

         how_many = how_many - 1;

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

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

   uirestore(state);
   set(fig,'units',fig_units);

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

end

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

fig = gcf;
current_char = [];

% Now wait for that buttonpress, and check for error conditions
waserr = 0;
try
  h=findall(fig,'type','uimenu','accel','C');   % Disabling ^C for edit menu so the only ^C is for
  set(h,'accel','');                            % interrupting the function.
  keydown = waitforbuttonpress;
  current_char = double(get(fig,'CurrentCharacter')); % Capturing the character.
  if~isempty(current_char) & (keydown == 1)           % If the character was generated by the
	  if(current_char == 3)                       % current keypress AND is ^C, set 'waserr'to 1
		  waserr = 1;                             % so that it errors out.
	  end
  end

  set(h,'accel','C');                                 % Set back the accelerator for edit menu.
catch
  waserr = 1;
end
drawnow;
if(waserr == 1)
   set(h,'accel','C');                                % Set back the accelerator if it errored out.
   error('Interrupted');
end

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

author: kristinbranson

repo: JAABA


text:
load('test_data.mat')
iou_ = compute_dense_overlap(ofx,ofy,stx,sty,vsx,vsy,dx1,dy1,dx2,dy2,gx1,gy1,gx2,gy2,1,1);
err = abs(iou_ - iou);
if all(err(:)) < 1e-12
    fprintf('Test for compute_dense_overlap [passed]\n');
else
    fprintf('Test for compute_dense_overlap [failed]\n');
end

author: peiyunh

repo: tiny

save_path: github-repos/MATLAB/peiyunh-tiny


text:
function p = legpoly(n, dom, normalize, method)
%LEGPOLY   Legendre polynomials.
%   P = LEGPOLY(N) computes a CHEBFUN of the Legendre polynomial of degree N on
%   the interval [-1,1]. N can be a vector of integers, in which case the output
%   is an array-valued CHEBFUN.
%
%   P = LEGPOLY(N, D) computes the Legendre polynomials as above, but on the
%   interval given by the domain D, which must be bounded.
%
%   P = LEGPOLY(N, D, 'norm') or P = LEGPOLY(N, 'norm') normalises so that
%   integral(P(:,j).*P(:,k)) = delta_{j,k}.
%
%   For N <= 1000 LEGPOLY uses a weighted QR factorisation of a 2*(N+1) x
%   2*(N+1) Chebyshev Vandermonde matrix. For scalar N > 1000 (or a short
%   vector) it uses the LEG2CHEB method and for a vector of N with any entry >
%   1000 it uses the standard recurrence relation. This default can be
%   overwritten by passing a fourth input LEGPOLY(N, D, NORM, METHOD), where
%   METHOD is 1, 2, or 3 respectively.
%
% See also CHEBPOLY, LEGPTS, and LEG2CHEB.

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

% Parse input:
methodIsSet = false;
if ( isempty(n) )
    p = chebfun;
    return
end
if ( nargin < 2 || isempty(dom) )
    dom = [-1, 1];
end
if ( nargin < 3 || isempty(normalize) )
    normalize = 0;
end
if ( ischar(dom) )
    if ( nargin == 3 )
        method = normalize;
        if ( ~isempty(method) )
            methodIsSet = true;
        end
    end
    normalize = dom;
    dom = [-1,1];
end
if ( strncmp(normalize, 'norm', 4) )
    normalize = 1;
elseif ( ischar(normalize) )
    normalize = 0;
end
if ( (nargin == 4) && ~isempty(method) )
    methodIsSet = true;
end

% Unbounded domains aren't supported/defined.
if ( any(isinf(dom)) )
    error('CHEBFUN:legpoly:infdomain', ...
        'Legendre polynomials are not defined over an unbounded domain.');
end

% Force a CHEBTECH basis.
defaultPref = chebfunpref();
pref = defaultPref;
tech = feval(pref.tech);
if ( ~isa(tech, 'chebtech') )
    pref.tech = @chebtech2;
end

% Useful values:
nMax = max(n);
nMax1 = nMax + 1;
domIn = dom;
dom = dom([1 end]);

% Determine which method:
if ( ~methodIsSet && nMax > 1000 )
    % Use LEG2CHEB():
    method = 3;
elseif ( ~methodIsSet )
    % Use QR orthogonalization of Chebyshev polynomials for moderate nmax.
    method = 2;
end

% If the user wants most of the Legendre polynomials then it is faster to use
% the recurrence: TODO: "most" is most than 20% of the [0:nMax]. Should this
% percentage vary with nMax?
if ( ~methodIsSet && (method == 3) && (nMax < 5000) && (numel(n) > nMax/5) )
    method = 1; % Do not go above 5000 with the recurrence.
end

switch method
    case 1 % Recurrence

        [aa, bb, cc] = unique(n);      %#ok<ASGLU>
        P = zeros(nMax1, length(n));   % Initialise storage
        x = chebpts(nMax1, 2);         % Chebyshev points
        L0 = ones(nMax1, 1); L1 = x;   % P_0 and P_1
        ind = 1;                       % Initialise counter
        for k = 2:nMax+2,              % The recurrence relation (k = degree)
            if ( aa(ind) == k-2 )
                if ( normalize )
                    invnrm = sqrt((2*k-3)/diff(dom));
                    P(:,ind) = L0*invnrm;
                else
                    P(:,ind) = L0;
                end
                ind = ind + 1;
            end
            tmp = L1;
            L1 = (2-1/k)*x.*L1 - (1-1/k)*L0;
            L0 = tmp;
        end
        C = chebtech2.vals2coeffs(P(:,cc));       % Convert to coefficients
    case 2 % QR

        pts = 2*nMax1;              % Expand on Chebyshev grid of twice the size
        [ignored, w] = chebpts(pts, 2);   % Grab the Clenshaw-Curtis weights
        theta = pi*(pts-1:-1:0)'/(pts-1);
        A = cos(theta*(0:nMax));                  % Vandemonde-type matrix
        D = spdiags(sqrt(w(:)), 0, pts, pts);     % C-C quad weights
        Dinv = spdiags(1./sqrt(w(:)), 0, pts, pts);
        [Q, ignored] = qr(D*A, 0);                % Weighted QR
        P = Dinv*Q;
        if ( normalize )
            PP = P(:,n+1) * diag(sqrt(2/diff(dom)) * sign(P(end,n+1)));
        else
            PP = P(:,n+1) * diag(1./P(end,n+1));
        end
        C = chebtech2.vals2coeffs(PP);            % Convert to coefficients
        C(nMax1+1:end,:) = [];                    % Trim coefficients > nMax

    case 3 % LEG2CHEB

        c_leg = zeros(nMax+1, numel(n));
        c_leg(n+1,:) = eye(numel(n));             % Legendre coefficients
        if ( normalize )
            C = leg2cheb(c_leg, 'norm');  % Chebyshev coefficients
        else
            C = leg2cheb(c_leg);          % Chebyshev coefficients
        end

end

% Construct CHEBFUN from coeffs:
p = chebfun(C, dom, pref, 'coeffs');

if ( numel(domIn) > 2 )
    p = restrict(p, domIn);
end

% Adjust orientation:
if ( size(n, 1) > 1 )
   p = p.';
end

end

repo: chebfun


text:
function s_factor_cell = generate_s_factor(omega, grid3d, m, R)
% For example, s_factor_cell{Axis.x, GT.dual} is the s-factor multiplied to Ex
% (or eps_ww).

chkarg(istypesizeof(omega, 'complex'), '"omega" should be complex.');
chkarg(istypesizeof(grid3d, 'Grid3d'), '"grid3d" should be instance of Grid3d.');

if nargin < 3  % no m
	m = 4;
end
chkarg(istypeof(m, 'real') && all(m(:) >= 1), 'element of "m" should be real and at least 1.0.');
chkarg(isexpandable2mat(m, Axis.count, Sign.count), ...
	'"m" should be scalar, length-%d vector, or %d-by-%d matrix.', Axis.count, Axis.count, Sign.count);
m = expand2mat(m, Axis.count, Sign.count);

if nargin < 4  % no R
	R = exp(-16);  % R = exp(-16) ~= 1e-7
end
chkarg(istypeof(R, 'real') && all(R(:) > 0) && all(R(:) <= 1), ...
	'element of "R" should be real and between 0.0 and 1.0.');
chkarg(isexpandable2mat(R, Axis.count, Sign.count), ...
	'"R" should be scalar, length-%d vector, or %d-by-%d matrix.', Axis.count, Axis.count, Sign.count);
R = expand2mat(R, Axis.count, Sign.count);
lnR = log(R);  % log() is the natural logarithm

s_factor_cell = cell(Axis.count, GT.count);
for w = Axis.elems
	Nw = grid3d.N(w);

	lpml = grid3d.lpml(w,:);  % locations of PML interfaces
	Lpml = grid3d.Lpml(w,:);  % thicknesses of PML

	for g = GT.elems
		l = grid3d.l{w, g};  % length(l) == Nw, rather than Nw+1.
		ind_pml = {l < lpml(Sign.n), l > lpml(Sign.p)};  % indices of locatiotns indise PML

		s_factor = ones(1, Nw);
		for s = Sign.elems
			s_factor(ind_pml{s}) = calc_s_factor(omega, abs(lpml(s) - l(ind_pml{s})), Lpml(s), m(w,s), lnR(w,s));
		end

		s_factor_cell{w,g} = s_factor;
	end
end

function s_factor = calc_s_factor(omega, depth, Lpml, m, lnR)
% assert(Lpml > 0);

sigma_max = -(m+1) * lnR/(2*Lpml);  % -(m+1) ln(R)/(2 eta Lpml), where eta = 1 in the unit of eta_0
sigma = sigma_max * (depth/Lpml).^m;

kappa_max = 1;
kappa = 1 + (kappa_max-1) * (depth/Lpml).^m;

ma = m;
amax = 0;
a = amax * (1 - depth/Lpml).^ma;

s_factor = kappa + sigma ./ (a + 1i*omega);  % s = kappa + sigma/(a + i omega eps), where eps = 1 in the unit of eps_0

author: wsshin

repo: maxwellfdfd

save_path: github-repos/MATLAB/wsshin-maxwellfdfd


text:
function y=getmondim(mnum)
% y = getmondim(monitornumber)  GET MONitor DIMensions.
%   "y" is the dimensions of the specified monitor represented as
%   [xstart,ystart,width,height] where values are in pixels.
%   (xstart,ystart) are the absolute coordinates of the lower left corner.
%   "monitornumber" is a number associated with each monitor (Default=1).
%       Monitor numbering is defined as the row number of the corresponding
%       monitor in the root property 'MonitorPositions'.
%   Examples:
%       y=getmondim; y=getmondim(2);

%   Copyright 2006 Mirtech, Inc.
%   created 09/26/2006  by Mirko Hrovat on Matlab Ver. 7.2
%   Mirtech, Inc.       email: mhrovat@email.com
%  Notes:
%       The primary monitor is always numbered 1 and always starts at (1,1).
%   'MonitorPositions' returns coordinates corresponding to the upper left
%   corner and the bottom right corner (Windows convention). Matlab 'Help'
%   on this property is incorrect. GETMONDIM converts these coordinates to
%   the Matlab convention where the lower left corner starts at (1,1).
%       There is a bug with the root property 'ScreenSize'. If the primary
%   monitor's resolution is adjusted after Matlab has started, then the
%   size parameters of 'ScreenSize' are not adjusted while the origin is
%   adjusted. GETMONDIM.M fixes this be using the 'ScreenSize' origin to
%   correct for the discrepancy. Note that on restarting Matlab the
%   'ScreenSize' property changes to the correct values!

if nargin==0,       mnum=[];    end
if isempty(mnum),   mnum=1;     end
savedunits=get(0,'Units');
set(0,'Units','pixels')         % make sure unit values are in pixels
mpos=get(0,'MonitorPositions'); % get list of monitor positions
% MonitorPositions are the coordinates for upper left & lower right corners.
primaryhght=mpos(1,4)-mpos(1,2)+1;
if any(mnum==1:size(mpos,1)),
    % need to convert MonitorPositions to Matlab convention.
    height=mpos(mnum,4)-mpos(mnum,2)+1;
    width=mpos(mnum,3)-mpos(mnum,1)+1;
    offset=1;       %#ok
% NOTE!
% Bugfix for the root property 'ScreenSize' in case the primary monitor's
% resolution has changed. If bug is ever fixed then delete these lines.
    screenxy=get(0,'ScreenSize');
    offset=2-screenxy(2);
% end bugfix

    % get lower left corner coordinates
    llpt=[mpos(mnum,1),primaryhght-mpos(mnum,4)+offset];
    y=[llpt,width,height];
else
    warning('  Monitor does not exist!');   %#ok
    y=[];
end
set(0,'Units',savedunits)

author: Sable

repo: mcbench-benchmarks


text:
%TIME_FAC_VS_FB
%
%   The purpose of this test is to compare the fastest of the factorization
%   routines against the fastest of the filter bank routines for different
%   problems.
%
%   Some results: The FB computation seems to be completely memory bound:
%   increasing the length of the window does not impact the running time.
%
%   For the common 2xoversampling case with a long signal, the
%   factorization routine falls thorugh and is more than 10 times slower,
%   for the other case it is 3 times slower.

Lr=[480000*sf^2,480000*sf^2,262144*sf^2,262144*sf^2,900*sf^2,600];
ar=[     600*sf,     600*sf,        512,        512,       2, 20];
Mr=[     800*sf,     800*sf,       1024,       1024,  600*sf, 30];
gr=[     800*sf,480000*sf^2,       1024,     40*512,  600*sf,400];
Wr=[          2,          2,          2,          2,       1,600];

for ii=1:length(Lr)

  L=Lr(ii);

  M=Mr(ii);
  a=ar(ii);

  gl=gr(ii);
  W=Wr(ii);

  N=L/a;
  c=gcd(a,M);
  p=a/c;
  q=M/c;
  d=N/q;

  disp(sprintf('L:%6i W:%2i a:%4i M:%4i gl:%4i c:%4i d:%4i',L,W,a,M,gl,c,d));

  f=rand(L,W);
  gfir=rand(gl,1);
  c2=mex_dgt_fb_2(f,gl,a,M,1);

  f=rand(L,W);
  gf=rand(p*q,c*d);
  c1=mex_dgt_fac_7(f,gf,a,M,1);

end;

repo: ltfat


text:
%-----------------------------------------------------------------------
% Job configuration created by cfg_util (rev $Rev: 4899 $)
%-----------------------------------------------------------------------
matlabbatch{1}.menu_cfg.menu_struct.conf_choice.type = 'cfg_choice';
matlabbatch{1}.menu_cfg.menu_struct.conf_choice.name = 'Dir Operations';
matlabbatch{1}.menu_cfg.menu_struct.conf_choice.tag = 'dir_ops';
matlabbatch{1}.menu_cfg.menu_struct.conf_choice.values = {
                                                          '<UNDEFINED>'
                                                          '<UNDEFINED>'
                                                          '<UNDEFINED>'
                                                          }';
matlabbatch{1}.menu_cfg.menu_struct.conf_choice.check = [];
matlabbatch{1}.menu_cfg.menu_struct.conf_choice.help = {''};

author: spm


text:
clc; clear variables; close all;
addpath(genpath(fileparts(mfilename('fullpath'))));

clear variables;

%% Main
filenames={
    'paramsTopOptProblem'
    };

for icases=1:size(filenames,1)
    clearvars -except filenames icases;
    close all;

    settings = SettingsTopOptProblem(filenames{icases});

    topOptProblem = TopOpt_Problem(settings);
    topOptProblem.computeVariables;
    topOptProblem.postProcess;
end
close all

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function [x, resvec, lsvec] = dense_overdetermined_lsqr(A, b, R, tol, maxit)
% x = dense_overdetermined_lsqr(A, b, R, tol, maxit)
% [x, resvec, lsvec] = dense_overdetermined_lsqr(A, b, R, tol, maxit)
%
% LSQR on a dense overdetermined matrix with a dense preconditioner R.
% Solves x = arg min (A*x - b)
% Inputs:
%   A, b
%   R - Upper triangular preconditioner. kappa(A * inv(R)) governs
%       the convergence.
%   tol - tolerance. Stop when
%      norm(inv(R') * A' * r) / (norm(A * inv(R), 'fro') * norm(r)) < tol
%   maxit - Maximum number of iterations
%
% Outputs:
%   x
%   optional: resvec, lsvec - values of norm(A' * r) and norm(r).
%
% 6-December 2009, Version 1.3
% Copyright (C) 2009, Haim Avron and Sivan Toledo.


text:
function [R, G, B] = Lab2RGB(L, a, b)
%LAB2RGB Convert an image from CIELAB to RGB
%
% function [R, G, B] = Lab2RGB(L, a, b)
% function [R, G, B] = Lab2RGB(I)
% function I = Lab2RGB(...)
%
% Lab2RGB takes L, a, and b double matrices, or an M x N x 3 double
% image, and returns an image in the RGB color space.  Values for L are in
% the range [0,100] while a* and b* are roughly in the range [-110,110].
% If 3 outputs are specified, the values will be returned as doubles in the
% range [0,1], otherwise the values will be uint8s in the range [0,255].
%
% This transform is based on ITU-R Recommendation BT.709 using the D65
% white point reference. The error in transforming RGB -> Lab -> RGB is
% approximately 10^-5.
%
% See also RGB2LAB.

% By Mark Ruzon from C code by Yossi Rubner, 23 September 1997.
% Updated for MATLAB 5 28 January 1998.
% Fixed a bug in conversion back to uint8 9 September 1999.
% Updated for MATLAB 7 30 March 2009.

if nargin == 1
  b = L(:,:,3);
  a = L(:,:,2);
  L = L(:,:,1);
end

if max(max(L)) < 1.1 || max(max(a)) < 1.1 || max(max(b)) < 1.1
  L = double(L) * 100;
  a = double(a) * 220;
  a = a - 110;
  b = double(b) * 220;
  b = b - 110;
end

% Thresholds
T1 = 0.008856;
T2 = 0.206893;

[M, N] = size(L);
s = M * N;
L = reshape(L, 1, s);
a = reshape(a, 1, s);
b = reshape(b, 1, s);

% Compute Y
fY = ((L + 16) / 116) .^ 3;
YT = fY > T1;
fY = (~YT) .* (L / 903.3) + YT .* fY;
Y = fY;

% Alter fY slightly for further calculations
fY = YT .* (fY .^ (1/3)) + (~YT) .* (7.787 .* fY + 16/116);

% Compute X
fX = a / 500 + fY;
XT = fX > T2;
X = (XT .* (fX .^ 3) + (~XT) .* ((fX - 16/116) / 7.787));

% Compute Z
fZ = fY - b / 200;
ZT = fZ > T2;
Z = (ZT .* (fZ .^ 3) + (~ZT) .* ((fZ - 16/116) / 7.787));

% Normalize for D65 white point
X = X * 0.950456;
Z = Z * 1.088754;

% XYZ to RGB
MAT = [ 3.240479 -1.537150 -0.498535;
       -0.969256  1.875992  0.041556;
        0.055648 -0.204043  1.057311];

RGB = max(min(MAT * [X; Y; Z], 1), 0);

R = reshape(RGB(1,:), M, N);
G = reshape(RGB(2,:), M, N);
B = reshape(RGB(3,:), M, N);

if nargout < 2
  R = uint8(round(cat(3,R,G,B) * 255));
end

author: happynear

repo: DeepVisualization


text:
function calpak_test585 ( )

%*****************************************************************************80
%
%% CALPAK_TEST585 tests YEAR_LENGTH_GREEK.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2013
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'CALPAK_TEST585\n' );
  fprintf ( 1, '  For the Greek calendar:\n' );
  fprintf ( 1, '  YEAR_LENGTH_GREEK determines the length of a year.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Year  Length\n' );
  fprintf ( 1, '\n' );

  for y = 1 : 10
    sy = y_to_s_greek ( y );
    fprintf ( 1, '  %10s  %d\n', sy, year_length_greek ( y ) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef prtBrvDpMixture < prtBrvMixture
    % ds = prtDataGenOldFaithful; x = ds.getX;
    %
    % mm = prtBrvDpMixture('components',repmat(prtBrvMvn,25,1));
    % mm.vbConvergenceThreshold = 1e-6;
    % mm.vbVerboseText = true;
    % mm.vbVerbosePlot = true;
    % mm.vbVerboseMovie = true;
    % [mmLearned, training] = mm.vb(x);

    methods

        function obj = prtBrvDpMixture(varargin)
            obj.mixing = prtBrvDiscreteStickBreaking;

            if nargin < 1
                return
            end
            obj = constructorInputParse(obj,varargin{:});
        end
    end
end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function func_GrayCSF
% Show Contrast Sensitive Function (CSF) via gray Sinwave stimulus
% This function is part of the toolbox "Basic introduction to HVS"
%
% Command line
% ----------------------
% func_GrayCSF
% input:  None
%
% output: CSF figure
%
% More information can be found in
% S. E. Palmer, "Vision Science: From Photons to Phenomenology," MIT Press,
% Cambridge, MA, 1999.
%
%
% Jing Tian Apr.24 2004
% Contact me : scuteejtian@hotmail.com
% Homepage : http://ikanchi.yeah.net
% This program is written in Apr.2003 during my postgraduate in
% NTU, Singapore.
% ----------------------

% parameter of Sinwave stimulus
% freq : frequency
% C    : Contrast
% low : step : high
freq = logspace(0.1, 0.9, 100)';
C = logspace(-2, 0, 100);
L = 100;

x = linspace(-1.5 * pi, 0.5 * pi, 100);
y = linspace(1, 100, 100);
[xx,yy] = meshgrid(x, y);
[newfreq , newC] = meshgrid(freq, C);

z = L .* (newC .* sin(2 .* pi .* newfreq .* xx) + 1);

imshow(z, []);
shading interp;
axis('off')


text:
function drawCube( origin, size,color,alpha,a)
% From
% http://www.mathworks.com/matlabcentral/newsreader/view_thread/235581
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri
%--------------------------------------------------------------------------
x=([0 1 1 0 0 0;1 1 0 0 1 1;1 1 0 0 1 1;0 1 1 0 0 0]-0.5)*size(1)+origin(1);
y=([0 0 1 1 0 0;0 1 1 0 0 0;0 1 1 0 1 1;0 0 1 1 1 1]-0.5)*size(2)+origin(2);
z=([0 0 0 0 0 1;0 0 0 0 0 1;1 1 1 1 0 1;1 1 1 1 0 1]-0.5)*size(3)+origin(3);
arad=a*pi/180;
R=[cos(arad) -sin(arad) 0; sin(arad) cos(arad) 0; 0 0 1];

for i=1:6
    h=patch(x(:,i),y(:,i),z(:,i),color);
    set(h,'facealpha',alpha)
    rotate(h,[0 0 1],a,[0 0 0]);
end

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
function net = LLRAddMasking(net, upsample_fac, bilinear_up, upsample_2x_per_layer)

if upsample_2x_per_layer
    post_name = '_bil_x2';
else
    post_name = '';
end
up_name = [num2str(upsample_fac) 'x'];
pre_up_name = [num2str(2 * upsample_fac) 'x'];

pool_size = (upsample_fac * 2) / bilinear_up + 1;
assert(pool_size == 9);
pad_size = floor(pool_size)/2;
% Dilation of class probabilities
net.addLayer(['prob_' pre_up_name], dagnn.SoftMax(), ['prediction_' pre_up_name post_name], ['prob_' pre_up_name], {});
net.addLayer(['prob_dilate_' pre_up_name], dagnn.Pooling('stride', [1 1], 'poolSize', [pool_size pool_size], ...
    'pad', [pad_size, pad_size, pad_size, pad_size]), ['prob_' pre_up_name], ['prob_' pre_up_name '_dilate'], {});

% Dilation of negative of class probabilities
net.addLayer(['neg_prob_' pre_up_name], Neg(), ['prob_' pre_up_name], ['neg_prob_' pre_up_name], {});
['neg_prob_' pre_up_name]
net.addLayer(['neg_prob_dilate_' pre_up_name], dagnn.Pooling('stride', [1 1], 'poolSize', [pool_size pool_size], ...
    'pad', [pad_size, pad_size, pad_size, pad_size]), ['neg_prob_' pre_up_name], ['neg_prob_' pre_up_name '_dilate'], {});

% Sum of two dilation
net.addLayer(['bound_mask' pre_up_name], dagnn.Sum(), {['prob_' pre_up_name '_dilate'], ['neg_prob_' pre_up_name '_dilate']}, ['bound_mask' pre_up_name]) ;
net.addLayer(['dot_prod_' up_name], DotProduct(), {['bound_mask' pre_up_name], ['prediction_' up_name '_add']}, ['pred_' up_name '_aft_DP']) ;
net.setLayerInputs(['sum' up_name], {['prediction_' pre_up_name post_name], ['pred_' up_name '_aft_DP']});

author: aimerykong


text:
%GETNAME High level routine for finding names of datasets and classifiers
%
%   NAME = getname(A,N)
%
% INPUT
%   A    Dataset, mapping or cell array of datasets or mappings
%   N    Number of characters in NAME (default: all)
%
% OUTPUT
%   NAME  Dataset name or cell array of names
%
% DESCRIPTION
% If N given, the return string has exactly N characters. This is done by
% truncation or by padding with blanks. This is useful for display purposes.
% SEE ALSO (PRTools Guide)
% DATASETS, MAPPINGS

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com

function name = getname(c,n)

  if nargin < 2, n = []; end
  if iscell(c)
    name = cell(1,numel(c));
    for j=1:numel(c)
      name{j} = getname(c{j},n);
    end
  elseif isa(c,'prdataset') || ismapping(c)
    % this will go to dataset or mapping getname
    name = getname(c,n);
  else
    error('Illegal data type')
  end

return

author: marianux

repo: ecg-kit


text:
% LOTKA-VOLTERRA system

clear all, close all, clc
figpath = '../FIGURES/'; mkdir(figpath)
datapath = '../DATA/';
addpath('../utils');

%% Load Models
ModelCollection = {'DMDc','SINDYc', 'NARX'};
Nmodels = length(ModelCollection);
InputSignalType = 'sphs'; % prbs; chirp; noise; sine2; sphs; mixed
TrainAlg = 'trainbr'; %trainlm,trainbr
dep_trainlength = 1;
dep_noise = 0;

% Ntrain_vec = [5:15,20:5:95,100:100:1000];%,1500:500:3000];
% eta_vec = 0.0;
% Nr = 1;

Ntrain_vec = [5:15,20:5:95,100:100:1000,1250,1500,2000,3000];%,1500:500:3000];
eta_vec = 0.05;
Nr = 1;

N_LENGTHS = length(Ntrain_vec);

%% Load data
ResultsALL(N_LENGTHS,1:Nmodels) = struct('x', [], 'u', [], 't', [], 'xref', [], 'J', [], 'eval_time', []);

ModelName = 'DMDc';
datapath1 = [datapath,'EX_LOTKA_Dependencies/',ModelName,'/'];
load(fullfile(datapath1,['EX_LOTKA_MPC_',ModelName,'_',InputSignalType,'_TrainLength','.mat']))
ResultsALL(:,1) = Results;

ModelName = 'SINDYc';
datapath1 = [datapath,'EX_LOTKA_Dependencies/',ModelName,'/'];
load(fullfile(datapath1,['EX_LOTKA_MPC_',ModelName,'_',InputSignalType,'_TrainLength','.mat']))
ResultsALL(:,2) = Results;

ModelName = 'NARX';
datapath1 = [datapath,'EX_LOTKA_Dependencies/',ModelName,'/',TrainAlg,'/'];
load(fullfile(datapath1,['EX_LOTKA_MPC_',ModelName,'_',InputSignalType,'_TrainLength','.mat']))
ResultsALL(:,3) = Results;

%% Correct cost
Q = [1,1]; R = 0.5; Ru = 0.5;
for iM = 1:Nmodels
    for iL = 1:N_LENGTHS
        ResultsALL(iL,iM).J = evalObjectiveFCN(ResultsALL(iL,iM).u,ResultsALL(iL,iM).x,ResultsALL(iL,iM).xref,diag(Q),R,Ru);
    end
end

%% Show trajectories
cb_hrzntl = 1;

cmap = jet(N_LENGTHS);
for jModel = 1:Nmodels
    figure, hold on, box on
    for iM = 1:N_LENGTHS
        plot(ResultsALL(iM,jModel).x(1,:),ResultsALL(iM,jModel).x(2,:),'-','Color',cmap(iM,:),'LineWidth',1)
    end
    plot(ResultsALL(iM,jModel).xref(1,1),ResultsALL(iM,jModel).xref(2,1),'ok','MarkerFaceColor','k')
    %     axis equal
    ylim([-20 100]), xlim([0 200])
    colormap(cmap),
    ch = colorbar;
    ch.Limits = [0 1];
    ch.Ticks = [0,0.25,0.5,0.75,1];
    ch.TickLabels = num2str( Ntrain_vec( [1, (N_LENGTHS-1)/4, (N_LENGTHS-1)/2+1, 3*(N_LENGTHS-1)/4+1, 4*(N_LENGTHS-1)/4+1] )' );
    if cb_hrzntl == 0
        ch.Position = [ch.Position(1),ch.Position(2)+0.14,ch.Position(3)-0.01,ch.Position(4)-0.1];
        set(gca,'Position',[0.08 0.2753 0.775 0.6497])
    elseif cb_hrzntl == 1
        ch.Location = 'southoutside';
        ch.Position = [0.2,0.1,0.75,ch.Position(4)];
        set(gca,'Position',[0.2 0.4 0.75 0.55])
    end
    xlabel('x1');
    ylabel('x2')
    set(gca,'LineWidth',1, 'FontSize',14)
    set(gcf,'Position',[100 100 300 200])
    set(gcf,'PaperPositionMode','auto')
    print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_CTRLPERF_',InputSignalType,'_Eta',sprintf('%03g',100*eta_vec),'_TrainingLength_State_',ModelCollection{jModel},'.eps']);
end

%% Show cost trajectory
cb_hrzntl = 1;

cmap = jet(N_LENGTHS);
for jModel = 1:Nmodels
    figure, hold on, box on
    for iM = 1:N_LENGTHS
        plot(ResultsALL(end,jModel).t,cumsum(ResultsALL(iM,jModel).J),'-','Color',cmap(iM,:),'LineWidth',1)
    end
    plot(ResultsALL(end,2).t,cumsum(ResultsALL(end,2).J),'--','Color','k','LineWidth',2)
    ylim([-10 10^6]), xlim([200 max(ResultsALL(end,jModel).t)])
    colormap(cmap),
    ch = colorbar;
    ch.Limits = [0 1];
    ch.Ticks = [0,0.25,0.5,0.75,1];
    ch.TickLabels = num2str( Ntrain_vec( [1, (N_LENGTHS-1)/4, (N_LENGTHS-1)/2+1, 3*(N_LENGTHS-1)/4+1, 4*(N_LENGTHS-1)/4+1] )' );
    if cb_hrzntl == 0
        ch.Position = [ch.Position(1),ch.Position(2)+0.14,ch.Position(3)-0.01,ch.Position(4)-0.1];
        set(gca,'Position',[0.08 0.2753 0.775 0.6497])
    elseif cb_hrzntl == 1
        ch.Location = 'southoutside';
        ch.Position = [0.2,0.1,0.75,ch.Position(4)];
        set(gca,'Position',[0.2 0.4 0.75 0.55])
    end
    ylim([10^3 10^6])
    set(gca,'yscale','log')
    xlabel('Time');
    ylabel('Cost')
    set(gca,'LineWidth',1, 'FontSize',14)
    set(gcf,'Position',[100 100 300 200])
    set(gcf,'PaperPositionMode','auto')
    print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_CTRLPERF_',InputSignalType,'_Eta',sprintf('%03g',100*eta_vec),'_TrainingLength_Cost_',ModelCollection{jModel},'.eps']);
end

%% Show execution time
clear ph
ccolors = [1 0 0; 0 1 0; 0.7,0.7,1];
lstyles = {'-', '--', '-.'};
symbols = {'o', 'd', 's'};

Tx = zeros(length(N_LENGTHS),Nmodels);
for iM = 1:Nmodels
    for iL = 1:N_LENGTHS
        Tx(iL,iM) = ResultsALL(iL,iM).eval_time;
    end
end
figure, hold on, box on
for iM = 1:Nmodels
    plot(Ntrain_vec,Tx(:,iM)./60,symbols{iM},'Color',ccolors(iM,:),'MarkerFaceColor',ccolors(iM,:))
end
xlim([4 1001])
ylim([0.5*10^-5 10^2])
set(gca,'ytick',[10.^([-5:2:2])])
set(gca,'XScale','log','YScale','log')
xlabel('Length of Training Data');
ylabel('Time [m]')
legend(ModelCollection,'Location','SouthEast')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_CTRLPERF_',InputSignalType,'_Eta',sprintf('%03g',100*eta_vec),'_TrainingLength_ExecTime','.eps']);
legend('off')
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_CTRLPERF_',InputSignalType,'_Eta',sprintf('%03g',100*eta_vec),'_TrainingLength_ExecTime','_noleg.eps']);

%% Add penalization
Q = 1*eye(2); %20
for iM = 1:Nmodels
    for iL = 1:N_LENGTHS
       if ResultsALL(iL,iM).t(end) == 0
           % Flag not converged / ended early
           ResultsALL(iL,iM).J(end) = 10^6;
       else
           % Flag if far from end state at end of time
           %ResultsALL(iL,iM).J(end) = ResultsALL(iL,iM).J(end) + abs(ResultsALL(iL,iM).xref(:,end) - ResultsALL(iL,iM).x(:,end))'*Q*abs(ResultsALL(iL,iM).xref(:,end) - ResultsALL(iL,iM).x(:,end));
       end
    end
end

%% Show Performance over training lengths
clear ph
ccolors = [1 0 0; 0 1 0; 0.7,0.7,1];
lstyles = {'-', '--', '-.'};
symbols = {'o', 'd', 's'};

Nstart = 20;
xx = Ntrain_vec(Nstart):2:1000;
yModel = zeros(length(xx),Nmodels);
JendAll = zeros(length(N_LENGTHS),Nmodels);
Jend = cell(Nmodels,1);
xJend = cell(Nmodels,1);

for iM = 1:Nmodels
%     Jend{iM} = zeros(length(N_LENGTHS),Nmodels);
    count = 0;
    for iL = 1:N_LENGTHS,
        tmp = cumsum(ResultsALL(iL,iM).J);
        if tmp(end)<10^10
            count = count + 1;
            xJend{iM}(count) = Ntrain_vec(iL);
            Jend{iM}(count) = tmp(end);
        end
        JendAll(iL,iM) = tmp(end);
    end
    yModel(:,iM) = spline(xJend{iM}(Nstart:end),Jend{iM}(Nstart:end),xx);
end

clear ph
figure, hold on, box on

if eta_vec == 0.05
    fillh1 = fill([4.5*10^0 2.85*10^3 2.85*10^3 4.5*10^0], [5.5*10^3 5.5*10^3 1.5*10^4 1.5*10^4],0.9*ones(1,3)); % /w noise, not performing well above
elseif  eta_vec == 0.0
    fillh1 = fill([4.5*10^0 2.85*10^3 2.85*10^3 4.5*10^0], [5.5*10^3 5.5*10^3 2.5*10^4 2.5*10^4],0.9*ones(1,3));
end
fillh1.EdgeColor = 0.9*ones(1,3); fillh1.FaceAlpha = 0.5;

[~,idx_best] = min(JendAll(:,2));
plot(Ntrain_vec,JendAll(idx_best,2)*ones(size(Ntrain_vec)),'-','Color',ccolors(2,:),'LineWidth',1)

for iM = 1:Nmodels
%     if iM <= Nmodels
%         plot(xx,yModel(:,iM),'-','Color',ccolors(iM,:), 'LineWidth',1, 'LineStyle',lstyles{iM}), hold on,
%     end
    ph(iM) = plot(Ntrain_vec,JendAll(:,iM),symbols{iM},'Color',ccolors(iM,:),'MarkerFaceColor',ccolors(iM,:),'MarkerSize',4);
end
if eta_vec == 0
    iM = 1; plot([6.5,6.5],[5*10^3 10^7],'Color',ccolors(iM,:),'LineStyle',lstyles{iM});
    iM = 2; plot([13.5,13.5],[5*10^3 10^7],'Color',ccolors(iM,:),'LineStyle',lstyles{iM});
    iM = 3; plot([100,100],[5*10^3 10^7],'Color',ccolors(iM,:),'LineStyle',lstyles{iM});
elseif eta_vec == 0.05
    iM = 1; plot([7.5,7.5],[5*10^3 10^7],'Color',ccolors(iM,:),'LineStyle',lstyles{iM});
    iM = 2; plot([70,70],[5*10^3 10^7],'Color',ccolors(iM,:),'LineStyle',lstyles{iM});
    iM = 3; plot([130,130],[5*10^3 10^7],'Color',ccolors(iM,:),'LineStyle',lstyles{iM});
end
xlim([4 3001])
% ylim([-100 5000])
% ylim([5*10^2 1.5*10^4])
% ylim([1*10^2 10^5])
ylim([5*10^3 10^7])
% set(gca,'YScale','log')
set(gca,'XScale','log','YScale','log', 'ytick', [10^3,10^4,10^5,10^6,10^7], 'xtick', [10^0,10^1,10^2,10^3]) %[10^3,5*10^3,10^4]
xlabel('Length of Training Data');
ylabel('Cost')
lh = legend(ph,ModelCollection,'Location','NorthEast');
lh.Position = [0.71 0.73 lh.Position(3)-0.05 lh.Position(4)];
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_CTRLPERF_',InputSignalType,'_Eta',sprintf('%03g',100*eta_vec),'_TrainingLength_Cost','.eps']);

lh.Location = 'NorthOutside';
lh.Orientation = 'horizontal';
lh.Box = 'off';
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_CTRLPERF_',InputSignalType,'_Eta',sprintf('%03g',100*eta_vec),'_TrainingLength_Cost','_legout.eps']);
legend('off')
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_CTRLPERF_',InputSignalType,'_Eta',sprintf('%03g',100*eta_vec),'_TrainingLength_Cost','_noleg.eps']);

%% Show best case for all models
% eta = 0.05: 8    38    21
% -->  Ntrain = 12        1250          65
% --> Jend = 9.0603e+03,  8.6549e+03, 2.0592e+04
% Jend for best eta==0: 8(12), 38(1250), 21(65)
% SINDY converged Ntrain = 200
% eta = 0.0: 12    25    37
% -->  Ntrain = 20          85        1000
% --> Jend = 9.0601e+03,  8.6595e+03, 8.9236e+03
% SINDY converged Ntrain = 200
BestModelIDX = zeros(1,Nmodels);
Jend = zeros(N_LENGTHS,Nmodels);
for iM = 1:Nmodels
   for iL = 1:N_LENGTHS
        tmp = cumsum(ResultsALL(iL,iM).J);
        Jend(iL,iM) = tmp(end);
   end
   [~,BestModelIDX(iM)] = min(Jend(:,iM));
end

% [val,idx1] = min(abs(Jend(:,1)-9.0601e+03))
% Ntrain_vec(idx1)
% [val,idx2] = min(abs(Jend(:,2)-8.6595e+03))
% Ntrain_vec(idx2)
% [val,idx3] = min(abs(Jend(:,3)-8.9236e+03))
% Ntrain_vec(idx3)

ccolors = [1 0 0; 0 1 0; 0.7,0.7,1];
lstyles = {'-', '--', '-.'};

figure,
hold on, box on
plot(ResultsALL(BestModelIDX(1),1).t,ResultsALL(1,1).xref(1,:),'-k')
for iM = 1:Nmodels
    plot(ResultsALL(BestModelIDX(iM),iM).t,ResultsALL(BestModelIDX(iM),iM).x(1,:),'Color',ccolors(iM,:),'LineStyle',lstyles{iM},'LineWidth',2)
end
ylim([65 110])
xlabel('Time'), ylabel('Prey')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 100])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_CTRLPERF_',InputSignalType,'_Eta',sprintf('%03g',100*eta_vec),'_TrainingLength_BestModels_x1','.eps']);

figure,
hold on, box on
plot(ResultsALL(BestModelIDX(1),1).t,ResultsALL(1,1).xref(2,:),'-k')
for iM = 1:Nmodels
    plot(ResultsALL(BestModelIDX(iM),iM).t,ResultsALL(BestModelIDX(iM),iM).x(2,:),'Color',ccolors(iM,:),'LineStyle',lstyles{iM},'LineWidth',2)
end
ylim([7 27])
xlabel('Time'), ylabel('Predator')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 100])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_CTRLPERF_',InputSignalType,'_Eta',sprintf('%03g',100*eta_vec),'_TrainingLength_BestModels_x2','.eps']);

figure,
hold on, box on
for iM = 1:Nmodels
    plot(ResultsALL(BestModelIDX(iM),iM).t,cumsum(ResultsALL(BestModelIDX(iM),iM).J),'Color',ccolors(iM,:),'LineStyle',lstyles{iM},'LineWidth',2)
end
if eta_vec == 0
    ylim([7*10^3 9.5*10^3])
elseif eta_vec == 0.05
    ylim([7*10^3 1*10^4])
end
set(gca,'yscale','linear','ytick',[7*10^3,8*10^3,9*10^3],'yticklabels',{'70','80','90'})
xlabel('Time'), ylabel('Cost')
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 100])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_CTRLPERF_',InputSignalType,'_Eta',sprintf('%03g',100*eta_vec),'_TrainingLength_BestModels_J','.eps']);

%% Show State Performance over training lengths
Q = 1*eye(2); %20
Jstate = zeros(length(ResultsALL(iL,iM).J),N_LENGTHS,Nmodels);
for iM = 1:Nmodels
    for iL = 1:N_LENGTHS
        for it = 1:length(ResultsALL(iL,iM).xref(1,:))
            Jstate(it,iL,iM) = abs(ResultsALL(iL,iM).xref(:,it) - ResultsALL(iL,iM).x(:,it))'*Q*abs(ResultsALL(iL,iM).xref(:,it) - ResultsALL(iL,iM).x(:,it));
        end
    end
end

clear ph
ccolors = [1 0 0; 0 1 0; 0.7,0.7,1];
lstyles = {'-', '--', '-.'};
symbols = {'o', 'd', 's'};

Nstart = 20;
xx = Ntrain_vec(Nstart):2:1000;
yModel = zeros(length(xx),Nmodels);
JendAll = zeros(length(N_LENGTHS),Nmodels);
Jend = cell(Nmodels,1);
xJend = cell(Nmodels,1);

for iM = 1:Nmodels
%     Jend{iM} = zeros(length(N_LENGTHS),Nmodels);
    count = 0;
    for iL = 1:N_LENGTHS,
        tmp = cumsum(Jstate(:,iL,iM));
        if tmp(end)<10^10
            count = count + 1;
            xJend{iM}(count) = Ntrain_vec(iL);
            Jend{iM}(count) = tmp(end);
        end
        JendAll(iL,iM) = tmp(end);
    end
    yModel(:,iM) = spline(xJend{iM}(Nstart:end),Jend{iM}(Nstart:end),xx);
end

figure, hold on, box on
for iM = 1:Nmodels
%     if iM < Nmodels
%         plot(xx,yModel(:,iM),'-','Color',ccolors(iM,:), 'LineWidth',1, 'LineStyle',lstyles{iM}), hold on,
%     end
    plot(Ntrain_vec,JendAll(:,iM),symbols{iM},'Color',ccolors(iM,:),'MarkerFaceColor',ccolors(iM,:),'MarkerSize',4)
end
xlim([4 1001])
ylim([5*10^3 5*10^6])
% set(gca,'YScale','log')
set(gca,'XScale','log','YScale','log', 'ytick', [10^4,10^5,10^6], 'xtick', [10^0,10^1,10^2,10^3])
xlabel('Length of Training Data');
ylabel('Cost')
lh = legend(ModelCollection,'Location','NorthEast');
lh.Position = [0.71 0.73 lh.Position(3)-0.05 lh.Position(4)];
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 200])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_CTRLPERF_',InputSignalType,'_Eta',sprintf('%03g',100*eta_vec),'_TrainingLength_CostState','.eps']);

lh.Location = 'NorthOutside';
lh.Orientation = 'horizontal';
lh.Box = 'off';
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_CTRLPERF_',InputSignalType,'_Eta',sprintf('%03g',100*eta_vec),'_TrainingLength_CostState','_legout.eps']);
legend('off')
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_CTRLPERF_',InputSignalType,'_Eta',sprintf('%03g',100*eta_vec),'_TrainingLength_CostState','_noleg.eps']);

return
%% Legend // Dummy plot
figure, hold on, box on
for iM = 1:Nmodels
    plot(Ntrain_vec,ones(size(Ntrain_vec)),'-','Marker',symbols{iM},'Color',ccolors(iM,:),'MarkerFaceColor',ccolors(iM,:))
end

% xlabel('Length of Training Data');
% ylabel('Cost')
axis off
lh = legend(ModelCollection,'Location','NorthEast');
lh.Location = 'NorthOutside';
lh.Orientation = 'horizontal';
lh.Box = 'off';
set(gca,'LineWidth',1, 'FontSize',14)
set(gcf,'Position',[100 100 300 100])
set(gcf,'PaperPositionMode','auto'),
print('-depsc2', '-painters', '-loose', '-cmyk', [figpath,'EX_LOTKA_CTRLPERF_',InputSignalType,'_Eta',sprintf('%03g',100*eta_vec),'_TrainingLength_Legend','.eps']);

return
%% TESTING

% [~,idx_best_DMDc] = min(JendAll(:,1));
% % plot(Ntrain_vec,JendAll(idx_best,2)*ones(size(Ntrain_vec)),'-','Color',ccolors(2,:),'LineWidth',1)
% iM = 1;
% iL = idx_best_DMDc;

%
iM = 3;
iL = find(Ntrain_vec == 85); % DMDc

% iM = 1;
% iL = find(Ntrain_vec == 10); % DMDc

% iM = 3;
% iL = find(Ntrain_vec == 3000); % NARX: 85, 100, 14

figure, hold on
plot(ResultsALL(end,2).x(1,:),ResultsALL(end,2).x(2,:),'-k','LineWidth',1)
plot(ResultsALL(iL,iM).x(1,:),ResultsALL(iL,iM).x(2,:),'--r','LineWidth',1)

figure, hold on
plot(ResultsALL(end,2).t,ResultsALL(end,2).x(2,:),'-k','LineWidth',1)
plot(ResultsALL(end,iM).t,ResultsALL(iL,iM).x(2,:),'--r','LineWidth',1)
plot(ResultsALL(end,2).t,ResultsALL(end,2).x(1,:),'-k','LineWidth',1)
plot(ResultsALL(end,iM).t,ResultsALL(iL,iM).x(1,:),'--r','LineWidth',1)

figure, hold on
plot(ResultsALL(end,2).t,cumsum(ResultsALL(end,2).J),'-k','LineWidth',1)
plot(ResultsALL(end,iM).t,cumsum(ResultsALL(iL,iM).J),'--r','LineWidth',1)

figure, hold on
plot(ResultsALL(end,2).t,ResultsALL(end,2).u,'-k','LineWidth',1)
plot(ResultsALL(end,iM).t,ResultsALL(iL,iM).u,'--r','LineWidth',1)

%% TESTING
figure,hold on
for iM = 1:Nmodels
    Jend = cumsum(ResultsALL(iM,1).J);
    semilogx(Ntrain_vec(iM),Jend(end),'or')
    Jend = cumsum(ResultsALL(iM,2).J);
    semilogx(Ntrain_vec(iM),Jend(end),'og')
    Jend = cumsum(ResultsALL(iM,3).J);
    semilogx(Ntrain_vec(iM),Jend(end),'ob')
end

%%
figure,hold on
for iM = 1:Nmodels
    Jend = cumsum(ResultsALL(iM,1).J);
    plot3(Ntrain_vec(iM).*ones(size(Jend)),ResultsALL(iM,1).t,Jend,'-r')
end

%%
Jend = zeros(Nmodels,1);
for iM = 1:Nmodels
    tmp = cumsum(ResultsALL(iM,1).J); Jend(iM) = tmp(end);
end

figure,plot(Jend)

figure,plot(ResultsALL(3,1).x(1,:))
figure,plot(ResultsALL(3,1).u)

author: eurika-kaiser

repo: SINDY-MPC


text:
function n = fitNormal(data, show_graph)
%FITNORMAL - Fit a plane to the set of coordinates
%
%For a passed list of points in (x,y,z) cartesian coordinates,
%find the plane that best fits the data, the unit vector
%normal to that plane with an initial point at the average
%of the x, y, and z values.
%
% :param data: Matrix composed of of N sets of (x,y,z) coordinates
%              with dimensions Nx3
% :type data: Nx3 matrix
%
% :param show_graph: Option to display plot the result (default false)
% :type show_graph: logical
%
% :return n: Unit vector that is normal to the fit plane
% :type n: 3x1 vector

	if nargin == 1
		show_graph = false;
	end

	for i = 1:3
		X = data;
		X(:,i) = 1;

		X_m = X' * X;
		if det(X_m) == 0
			can_solve(i) = 0;
			continue
		end
		can_solve(i) = 1;

		% Construct and normalize the normal vector
		coeff = (X_m)^-1 * X' * data(:,i);
		c_neg = -coeff;
		c_neg(i) = 1;
		coeff(i) = 1;
		n(:,i) = c_neg / norm(coeff);

	end

	if sum(can_solve) == 0
		error('Planar fit to the data caused a singular matrix.')
		return
	end

	% Calculating residuals for each fit
	center = mean(data);
	off_center = [data(:,1)-center(1) data(:,2)-center(2) data(:,3)-center(3)];
	for i = 1:3
		if can_solve(i) == 0
			residual_sum(i) = NaN;
			continue
		end

		residuals = off_center * n(:,i);
		residual_sum(i) = sum(residuals .* residuals);

	end

	% Find the lowest residual index
	best_fit = find(residual_sum == min(residual_sum));

	% Possible that equal mins so just use the first index found
	n = n(:,best_fit(1));

	if ~show_graph
		return
	end

	range = max(max(data) - min(data)) / 2;
	mid_pt = (max(data) - min(data)) / 2 + min(data);
	xlim = [-1 1]*range + mid_pt(1);
	ylim = [-1 1]*range + mid_pt(2);
	zlim = [-1 1]*range + mid_pt(3);

	L=plot3(data(:,1),data(:,2),data(:,3),'ro','Markerfacecolor','r'); % Plot the original data points
	hold on;
	set(get(L, 'Parent'),'DataAspectRatio',[1 1 1],'XLim',xlim,'YLim',ylim,'ZLim',zlim);

	norm_data = [mean(data); mean(data) + (n' * range)];

	% Plot the original data points
	L=plot3(norm_data(:,1),norm_data(:,2),norm_data(:,3),'b-','LineWidth',3);
	set(get(get(L,'parent'),'XLabel'),'String','x','FontSize',14,'FontWeight','bold')
	set(get(get(L,'parent'),'YLabel'),'String','y','FontSize',14,'FontWeight','bold')
	set(get(get(L,'parent'),'ZLabel'),'String','z','FontSize',14,'FontWeight','bold')
	title(sprintf('Normal Vector: <%0.3f, %0.3f, %0.3f>',n),'FontWeight','bold','FontSize',14)
	grid on;
	axis square;
	hold off;
end