text:
function a = bivand2 ( n, alpha, beta )

%*****************************************************************************80
%
%% BIVAND2 returns a bidimensional Vandermonde2 matrix.
%
%  Discussion:
%
%    N = 3, ALPHA = ( 1, 2, 3 ), BETA = ( 10, 20, 30 )
%
%    (x,y)   | (1,10) (2,10) (3,10) (1,20) (2,20) (3,20) (1,30) (2,30) (3,30)
%    --------+---------------------------------------------------------------
%    1       |     1      1      1      1      1      1      1      1      1
%    x       |     1      2      3      1      2      1      1      2      3
%    x^2     |     1      4      9      1      4      1      1      4      9
%       y    |    10     10     10     20     20     20     30     30     30
%    x  y    |    10     20     30     20     40     60     30     60     90
%    x^2y    |    10     40     90     20     80    180     30    120    270
%       y^2  |   100    100    100    400    400    400    900    900    900
%    x  y^2  |   100    200    300    400    800   1200    900   1800   2700
%    x^2y^2  |   100    400    900    400   1600   3600    900   3600   8100
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 May 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of the data vectors.
%
%    Input, real ALPHA(N), BETA(N), the values that define A.
%
%    Output, real A(N^2,N^2), the matrix.
%
  nn = n * n;
  a = zeros ( nn, nn );

  i = 0
  for iy = 1 : n
    for ix = 1 : n
      i = i + 1;
      j = 0;
      for jy = 1 : n
        for jx = 1 : n
          j = j + 1;
          a(i,j) = alpha(jx) ^ ( ix - 1 ) * beta(jy) ^ ( iy - 1 );
        end
      end
    end
  end

  return
end


text:
function plot_StressVarBeta(vVariance,vBeta)
% function plot_StressVarBeta(vVariance,vBeta)
% --------------------------------------------
% Plot linearly the variance of the stress tensor versus the angle beta,
% representing the difference between the assumed direction of traction and the inverted
%
% Incoming variables:
% vVariance : Vector of Variance values
% vBeta     : Vector of Beta values
%
% j.woessner@sed.ethz.ch
% last update: 24.02.2005

% Deselect nans
vSel = (~isnan(vVariance) & ~isnan(vBeta));
vVariance = vVariance(vSel,:);
vBeta = vBeta.subset(vSel);

% Plot figure
figure
plot(vBeta,vVariance,'k^')
xlabel('\beta [degree]','FontSize',12,'Fontweight','bold')
ylabel('Variance','FontSize',12,'Fontweight','bold')
set(gca,'Linewidth',2,'FontSize',12,'Fontweight','bold')

author: CelsoReyes


text:
function fx = p18_fx ( x )

%*****************************************************************************80
%
%% P18_FX evaluates the function for problem 18.
%
%  Discussion:
%
%    F(X) = 10^14 * (x-1)^7, but is written in term by term form.
%
%    This polynomial becomes difficult to evaluate accurately when
%    written this way.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 October 2011
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Cleve Moler,
%    Numerical Computing with MATLAB,
%    SIAM, 2004,
%    ISBN13: 978-0-898716-60-3,
%    LC: QA297.M625.
%
%  Parameters:
%
%    Input, real X(*), the point at which F is to be evaluated.
%
%    Output, real FX(*), the value of the function at X.
%
  fx = 10.0^14 * ( ...
            x.^7 ...
     -  7 * x.^6 ...
     + 21 * x.^5 ...
     - 35 * x.^4 ...
     + 35 * x.^3 ...
     - 21 * x.^2 ...
     +  7 * x ...
     -  1 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function mesh = mshCube(N,L)
%+========================================================================+
%|                                                                        |
%|                 OPENMSH - LIBRARY FOR MESH MANAGEMENT                  |
%|           openMsh is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : mshCube.m                                     |
%|    #    |   VERSION    : 0.40                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2017                                    |
%|  / 0 \  |   LAST MODIF : 14.03.2018                                    |
%| ( === ) |   SYNOPSIS   : Build uniform mesh for a cube                 |
%|  `---'  |                                                              |
%+========================================================================+

% Optimal number of point for each dimension
n = L/min(L);
n = round( n * (N/prod(n))^(1/3) );

% Delaunay mesh
x = -0.5*L(1) : L(1)/n(1) : 0.5*L(1);
y = -0.5*L(2) : L(2)/n(2) : 0.5*L(2);
z = -0.5*L(3) : L(3)/n(3) : 0.5*L(3);
[x,y,z] = meshgrid(x,y,z);
DT      = delaunayTriangulation([x(:) y(:) z(:)]);

% Build mesh
mesh = msh(DT.Points,DT.ConnectivityList);
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function y = vl_dsigmoid(x)
% VL_DSIGMOID  Derivative of the sigmoid function
%   Y = VL_DSIGMOID(X) returns the derivative of VL_SIGMOID(X). This is
%   calculated as - VL_SIGMOID(X) * (1 - VL_SIGMOID(X)).
%
%   See also: VL_SIGMOID(X), VL_HELP().

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

t = vl_sigmoid(x) ;
y = t .* (1 - t) ;

author: yihui-he

repo: panorama

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


text:
function ui = ic_spike ( x )

%*****************************************************************************80
%
%% IC_SPIKE evaluates the initial condition for a spike function.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(*), the node coordinates.
%
%    Output, real UI(*), the value of the initial condition at each node.
%
  ui = max ( 1.0 - 3.0 * abs ( x ), 0.0 );

  return
end


text:
function h = scimat_imagesc(scimat, varargin)
% SCIMAT_IMAGESC  Display 2D image in real world coordinate axes.
%
% SCIMAT_IMAGESC displays the 2D image in SCIMAT.data scaling the axes to
% real world coordinates using the metadata in SCIMAT.axis.
%
% This function is a thin interface to imagesc(), with the same syntax,
% except that it saves you from having to compute X, Y for
% imagesc(X, Y, SCIMAT.data).
%
% SCIMAT_IMAGESC(SCIMAT)
%
%   SCIMAT is a struct with an image and metadata (see "help scimat").
%
% SCIMAT_IMAGESC(SCIMAT, ...)
%
%   See "help image" for extra options.
%
% See also: image, imagesc.

% Author: Ramon Casero <rcasero@gmail.com>
% Copyright © 2015 University of Oxford
% Version: 0.1.0
%
% University of Oxford means the Chancellor, Masters and Scholars of
% the University of Oxford, having an administrative office at
% Wellington Square, Oxford OX1 2JD, UK.
%
% This file is part of Gerardus.
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. The offer of this
% program under the terms of the License is subject to the License
% being interpreted in accordance with English Law and subject to any
% action against the University of Oxford being under the jurisdiction
% of the English Courts.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% check arguments
nargoutchk(0, 1);

% image size
N = size(scimat.data, 1);
M = size(scimat.data, 2);

if (size(scimat.data, 3) > 1 || size(scimat.data, 4) > 1)
    error('SCIMAT must contain a 2D image, and no more than one frame')
end

% real world coordinates of the first and last voxels
xmin = scimat_index2world([1, 1], scimat);
xmax = scimat_index2world([N M], scimat);

% display image
if (nargout == 0)
    imagesc(...
        linspace(xmin(1), xmax(1), M), ...
        linspace(xmin(2), xmax(2), N), ...
        squeeze(scimat.data), ...
        varargin{:} ...
        )
else
        h = imagesc(...
        linspace(xmin(1), xmax(1), M), ...
        linspace(xmin(2), xmax(2), N), ...
        squeeze(scimat.data), ...
        varargin{:} ...
        );
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [steep] = steep(B,V,C,s,o,d)
%Steepest ascent/descent is a simple and efficient optimization method based on
% statistical design of experiments.
%  In order to optimize a response with respect to a set of quantitative variables,
%  often by using sequential experimentation, we use polynomial models to approximate
%  the response surface. The process often begins by using a 2^p factorial design with
%  additional center points, then fitting a first-order model. If curvature is detected
%  (lack of fit to the first-order model), then we conclude that we are near the optimum
%  and add points to obtain a second-order design. If curvature is not detected, we are
%  far from the optimum and use the method of steepest ascent/descent to get points
%  headed (hopefully) toward to optimum. As we recognize that we are near the optimum,
%  we use a second-order design.
%  If there is no lack-of-fit then the linear model seems to approximate the response
%  surface in this region, which means we are far from an optimum. In this case we want
%  to obtain additional responses (design points) that go closer to the optimum. To do
%  this we use the method of steepest ascent/descent, which just means that future design
%  points should increase/decrease the levels of IV's (factors) in proportion to b's until
%  we are near the optimum response (maximum or minimum).
%  We know that to maximize/minimize a response, the movement of the design center must
%  be in the direction of the directional derivatives of the response function, that is,
%  in the direction of
%
%            df/dx = df/dx_1, . . .,df/dx_p .
%
%  We then multiply by a constant A so that
%
%            Dx = A*df/dx,
%  where
%
%            A = R/sqrt(Sum(df/dx_i)^2) .
%
%  Thus,
%
%            Sum(Dx_i^2) = R^2 .
%
%  For a first order model,
%
%            df/dx_i = b_i .
%
%  So,
%
%            Dx_i = A*b_i    and     A = R/sqrt(Sum(b_i^2)) .
%
%  From this we see that the movement of x_i up the path of steepest ascent/descent
%  is proportional to b_i. Since this is the case it is easier not to pick particular
%  values of R but rather fix a value of b_i and make the other changes proportional
%  to it.
%
%                               Gradient vector
%                  |                ^
%                  |    \    \    \/   O
%                  |\    \    \   /\ New\
%                  | \    \    \ /  O trials
%              x_2 |  \    \    /    \    \
%                  |   \ O--\----\O   \    \
%                  |    \| Start  |    \    \
%                  |     | design |\    \    \
%                  |     O\----\--O \    \    \
%                  |_ _ _ _\_ _ \_ _ \_ _ \_ _ \_ _
%                                 x_1
%
%      First-order response surface and path of steepest ascent.
%
%  Thus, the goal is to optimize the response variable Y. It is assumed that the
%  factors are continuous and controllable by the experimenter with negligible error.
%  If b_1, b_2,...,b_p are the parameters and observed response is Y then the parameter
%  of the polynomial is estimated by method of least squares.
%  The method of steepest ascent/descent is followed to reach the optimun point where
%  best surface is obtained. That is direction of the maximum increase/minimum decrease
%  of the response in the case of surface finish cosiderations. Direction of steep Y is
%  increasing/decreasing. Experiments are conducted along the path of steepest
%  ascent/descent until no further decrease in response is observed. Then a new
%  first-order model may be fit, a new path of steepest ascent/descent determined the
%  procedure continues. Eventually the vicinity of optimum is arrived.
%
%  Syntax: steep(B,V,C,s,o,d)
%
%  Inputs:
%       B - vector of regression coefficients (parameters) of the fitted linear equation
%           (first order): [intercept (1), linear (p)]. It must to be enter in that
%           strictly order.
%       V - vector of the factor values at starting point (current factor setting).
%       C - vector of units increment chosen between the low and medium levels and
%           between the medium and high levels for each parameter.
%       s - number of interested steps.
%       o - change relative to one unit change regarding to the maximum absolute
%           parameter [o = 1 (default)] or parameter choice by experimenter (o = 2).
%       d - steepest ascent [d = 1 (default)] or descent (d = 2) procedure.
%
%  Outputs:
%       A complete summary of the estimation of the selected direction of the
%       selected steepest method.
%
%  To see the procedure consider the following example taken from Montgomery DOX 5E
%  (Example 11-1, pg. 431). Based on the fitted first-order model,
%                     y = 40.44 + 0.775*x1 + 0.325*x2
%  with actual factor-values on the origin of 35 and 155 for time-reaction (min) and
%  temperature (oF) and step size of 5 and 5, respectively. We are interested to obtain
%  the predicted response by the steepest ascent metod for 12 steps.
%
%  Input data must be:
%  B=[40.44 0.775 0.325];
%  V=[35 155];
%  C=[5 2];
%  s = 12;
%  o = 1;
%  d = 1;
%
%  Calling on Matlab the function:
%             steep(B,V,C,s,o,d)
%
%  Answer is:
%
%  Estimation of the selected direction of steepest ascent.
%  The interested steps were 12
%  Leading from the center of the design region the operating point for the IV, 2
%  --------------------------------------------------------------------------------------------------------------
%   Run          Operating point        Predicted response
%  --------------------------------------------------------------------------------------------------------------
%     0        35.000      155.000            40.440
%     1        40.000      157.097            41.351
%     2        45.000      159.194            42.263
%     3        50.000      161.290            43.174
%     4        55.000      163.387            44.085
%     5        60.000      165.484            44.996
%     6        65.000      167.581            45.908
%     7        70.000      169.677            46.819
%     8        75.000      171.774            47.730
%     9        80.000      173.871            48.642
%    10        85.000      175.968            49.553
%    11        90.000      178.065            50.464
%    12        95.000      180.161            51.375
%  --------------------------------------------------------------------------------------------------------------
%
%  Created by A. Trujillo-Ortiz, R. Hernandez-Walls, A. Castro-Perez and
%             F.J. Marquez-Rocha
%             Facultad de Ciencias Marinas
%             Universidad Autonoma de Baja California
%             Apdo. Postal 453
%             Ensenada, Baja California
%             Mexico.
%             atrujo@uabc.mx
%             And the special collaboration of the post-graduate students of the 2005:1
%             Multivariate Statistics Course: Cesar Orlando Chavira-Ortega, Alfredo Frias-Velasco,
%             and Herlinda Gomez-Villa.
%  Copyright (C) May 20, 2005.
%
%  To cite this file, this would be an appropriate format:
%  Trujillo-Ortiz, A., R. Hernandez-Walls, A. Castro-Perez, F.J Marquez-Rocha,
%    C.O. Chavira-Ortega, A. Frias-Velasco and H. Gomez-Villa. (2005). steep:
%    Steepest ascent/descent as a simple and efficient optimization method based
%    on statistical design of experiments. A MATLAB file. [WWW document]. URL
%    http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=7737
%
%  References:
%
%  Box, G. E. P. and Draper, N. R. (1987), Empirical Model-Building and Response Surfaces.
%          John Wiley & Sons: New York.
%  Box, G. E. P. and Wilson, K. B. (1951), On the Experimental Attainment of Optimum
%          Conditions. J. Royal Stat. Soc. Ser. B., 13:1-45.
%  Montgomery, D. C. (2001), Design and Analysis of Experiments, 5th Ed., John Wiley & Sons: New York.
%  Myers, R. H. and Montgomery, D. C. (2002),  Response Surface Methodology: Process and
%          Product Optimization Using Designed Experiments. 2nd. Ed., John Wiley & Sons: New York.
%

if nargin < 6,
   d = 1 %(deafault);
end;

if nargin < 5,
   o = 1 %(deafault);
end;

b = B(2:end);
b = b(:);  %vector of  from coded regression model
V = V(:);  %base level (center point)
C = C(:);  %unit change
s = s;     %
v = length(B)-1; %number of independent variables

if o == 1,
   f = max(abs(b));
else (o == 2),
   l = input('Give me the interested factor: ');
   f = b(l);
end;

S = [];
for i = 0:s,
    if d == 1;      %steepest ascent
        x = [V + i*diag(b./f*C')];
    else (d == 2);  %steepest descent
        x = [V - i*diag(b./f*C')];
    end;
    S = [S x];
end;

S = S';
V = V';
[m n] = size(S);
O  = (S - repmat(V,m,1))./repmat(C',m,1);

Y = B(1)+O*b;

T = [];
for i = 0:s,
   T = [T;i];
end;

if d == 1;
    d = 'ascent.';
else (d == 2);
    d = 'descent.';
end;

[nul,fs] = size(S);
form = ['%12.3f '];
form2 = [' Run '];
form3 = [' Predicted response '];
for k = 2:fs,
    form = [form '%12.3f '];
    form2 = [form2 '       '];
    form3 = ['      ' form3];
end;
RR = [T S Y];

disp(' ')
fprintf('Estimation of the selected direction of steepest %s\n',d);
fprintf('The interested steps were %i\n',s);
fprintf('Leading from the center of the design region the operating point for the IV, %i\n',v);
fprintf('--------------------------------------------------------------------------------------------------------------\n');
disp([form2 '  Operating point '  form3])
fprintf('--------------------------------------------------------------------------------------------------------------\n');
ffa=['%4i  ' form  ' %16.3f\n'];
fprintf(ffa,RR');
fprintf('--------------------------------------------------------------------------------------------------------------\n');
disp(' ')

return,

author: Sable

repo: mcbench-benchmarks


text:
function test_issue1952

% WALLTIME 00:20:00
% MEM 4gb

% DEPENDENCY ft_denoise_ssp

%% PATHS

% path = '/home/lau/Dokumenter/wakeman_henson_face_data/ds117/sub001/MEG';
% dataset = fullfile(path, 'run_01_raw.fif');

dataset = dccnpath('/home/common/matlab/fieldtrip/data/ftp/tutorial/epilepsy/raw/case1/neuromag/case1_cHPI_raw.fif');
if ~exist(dataset, 'file')
  % probably not on filesystem at Donders
  datadir = tempdir;
  s = ftp('ftp.fieldtriptoolbox.org');
  cd(s,'pub/fieldtrip/tutorial/epilepsy/raw/case1/neuromag');
  mget(s, 'case1_cHPI_raw.fif', datadir);
  close(s);
  dataset = fullfile(datadir, 'case1_cHPI_raw.fif');
end

%% SEGMENT

cfg                         = [];
cfg.dataset                 = dataset;

hdr = ft_read_header(dataset, 'checkmaxfilter', false);
fs  = hdr.Fs;

cfg.trl = [(fs:fs:100*fs)+1;(2*fs:fs:101*fs)]';
cfg.trl(:,3) = 0;

% cfg.trialdef.eventtype      = 'STI101';
% cfg.trialdef.eventvalue     = 5;
% cfg.trialdef.prestim        = 0.200;
% cfg.trialdef.poststim       = 0.600;
% cfg.continuous              = 'yes';
% cfg = ft_definetrial(cfg);

% preprocess the data
cfg.channel                 = {'MEG'};
cfg.demean                  = 'yes';     % apply baselinecorrection
%cfg.baselinewindow          = [-0.200 0]; % on basis of mean signal between -0.05 and 0
cfg.checkmaxfilter = false;
data_fif = ft_preprocessing(cfg);

%% TIMELOCK

cfg = [];

timelock = ft_timelockanalysis(cfg, data_fif);

%% PLOT TIMELOCK

cfg = [];
cfg.layout = 'neuromag306mag_helmet.mat';
ft_multiplotER(cfg, timelock);

%% DENOISE WITH SSP

cfg = [];
cfg.ssp = 'all';

denoised_timelock = ft_denoise_ssp(cfg, timelock);

%% DENOISE WITH SSP CELL ARRAY

cfg = [];
% cfg.ssp = {'mag_ssp_upright_fif___pca_mags_v1', ...
%            'mag_ssp_upright_fif___pca_mags_v2', ...
%            'mag_ssp_upright_fif___pca_mags_v3'};
cfg.ssp = {'magssp68iason_fif___pca_v1', ...
  'magssp68iason_fif___pca_v2'};

denoised_timelock_specific = ft_denoise_ssp(cfg, timelock);

%% PLOT SINGLE PLOT

cfg = [];
cfg.channel = 'MEG2611';
cfg.layout = 'neuromag306mag_helmet.mat';

ft_singleplotER(cfg, timelock, denoised_timelock, denoised_timelock_specific);

%% PLOT DENOISED TIMELOCK

cfg = [];
cfg.layout = 'neuromag306mag_helmet.mat';
ft_multiplotER(cfg, timelock, denoised_timelock, denoised_timelock_specific); %% fails

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [node,edge] = ...
        getiso2(xpos,ypos,zdat,ilev,filt)
%GETISO2 extract an iso-contour from a structured two-dimen-
%sional data-set.
%   [NODE,EDGE] = GETISO2(XPOS,YPOS,ZFUN,ZLEV) returns the
%   contour ZFUN(XPOS,YPOS) = ZLEV as a PSLG, by post-proce-
%   ssing the output of the CONTOUR function. The arguments
%   XPOS, YPOS and ZFUN must all be N-by-M arrays, and ZLEV
%   a scalar contouring value.
%
%   See also GETNAN2, FIXGEO2, BFSGEO2, REFINE2

%-----------------------------------------------------------
%   Darren Engwirda : 2018 --
%   Email           : de2363@columbia.edu
%   Last updated    : 25/03/2018
%-----------------------------------------------------------

    if (nargin < +5), filt = +0. ; end

%---------------------------------------------- basic checks
    if ( ~isnumeric(xpos) || ...
         ~isnumeric(ypos) || ...
         ~isnumeric(zdat) || ...
         ~isnumeric(ilev) || ...
         ~isnumeric(filt) )
        error('getiso2:incorrectInputClass' , ...
            'Incorrect input class.') ;
    end

%---------------------------------------------- basic checks
    if (ndims(xpos) ~= +2 || ...
        ndims(ypos) ~= +2 || ...
        ndims(zdat) ~= +2 )
        error('getiso2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

    if (isvector(xpos))
        xnum = length(xpos);
    else
        xnum = size(xpos,2);
    end

    if (isvector(ypos))
        ynum = length(ypos);
    else
        ynum = size(ypos,1);
    end

    if (xnum ~= size(zdat,2) || ...
        ynum ~= size(zdat,1) || ...
        numel(ilev) ~= +1 || ...
        numel(filt) ~= +1 )
        error('getiso2:incorrectDimensions' , ...
            'Incorrect input dimensions.');
    end

%------------------------------------ compute the isocontour
    cmat = contourc( ...
        xpos,ypos,zdat,[ilev,ilev]) ;

%------------------------------------ "walk" contour segment
    node = [] ; edge = [] ; ipos = +1 ;

    while (ipos < size(cmat,2))

        numc = cmat(2,ipos);
        ppts =[cmat(1,ipos+1:ipos+numc)', ...
               cmat(2,ipos+1:ipos+numc)'
              ] ;

        pmin = min(ppts,[],1);
        pmax = max(ppts,[],1);

        pdel = pmax - pmin ;

        if (min(pdel)>=filt)

            if all(ppts(1,:) == ppts(end,:))

    %-------------------------------- closed - back to start
            enew = ...
           [(1:numc-1)',(2:numc-0)'; numc,1] ;

            else

    %-------------------------------- open - dangling endpts
            enew = ...
           [(1:numc-1)',(2:numc-0)'] ;

            end

            enew = ...
            enew + size(node,1);

            node = [node; ppts];
            edge = [edge; enew];

        end

        ipos = ipos + numc + 1 ;

    end

end

author: CHLNDDEV


text:
function frac=dtiClusteringReability(E, reps, nClusters)
%O'Donnell PhD Thesis (2006): measure clustering solution reliability
%We quantified this clustering variability by testing whether pairs
%of data points (paths from tractography) clustered always in the same cluster, never
%in the same cluster, or sometimes in the same and sometimes in different clusters.
%The fraction of pairs that sometimes were in the same cluster, and sometimes in
%different clusters, was measured for various numbers of clusters -- thus
%producing evaluation graph

%"In any case, the fraction of pairs that are inconsistent is less than 5 percent for any number of clusters greater
%than about 75."
nFibers=size(E, 1);
nPairs=nFibers*(nFibers-1)/2;
reps=500; nClusters=23;

for  rep=1:reps
rep
[IDX(:, rep), C(:, : , rep)]=kmeans(E, nClusters);
end
same=zeros(1, reps);
for fiber1=1:size(E, 1)
fiber1
for fiber2=fiber1+1:size(E, 1)
%For those that are sometimes in same, sometimes in differnt
%clusters == "clusters that were not always in the same or always in
%different clusters" . The lowest values show best performance

if sum(IDX(fiber1, :)==IDX(fiber2, :))>0 & sum(IDX(fiber1, :)==IDX(fiber2, :))<reps
same(sum(IDX(fiber1, :)==IDX(fiber2, :)))=            same(sum(IDX(fiber1, :)==IDX(fiber2, :)))+1;
end
end
end

%Fraction of pairs of fibers always in different or always (out of reps) in same
%clusters. "consistency" index.
frac=1-sum(same(:))/nPairs;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function v=v_dlyapsq(a,b)
%V_DLYAPSQ Solves the discrete Lyapunov equation AV'VA' - V'V + BB' = 0
% V is upper triangular with real non-negative diagonal entries
% this is equivalent to v=chol(dlyap(a,b*b')) but better conditioned numerically

%      Copyright (C) Mike Brookes 2002
%      Version: $Id: v_dlyapsq.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[q,s]=schur(a');
[q,s]=rsf2csf(q,s);
[qd,r]=qr(b'*q,0);
% save r for testing
r0=r;
[m,n]=size(r);
u=zeros(n,n);
if m==1
   for i=1:n-1
      in=i+1:n;
      si=s(i,i);
      aa=sqrt(1-si*si');
      u(i,i)=r(1)/aa;
      u(i,in)=(u(i,i)*si'*s(i,in)+aa*r(2:end))/(eye(n-i)-si'*s(in,in));
      r=aa*(u(i,i)*s(i,in)+u(i,in)*s(in,in))-si*r(2:end);
   end
   u(n,n)=r/sqrt(1-s(n,n)*s(n,n)');

else
   w=zeros(m,1); w(m)=1;
   em=eye(m);
   for i=1:n-m
      in=i+1:n;
      si=s(i,i);
      aa=sqrt(1-si*si');
      u(i,i)=r(1,1)/aa;
      u(i,in)=(u(i,i)*si'*s(i,in)+aa*r(1,2:end))/(eye(n-i)-si'*s(in,in));
      vv=aa*(u(i,i)*s(i,in)+u(i,in)*s(in,in))-si*r(1,2:end);
      rr=zeros(m,n-i);
      rr(1:m-1,:)=r(2:end,2:end);
      [qq,r]=qrupdate(em,rr,w,vv');
   end
   for i=n-m+1:n-1
      in=i+1:n;
      si=s(i,i);
      aa=sqrt(1-si*si');
      u(i,i)=r(1,1)/aa;
      u(i,in)=(u(i,i)*si'*s(i,in)+aa*r(1,2:end))/(eye(n-i)-si'*s(in,in));
      vv=aa*(u(i,i)*s(i,in)+u(i,in)*s(in,in))-si*r(1,2:end);
      rr=zeros(n-i+1,n-i);
      rr(1:n-i,:)=r(2:end,2:end);
      [qq,rr]=qrupdate(eye(n-i+1),rr,w(m-n+i:end),vv');
      r=rr(1:n-i,:);
   end

   u(n,n)=r/sqrt(1-s(n,n)*s(n,n)');

end

v=triu(qr(u*q'));
dv=diag(v);
ix=dv~=0;
v(ix,:)=diag(abs(dv(ix))./dv(ix))*v(ix,:);
if isreal(a) & isreal(b)
   v=real(v);
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
%% viewFourthOrderTensor
% Below is a demonstration of the features of the |viewFourthOrderTensor2D| function

%%
clear; close all; clc;

%% Syntax
% |viewFourthOrderTensor(C,numDigits,fontSizeIm,fontSize);|

%% Description
% This function creates the 9x9 and the 6x6 (aka Voigt) array mappings for
% 3D 4th order input tensor C. The function can take up to 4 inputs:
%%
% *  |C|, a fourth order tensor (3x3x3x3)
% *  |numDigits|, the number of decimal places to use for numerical values (default 5)
% *  |fontSizeIm|, the font size in the image (default 15)
% *  |fontSize|, the font size of the axis window (default 15)

%% Examples

%% Viewing fourth-order stiffness tensors

%%
% Creating the stiffness tensor for Hooke's law of linear elasticity

%Constructing 4th order base tensor set
I=eye(3,3); %The 2nd order identity tensor
II1=dyadicProduct(I,I,1); %4th order base tensor 1
II3=dyadicProduct(I,I,3); %4th order base tensor 3

%Parameters for Hooke's law
mu=1; %The shear modulus
lambda=5; %The lambda lame parameter
C=lambda.*II1+2.*mu.*II3; %Construct 4th order stiffness tensor

%%
% Visualizing the tensor using |viewFourthOrderTensor|

viewFourthOrderTensor(C); %Visualize tensor C

%% Viewing fourth-order stiffness tensors with symbolic variables

try
    %%

    syms mu lambda; %Create symbolic Lame parameters

    %Parameters for Hooke's law
    C=lambda.*II1+2.*mu.*II3; %Construct 4th order stiffness tensor

    %%
    % Visualizing the tensor using |viewFourthOrderTensor|

    numDigits=0;
    fontSizeIm=15;
    fontSize=25;
    viewFourthOrderTensor(C,numDigits,fontSizeIm,fontSize); %Visualize tensor C

    %%
    syms mu k; %Create symbolic Lame parameters

    %Construct 4th order stiffness tensor
    k=lambda+2/3*mu; %Bulk modulus

    %Construct 4th order stiffness tensor
    C=(k-2/3*mu).*II1+2.*mu.*II3; %Construct 4th order stiffness tensor

    %%
    % Visualizing the tensor using |viewFourthOrderTensor|

    numDigits=0;
    fontSizeIm=15;
    fontSize=25;
    viewFourthOrderTensor(C,numDigits,fontSizeIm,fontSize); %Visualize tensor C

end

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function m = multiplicity(mori,varargin)
% multiplicity of a misorientation
%
% Syntax
%   m = multiplicity(mori)
%
% Input
%  mori - mis@orientation
%
% Output
%  m - integer
%

m = numSym(mori.CS) * numSym(mori.SS) ./ ...
  length(unique(symmetrise(mori),'noSymmetry'));

author: mtex-toolbox

repo: mtex

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


text:
function q = vgg_quat_from_rotation_matrix( R )
% vgg_quat_from_rotation_matrix Generates quaternion from rotation matrix
%            q = vgg_quat_from_rotation_matrix(R)

q = [	(1 + R(1,1) + R(2,2) + R(3,3))
	(1 + R(1,1) - R(2,2) - R(3,3))
	(1 - R(1,1) + R(2,2) - R(3,3))
	(1 - R(1,1) - R(2,2) + R(3,3)) ];

if ~issym(q)
  % Pivot to avoid division by small numbers
  [b I] = max(abs(q));
else
  % For symbolic quats, just make sure we're nonzero
  for k=1:4
    if q(k) ~= 0
      I = k;
      break
    end
  end
end

q(I) = sqrt(q(I)) / 2 ;

if I == 1
	q(2) = (R(3,2) - R(2,3)) / (4*q(I));
	q(3) = (R(1,3) - R(3,1)) / (4*q(I));
	q(4) = (R(2,1) - R(1,2)) / (4*q(I));
elseif I==2
	q(1) = (R(3,2) - R(2,3)) / (4*q(I));
	q(3) = (R(2,1) + R(1,2)) / (4*q(I));
	q(4) = (R(1,3) + R(3,1)) / (4*q(I));
elseif I==3
	q(1) = (R(1,3) - R(3,1)) / (4*q(I));
	q(2) = (R(2,1) + R(1,2)) / (4*q(I));
	q(4) = (R(3,2) + R(2,3)) / (4*q(I));
elseif I==4
	q(1) = (R(2,1) - R(1,2)) / (4*q(I));
	q(2) = (R(1,3) + R(3,1)) / (4*q(I));
	q(3) = (R(3,2) + R(2,3)) / (4*q(I));
end

author: jmmanley

repo: VGG-Multiple-View-Geometry


text:
function CERRStatusString(str, dispFlag)
%"CERRStatusString"
%   Display a status string in the console and/or in the CERR Gui.
%
%LM:  JOD, 3 Jan 02.
%     JOD, 7 May 03, Changed CERRStatusString format.
%
%   If flag = 'both' or is undefined, status goes to console and CERR
%   GUI.  If flag = 'gui', status goes only to gui.  If flag = 'console'
%   status goes only to console.
%
%Usage:
%   CERRStatusString(str, dispFlag)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global stateS

%Check options to see if StatusString should be muted, only display to the GUI.
if isempty(stateS)
    optS = opts4Exe([getCERRPath,'CERROptions.json']);
elseif isfield(stateS, 'optS')
    optS = stateS.optS;
else
    optS = [];
end

if ~isempty(optS) && isfield(optS, 'CERRStatusStringEnabled')
    if ~optS.CERRStatusStringEnabled
        dispFlag = 'gui';
    end
end

if ~exist('dispFlag','var')
    dispFlag = 'both';
end

if ishandle(1)
    g = gcf;
end

%Also write non-empty strings to the matlab console:
%if strcmpi(dispFlag, 'both') || strcmpi(dispFlag, 'console')
if any(ismember(lower(dispFlag),{'both','console'}))
    if ~strcmp(str,'')
        disp(['CERR>>  ' str])
    end
end

h = [];
%try
if ~isempty(stateS) && isfield(stateS,'handle') && ...
        isfield(stateS.handle,'CERRSliceViewer') && ishandle(stateS.handle.CERRSliceViewer)
    h = stateS.handle.CERRSliceViewer;
end
%end

if ishandle(h)
    set(0,'CurrentFigure',h);

    %if strcmpi(dispFlag, 'both') || strcmpi(dispFlag, 'gui')
    if any(ismember(lower(dispFlag),{'both','gui'}))
        %try
            set(stateS.handle.CERRStatus,'string',str)
        %    drawnow
        %end

        %return to old figure
        if ishandle(1)
            set(0,'CurrentFigure',g);
        end
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%% Machine Learning Online Class - Exercise 1: Linear Regression

%  Instructions
%  ------------
%
%  This file contains code that helps you get started on the
%  linear exercise. You will need to complete the following functions
%  in this exericse:
%
%     warmUpExercise.m
%     plotData.m
%     gradientDescent.m
%     computeCost.m
%     gradientDescentMulti.m
%     computeCostMulti.m
%     featureNormalize.m
%     normalEqn.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%
% x refers to the population size in 10,000s
% y refers to the profit in $10,000s
%

%% Initialization
clear ; close all; clc

%% ==================== Part 1: Basic Function ====================
% Complete warmUpExercise.m
fprintf('Running warmUpExercise ... \n');
fprintf('5x5 Identity Matrix: \n');
warmUpExercise()

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

%% ======================= Part 2: Plotting =======================
fprintf('Plotting Data ...\n')
data = load('ex1data1.txt'); % read comma separated data
X = data(:, 1); y = data(:, 2);
m = length(y); % number of training examples

% Plot Data
% Note: You have to complete the code in plotData.m
plotData(X, y);

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

%% =================== Part 3: Gradient descent ===================
fprintf('Running Gradient Descent ...\n')

X = [ones(m, 1), data(:,1)]; % Add a column of ones to x
theta = zeros(2, 1); % initialize fitting parameters

% Some gradient descent settings
iterations = 1500;
alpha = 0.01;

% compute and display initial cost
computeCost(X, y, theta)

% run gradient descent
theta = gradientDescent(X, y, theta, alpha, iterations);

% print theta to screen
fprintf('Theta found by gradient descent: ');
fprintf('%f %f \n', theta(1), theta(2));

% Plot the linear fit
hold on; % keep previous plot visible
plot(X(:,2), X*theta, '-')
legend('Training data', 'Linear regression')
hold off % don't overlay any more plots on this figure

% Predict values for population sizes of 35,000 and 70,000
predict1 = [1, 1.8] *theta;
fprintf('For population = 18,000, we predict a profit of %f\n',...
    predict1*10000);
predict2 = [1, 2.5] * theta;
fprintf('For population = 25,000, we predict a profit of %f\n',...
    predict2*10000);

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

%% ============= Part 4: Visualizing J(theta_0, theta_1) =============
fprintf('Visualizing J(theta_0, theta_1) ...\n')

% Grid over which we will calculate J
theta0_vals = linspace(-10, 10, 100);
theta1_vals = linspace(-1, 4, 100);

% initialize J_vals to a matrix of 0's
J_vals = zeros(length(theta0_vals), length(theta1_vals));

% Fill out J_vals
for i = 1:length(theta0_vals)
    for j = 1:length(theta1_vals)
	  t = [theta0_vals(i); theta1_vals(j)];
	  J_vals(i,j) = computeCost(X, y, t);
    end
end

% Because of the way meshgrids work in the surf command, we need to
% transpose J_vals before calling surf, or else the axes will be flipped
J_vals = J_vals';
% Surface plot
figure;
surf(theta0_vals, theta1_vals, J_vals)
xlabel('\theta_0'); ylabel('\theta_1');

% Contour plot
figure;
% Plot J_vals as 15 contours spaced logarithmically between 0.01 and 100
contour(theta0_vals, theta1_vals, J_vals, logspace(-2, 3, 20))
xlabel('\theta_0'); ylabel('\theta_1');
hold on;
plot(theta(1), theta(2), 'rx', 'MarkerSize', 10, 'LineWidth', 2);

author: atinesh-s

repo: Coursera-Machine-Learning-Stanford


text:
function y = mldivide(X,Y)
%MLDIVIDE (overloaded)

if (isa(X,'sdpvar'))
  error('Division of matrix variables not possible.')
end

try
  lmi_variables = getvariables(Y);
  nv = length(lmi_variables);
  y  = Y;
  n = Y.dim(1);
  m = Y.dim(2);
  if m==1
    y.basis = X\Y.basis;
    y.dim(1) = size(y.basis,1);
    y.dim(2) = 1;
  else % FIX : VECTORIZE THIS...
    [L,U] = qr(X);
    temp = U\(L\reshape(Y.basis(:,1),n,m));
    y.basis = temp(:);
    for i = 1:nv
        temp = U\(L\reshape(Y.basis(:,i+1),n,m));
        y.basis(:,i+1) = temp(:);
    end;
    y.dim(1) = size(temp,1);
    y.dim(2) = size(temp,2);
  end
catch
  error(lasterr);
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
 function [arg, ii] = arg_get(list, what, default)
%function [arg, ii] = arg_get(list, what, default)
% extract an argument (e.g., nx or pixel_size, etc.)
% from a list (vertically concatenated strings)
% e.g., as generated by arg_pair() or strvcat()
% the *first* match is used.
% Copyright 2005, Jeff Fessler, The University of Michigan

if nargin < 1, ir_usage, end
if nargin == 1 && streq(list, 'test'), arg_get_test, return, end

list(list == '	') = ' ';	% replace tabs with spaces
what = [what ' '];
for ii=1:size(list,1)
	if streq(list(ii,:), what, length(what))
		arg = list(ii,(1+length(what)):end);
		return
	end
end
if isvar('default') && ~isempty(default)
	if isnumeric(default)
		arg = sprintf('%d', default);
	else
		arg = default;
	end
	ii = 0;
else
	fail('cannot find "%s"', what)
end

function arg_get_test
tmp = strvcat('na 9', 'nb 8', 'nxyz 7');
tmp = arg_get(tmp, 'na');
tmp = str2num(tmp);
jf_equal(tmp, 9)

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
clear all; close all; clc

addpath('apm')

s = 'http://byu.apmonitor.com';
c = 'mpc';

% Connect to Arduino
tclab;

% Run time in minutes
run_time = 5.0;

% Number of cycles (1 cycle per second)
loops = round(60*run_time);

% Temperature (degC)
T1 = ones(1,loops) * T1C(); % measured T
Tsp1 = ones(1,loops) * 30;
T2 = ones(1,loops) * T2C(); % measured T
Tsp2 = ones(1,loops) * 23;
time = zeros(1,loops);

% Changes in set point
Tsp1(50:200) = 40;
Tsp1(201:300) = 35;

% milli-volts input
Q1 = zeros(1,loops);
Q2 = zeros(1,loops);

% model predictive control initialization
mpc_init(s,c);

start_time = clock;
prev_time = start_time;

% dynamic plot (note: subplots needs to be declared here first)
figure(1)
subplot(3,1,1)
hold on, grid on
anexp1 = animatedline('LineStyle','-', 'Color', 'k', 'LineWidth', 2);
ansp1 = animatedline('LineStyle','--','Color', 'r','LineWidth', 2);
ylabel('Temperature \circC')
legend('T_1 Measured', 'T_1 Set Point', ...
    'Location', 'northwest')
subplot(3,1,2)
hold on, grid on
anexp2 = animatedline('LineStyle','-', 'Color', 'k', 'LineWidth', 2);
ansp2 = animatedline('LineStyle','--','Color', 'r','LineWidth', 2);
ylabel('Temperature \circC')
legend('T_2 Measured', 'T_2 Set Point', ...
    'Location', 'northwest')
subplot(3,1,3)
hold on, grid on
anQ1 = animatedline('LineStyle','-', 'Color', 'k', 'LineWidth', 2);
anQ2 = animatedline('LineStyle','--', 'Color', 'b', 'LineWidth', 2);
ylabel('Power Level Q (%)')
legend('Q_1', 'Q_2', 'Location', 'northwest')
xlabel('Time (sec)')

for ii = 1:loops
    % Pause Sleep time
    pause_max = 1.0;
    pause_time = pause_max - etime(clock,prev_time);
    if pause_time >= 0.0
        pause(pause_time - 0.01)
    else
        pause(0.01)
    end

    % Record time and change in time
    t = clock;
    dt = etime(t,prev_time);
    if ii>=2
        time(ii) = time(ii-1) + dt;
    end
    prev_time = t;

    % read and record from temperature controller
    T1(ii) = T1C();
    T2(ii) = T2C();

    % model predictive control
    Q1(ii) = mpc2(T1(ii),Tsp1(ii));

    % adjust power level
    h1(Q1(ii));
    h2(Q2(ii));

    % plot
    addpoints(anexp1,time(ii),T1(ii))
    addpoints(ansp1,time(ii),Tsp1(ii))
    addpoints(anexp2,time(ii),T2(ii))
    addpoints(ansp2,time(ii),Tsp2(ii))
    addpoints(anQ1,time(ii),Q1(ii))
    addpoints(anQ2,time(ii),Q2(ii))
    drawnow

    if ii==10
        apm_web(s,c);
    end
end

h1(0);
h2(0);
disp('Heaters off')
% turn off heater but keep LED on if T > 50
if (T1C() || T2C()) > 50
    led(1)
    disp(['Warning, heater temperature 1 =', num2str(T1C())])
    disp(['Warning, heater temperature 2 =', num2str(T2C())])
else
    led(0)
end

% save txt file with data
data = [time',Q1',Q2',T1',T2',Tsp1',Tsp2'];
csvwrite('data.txt',data);

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function fnew=removeisolatedsurf(v,f,maxdiameter)
%
% fnew=removeisolatedsurf(v,f,maxdiameter)
%
% remove disjointed surface fragment filtered by using mesh diameter
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%    v: list of nodes of the input surface
%    f: list of triangles of the input surface
%    maxdiameter: maximum bounding box size for surface removal
%
% ouput:
%    fnew: new face list after removing the components smaller than
%          maxdiameter
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

fc=finddisconnsurf(f);
for i=1:length(fc)
    xdia=v(fc{i},1);
    xdia=max(xdia(:))-min(xdia(:));
    if(xdia<=maxdiameter) fc{i}=[]; continue; end
    ydia=v(fc{i},2);
    ydia=max(ydia(:))-min(ydia(:));
    if(ydia<=maxdiameter) fc{i}=[]; continue; end
    zdia=v(fc{i},3);
    zdia=max(zdia(:))-min(zdia(:));
    if(zdia<=maxdiameter) fc{i}=[]; continue; end
end
fnew=[];
for i=1:length(fc)
	if(length(fc{i})) fnew=[fnew;fc{i}]; end
end
if(size(fnew,1)~=size(f,1))
   fprintf(1,'removed %d elements of small isolated surfaces',size(f,1)-size(fnew,1));
end


text:
%% patchCentre
% Below is a demonstration of the features of the |patchCentre| function

%%
clear; close all; clc;

%% Syntax
% |[Vm]=patchCentre(F,V);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.


text:
% eml_sps_os_test.m
% compare aspire and matlab versions of E-ML-SPS-OS
% Copyright Apr 2000, Jeff Fessler, The University of Michigan

%
% generate data
%
if ~isvar('yi'), printm 'data'
	if has_aspire
		f.dir	= test_dir;
		f.wtf	= [f.dir 't,g.wtf'];
		f.wtr	= strrep(f.wtf, 'wtf', 'wtr');
		f.yi	= [f.dir 'yi.fld'];
		f.ci	= [f.dir 'ci.fld'];
		f.ri	= [f.dir 'ri.fld'];
		f.mask	= [f.dir 'mask.fld'];
	end
	em_test_setup
%	em3_test_setup, n.a = n.n2; % kludge
prompt
end

if ~isvar('Gb'), printm 'make Gb'
	f.nblock = 5;
	Gb{1} = Gblock(G, 1);
	Gb{2} = Gblock(G, f.nblock);
prompt
end

%
% matlab iterations
%
if ~isvar('xmat'), printm 'matlab E-ML-SPS-OS'
	f.niter = 9;
	f.pixmax = 6;

	f.curvs = {'oc', 'pc'};

	xinit = ones(size(xtrue));
	for ic = 1:length(f.curvs)
		tmp = eql_sps_os(xinit(ig.mask), Gb{ic}, yi, ci, ri, [], ...
			f.niter, f.pixmax, f.curvs{ic});
		xmat{ic} = ig.embed(tmp);
	end
	im clf, im(xmat{1}, 'Matlab E-ML-SPS-OS iterations')
prompt
end

if ~has_aspire, return, end

%
% aspire iterations
%
if ~isvar('xasp'), printm 'aspire E-ML-SPS-OS'

	f.init	= [f.dir 'init.fld'];
	f.out	= [f.dir 'out.fld'];
	fld_write(f.init, xinit, 'check', 0)

	f.saver	= 'stack,1';
	f.fitype = ['2z@' f.wtr '@-'];

	for ic = 1:length(f.curvs)
		if exist(f.out, 'file'), delete(f.out), end
		f.alg = sprintf('ospsc,%s,%d,%d,1,0', f.curvs{ic}, Gb{ic}.nblock, sg.na);
		f.penal	= sprintf('%g,quad,0,-', -100);
		f.method = sprintf('@%d@%s@%s', f.niter-1, f.alg, f.penal);
		f.com = sprintf(['i -chat 0 empl3 %s %s  %s %s 1 %s 1 %s -' ...
				' %s %s 0 1 %g 0 -'], ...
			f.out, f.init, f.yi, f.ci, f.ri, f.fitype, ...
			f.method, f.saver, f.pixmax);
		os_run(f.com)

		xasp{ic} = double(fld_read(f.out));
	end
end

if 1
	for ic = 1:length(f.curvs)
		t = vcorrcoef(xasp{ic}, xmat{ic});
		printf('corr. [curv=%s block=%d] %g,%g', ...
			f.curvs{ic}, Gb{ic}.nblock, t, t-1)
	end
	ic = 2;

	im clf, im(221, xmat{ic}, 'xhat matlab')
	im(222, xasp{ic}, 'xhat aspire')
	im(223, (xasp{ic}-xmat{ic})/max(col(xmat{ic})), 'aspire-matlab'), cbar

	t1 = eql_obj(xmat{ic}, G, yi(:), ci(:), ri(:), [], ig.mask);
	t2 = eql_obj(xasp{ic}, G, yi(:), ci(:), ri(:), [], ig.mask);

	if im
		subplot(224)
		plot(0:f.niter-1, t1-t1(1), '-o', 0:f.niter-1, t2-t1(1), '-x')
		xlabel iteration, ylabel '\Phi change', legend('mat', 'asp', 4)
		title(sprintf('E-ML-SPS-OS, Nsubset=%d', Gb{ic}.nblock))
	end
end


text:
function h = set(h,varargin)
%SET: Set drumplot properties
%   val = get(drumplot,prop_name,prop_val)
%
%   Valid property names:
%       WAVE, CATALOG, MPL, TRACE_COLOR, EVENT_COLOR, SCALE
%
%   See also DRUMPLOT, DRUMPLOT/SET
%
% Author: Dane Ketner, Alaska Volcano Observatory
% Modified: Glenn Thompson 2016/04/20
% $Date$
% $Revision$

if ~isa(h,'drumplot')
   error('DRUMPLOT/SET: Not a valid drumplot object')
end

if ~rem(nargin-1,2) == 0
      error(['DRUMPLOT/SET: Arguments after h must appear in ',...
             'property name/val pairs'])
end

for n = 1:2:nargin-2
   nw = numel(h.wave); % Number of waveforms in h
   pn = varargin{n};   % Property name
   pv = varargin{n+1}; % Property value
   nv = numel(pv);     % Number of elements in property value
   sv = size(pv);      % Size of property value

   switch lower(pn)
      case 'wave'
         if ~isa(pv,'waveform')
            error('DRUMPLOT/SET: wave property is not a waveform')
         end
         if numel(pv)~=nw
            error(['DRUMPLOT/SET: new waveform dimensions must match',...
                   ' existing waveform dimensions'])
         end
         h.wave = pv;

      case 'catalog'
         if isa(pv,'Catalog')
             h.catalog = pv;
         elseif isempty(pv)
            h.catalog = Catalog();
         else
            error('DRUMPLOT/SET: wrong Catalog format')
         end

      case 'mpl'
         if isnumeric(pv) && nv==1
            h.mpl = pv;
         else
            error('DRUMPLOT/SET: wrong mpl format')
         end

      case 'trace_color'
         if sv(2)==3 && isnumeric(pv)
            pv = {pv};
            nv = 1;
         end
         if (nw==nv) && iscell(pv)
            for m = 1:nv
               if size(pv{m},2)==3
                  for mm = 1:numel(pv{m})
                     if ~(pv{m}(mm)<=1 && pv{m}(mm)>=0)
                        error(['DRUMPLOT/SET: trace_color RGB values must ',...
                           'be between 0 and 1'])
                     end
                  end
               end
            end
            h.trace_color = pv;
         else
            error('DRUMPLOT/SET: wrong trace_color format')
         end

      case 'event_color'
         h.event_color = pv;

      case 'scale'
         if isnumeric(pv) && nv==nw
            h.scale = pv;
         elseif isnumeric(pv) && nv==1
            h.scale = pv*ones(1,nw);
         else
            error('DRUMPLOT/SET: wrong scale format')
         end

      otherwise
         error('DRUMPLOT/GET: Not a valid drumplot property')
    end
end

author: geoscience-community-codes

repo: GISMO


text:
classdef CarlaEnvironment < matlab.System & matlab.system.mixin.Propagates
    % Untitled Add summary here
    %
    % This template includes the minimum set of functions required
    % to define a System object with discrete state.

    % Public, tunable properties
    properties

    end

    properties(DiscreteState)

    end

    % Pre-computed constants
    properties(Access = private)
        % Semantic segmentation resolution
        width = 960;
        height = 540;

        % RGB resolution
        rgb_width = 1920;
        rgb_height = 1080;

        tesla;
        last_position;

        ss_rgb;
        module_ss_rgb;

        rgb;
        module_rgb;

        last_update_time;
        history_array = zeros(2, 8);
        steer_spline;

        steer_window_avg = 0;
        push_back_time = 0;
        frames_since_last_update = 0;

        bias = -0.033;
    end

    methods(Access = protected)
        function setupImpl(obj)
            % Perform one-time calculations, such as computing constants
            port = int16(2000);
            client = py.carla.Client('localhost', port);
            client.set_timeout(10.0);
            world = client.get_world();

            % Spawn Vehicle
            blueprint_library = world.get_blueprint_library();
            car_list = py.list(blueprint_library.filter("model3"));
            car_bp = car_list{1};
            spawn_point = py.random.choice(world.get_map().get_spawn_points());
%             spawn_point.location.x = -88.8;
%             spawn_point.location.y = 149.0;
%             spawn_point.rotation.yaw = 90;

            obj.tesla = world.spawn_actor(car_bp, spawn_point);
            obj.tesla.set_autopilot(false);

            control = obj.tesla.get_control();
            control.throttle = 0.6;
            control.steer = 0;
            obj.tesla.apply_control(control);

            % Semantic Segmentation
            blueprint = world.get_blueprint_library().find('sensor.camera.semantic_segmentation');
            blueprint.set_attribute('image_size_x', num2str(obj.width))
            blueprint.set_attribute('image_size_y', num2str(obj.height))

            transform = py.carla.Transform(py.carla.Location(pyargs('x',0.8, 'z',1.7)));
            obj.ss_rgb = world.spawn_actor(blueprint, transform, pyargs('attach_to',obj.tesla));

            obj.module_ss_rgb = sensorBind(obj.ss_rgb, "ss_rgb", "semantic_segmentation_rgb", "array");

            % RGB
            blueprint = world.get_blueprint_library().find('sensor.camera.rgb');
            blueprint.set_attribute('image_size_x', num2str(obj.rgb_width))
            blueprint.set_attribute('image_size_y', num2str(obj.rgb_height))

            transform = py.carla.Transform(py.carla.Location(pyargs('x',-5.0, 'z',2.5)));
            obj.rgb = world.spawn_actor(blueprint, transform, pyargs('attach_to',obj.tesla));

            obj.module_rgb = sensorBind(obj.rgb, "rgb", "rgb", "array");

            obj.last_position = obj.tesla.get_location();

            obj.last_update_time = cputime;
        end

        function [SEMANTIC_SEGMENTATION_RGB, RGB] = stepImpl(obj)

            % Semantic Segmentation
            SEMANTIC_SEGMENTATION_RGB = uint8(py.getattr(obj.module_ss_rgb, 'array'));

            %% Lane Detection
            control = obj.tesla.get_control();

            % Crop the image to region of interest
            horizontal_crop = 0.0;
            hor_start_pos = ceil(obj.width/2);
            hor_end_pos = floor(obj.width * (1 - horizontal_crop / 2));

            vertical_crop = .5;
            vert_start_pos = 1 + ceil(obj.height * vertical_crop);

            filtered_lanes = SEMANTIC_SEGMENTATION_RGB(vert_start_pos:obj.height,hor_start_pos:hor_end_pos,1) == 157;

            % Apply hough transformations to get the lines from the sensor
            % image
            BW = filtered_lanes;
            [H,T,R] = hough(BW);
            P  = houghpeaks(H,5,'threshold',ceil(0.0*max(H(:))));
            lines = houghlines(BW,T,R,P,'FillGap',50,'MinLength',20);

            % Calculate the angles of the lines based on the x-y
            % coordinates
            angles = zeros(1, length(lines));
            for k = 1:length(lines)
                angles(k) = rad2deg(atan((lines(k).point1(2) - lines(k).point2(2))/(lines(k).point1(1) - lines(k).point2(1))));
            end

            % Normally from the camera position, the angle between the road
            % and the car if the car is moving straight
            ref_angle = 57;
            max_angle = max(angles);

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Loop to draw the nearest right lanes detected
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            imshow(filtered_lanes), hold on

            for k = 1:length(lines)
               angle = rad2deg(atan((lines(k).point1(2) - lines(k).point2(2))/(lines(k).point1(1) - lines(k).point2(1))));
               if ~(angle == max_angle)
                  continue
               end

               xy = [lines(k).point1; lines(k).point2];
               plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');

               % Plot beginnings and ends of lines
               plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');
               plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red');

            end

            set(gcf, 'Visible', 'on');
            hold off;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            % If the vehicle has just spawned and not started moving, dont
            % update steer history
            started = false;
            current_location = obj.tesla.get_location();
            thres_loc = 0.001;

            if abs(current_location.x - obj.last_position.x) > thres_loc...
                    && abs(current_location.y - obj.last_position.y) > thres_loc
                started = true;
            end

            % Calculate the time it took to process one frame. Really
            % necessary to make the calculations independent of the frame
            % rate. Otherwise, running simulation on a pc that can run it
            % on 20 FPS will steer 20 times vs 2 FPS steers 2 times, dt
            % scales that
            dt = cputime - obj.last_update_time;
            fac = 0.5 * tanh(0.1 * abs(ref_angle - max_angle)) * dt;

            obj.last_update_time = cputime;

            % How much effect the last steering has on the future steer
            mem = 0.0025;

            % Filter out lines that are almost horizontal
            if ~isempty(fac) && ((max_angle < 85) && (max_angle > 25)) && started
                if max_angle < ref_angle
                    steer = mem * control.steer + fac;
                else
                    steer = mem * control.steer - fac;
                end
            else
                steer = obj.bias;
            end

            % Update the spline history ~ every .2 seconds. A very high FPS
            % would polute the spline fitting
            if cputime - obj.push_back_time > .2

                % Delete the oldest element in the array
                obj.history_array(:,1) = [];

                % Most recent history weight
                weight = 1.25;

                % Average of the steer since the last update
                steer = (obj.steer_window_avg + steer)/(obj.frames_since_last_update + 1);

                % Add to the histroy
                obj.history_array = [obj.history_array, [cputime;steer]];

                % Pre-processing for the curve fitting
                [timeData, steerData] = prepareCurveData( obj.history_array(1,:),...
                                                          obj.history_array(2,:));

                obj.history_array(end) = obj.history_array(end) * weight;

                % Fit a smoothing spline to the data
                obj.steer_spline = fit( timeData, steerData, 'smoothingspline', ...
                                        'Normalize', 'on', 'SmoothingParam',0.95);

                obj.history_array(end) = obj.history_array(end) / weight;

                obj.steer_window_avg = 0;
                obj.frames_since_last_update = 0;

                obj.push_back_time = cputime;

            else
                % If not yet time to update, keep track of the steer values
                % in this time slot
                obj.steer_window_avg = obj.steer_window_avg + steer;
                obj.frames_since_last_update = obj.frames_since_last_update + 1;
            end

            % Predict and steer the car
            control.steer = obj.steer_spline(cputime + dt/2);

            obj.tesla.apply_control(control);
            obj.last_position = current_location;

            %% RGB sensor
            RGB = uint8(py.getattr(obj.module_rgb, 'array'));
        end

        function [SEMANTIC_SEGMENTATION_RGB, RGB] = isOutputComplexImpl(~)
            SEMANTIC_SEGMENTATION_RGB = false;
            RGB = false;
        end

        function [SEMANTIC_SEGMENTATION_RGB, RGB] = getOutputSizeImpl(obj)
            SEMANTIC_SEGMENTATION_RGB = [obj.height obj.width 3];
            RGB = [obj.rgb_height obj.rgb_width 3];
        end

        function [SEMANTIC_SEGMENTATION_RGB, RGB] = getOutputDataTypeImpl(~)
            SEMANTIC_SEGMENTATION_RGB = 'uint8';
            RGB = 'uint8';
        end

        function [SEMANTIC_SEGMENTATION_RGB, RGB] = isOutputFixedSizeImpl(~)
            SEMANTIC_SEGMENTATION_RGB = true;
            RGB = true;
        end

        function resetImpl(~)
            % Initialize / reset discrete-state properties
        end
    end

    methods(Access= public)
        function delete(obj)
            % Delete the car from the Carla world
            if ~isempty(obj.tesla)
                obj.tesla.destroy();
            end

            % Semantic segmentation sensor
            if ~isempty(obj.ss_rgb)
                obj.ss_rgb.destroy();
            end

            % RGB camera
            if ~isempty(obj.rgb)
                obj.rgb.destroy();
            end
        end
    end
end

author: darkscyla

repo: MATLAB-Carla-Interface


text:
function h = supportFunction(polygon, varargin)
%SUPPORTFUNCTION Compute support function of a polygon
%
%   H = supportFunction(POLYGON, N)
%   uses N points for suport function approximation
%
%   H = supportFunction(POLYGON)
%   assume 24 points for approximation
%
%   H = supportFunction(POLYGON, V)
%   where V is a vector, uses vector V of angles to compute support
%   function.
%
%   See also:
%   polygons2d, convexification
%
%   ---------
%   author : David Legland
%   INRA - TPV URPOI - BIA IMASTE
%   created the 20/12/2004.
%

N = 24;
u = 0:2*pi/N:2*pi*(1-1/N);

if length(varargin)==1
    var = varargin{1};
    if length(var)==1
        N = var;
        u = 0:2*pi/N:2*pi*(1-1/N);
    else
        u = var;
    end
end

% ensure u vertical vector
if size(u, 1)==1
    u=u';
end

h = zeros(size(u));

for i=1:length(u)

    v = repmat([cos(u(i)) sin(u(i))], [size(polygon, 1), 1]);

    h(i) = max(dot(polygon, v, 2));
end

author: rpng

repo: lips

save_path: github-repos/MATLAB/rpng-lips


text:
disp('================================== transform ==============================')
transform

disp('================================== trajectory ==============================')
trajectory

disp('================================== jacobian ==============================')
jacobian

disp('================================== kinematics ==============================')
kinematics

disp('================================== dynamics ==============================')
dynamics

disp('================================== demos ==============================')
pause off
rttrdemo
rtfkdemo
rtikdemo
rttgdemo
rtandemo
rtjademo
rtfddemo
rtidemo

author: petercorke

repo: robotics-toolbox-matlab


text:
function square_exactness_test01 ( )

%*****************************************************************************80
%
%% SQUARE_EXACTNESS_TEST01: product Gauss-Legendre rules for Legendre 2D integral.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 May 2014
%
%  Author:
%
%    John Burkardt
%
  a(1:2) = -1.0;
  b(1:2) = +1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SQUARE_EXACTNESS_TEST01\n' );
  fprintf ( 1, '  Product Gauss-Legendre rules for the 2D Legendre integral.\n' );
  fprintf ( 1, '  Density function rho(x) = 1.\n' );
  fprintf ( 1, '  Region: -1 <= x <= +1.\n' );
  fprintf ( 1, '  Region: -1 <= y <= +1.\n' );
  fprintf ( 1, '  Level: L\n' );
  fprintf ( 1, '  Exactness: 2*L+1\n' );
  fprintf ( 1, '  Order: N = (L+1)*(L+1)\n' );

  for l = 0 : 5

    n_1d = l + 1;
    n = n_1d * n_1d;
    t = 2 * l + 1;

    [ w, x, y ] = legendre_2d_set ( a, b, n_1d, n_1d );

    p_max = t + 1;
    legendre_2d_exactness ( a, b, n, x, y, w, p_max );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% In this example we consider a 127 dimensional space.
% We choose vectors in this space by removing a row from a 128 x 128 Hadamard matrix.
% We then construct 2 d-dimensional subspaces by choosing vectors from this set
% randomly.
% We then measure the principal angle between such a pair of subspaces.
% We repeat it for t trials of random d-dimensional subspaces.
% We vary d from 1 to 64.
% Finally we plot the variation of principal angle with number of dimensions.
% We note that as the subspace dimension increases, the principal angle between
% two randomly chosen subspaces [with basis vectors from the Hadamard matrix]
% keeps decreasing.

clc; clear all;close all;
% Construct a Hadamard matrix
n = 128;
h = hadamard(n);
% drop one row
h2 = h(1:end-1, :);
% normalize the vectors
h2 = spx.norm.normalize_l2(h2);
% compute it's Gram matrix
g = h2' * h2;
% verify that absolute value of all the non-zero elements is 1/(n -1).

trials = 2;

least_angles = [];
ds = [];

% The size of subspace
for d = 1:n/2
    fprintf('Choosing subspace dimension: %d\n', d);
    for t=1:trials
        % select two independent subspaces
        indices = randperm(n, 2*d);
        A = h2(:, indices(1:d));
        B = h2(:, indices(d+1:end));
        angles = spx.la.spaces.principal_angles_cos(A, B);
        least_angle = angles(1);
        least_angle_deg = rad2deg(acos(least_angle));

        fprintf('Trial: %d, smallest principal angle: %.4f, %.2f degrees\n', t, ...
            least_angle, least_angle_deg);
    end
    ds = [ds d];
    least_angles = [least_angles least_angle_deg];
end
plot(ds, least_angles);
xlabel('D');
ylabel('Angle (degrees)');
grid on;

author: indigits

repo: sparse-plex


text:
function [A,B,D,E,I]=twt(X);
% subprogram for Twin T circuit
%
N=3; % N = no. of L's & C's
M=1; % M = no. of indep inputs
U=2; % U = no. of dep nodes in dc equivalent circuit
Y=2; % Y = output node in dc equiv ckt = V3
%
% X = [R1 R3 R5 R7 C2 C4 C6];
%      1  2  3  4  5  6  7
%
R1=X(1);R3=X(2);R5=X(3);R7=X(4);C2=X(5);C4=X(6);C6=X(7);
% Create array space
A1=zeros(U+N);B2=zeros(U+N,N+M);V=zeros(U+N,N+M);H=zeros(N,M+N);
%
% Build A1 matrix.
%
A1=[1/R5 0 0 -1 -1;
   0 -1/R3  1 0 0;
   0 1/R3+1/R7 0 0 1;
   1 0 0 0 0;
   -1 1 0 0 0];
%
% Fill in B2 array
E2=1;E4=1;E6=1;Ein=1;
%
B2=[0 0 0 0;
   -E2*(1/R1+1/R3) 0 0 Ein/R1;
   E2/R3 0 0 0;
   0 -E4 0 Ein;
   0 0 E6 0];
%
P=diag([C2 C4 C6]);
%
% As stated previously, the following code is
% the same for every circuit.
%
V=A1\B2;H=V(U+1:U+N,1:N+M);I=eye(N);
AB=P\H;A=AB(1:N,1:N);B=AB(1:N,N+1:N+M);
D=V(Y:Y,1:N);E=V(Y:Y,N+1:N+M);

author: Sable

repo: mcbench-benchmarks


text:
function d = delimiter()

if ispc
    d = ';';
else
    d = ':';
end

author: BUNPC


text:
function title = p15_title ( )

%*****************************************************************************80
%
%% P15_TITLE returns a title for problem 15.
%
%  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:
%
%    Output, string TITLE, a title for the problem.
%
  title = '#15: Sangbum Kim''s forward step.';

  return
end


text:
function out = DN_Burstiness(y)
% DN_Burstiness     Burstiness statistic of a time series
%
% Returns the 'burstiness' statistic from
%
% Goh and Barabasi, 'Burstiness and memory in complex systems' Europhys. Lett.
% 81, 48002 (2008).
%
%---INPUT:
% y, the input time series
%
%---OUTPUT:
% The burstiness statistic, B.

% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation, either version 3 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program. If not, see <http://www.gnu.org/licenses/>.
% ------------------------------------------------------------------------------

r = std(y)/mean(y); % coefficient of variation

%-------------------------------------------------------------------------------
% Original Goh and Barabasi burstiness statistic, B:
out.B = (r - 1)/(r + 1);
% B = (std(y) - mean(y))/(std(y) + mean(y));

%-------------------------------------------------------------------------------
% Improved burstiness statistic, accounting for scaling for finite time series
% Kim and Jo, 2016, http://arxiv.org/pdf/1604.01125v1.pdf
N = length(y);
out.B_Kim = (sqrt(N+1)*r - sqrt(N-1))/((sqrt(N+1)-2)*r + sqrt(N-1));

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
% Solves the proximity problem associated with the 1- or 2-dimensional Total Variation Lp norm.
% Depending on the dimension and norm of choice, a different algorithm is used for the
% optimization.
% Currently only p=1,2 are accepted
% Note p=inf is slow.
%
% Inputs:
%   - y: input of the proximity operator.
%   - lambda: premultiplier of the norm.
%   - p: norm.
%   - [threads]: number of threads (default 1). Used only for 2-D signals.
%
% Outputs:
%   - x: solution of the proximity problem.
%   - info: statistical info of the run algorithm:
%       info.iters: number of iterations run (major iterations for the 2D case)
%       info.stop: value of the stopping criterion.
function [x,info] = prox_TVLp(y,lambda,p,threads)
    % Check inputs
    if nargin < 4, threads=1; end;

    % Choose an algorithm depending on the norm and the dimension of the input
    if isvector(y)
        switch p
            % L1: Projected Newton with Armijo step
            case 1
                [x,in] = solveTV1_PNc(y,lambda);
                info.iters = in(1);
                info.stop = in(2);
            % L2: Hybrid More-Sorensen + Projected Gradient
            case 2
                [x,in] = solveTV2_morec2(y,lambda);
                info.iters = in(1);
                info.stop = in(2);
            otherwise
                fprintf(1,'ERROR (prox_TVLp): unacceptable norm p=%d, returning input.\n',p);
                x = y;
        end
    else if length(size(y)) == 2
        [x,in] = solveTVgen_PDykstrac(y,lambda*[1 1],[1 2],[p p],threads);
        info.iters = in(1);
        info.stop = in(2);
    else
        fprintf(1,'ERROR (prox_TVLp): for (N>2)-dimensional inputs use prox_TVgen function. Returning input.\n');
        x = y;
    end
end

author: beckel

repo: nilm-eval

save_path: github-repos/MATLAB/beckel-nilm-eval


text:
function diagnostic = solvesdp(varargin)
%SOLVESDP Obsolete command, please use OPTIMIZE

persistent CACHED_SOLVERS
persistent allsolvers
persistent EXISTTIME
persistent NCHECKS

yalmiptime = clock; % Let us see how much time we spend

% *********************************
% CHECK INPUT
% *********************************
nargin = length(varargin);

% First check of objective for early transfer to multiple solves
if nargin>=2
    if isa(varargin{2},'double')
        varargin{2} = [];
    elseif isa(varargin{2},'sdpvar') && numel(varargin{2})>1
        % Several objectives
        diagnostic = solvesdp_multiple(varargin{:});
        return
    end
end

% Early jump to bisection solver
if nargin >= 3
    if isa(varargin{3},'struct')
        if isfield(varargin{3},'solver')
            if isequal(varargin{3}.solver,'bisection')
                if any(is(varargin{1},'sos'))
                    [F_sos,h_sos] = compilesos(varargin{1},varargin{2},varargin{3});
                    varargin{1} = F_sos;
                    varargin{2} = h_sos;
                end
                diagnostic = bisection(varargin{:});
                return
            end
        end
    end
end

if nargin<1
    help solvesdp
    return
else
    F = varargin{1};
    if isa(F,'constraint')
        F = lmi(F);
    end
    if isa(F,'lmi')
        F = flatten(F);
    end

    if isa(F,'sdpvar')
        % We do allow sloppy coding of logic constraints, i.e writing a
        % constraints as [a|b true(a)]
        Fnew = [];
        for i = 1:length(F)
            if length(getvariables(F(i)))>1
                Fnew = nan;
                break
            end
            operator = yalmip('extstruct',getvariables(F(i)));
            if isempty(operator)
                Fnew = nan;
                break
            end
            if length(operator)>1
                Fnew = nan;
                break
            end
            if ~strcmp(operator.fcn,'or')
                Fnew = nan;
                break
            end
            Fnew = Fnew + (true(F(i)));
        end
        if isnan(Fnew)
            error('First argument (F) should be a constraint object.');
        else
            F = Fnew;
        end
    elseif isempty(F)
        F = lmi([]);
    elseif ~isa(F,'lmi')
        error('First argument (F) should be a constraint object.');
    end
end

if nargin>=2
    h = varargin{2};
    if isa(h,'double')
        h = [];
    end
    if ~(isempty(h) | isa(h,'sdpvar') | isa(h,'logdet') |  isa(h,'ncvar'))
        if isa(h,'struct')
            error('Second argument (the objective function h) should be an sdpvar or logdet object (or empty). It appears as if you sent an options structure in the second argument.');
        else
            error('Second argument (the objective function h) should be an sdpvar or logdet object (or empty).');
        end
    end
    if isa(h,'logdet')
        logdetStruct.P     = getP(h);
        logdetStruct.gain  = getgain(h);
        h = getcx(h);
        if isempty(F)
            F = ([]);
        end
    else
        logdetStruct = [];
    end
else
    logdetStruct = [];
    h = [];
end

if ~isempty(F)
    if any(is(F,'sos'))
        diagnostic = solvesos(varargin{:});
        return
    end
end

if isa(h,'sdpvar')
    if is(h,'complex')
        error('Complex valued objective does not make sense.');
    end
end

if nargin>=3
    options = varargin{3};
    if ~(isempty(options) | isa(options,'struct'))
        error('Third argument (options) should be an sdpsettings struct (or empty).');
    end
    if isempty(options)
        options = sdpsettings;
    end
else
    options = sdpsettings;
end
options.solver = lower(options.solver);

% If user has logdet term, but no preference on solver, we try to hook up
% with SDPT3 if possible.
if ~isempty(logdetStruct)
    if strcmp(options.solver,'')
     %   options.solver = 'sdpt3,*';
    end
end

% Call chance solver?
% if length(F) > 0
%     rand_declarations = is(F,'random');
%     if any(rand_declarations)
%     %    diagnostic = solverandom(F(find(~rand_declarations)),h,options,recover(getvariables(sdpvar(F(find(unc_declarations))))));
%         return
%     end
% end

% Call robust solver?
if length(F) > 0
    unc_declarations = is(F,'uncertain');
    if any(unc_declarations)
        try
            diagnostic = solverobust(F(find(~unc_declarations)),h,options,recover(getvariables(sdpvar(F(find(unc_declarations))))));
            return
        catch
            if strfind(lasterr,'Undefined function ''solverobust'' ')
                display('***');
                display('It looks like you have failed to add yalmip/modules/robust to your path');
                display('Information on required paths https://yalmip.github.io/tutorial/installation/');
                display(['For complete path, use addpath(genpath(''' fileparts(which('yalmiptest.m')) '''))']);
                display('***');
                error(lasterr)
            else
                error(lasterr)
            end
        end
    end
end

if isequal(options.solver,'mpt') | nargin>=4
    solving_parametric = 1;
else
    solving_parametric = 0;
end

% Just for safety
if isempty(F) & isempty(logdetStruct)
    F = lmi;
end

if any(is(F,'sos'))
    error('You have SOS constraints. Perhaps you meant to call SOLVESOS.');
end

% Super stupido
if length(F) == 0 & isempty(h) & isempty(logdetStruct)
   diagnostic.yalmiptime = 0;
   diagnostic.solvertime = 0;
   diagnostic.info = 'No problems detected (YALMIP)';
   diagnostic.problem = 0;
   diagnostic.dimacs = [NaN NaN NaN NaN NaN NaN];
   return
end

% *************************************************************************
%% LOOK FOR AVAILABLE SOLVERS
% Finding solvers can be very slow on some systems. To alleviate this
% problem, YALMIP can cache the list of available solvers.
% *************************************************************************
if (options.cachesolvers==0) | isempty(CACHED_SOLVERS)
    getsolvertime = clock;
    solvers = getavailablesolvers(0,options);
    getsolvertime = etime(clock,getsolvertime);
    % CODE TO INFORM USERS ABOUT SLOW NETWORKS!
    if isempty(EXISTTIME)
        EXISTTIME = getsolvertime;
        NCHECKS = 1;
    else
        EXISTTIME = [EXISTTIME getsolvertime];
        NCHECKS = NCHECKS + 1;
    end
    if (options.cachesolvers==0)
        if ((NCHECKS >= 3 & (sum(EXISTTIME)/NCHECKS > 1)) | EXISTTIME(end)>2)
            if warningon
                info = 'Warning: YALMIP has detected that your drive or network is unusually slow.\nThis causes a severe delay in OPTIMIZE when I try to find available solvers.\nTo avoid this, use the options CACHESOLVERS in SDPSETTINGS.\nSee the FAQ for more information.\n';
                fprintf(info);
            end
        end
    end
    if length(EXISTTIME) > 5
        EXISTTIME = EXISTTIME(end-4:end);
        NCHECKS = 5;
    end
    CACHED_SOLVERS = solvers;
else
    solvers = CACHED_SOLVERS;
end

%here we are not actually selecting a solver, just checking whether it supports complex numbers
solverindex = min(find(strcmpi(lower({solvers.tag}),lower(options.solver))));
if isempty(solverindex)
	complexsolver = 0;
else
	solver = solvers(solverindex);
	complexsolver = solver.complex;
end

% Dualize the problem?
if ~isempty(F)
    if options.dualize == -1
        sdp = find(is(F,'sdp'));
        if ~isempty(sdp)
            if all(is(F(sdp),'sdpcone'))
                options.dualize = 1;
            end
        end
    end
end
if options.dualize == 1
    [Fd,objd,aux1,aux2,aux3,complexInfo] = dualize(F,h,[],[],[],options,complexsolver);
    options.dualize = 0;
    diagnostic = solvesdp(Fd,-objd,options);
    if ~isempty(complexInfo)
        for i = 1:length(complexInfo.replaced)
            n = size(complexInfo.replaced{i},1);
            re = 2*double(complexInfo.new{i}(1:n,1:n));
            im = 2*double(complexInfo.new{i}(1:n,n+1:end));
            im=triu((im-im')/2)-(triu((im-im')/2))';
            assign(complexInfo.replaced{i},re + sqrt(-1)*im);
        end
    end
    return
end

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
% DID WE SELECT THE MOMENT SOLVER
% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
if isequal(options.solver,'moment')
    if ~isempty(logdetStruct)
        error('Cannot dualize problems with logaritmic objective')
    end
    options.solver = options.moment.solver;
    [diagnostic,x,momentdata] = solvemoment(F,h,options,options.moment.order);
    diagnostic.momentdata = momentdata;
    diagnostic.xoptimal = x;
    return
end

% ******************************************
% COMPILE IN GENERALIZED YALMIP FORMAT
% ******************************************
[interfacedata,recoverdata,solver,diagnostic,F,Fremoved,ForiginalQuadratics] = compileinterfacedata(F,[],logdetStruct,h,options,0,solving_parametric,solvers);

% ******************************************
% FAILURE?
% ******************************************
if ~isempty(diagnostic)
    diagnostic.yalmiptime = etime(clock,yalmiptime);
    return
end

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
% DID WE SELECT THE LMILAB SOLVER WITH A KYP
% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
if  strcmpi(solver.tag,'lmilab') & any(is(F,'kyp'))
    [diagnostic,failed] = calllmilabstructure(F,h,options);
    if ~failed % Did this problem pass (otherwise solve using unstructured call)
        diagnostic.yalmiptime = etime(clock,yalmiptime)-diagnostic.solvertime;
        return
    end
end

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
% DID WE SELECT THE KYPD SOLVER
% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
if strcmpi(solver.tag,'kypd')
    diagnostic = callkypd(F,h,options);
    diagnostic.yalmiptime = etime(clock,yalmiptime)-diagnostic.solvertime;
    return
end

% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
% DID WE SELECT THE STRUL SOLVER
% ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
if strfind(solver.tag,'STRUL')
    diagnostic = callstrul(F,h,options);
    diagnostic.yalmiptime = etime(clock,yalmiptime)-diagnostic.solvertime;
    return
end

%******************************************
% DID WE SELECT THE MPT solver (backwards comb)
%******************************************
actually_save_output = interfacedata.options.savesolveroutput;
if strcmpi(solver.tag,'mpt') | strcmpi(solver.tag,'mpt-2') | strcmpi(solver.tag,'mpt-3') | strcmpi(solver.tag,'mpcvx') | strcmpi(solver.tag,'mplcp') | strcmpi(solver.tag,'pop')
    interfacedata.options.savesolveroutput = 1;
    if isempty(interfacedata.parametric_variables)
        if (nargin < 4 | ~isa(varargin{4},'sdpvar'))
            error('You must specify parametric variables.')
        else
            varargin{4} = reshape(varargin{4},[],1);
            interfacedata.parametric_variables = [];
            for i = 1:length(varargin{4})
                  interfacedata.parametric_variables = [interfacedata.parametric_variables;find(ismember(recoverdata.used_variables,getvariables(varargin{4}(i))))];
            end
            if isempty(varargin{5})
                interfacedata.requested_variables = [];
            else
                interfacedata.requested_variables = [];
                for i = 1:length(varargin{5})
                    interfacedata.requested_variables = [interfacedata.requested_variables;find(ismember(recoverdata.used_variables,getvariables(varargin{5}(i))))];
                end
            end
        end
    end
end

% *************************************************************************
% Just return the YALMIP model. Used when solving multiple objectives
% *************************************************************************
if isfield(options,'pureexport')
    interfacedata.recoverdata = recoverdata;
    diagnostic = interfacedata;
    return
end

if strcmpi(solver.version,'geometric') || (strcmpi(solver.tag,'bnb') && strcmpi(solver.lower.version,'geometric'))
    % Actual linear user variables
    if options.assertgpnonnegativity
        check = find(interfacedata.variabletype==0);
        check = setdiff(check,interfacedata.aux_variables);
        check = setdiff(check,interfacedata.evalVariables);
        check = setdiff(check,interfacedata.extended_variables);
        [lb,ub] = find_lp_bounds(interfacedata.F_struc,interfacedata.K);
        if ~all(lb(check)>=0)
            % User appears to have explictly selected a GP solver
            userdirect = ~isempty(strfind(options.solver,'geometric')) || ~isempty(strfind(options.solver,'mosek')) || ~isempty(strfind(options.solver,'gpposy'));
            userindirect = strcmpi(solver.tag,'bnb') && strcmpi(solver.lower.version,'geometric');
            if userdirect || userindirect
                % There are missing non-negativity bounds
                output = createOutputStructure(zeros(length(interfacedata.c),1)+NaN,[],[],18,yalmiperror(18,''),[],[],nan);
                diagnostic.yalmiptime = etime(clock,yalmiptime);
                diagnostic.solvertime = output.solvertime;
                try
                    diagnostic.info = output.infostr;
                catch
                    diagnostic.info = yalmiperror(output.problem,solver.tag);
                end
                diagnostic.problem = output.problem;
                if options.dimacs
                    diagnostic.dimacs = dimacs;
                end
                return
            else
                % YALMIP selected solver and picked a GP solver. As this is
                % no GP, we call again, but this time explicitly tell
                % YALMIP that it isn't a GP
                options.thisisnotagp = 1;
                varargin{3} = options;
                diagnostic = solvesdp(varargin{:});
                return
            end
        end
    end
end

% *************************************************************************
% TRY TO SOLVE PROBLEM
% *************************************************************************
if options.debug
    eval(['output = ' solver.call '(interfacedata);']);
else
    try
        eval(['output = ' solver.call '(interfacedata);']);
    catch
        output = createOutputStructure(zeros(length(interfacedata.c),1)+NaN,[],[],9,yalmiperror(9,lasterr),[],[],nan);
    end
end

if options.dimacs
    try
        b = -interfacedata.c;
        c = interfacedata.F_struc(:,1);
        A = -interfacedata.F_struc(:,2:end)';
        x = output.Dual;
        y = output.Primal;
        % FIX this nonlinear crap (return variable type in
        % compileinterfacedata)
        if options.relax == 0 & any(full(sum(interfacedata.monomtable,2)~=0))
            if ~isempty(find(sum(interfacedata.monomtable | interfacedata.monomtable,2)>1))
                z=real(exp(interfacedata.monomtable*log(y+eps)));
                y = z;
            end
        end

        if isfield(output,'Slack')
            s = output.Slack;
        else
            s = [];
        end

        dimacs = computedimacs(b,c,A,x,y,s,interfacedata.K);
    catch
        dimacs = [nan nan nan nan nan nan];
    end
else
    dimacs = [nan nan nan nan nan nan];
end

% ********************************
% ORIGINAL COORDINATES
% ********************************
output.Primal = recoverdata.x_equ+recoverdata.H*output.Primal;

% ********************************
% OUTPUT
% ********************************
diagnostic.yalmipversion = yalmip('ver');
diagnostic.matlabversion = version;
diagnostic.yalmiptime = etime(clock,yalmiptime)-output.solvertime;
diagnostic.solvertime = output.solvertime;
try
    diagnostic.info = output.infostr;
catch
    diagnostic.info = yalmiperror(output.problem,solver.tag);
end
diagnostic.problem = output.problem;
if options.dimacs
    diagnostic.dimacs = dimacs;
end

% Some more info is saved internally
solution_internal = diagnostic;
solution_internal.variables = recoverdata.used_variables(:);
solution_internal.optvar = output.Primal;

if ~isempty(interfacedata.parametric_variables)
    diagnostic.mpsol = output.solveroutput;
    options.savesolveroutput = actually_save_output;
end;

if interfacedata.options.savesolveroutput
    diagnostic.solveroutput = output.solveroutput;
end
if interfacedata.options.savesolverinput
    diagnostic.solverinput = output.solverinput;
end
if interfacedata.options.saveyalmipmodel
    diagnostic.yalmipmodel = interfacedata;
end

if options.warning & warningon & isempty(strfind(diagnostic.info,'No problems detected'))
    disp(['Warning: ' output.infostr]);
end

if ismember(output.problem,options.beeponproblem)
    try
        beep; % does not exist on all ML versions
    catch
    end
end

% And we are done! Save the result
if ~isempty(output.Primal)
    if size(output.Primal,2)>1
        for j = 1:size(output.Primal,2)
            temp = solution_internal;
            temp.optvar = temp.optvar(:,j);
            yalmip('setsolution',temp,j);
        end
    else
        yalmip('setsolution',solution_internal);
    end
end
if interfacedata.options.saveduals & solver.dual
    if isempty(interfacedata.Fremoved) | (nnz(interfacedata.Q)>0)
        try
            setduals(F,output.Dual,interfacedata.K);
        catch
        end
    else
        try
            % Duals related to equality constraints/free variables
            % have to be recovered b-A*x-Ht == 0
            b = -interfacedata.oldc;
            A = -interfacedata.oldF_struc(1+interfacedata.oldK.f:end,2:end)';
            H = -interfacedata.oldF_struc(1:interfacedata.oldK.f,2:end)';
            x = output.Dual;
            b_equ = b-A*x;
            newdual =  H\b_equ;
            setduals(interfacedata.Fremoved + F,[newdual;output.Dual],interfacedata.oldK);
        catch
             % this is a new feature...
            disp('Dual recovery failed. Please report this issue.');
        end
    end
end
% Hack to recover original QCQP duals from gurobi
if strcmp(solver.tag,'GUROBI-GUROBI')
    if length(ForiginalQuadratics) > 0
        if isfield(output,'qcDual')
            if length(output.qcDual) == length(ForiginalQuadratics)
             %   Ktemp.l = length(output.qcDual);
             %   Ktemp.f = 0;
             %   Ktemp.q = 0;
             %   Ktemp.s = 0;
             %   Ktemp.r = 0;
                Ftemp = F + ForiginalQuadratics;
                K = interfacedata.K;
                Ktemp = K;
                Ktemp.l = Ktemp.l + length(ForiginalQuadratics);
                tempdual = output.Dual;
                tempdual = [tempdual(1:K.f + K.l);-output.qcDual;tempdual(1+K.f+K.l:end)];
                setduals(Ftemp,tempdual,Ktemp);
%                setduals(ForiginalQuadratics,-output.qcDual,Ktemp);
            end
        end
    end
end

function yesno = warningon

s = warning;
yesno = isequal(s,'on');


text:
function [ ] = ml_plot_gmm_pdf(X, Priors, Mu, Sigma )
%ML_PLOT_GMM_PDF Summary of this function goes here
%   Detailed explanation goes here

% Auxiliary Variables
K = length(Priors);

% Contour Function for GMM
xplot = linspace(min(X(1,:)), max(X(1,:)), 100)';
yplot = linspace(min(X(2,:)), max(X(2,:)), 100)';
[X_, Y_] = meshgrid(xplot, yplot);
vals = zeros(size(X_));

for i = 1:length(xplot)
    for j = 1:size(yplot)
        x = [xplot(i);yplot(j)];
        vals(j,i) =  ml_gmm_pdf(x, Priors, Mu, Sigma);
    end
end

% Plot the GMM PDF Contours
figure;
% contour(X_,Y_, vals, 50, 'LineStyle', 'none');hold on;
contour(X_,Y_, vals, 100);hold on;
% surf(X_,Y_, vals, 100);hold on;
colormap jet
colorbar
colors = hsv(K);
ml_plot_centroid(Mu',colors);hold on;
options.class_names = {};
options.labels      = [];
options.points_size = 30;
options.colors      = [0 1 0.5];
options.plot_figure = true;
ml_plot_data(X',options);
title ('GMM-PDF Contours');

end

author: epfl-lasa


text:
% TRUNCATE_FILTER Truncates the Fourier transform of a filter
%
% Usage
%    filter = TRUNCATE_FILTER(filter_f, threshold, lowpass)
%
% Input
%    filter_f (numeric): The Fourier representation of the filter.
%
% Output
%    filter (struct): The truncated representation of the filter. See descrip-
%       tion for more details.
%
% Description
%    By extracting and storing only the Fourier transform coefficients whose
%    are above a certain threshold relative to the maximum value, storage
%    requirements for the filters are lessened and computation is sped up
%    since only non-zero coefficients are multiplied during convolution.
%
%    The support of the Fourier transform is defined so that all coefficients
%    with a magnitude above threshold*fmax are kept, where fmax is the maxi-
%    mum magnitude, and so that length(filter_f) divided by the size of the
%    support is a power of 2.
%
%    The output filter contains the fields:
%       filter.type (char): Fixed to 'fourier_truncated'.
%       filter.N (int): The original size of the filter.
%       filter.recenter (boolean): Indicates whether the Fourier transform
%          should be recentered after convolution. This is always true.
%       filter.start (int): The frequency index where the support of the
%          Fourier transform starts.
%       filter.coefft (numeric): The values of the Fourier coefficients on the
%          support.
%
% See also
%    OPTIMIZE_FILTER, PERIODIZE_FILTER

function filter = truncate_filter(filter_f, threshold, lowpass)
	N = length(filter_f);

	filter.type = 'fourier_truncated';
	filter.N = N;

	% Could have filter.recenter = lowpass, but since we don't know if we're
	% taking the modulus or not, we always need to recenter.
	filter.recenter = 1;

	[temp,ind_max] = max(filter_f);
	filter_f = circshift(filter_f,N/2-ind_max);
	ind1 = find(abs(filter_f)>(max(abs(filter_f))*threshold),1);
	ind2 = find(abs(filter_f)>(max(abs(filter_f))*threshold),1,'last');

	len = ind2-ind1+1;
	len = filter.N/2^(floor(log2(filter.N/len)));

	ind1 = round(round((ind1+ind2)/2)-len/2);
	ind2 = ind1+len-1;

	filter_f = filter_f(mod([ind1:ind2]-1,filter.N)+1);

	filter.coefft = filter_f;
	filter.start = ind1-(N/2-ind_max);
end

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function varargout=polygonwind(varargin)
%
% For calling details please see v_polygonwind.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_polygonwind.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: polygonwind.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_polygonwind(varargin{:});
else
    v_polygonwind(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function uniform_run ( )

%*****************************************************************************80
%
%% UNIFORM_RUN varies the parameter DELTA in a log uniform way.
%
%  Discussion:
%
%    Our base value for DELTA_BASE is 0.01.
%
%    Our quantity of interest Q is the time at which the solution achieves
%    the value 0.99.
%
%    Our parameter U is uniformly distributed in [-1,+1], and determines
%    our actual value of DELTA by
%
%      DELTA = 2^U * DELTA_BASE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
  delta_base = 0.01;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'UNIFORM_RUN\n' );
  fprintf ( 1, '  Vary the value of delta using a log-uniformly\n' );
  fprintf ( 1, '  distributed factor between 1/2 and 2.\n' );
  fprintf ( 1, '  Our quantity of interest Q is the time at which the\n' );
  fprintf ( 1, '  solution reaches 0.99.\n' );
  fprintf ( 1, '  Plot the solution curves, and plot Delta versus Q\n' );
  fprintf ( 1, '\n' );
%
%  Set up an option to have ODE45 return when the solution reached 0.99.
%
  options = odeset ( 'EVENTS', @event_function );

  figure ( 1 )
  clf
  hold on

  for trial = 1 : 10
%
%  Uniformly select
%    U, the logarithm base 2, between -1 and +1, of
%    F, a factor between 1/2 and 2, which multiplies
%    DELTA_BASE, the base value, which gives us
%    DELTA, the initial size of the flame.
%
    u = 2 * rand ( 1, 1 ) - 1.0;
    f = 2^u;
    delta = f * delta_base;
    d_plot(trial) = delta;
    fprintf ( 1, '  U = %g, factor = %g, DELTA = %g\n', u, f, delta );
%
%  Set the starting point.
%
    t_start = 0.0;
    y_start(1,1) = delta;
%
%  Get the stopping point.
%
    t_stop = 250.0;
%
%  Call ODE45 to solve the problem, stopping immediately if the event
%  (y=0.99) is observed.
%
    t_span = [ t_start, t_stop ];

    [ t, y ] = ode45 ( @flame_fun, t_span, y_start, options );
    fprintf ( '  Y(T) = 0.99 at T = %g\n', t(end) );
    q_plot(trial) = t(end);

    c = [ trial, 0, 11 - trial ] / 10;

    plot ( t, y(:,1), 'b-', 'Linewidth', 3, 'Color', c );

  end

  title ( 'Shampine Flame, Multiple values of DELTA', 'Fontsize', 24 );
  grid on
  xlabel ( '<--- T --->', 'Fontsize', 24 );
  ylabel ( '<--- X(T) --->', 'Fontsize', 24 );

  hold off

  filename = 'uniform_run.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Multiple solutions plotted in file "%s".\n', filename );
%
%  Since our deltas were computed randomly, they aren't in order.
%  We need to sort them, and their corresponding function values,
%  before plotting.
%
  [ d_plot, i ] = sort ( d_plot );
  q_plot = q_plot(i);
%
%  Plot the observed values of delta versus F, the time at
%  which the solution reached 0.99.
%
  figure ( 2 )
  plot ( d_plot, q_plot, 'r.-', 'Linewidth', 3, 'Markersize', 25 );
  grid on
  xlabel ( '<-- Parameter Delta -->', 'Fontsize', 24 )
  ylabel ( '<-- Ignition time Q(Delta) -->', 'Fontsize', 24 )
  title ( 'Ignition time as a function of parameter Delta', 'Fontsize', 24 )

  filename = 'uniform_qoi.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Quantity of Interest plotted in file "%s".\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'UNIFORM_RUN\n' );
  fprintf ( 1, '  Normal end of execution\n' );

  return
end


text:
function encoder_save(encoder, filePath)
if isfield(encoder, 'net')
  encoder.net = vl_simplenn_move(encoder.net, 'cpu') ;
  encoder.net.useGpu = false ;
end
save(filePath, '-struct', 'encoder') ;

author: mcimpoi

repo: deep-fbanks

save_path: github-repos/MATLAB/mcimpoi-deep-fbanks


text:
function [dynpcm2] = inH2O2dynpcm2(inH2O)
% Convert pressure from inches of water column at 4 degrees to dyne/cm2
% Chad Greene 2012
dynpcm2 = inH2O*2490.89;