text:
function CA184rule()
n=1000;
X=rand(1,n)>0.9;
X(n+1)=X(1);
Timestep=500;
data=zeros(Timestep,n);
Y=zeros(1,n);
for j=1:Timestep
    for i=1:n
        if i~=1
            if X(i-1)==1&&X(i)==0
                Y(i)=1;
           % elseif X(i)==1&&X(i+1)==0
              %  Y(i)=0;
            elseif X(i)==1&&X(i+1)==1
                Y(i)=1;
            else Y(i)=0;
            end
        else
            if X(n)==1&&X(i)==0
                Y(i)=1;
            %elseif X(i)==1&&X(i+1)==0
               % Y(i)=0;
            elseif X(i)==1&&X(i+1)==1
                Y(i)=1;
            else Y(i)=0;
            end
        end
    end
   X=Y;
   X(n+1)=Y(1);
   data(j,:)=Y;
end
imshow(data);
axis on;

author: Eurus-Holmes


text:
function [Btuph] = hp2Btuph(hp)
% Convert power from mechanical horsepower to British
% thermal units per hour.
% Chad A. Greene 2012
Btuph = hp*2544.433577644;

author: Sable

repo: mcbench-benchmarks


text:
function X = arg_setdirect(X,value)
% Recursively set the arg_direct flag in a data structure to the given value.
% [Result] = arg_setdirect(Data,Value)
%
% This function treats occurrences of the string arg_direct in cell arrays as a part of a
% name-value pair to update, and fields named arg_direct in structs or struct arrays as
% fields to update. This function does not append the flag if none is currently present.
%
% In:
%	Data : the data structure to update
%
%	Value : value of the arg_direct flag (true or false)
%
% Out:
%	Result : data structure with arg_direct flag updated;
%
% Notes:
%	This function is not yet particularly fast for large data structures.
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2012-02-20

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

if nargin < 2
    error('please supply a value.'); end

if isempty(X), return; end

if iscell(X)
    for k=1:numel(X)
         X{k} = arg_setdirect(X{k},value);
        if ischar(X{k}) && strcmpi(X{k},'arg_direct')
            X{k+1} = value; end
    end
elseif isstruct(X)
    if length(X) > 1
        for k=1:numel(X)
            X(k) = arg_setdirect(X(k),value); end
    else
        for fn=fieldnames(X)'
            X.(fn{1}) = arg_setdirect(X.(fn{1}),value);
            if strcmp(fn{1},'arg_direct')
                X.(fn{1}) = value; end
        end
    end
end

repo: REST


text:
function [P_music, est_dirs] = sphMUSIC(sphCOV, grid_dirs, nSrc)
%SPHMUSIC DoA estimation using MUSIC in the SHD
%
%   This routine computes the MUSIC pseudo-spectrum directly in the SHD.
%   Subspace approaches such as MUSIC can offer higher spatial resolution
%   than beamforming approaches such as the steered-response power, as long
%   as the source signals are not correlated between them and with the
%   reverberant/diffuse sound.
%
%   Inputs:
%       sphCOV: (order+1)^2x(order+1)^2 covariance/correlation matrix of
%           SH signals
%       grid_dirs:  Kx2 directions of [azi elev] in rads that define a grid
%           for the power map. For easy plotting of the directional maps,
%           use grid2dirs.m to generate the directions
%       nSrc:   (optional) number of peaks to try to find, as potential DoA
%           estimates (Von-Mises peak-finding contributed by Dr. Sakari Tervo)
%
%   Outputs:
%       P_pwd:      Kx1 vector of output powers, evaluated at grid directions
%       est_dirs:   nSrcx2 [azi elev] of estimated directions from
%           peak-finding
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SPHMUSIC.M - 5/10/2016
% Archontis Politis, archontis.politis@aalto.fi
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nSH = size(sphCOV,1);
order = sqrt(nSH)-1;
% sorted eigenvalue decomposition
V = sorted_eig(sphCOV, 'descend');
% noise subspace
Vn = V(:,nSrc+1:end);

% get steering vectors for grid points
nGrid = size(grid_dirs,1);
grid_xyz = unitSph2cart(grid_dirs);
grid_dirs2 = [grid_dirs(:,1) pi/2-grid_dirs(:,2)]; % go from azi-elevation to azi-inclination
Y_grid = getSH(order, grid_dirs2, 'real');

% get MUSIC spectrum
P_music = zeros(nGrid,1);
for ng = 1:nGrid
    stVec = Y_grid(ng,:).';
    P_music(ng) = 1 / (stVec' * Vn * Vn' * stVec);
end

% peak finding, if asked
if nargout==2 && nargin==3

    kappa  = 50; % Von-Mises concentration factor
    P_minus_peak = P_music;
    est_dirs = zeros(nSrc, 2);
    for k = 1:nSrc
        [~, peak_idx] = max(P_minus_peak);
        est_dirs(k,:) = grid_dirs(peak_idx,:);
        VM_mean = grid_xyz(peak_idx,:); % orientation of VM distribution
        VM_mask = kappa/(2*pi*exp(kappa)-exp(-kappa)) * exp(kappa*grid_xyz*VM_mean'); % VM distribution
        VM_mask = 1./(0.00001+VM_mask); % inverse VM distribution
        P_minus_peak = P_minus_peak.*VM_mask;
    end
end

end

author: polarch

repo: Spherical-Array-Processing


text:
function x=SeamPut(x,SeamVector)
% SEAMPUT takes as input a RGB or grayscale image and SeamVector array to
% find the pixels contained in the seam, and to adds a seam of interpolated
% pixels to the image. Each column of SeamVector must be a single seam.
%
% Author: Danny Luong
%         http://danluong.com
%
% Last updated: 12/20/07

[rows cols dim]=size(x);
[SVrows SVcols SVdim]=size(SeamVector);

if rows~=SVrows
    error('SeamVector and image dimension mismatch');
end

for k=1:SVcols
    for i=1:dim
        for j=1:rows
            if SeamVector(j,k)==1
                PutImg(j,:,i)=[x(j,1,i) x(j,1:cols,i)];
            elseif SeamVector(j,k)==cols
                PutImg(j,:,i)=[x(j,1:cols,i) x(j,cols,i)];
            else
                PutImg(j,:,i)=[x(j,1:SeamVector(j,k),i) 1/2*(x(j,SeamVector(j,k),i)+x(j,SeamVector(j,k)+1,i)) x(j,SeamVector(j,k)+1:cols,i)];
            end
        end
    end
    x=PutImg;
    clear PutImg
    [rows cols dim]=size(x);
end


text:
%% Creating one filter
% In this demo, we will show how to build some filters in spatial domain or
% Fourier domain. Gabor, Morlet and Gaussian filters are pointed out.
%
%% Reviewed functions
% The following function will be demonstrated:
%%
%
% # gabor_2d
% # gaussian_2d
% # morlet_2d_noDC
% # morlet_2d_pyramid
%
%% Demonstration
%
% Here are shown 4 examples of filters:

N=50;
M=50;
sigma=10;
slant=3;
theta=pi/3;
xi=2;

gab = gabor_2d(N, M, sigma, slant, xi, theta);
morlet = morlet_2d_noDC(N, M, sigma, slant, xi, theta);
morlet2 = morlet_2d_pyramid(N, M, sigma, slant, xi, theta);
gaussian = gaussian_2d(N, M, sigma);
subplot(2,2,1)
imagesc(real(gab));
axis off
title('Real part of Gabor wavelet in Fourier domain')
subplot(2,2,2)
imagesc(real(morlet));
axis off
title('Real part of Morlet wavelet in Fourier domain')
subplot(2,2,3)
imagesc(real(morlet2));
axis off
title('Real part of Morlet wavelet in Spatial domain')
subplot(2,2,4)
imagesc(gaussian);
axis off
title('Gaussian wavelet')

author: scatnet

repo: scatnet

save_path: github-repos/MATLAB/scatnet-scatnet


text:
function fem2d_bvp_linear_test01 ( )

%*****************************************************************************80
%
%% FEM2D_BVP_LINEAR_TEST01 carries out test case #1.
%
%  Discussion:
%
%    Use A1, C1, F1, EXACT1, EXACT_UX1, EXACT_UY1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 June 2014
%
%  Author:
%
%    John Burkardt
%
  nx = 3;
  ny = 3;

% nx = 5;
% ny = 5;

% nx = 9;
% ny = 9;

% nx = 17;
% ny = 17;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'FEM2D_BVP_LINEAR_TEST01\n' );
  fprintf ( 1, '  Solve - del ( A del U ) + C U = F \n' );
  fprintf ( 1, '  on the unit square with zero boundary conditions.\n' );
  fprintf ( 1, '  A1(X,Y) = 1.0\n' );
  fprintf ( 1, '  C1(X,Y) = 0.0\n' );
  fprintf ( 1, '  F1(X,Y) = 2*X*(1-X)+2*Y*(1-Y).\n' );
  fprintf ( 1, '  U1(X,Y) = X * ( 1 - X ) * Y * ( 1 - Y )\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The grid uses %d by %d nodes.\n', nx, ny );
%
%  Geometry definitions.
%
  x = linspace ( 0.0, 1.0, nx );
  y = linspace ( 0.0, 1.0, ny );

  u = fem2d_bvp_linear ( nx, ny, @a1, @c1, @f1, x, y );

  if ( 0 )
    [ X, Y ] = meshgrid ( x, y );
    surf ( X, Y, u )
  end

  if ( nx * ny <= 25 )
    fprintf ( 1, '\n' );
    fprintf ( 1, '     I     J    X         Y         U         Uexact    Error\n' );
    fprintf ( 1, '\n' );

    for j = 1 : ny
      for i = 1 : nx
        uexact = exact1 ( x(i), y(j) );
        fprintf ( 1, '  %4d  %4d  %8f  %8f  %8f  %8f  %8e\n', ...
          i, j, x(i), y(j), u(i,j), uexact, abs ( u(i,j) - uexact ) );
      end
    end

  end

  e1 = fem2d_l1_error ( nx, ny, x, y, u, @exact1 );
  e2 = fem2d_l2_error_linear ( nx, ny, x, y, u, @exact1 );
  h1s = fem2d_h1s_error_linear ( nx, ny, x, y, u, @exact_ux1, @exact_uy1 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  l1 error   = %g\n', e1 );
  fprintf ( 1, '  L2 error   = %g\n', e2 );
  fprintf ( 1, '  H1S error  = %g\n', h1s );

  return
end
function value = a1 ( x, y )

%*****************************************************************************80
%
%% A1 evaluates A function #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the evaluation point.
%
%    Output, real VALUE, the value of A(X).
%
  value = 1.0;

  return
end
function value = c1 ( x, y )

%*****************************************************************************80
%
%% C1 evaluates C function #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the evaluation point.
%
%    Output, real VALUE, the value of C(X).
%
  value = 0.0;

  return
end
function value = exact1 ( x, y )

%*****************************************************************************80
%
%% EXACT1 evaluates exact solution #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the evaluation point.
%
%    Output, real VALUE, the value of the solution.
%
  value = x .* ( 1.0 - x ) .* y .* ( 1.0 - y );

  return
end
function value = exact_ux1 ( x, y )

%*****************************************************************************80
%
%% EXACT_UX1 evaluates the derivative dUdX of exact solution #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the evaluation point.
%
%    Output, real VALUE, the value of dUdX.
%
  value = ( 1.0 - 2.0 * x ) .* ( y - y .* y );

  return
end
function value = exact_uy1 ( x, y )

%*****************************************************************************80
%
%% EXACT_UY1 evaluates the derivative dUdY of exact solution #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the evaluation point.
%
%    Output, real VALUE, the value of dUdX.
%
  value = ( x - x .* x ) .* ( 1.0 - 2.0 * y );

  return
end
function value = f1 ( x, y )

%*****************************************************************************80
%
%% F1 evaluates right hand side function #1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 June 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X, Y, the evaluation point.
%
%    Output, real VALUE, the value of the right hand side.
%
  value = 2.0 * x .* ( 1.0 - x ) ...
        + 2.0 * y .* ( 1.0 - y );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [handleImageArray, nameImageArray] = get_all_image_objects(this, type)
% Returns cell of handles & names for all/selected image objects in MrSeries
%
%   Y = MrSeries()
%  [handleImageArray, nameImageArray] = Y.get_all_image_objects(inputs)
%
% This is a method of class MrSeries.
%
% IN
%
% OUT
%   handleImageArray    cell(nImages, 1) of handles to MrImages
%   nameImageArray      cell(nImages, 1) of names of properties of MrSeries
%                       which are MrImages
%
% EXAMPLE
%   get_all_image_objects
%
%   See also MrSeries

% Author:   Saskia Klein & Lars Kasper
% Created:  2014-07-09
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.

%handleImageArray = findall(this, 'MrImage');

if nargin < 2
    type = 'all';
end

nameStatImageArray = {'mean'; 'sd'; 'snr'; 'coeffVar'; 'diffLastFirst'; 'diffOddEven'};
nameDataImageArray = {'data'};
nameAnatomy = {'anatomy'};

switch lower(type)
    case 'data'
        nameImageArray = nameDataImageArray;
    case {'stat', 'stats', 'statistical'}
        nameImageArray = nameStatImageArray;
    case 'anatomy'
        nameImageArray = nameAnatomy;
    case 'all'
        nameImageArray = [nameDataImageArray; nameStatImageArray; ...
            nameAnatomy];
end

nImages = numel(nameImageArray);

for iImage = 1:nImages
    handleImageArray{iImage} = this.(nameImageArray{iImage});
end

author: translationalneuromodeling

repo: tapas


text:
classdef ImageNet < datasets.ImageFolder
%IMAGENET ImageNet dataset
%   Encapsulates the ImageNet dataset for training.
%
%   D = datasets.ImageNet('dataDir', '/data/imagenet') loads the ImageNet
%   dataset from the directory '/data/imagenet'.
%   Detailed instructions on how to set it up are below.
%
%   D.train() returns a cell array with mini-batches, from the shuffled
%   training set. Each mini-batch consists of a set of indexes.
%
%   D.val() returns a cell array with mini-batches, from the validation set
%   (without shuffling). Each mini-batch consists of a set of indexes.
%
%   [IMAGES, LABELS] = D.get(BATCH) returns a tensor of images and the
%   corresponding labels for the given mini-batch BATCH.
%
%   datasets.ImageNet(...,'option', value, ...) sets several properties,
%   which are all inherited from the datasets.ImageFolder class.
%   See 'help datasets.ImageFolder' for more information.
%
%   See 'autonn/examples/cnn/imagenet_example.m' for a full example.
%
%
%   * ImageNet data setup:
%
%   The ImageNet/ILSVRC data ships in several TAR archives that can be
%   obtained from the ILSVRC challenge website. You will need:
%
%   ILSVRC2012_img_train.tar
%   ILSVRC2012_img_val.tar
%   ILSVRC2012_img_test.tar
%   ILSVRC2012_devkit.tar
%
%   Note that images in the CLS-LOC challenge are the same for the 2012,
%   2013, and 2014 edition of ILSVRC, but that the development kit is
%   different. However, all devkit versions should work.
%
%   Create the following hierarchy:
%
%   <dataDir>/images/train/ : content of ILSVRC2012_img_train.tar
%   <dataDir>/images/val/ : content of ILSVRC2012_img_val.tar
%   <dataDir>/images/test/ : content of ILSVRC2012_img_test.tar
%   <dataDir>/ILSVRC2012_devkit : content of ILSVRC2012_devkit.tar
%
%   In order to speedup training and testing, it is a good idea to
%   preprocess the images to have a fixed size (e.g. 256 pixels high).
%   There is a script in MatConvNet (utils/preprocess-imagenet.sh) that
%   achieves this. It may also be necessary to store the images in RAM
%   disk. Reading images off disk with a sufficient speed is crucial for
%   fast training.

% Copyright (C) 2018 Joao F. Henriques, Andrea Vedaldi.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

  properties
    labels
    sets
    classes
  end

  methods
    function o = ImageNet(varargin)
      % call parent class's constructor, but delay initialization
      o = o@datasets.ImageFolder('skipInitialization', true) ;

      % parse generic ImageFolder arguments
      varargin = o.parseGenericArgs(varargin) ;
      assert(isempty(varargin), 'Unknown arguments.') ;

      assert(~isempty(o.dataDir), 'Must specify dataset.dataDir.') ;

      cache = [o.dataDir '/imagenet.mat'] ;
      if exist(cache, 'file')
        % load from cache file
        s = load(cache) ;
        o.filenames = s.filenames ;
        o.labels = s.labels ;
        o.sets = s.sets ;
        o.classes = s.classes ;
      else
        % use matconvnet's example imagenet imdb as-is (see below)
        imdb = cnn_imagenet_setup_data('dataDir', o.dataDir) ;

        % read info from imdb
        o.filenames = strcat(['images' filesep], imdb.images.name) ;
        o.labels = imdb.images.label ;
        o.sets = imdb.images.set ;
        o.classes = imdb.classes ;

        % save to cache file
        [filenames, labels, sets, classes] = deal(o.filenames, o.labels, o.sets, o.classes) ;  %#ok<ASGLU>
        save(cache, 'filenames', 'labels', 'sets', 'classes', '-v6') ;  % no compression, < 2 GB
      end

      o.trainSet = find(o.sets == 1) ;
      o.valSet = find(o.sets == 2) ;

      o.initialize() ;
    end

    function [images, labels] = get(o, batch)
      % return a single batch of images, with optional prefetching
      [images, idx] = o.get@datasets.ImageFolder(batch) ;
      labels = o.labels(idx) ;
    end
  end
end

function imdb = cnn_imagenet_setup_data(varargin)
% CNN_IMAGENET_SETUP_DATA  Initialize ImageNet ILSVRC CLS-LOC challenge data
%    This function creates an IMDB structure pointing to a local copy
%    of the ILSVRC CLS-LOC data.

opts.dataDir = fullfile('data','imagenet12') ;
opts = vl_argparse(opts, varargin) ;

% -------------------------------------------------------------------------
%                                                  Load categories metadata
% -------------------------------------------------------------------------

d = dir(fullfile(opts.dataDir, '*devkit*')) ;
if numel(d) == 0
  error('Make sure that ILSVRC data is correctly installed in %s', ...
    opts.dataDir) ;
end
devkitPath = fullfile(opts.dataDir, d(1).name) ;

% find metadata
mt = dir(fullfile(devkitPath, 'data', 'meta_clsloc.mat')) ;
if numel(mt) == 0
  mt = dir(fullfile(devkitPath, 'data', 'meta.mat')) ;
end
metaPath = fullfile(devkitPath, 'data', mt(1).name) ;

% find validation images labels
tmp = dir(fullfile(devkitPath, 'data', '*_validation_ground_truth*')) ;
valLabelsPath = fullfile(devkitPath, 'data', tmp(1).name) ;

% find validation images blacklist
tmp = dir(fullfile(devkitPath, 'data', '*_validation_blacklist*')) ;
if numel(tmp) > 0
  valBlacklistPath = fullfile(devkitPath, 'data', tmp(1).name) ;
else
  valBlacklistPath = [] ;
  warning('Could not find validation images blacklist file');
end

fprintf('using devkit %s\n', devkitPath) ;
fprintf('using metadata %s\n', metaPath) ;
fprintf('using validation labels %s\n', valLabelsPath) ;
fprintf('using validation blacklist %s\n', valBlacklistPath) ;

meta = load(metaPath) ;
cats = {meta.synsets(1:1000).WNID} ;
descrs = {meta.synsets(1:1000).words} ;

imdb.classes.name = cats ;
imdb.classes.description = descrs ;
imdb.imageDir = fullfile(opts.dataDir, 'images') ;

% -------------------------------------------------------------------------
%                                                           Training images
% -------------------------------------------------------------------------

fprintf('searching training images ...\n') ;
names = {} ;
labels = {} ;
for d = dir(fullfile(opts.dataDir, 'images', 'train', 'n*'))'
  [~,lab] = ismember(d.name, cats) ;
  ims = dir(fullfile(opts.dataDir, 'images', 'train', d.name, '*.JPEG')) ;
  names{end+1} = strcat([d.name, filesep], {ims.name}) ;
  labels{end+1} = ones(1, numel(ims)) * lab ;
  fprintf('.') ;
  if mod(numel(names), 50) == 0, fprintf('\n') ; end
  %fprintf('found %s with %d images\n', d.name, numel(ims)) ;
end
names = horzcat(names{:}) ;
labels = horzcat(labels{:}) ;

if numel(names) ~= 1281167
  warning('Found %d training images instead of 1,281,167. Dropping training set.', numel(names)) ;
  names = {} ;
  labels =[] ;
end

names = strcat(['train' filesep], names) ;

imdb.images.id = 1:numel(names) ;
imdb.images.name = names ;
imdb.images.set = ones(1, numel(names)) ;
imdb.images.label = labels ;

% -------------------------------------------------------------------------
%                                                         Validation images
% -------------------------------------------------------------------------

ims = dir(fullfile(opts.dataDir, 'images', 'val', '*.JPEG')) ;
names = sort({ims.name}) ;
labels = textread(valLabelsPath, '%d') ;

if numel(ims) ~= 50e3
  warning('Found %d instead of 50,000 validation images. Dropping validation set.', numel(ims))
  names = {} ;
  labels =[] ;
else
  if ~isempty(valBlacklistPath)
    black = textread(valBlacklistPath, '%d') ;
    fprintf('blacklisting %d validation images\n', numel(black)) ;
    keep = setdiff(1:numel(names), black) ;
    names = names(keep) ;
    labels = labels(keep) ;
  end
end

names = strcat(['val' filesep], names) ;

imdb.images.id = horzcat(imdb.images.id, (1:numel(names)) + 1e7 - 1) ;
imdb.images.name = horzcat(imdb.images.name, names) ;
imdb.images.set = horzcat(imdb.images.set, 2*ones(1,numel(names))) ;
imdb.images.label = horzcat(imdb.images.label, labels') ;

% -------------------------------------------------------------------------
%                                                               Test images
% -------------------------------------------------------------------------

ims = dir(fullfile(opts.dataDir, 'images', 'test', '*.JPEG')) ;
names = sort({ims.name}) ;
labels = zeros(1, numel(names)) ;

if numel(labels) ~= 100e3
  warning('Found %d instead of 100,000 test images', numel(labels))
end

names = strcat(['test' filesep], names) ;

imdb.images.id = horzcat(imdb.images.id, (1:numel(names)) + 2e7 - 1) ;
imdb.images.name = horzcat(imdb.images.name, names) ;
imdb.images.set = horzcat(imdb.images.set, 3*ones(1,numel(names))) ;
imdb.images.label = horzcat(imdb.images.label, labels) ;

% -------------------------------------------------------------------------
%                                                            Postprocessing
% -------------------------------------------------------------------------

% sort categories by WNID (to be compatible with other implementations)
[imdb.classes.name,perm] = sort(imdb.classes.name) ;
imdb.classes.description = imdb.classes.description(perm) ;
relabel(perm) = 1:numel(imdb.classes.name) ;
ok = imdb.images.label >  0 ;
imdb.images.label(ok) = relabel(imdb.images.label(ok)) ;

end

repo: MFT


text:
classdef DWU < ALGORITHM
% <multi> <real/integer/label/binary/permutation>
% Dominance-weighted uniformity multi-objective evolutionary algorithm

%------------------------------- Reference --------------------------------
% G. Moreira and L. Paquete, Guiding under uniformity measure in the
% decision space, Proceedings of the IEEE Latin American Conference on
% Computational Intelligence, 2019.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Gladston Moreira

    methods
        function main(Algorithm,Problem)
            %% Generate random population
            Population = Problem.Initialization();
            [Population,FrontNo] = EnvironmentalSelection(Population,Problem.N);

            %% Optimization
            while Algorithm.NotTerminated(Population)
                MatingPool = TournamentSelection(2,Problem.N,FrontNo);
                Offspring  = OperatorGA(Problem,Population(MatingPool));
                [Population,FrontNo,] = EnvironmentalSelection([Population,Offspring],Problem.N);
            end
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function[f]=fourier(varargin)
%FOURIER  Returns the Fourier frequencies for a given length time series.
%
%   F=FOURIER(M) returns the one-sided (or positive) Fourier frequencies
%   for a time series of length M.
%
%   F is a radian or angular frequency so that the Nyquist is at PI.
%
%   F=FOURIER(M,'two') instead returns the two-sided Fourier frequencies.
%   The default behavior is equivalent to F=FOURIER(M,'one').
%
%   F=FOURIER(DT,M) uses sample rate DT in calculating the frequencies, so
%   that the Nyquist will be at PI/DT.
%
%   Note that the highest resolved frequency, MAX(F), differs for even or
%   odd M.  For even M, it is the Nyquist PI/DT, but for odd M the Nyquist
%   is not resolved and the highest resolved frequency is (M-1)/M * PI/DT.
%
%   For the one-sided option, F has length FLOOR(M/2)+1, or M/2+1 for even
%   M, and (M+1)/2 for odd M.  For the two-sided option, F has length M.
%   __________________________________________________________________
%
%   Array input
%
%   F=FOURIER(M) also works if M an array instead of a scalar.  In this
%   case, F is a cell array with LENGTH(M) elements.
%
%   F=FOURIER(DT,M) with M being an array works provided DT is either a
%   scalar or an array of the same length as M.
%   __________________________________________________________________
%
%   See also CENTEREDTIMES.
%
%   Usage: f=fourier(M);
%          f=fourier(dt,M);
%          f=fourier(dt,M,'two');
%   __________________________________________________________________
%   This is part of JLAB --- type 'help jlab' for more information
%   (C) 2011--2019 J.M. Lilly --- type 'help jlab_license' for details

if strcmpi(varargin{1}, '--t')
    fourier_test,return
end

str='one';
if ischar(varargin{end})
    str=varargin{end}(1:3);
    varargin=varargin(1:end-1);
end

if length(varargin)==1
    N=varargin{1};
    dt=1;
else
    dt=varargin{1};
    N=varargin{2};
end

if length(N)==1
    f=fourier_one(dt,N,str);
else
    if length(dt)==1
        dt=dt+zeros(size(N));
    elseif length(dt)~=length(N)
        error('If N is an array, DT must either be a scalar or an array of the same length.')
    end
    for i=1:length(N)
        f{i,1}=fourier_one(dt,N(i),str);
    end
end

function[f]=fourier_one(dt,N,str)

if strcmpi(str,'one')
    f=2*pi*(0:floor(N/2))'./N;
else
    f=2*pi*(0:N-1)'./N;
    index=find(f>pi);
    f(index)=f(index)-2*pi;
end

% if iseven(N)
%     f=(0:1./N:1/2)';
% elseif isodd(N)
%     f=(0:1./N:1/2*frac(N-1,N))';
% end

f=f./dt;

function[]=fourier_test

reporttest('FOURIER is length N/2+1 for even N',length(fourier(10))==6);
reporttest('FOURIER is length (N+1)/2 for odd N',length(fourier(11))==6);

%Even and odd length frequencies behave differently, which you can see with
%abs(fft(cos(pi*[1:10]'))).^2
%abs(fft(cos(pi*[1:11]'))).^2
%abs(fft(cos((10/11)*pi*[1:11]'))).^2

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


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

BV = [];
BV3 = [];
mag = [];
me = [];
av2=[];
Nmin = 50;
def = {'10','5'};
tit ='b with depth input parameters';
prompt={ 'Depth range sampled in km',...
    'move window in km',...
    };

ni2 = inputdlg(prompt,tit,1,def);
l = ni2{1};
drkm = str2double(l);
l = ni2{2};
mokm = str2double(l);

ButtonName=questdlg('Mc determination?', ...
    ' Question', ...
    'Automatic','Fixed Mc=Mmin','Money');

think
[s,is] = sort(newt2.Depth);
newt1 = newt2(is(:,1),:) ;
watchon;

for t = min(newt2.Depth):mokm:max(newt2.Depth)
    % calculate b-value based an weighted LS
    l = newt2.Depth >= t & newt2.Depth < t+drkm;
    b = newt2(l,:);

    switch ButtonName
        case 'Automatic'
            % mcperc_ca3;
            %if isnan(Mc95) == 0
            %    magco = Mc95;
            %elseif isnan(Mc90) == 0
            %    magco = Mc90;
            %else
            [bv magco stan av me mer me2,  pr] =  bvalca3(b,1,1);
            %end
        case 'Fixed Mc=Mmin'
            magco = min(newt1(:,6))
    end

    l = b(:,6) >= magco-0.05;
    if length(b(l,:)) >= Nmin
        %[bv magco0 stan av me mer me2,  pr] =  bvalca3(b(l,:),2,2);
        [mea bv stan,  av] =  bmemag(b(l,:));
    else
        bv = nan; bv2 = nan, magco = nan; av = nan; av2 = nan;
    end
    BV = [BV ; bv min(b(:,7)) ; bv max(b(:,7)) ; inf inf];
    BV3 = [BV3 ; bv mean(b(:,7)) stan ];
    %mag = [mag ; av newt1(t+round(ni/2),7)];

    % calculate b-value based on maximum likelihood
    %av2 = [av2 ;   av  newt1(t+round(ni/2),7) stan bv];

end

watchoff

% Find out of figure already exists
%
[existFlag,figNumber]=figure_exists('b-value with depth',1);
newdepWindowFlag=~existFlag;
bdep= figNumber;

% Set up the Cumulative Number window

if newdepWindowFlag
    bdep = figure_w_normalized_uicontrolunits( ...
        'Name','b-value with depth',...
        'NumberTitle','off', ...
        'MenuBar','none', ...
        'NextPlot','add', ...
        'backingstore','on',...
        'Visible','on', ...
        'Position',[ 150 150 (ZmapGlobal.Data.map_len - [50 20]) ]);

    uicontrol('Units','normal',...
        'Position',[.0 .85 .08 .06],'String','Info ',...
         'Callback','infoz(1)');

    matdraw
end

figure_w_normalized_uicontrolunits(bdep)
delete(gca)
delete(gca)
delete(gca)
delete(gca)
hold off

axis off
hold on
orient tall
%rect = [0.15 0.65 0.7 0.25];
rect = [0.25 0.15 0.5 0.75];
axes('position',rect)
ple = errorbar(BV3(:,2),BV3(:,1),BV3(:,3),BV3(:,3),'k')
set(ple(1),'color',[0.5 0.5 0.5]);

hold on
pl = plot(BV(:,2),BV(:,1),'color',[0.5 0.5 0.5]);

pl = plot(BV3(:,2),BV3(:,1),'sk')

set(pl,'LineWidth',1.0,'MarkerSize',4,...
    'MarkerFaceColor','w','MarkerEdgeColor','k','Marker','s');

set(gca,'box','on',...
    'SortMethod','childorder','TickDir','out','FontWeight',...
    'bold','FontSize',ZmapGlobal.Data.fontsz.m,'Linewidth',1.,'Ticklength',[ 0.02 0.02])

bax = gca;
strib = [name ', ni = ' num2str(ni), ', Mmin = ' num2str(min(newt2.Magnitude)) ];
ylabel('b-value')
xlabel('Depth [km]')
title2(strib,'FontWeight','bold',...
    'FontSize',ZmapGlobal.Data.fontsz.m,...
    'Color','k')

xl = get(gca,'Xlim');
view([90 90])

author: CelsoReyes


text:
function amax_index = r8vec_amax_index ( n, a )

%*****************************************************************************80
%
%% R8VEC_AMAX_INDEX returns the index of the maximum absolute value in an R8VEC.
%
%  Discussion:
%
%    An R8VEC is a vector of R8 values.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the array.
%
%    Input, real A(N), the array.
%
%    Output, integer AMAX_INDEX, the index of the entry of largest magnitude.
%
  if ( n <= 0 )

    amax_index = -1;

  else

    amax_index = 1;
    amax = abs ( a(1) );

    for i = 2 : n
      if ( amax < abs ( a(i) ) )
        amax_index = i;
        amax = abs ( a(i) );
      end
    end

  end

  return
end


text:
function [ n_data, z, s, a, fx ] = lerch_values ( n_data )

%*****************************************************************************80
%
%% LERCH_VALUES returns some values of the Lerch transcendent function.
%
%  Discussion:
%
%    The Lerch function is defined as
%
%      Phi(z,s,a) = Sum ( 0 <= k < Infinity ) z^k / ( a + k )^s
%
%    omitting any terms with ( a + k ) = 0.
%
%    In Mathematica, the function can be evaluated by:
%
%      LerchPhi[z,s,a]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, real Z, the parameters of the function.
%
%    Output, integer S, the parameters of the function.
%
%    Output, real A, the parameters of the function.
%
%    Output, real FX, the value of the function.
%
  n_max = 12;

  a_vec = [ ...
     0.0, ...
     0.0, ...
     0.0, ...
     1.0, ...
     1.0, ...
     1.0, ...
     2.0, ...
     2.0, ...
     2.0, ...
     3.0, ...
     3.0, ...
     3.0 ];

  fx_vec = [ ...
     0.1644934066848226E+01, ...
     0.1202056903159594E+01, ...
     0.1000994575127818E+01, ...
     0.1164481052930025E+01, ...
     0.1074426387216080E+01, ...
     0.1000492641212014E+01, ...
     0.2959190697935714E+00, ...
     0.1394507503935608E+00, ...
     0.9823175058446061E-03, ...
     0.1177910993911311E+00, ...
     0.3868447922298962E-01, ...
     0.1703149614186634E-04 ];

  s_vec = [ ...
     2, 3, 10, ...
     2, 3, 10, ...
     2, 3, 10, ...
     2, 3, 10 ];

  z_vec = [ ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.1000000000000000E+01, ...
     0.5000000000000000E+00, ...
     0.5000000000000000E+00, ...
     0.5000000000000000E+00, ...
     0.3333333333333333E+00, ...
     0.3333333333333333E+00, ...
     0.3333333333333333E+00, ...
     0.1000000000000000E+00, ...
     0.1000000000000000E+00, ...
     0.1000000000000000E+00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    z = 0.0;
    s = 0;
    a = 0.0;
    fx = 0.0;
  else
    z = z_vec(n_data);
    s = s_vec(n_data);
    a = a_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
%% check Clebsch Gordan Tensor

%% the reference

% some arbitrary rotation
g = rotation.byEuler(-72*degree,-88*degree,-134*degree);

% we want to express the product of two wigner D functions
D2 = WignerD(g,'order',2);
D1D2_ref = D2(:) * D2(:).';

%% expansion into Wigner functions of lower order

% zero order component
D0 = WignerD(g,'order',0);
CG0 = ClebschGordanTensor(2,2,0);
C0 = EinsteinSum(CG0,[1 3 -1],D0,[-1,-2],CG0,[2 4 -2])

% first order component
D1 = WignerD(g,'order',1);
CG1 = ClebschGordanTensor(2,2,1);
C1 = EinsteinSum(CG1,[1 3 -1],D1,[-1 -2],CG1,[2 4 -2])

% second order component
D2 = WignerD(g,'order',2);
CG2 = ClebschGordanTensor(2,2,2);
C2 = EinsteinSum(CG2,[1 3 -1],D2,[-1 -2],CG2,[2 4 -2])

% third order component
D3 = WignerD(g,'order',3);
CG3 = ClebschGordanTensor(2,2,3);
C3 = EinsteinSum(CG3,[1 3 -1],D3,[-1 -2],CG3,[2 4 -2])

% fourth order component
D4 = WignerD(g,'order',4);
CG4 = ClebschGordanTensor(2,2,4);
C4 = EinsteinSum(CG4,[1 3 -1],D4,[-1 -2],CG4,[2 4 -2])

a = reshape(matrix(C0 + C1 + C2 + C3 + C4),[25,25]) ./ D1D2_ref;

imagesc(real(a))
mtexColorMap white2black

%% next we expand D1 * D1 * D1 * D1 in the same way

%% the reference

% some arbitrary rotation
g = rotation.byEuler(-72*degree,-88*degree,-134*degree);

% we want to express the product of four wigner D functions
D1 = WignerD(g,'order',1);

T2D1_ref = D1(:) * D1(:).';
T4D1_ref = T2D1_ref(:) * T2D1_ref(:).';

%% expansion into Wigner functions of lower order

%% zero order component

T4D1 = tensor(zeros(repmat(3,1,8)));
for J = 0:4

  DJ = WignerD(g,'order',J);

  CGJ = tensor(zeros([repmat(3,1,8),2*J+1,2*J+1]),'rank',10);
  for j1 = 0:2
    for j2 = 0:2
      CGj1 = ClebschGordanTensor(1,1,j1);
      CGj2 = ClebschGordanTensor(1,1,j2);
      CGj1j2J = ClebschGordanTensor(j1,j2,J);
      C = EinsteinSum(...
        CGj1,[1 3 -1],...
        CGj1,[2 4 -2],...
        CGj2,[5 7 -3],...
        CGj2,[6 8 -4],...
        CGj1j2J,[-1 -3 9],...
        CGj1j2J,[-2 -4 10]);
      CGJ = CGJ + C;
    end
  end
  T4D1 = T4D1 + EinsteinSum(CGJ,[1:8 -1 -2],DJ,[-1 -2])

end

%%

a = reshape(double(T4D1),[3*3*3*3,3*3*3*3]) ./ T4D1_ref;

imagesc(real(a))

author: mtex-toolbox

repo: mtex

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


text:
function Frames=calib_towncenter(Frames,namefile)

fid=fopen(namefile);
C=textscan(fid,'%s');
c=C{1};

fx=str2double(cell2mat(c(3)));
fy=str2double(cell2mat(c(6)));
px=str2double(cell2mat(c(9)));
py=str2double(cell2mat(c(12)));
sk=str2double(cell2mat(c(15)));
tx=str2double(cell2mat(c(18)));
ty=str2double(cell2mat(c(21)));
tz=str2double(cell2mat(c(24)));
rx=str2double(cell2mat(c(27)));
ry=str2double(cell2mat(c(30)));
rz=str2double(cell2mat(c(33)));
rw=str2double(cell2mat(c(36)));
k1=str2double(cell2mat(c(39)));
k2=str2double(cell2mat(c(42)));
p1=str2double(cell2mat(c(45)));
p2=str2double(cell2mat(c(48)));

R=[1-2*ry^2-2*rz^2,2*rx*ry-2*rz*rw,2*rx*rz+2*ry*rw;2*rx*ry+2*rz*rw,1-2*rx^2-2*rz^2,2*ry*rz-2*rx*rw;2*rx*rz-2*ry*rw,2*ry*rz+2*rx*rw,1-2*rx^2-2*ry^2];

R=[R,[tx;ty;tz]];

K=[fx,0,px;0,fy,py;0,0,1];

P=K*R;
P(:,3)=[];

R(:,3)=[];

for fr=1:numel(Frames)
    for p=1:numel(Frames(fr).id)

            u=double(Frames(fr).ximg(p));
            v=double(Frames(fr).yimg(p));

            x2=(u-px)/fx;
            y2=(v-py)/fy;

            Pu=undistort(k1,k2,p1,p2,[x2;y2]);
            xu=Pu(1);
            yu=Pu(2);

            Pd2=homotrans(inv(R),[Pu;1]);

           % Pd=homotrans(inv(P),[u;v;1]);

            Frames(fr).x(p)=Pd2(1);
            Frames(fr).y(p)=Pd2(2);
            Frames(fr).z(p)=0;
            Frames(fr).vx(p)=0;
            Frames(fr).vy(p)=0;
            Frames(fr).vz(p)=0;
            Frames(fr).ximg(p)=round(Frames(fr).ximg(p));
            Frames(fr).yimg(p)=round(Frames(fr).yimg(p));
    end
end

end

function x=undistort(k1,k2,p1,p2,xd)

    xd=double(xd);
    x = xd;                             % initial guess

    for kk=1:20,

        r_2 = sum(x.^2);
        k_radial =  double(1 + k1 * r_2 + k2 * r_2.^2 );
        delta_x = double([2*p1*x(1,:).*x(2,:) + p2*(r_2 + 2*x(1,:).^2);
        p1 * (r_2 + 2*x(2,:).^2)+2*p2*x(1,:).*x(2,:)]);
        x = (xd - delta_x)./(ones(2,1)*k_radial);

    end;

end

function t = homotrans(P,v)

[dim,npts] = size(v);

if ~all(size(P)==dim)
    error('Transformation matrix and point dimensions do not match');
end

t = P*v;  % Transform

for r = 1:dim-1     %  Now normalise
    t(r,:) = t(r,:)./t(end,:);
end

t(end,:) = ones(1,npts);

end

author: VisDrone

repo: DroneCrowd

save_path: github-repos/MATLAB/VisDrone-DroneCrowd


text:
function test_vbgppcamv(vbpca, ppca)

rand('state', 6);
randn('state', 6);

if nargin < 1
  vbpca = false;
end
if nargin < 2
  ppca = false;
end

inW = [120*rand(1,30); 80*rand(1,30)];
inX = 1:300;

M = length(inW);
N = length(inX);
D = 2;
X = zeros(D,N);
W = zeros(M,D);

% Covariance matrices
logthetaW = log(3e3);
logthetaX = log(10);

% Covariance functions
gpcovW = {@gpcovScale, @(logtheta,x1,x2) gpcov(logtheta,x1,x2,@sqdistEarth)};
logthetaW = [3; logthetaW(:)];
pseudoW = cell(D,1);%zeros([dimW,Mp,D]);
Mp = ceil(0.2*M);
for d=1:D
  permM = randperm(M);
  pseudoW{d} = inW(:,permM(1:Mp));% + randn(2,Mp);
end

gpcovX = @gpcov;

% Generate latent variables
for d=1:D
% $$$   W(:,d) = mvnrnd(zeros(1,M),CovW)';
% $$$   X(d,:) = mvnrnd(zeros(1,N),CovX);

  W(:,d) = gprnd(inW, logthetaW, gpcovW);

  lsX = 10*ceil(exp(logthetaX(end)));
  X(d,:) = filter(hamming(lsX), 1, randn(N,1));
% $$$   W(:,d) = filter(ones(lsW,1)/lsW, 1, randn(M,1));
% $$$   X(d,:) = filter(ones(lsX,1)/lsX, 1, randn(N,1));
  covfunW{d} = gpcovW;
  covfunX{d} = gpcovX;
  initthetaW{d} = logthetaW-1;
  initthetaX{d} = logthetaX-1;
end

W = W * diag( sqrt(M./rowsum(W.^2)) );
X = diag( sqrt(N./colsum(X.^2)) ) * X;

figure
plot(X');
return

% $$$ tsplot(X)
% $$$ tsplot(W')
% $$$ return

% Generate data
Y = W*X;

% Noise level
s2 = 0.1;

% Make noisy observations
Yn = Y + sqrt(s2)*randn(M,N);  % noise
Ynm = Yn; Ynm(rand(M,N)<0.5) = nan; % missing values
Ytest = Yn; Ytest(~isnan(Ynm)) = nan;

initthetaW

% Learn GP PCA
maxiter = 3e1;
Q = vbgppcamv(Ynm,D,inW,inX,covfunW,initthetaW,covfunX,initthetaX, ...
              'maxiter',maxiter, 'pseudodensityx', 0.4, 'pseudodensityw', ...
              0.4, 'loglikelihood', true, 'updatehyper', 5, 'maxsearchx', ...
              3, 'maxsearchw', 200, 'initpseudow', {pseudoW});

est_noise = 1/Q.tau
real_noise = s2

% Learn other models
if ppca
  Qppca = pca_full(Ynm,D,'maxiters',maxiter,'rotate2pca',true, ...
                   'algorithm','ppca');
  Yh_ppca = Qppca.A * Qppca.S + repmat(Qppca.Mu,1,N);
  testrmse_ppca = rmse(Ytest, Yh_ppca)
  noiselessrmse_ppca = rmse(Y, Yh_ppca)
end
if vbpca
  Qvbpca = vbpcamv(Ynm,M-1,'maxiters',maxiter);
  Yh_vbpca = Qvbpca.W * Qvbpca.X + repmat(Qvbpca.mu,1,N);
  testrmse_vbpca = rmse(Ytest, Yh_vbpca)
  noiselessrmse_vbpca = rmse(Y, Yh_vbpca)
end

Yh_gppca = Q.W * Q.X;
testrmse_gppca = rmse(Ytest, Yh_gppca)
noiselessrmse_gppca = rmse(Y, Yh_gppca)
figure
subplot(5,1,1);
plot(Y');
title('Original noiseless data')
subplot(5,1,2);
plot(Ynm');
title('Observed data')
subplot(5,1,3);
plot(Yh_gppca');
title('Reconstruction of GP VB PCA')
if vbpca
  subplot(5,1,4);
  plot(Yh_vbpca');
  title('Reconstruction of VB PCA')
end
if ppca
  subplot(5,1,5);
  plot(Yh_ppca');
  title('Reconstruction of PPCA')
end

vX = Q.varX;%diag(Q.CovX);
eX = 2*sqrt(vX);%sqrt( vX(reshape(1:(N*D), D, N)) );
vW = Q.varW;%diag(Q.CovW);
eW = 2*sqrt(vW);%sqrt( vW(reshape(1:(M*D), M, D)) );
tsgpplot(inX, Q.X', eX', 'pseudoinputs', {Q.pseudoX});

%tsgpplot(inW, Q.W, eW, 'pseudoinputs', {Q.pseudoW});
figure
for d=1:D
  subplot(D,1,d);
  mapproj('global-ellipse');
  mapplot(inW(1,:),inW(2,:),'ro');
  hold on
  mapplot(Q.pseudoW{d}(1,:),Q.pseudoW{d}(2,:),'k+');
  mapcoast
end

figure
mapproj('global-ellipse');
gpmapcolor(Q.pseudoW, Q.Wp, Q.CovWp, -180:5:180, -90:5:90, Q.logthetaW, covfunW);
% $$$ gpmapcolor({Q.pseudoW}, {Q.Wp}, {Q.CovWp}, 0:10:360, 0:10:180, {Q.logthetaW}, ...
% $$$            covfunW);

% $$$ exp(Q.logthetaW{1})
% $$$ exp(Q.logthetaW{2})
% $$$
% $$$ exp(Q.logthetaX{1})
% $$$ exp(Q.logthetaX{2})

% $$$ % DEBUG: (test the effect of CovXp)
% $$$ X = zeros(D,N);
% $$$ varX = zeros(D,N);
% $$$ for d=1:D
% $$$   [X(d,:), varX(d,:)] = gppred(Q.pseudoX{d}, Q.Xp{d}, 0, inX, covfunX{d}, ...
% $$$                                Q.logthetaX{d});
% $$$ end
% $$$ varX;
% $$$ tsgpplot(inX, X', 2*sqrt(varX'), 'pseudoinputs', {Q.pseudoX});
% $$$ exp(Q.logthetaX{1})

% $$$ figure
% $$$ pcolor(Q.CovXp{1});
% $$$ figure
% $$$ pcolor(Q.CovWp{1});

return

%inv(funcK_noiseless(D,[10;10;1;10;10]));return
[p, loglike] = gplearn(ymv, @(p) funcK_noisy(D,p), [10;10;1;300;300;0.1])
%[p, loglike] = gplearn(ymv, @(p) funcK_noisy(D,p), [300;300;10;1;1])
%[p, loglike] = gplearn(ymv, @(p) funcK_noisy(D,p), [300;300;10;5;5;10;0.1])
%[p, loglike] = gplearn(ymv, @(p) funcK_noisy(D,p), [1000;100;0.1;3;100;0.059;1;0.1])
%[p, loglike] = gplearn(ymv, @(p) funcK_noisy(D,p), [100;100;3;100;0.059;1;0.1])

%p(1:4) = [2000;30;0.1;10]
%p(3) = 0.1;
[mu,Cov] = gppred(th,tmv,ymv, @(x1,x2) funcK_noiseless(gpdist(x1,x2),p(1:(end-1))), p(end));
yh = mnorm_rnd(mu, Cov, 10);

figure
clf
plot(th, yh, 'r')
figure
clf
gpplot(th,mu,Cov);
hold on
plot(tmv,ymv,'k+')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function gpmapplot(coord, y)
x = coord';
% Use haversine distance measure
mydist = @(x1,x2) gpdist(x1,x2,@dist_coord);
D = mydist(x,x);
funcK = @(p) gpK(@() gpK_ratquad(D,p(1),p(2),p(3)), ...
                 @() gpK_noise(length(D),p(4)));
[p, loglike] = gplearn(y, funcK, [10;5;1;1e-5])
funcK_noiseless = @(D) gpK(@() gpK_ratquad(D, p(1),p(2),p(3)));
[LONI, LATI] = get_grid(40, 40);
xh = [LONI(:)';LATI(:)'];
[mu,Cov] = gppred(xh,x,y, @(x1,x2) funcK_noiseless(mydist(x1,x2)), p(4));
yh = mu;
ZI = reshape(yh, size(LONI));
plot_map
hold on
m_pcolor(LONI, LATI, ZI);

% Set nice colormap
colormap(climcolmap);
shading flat;
lim = max( -min(yh), max(yh) );
set(gca, 'clim', [-lim lim]);

return

% $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ function D = gpdist(x1, x2, funcDist)
% $$$ if isvector(x1)
% $$$   x1 = x1(:)';
% $$$   x2 = x2(:)';
% $$$ end
% $$$ if nargin < 3
% $$$   funcDist = @(z1,z2) abs(z1-z2);
% $$$ end
% $$$ %[X2,X1] = meshgrid(x2,x1);
% $$$ n1 = size(x1,2);
% $$$ n2 = size(x2,2);
% $$$ D = zeros(n1,n2);
% $$$ for i=1:n1
% $$$   for j=1:n2
% $$$     D(i,j) = funcDist(x1(:,i),x2(:,j));
% $$$   end
% $$$ end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = dist_coord(coord1, coord2, varargin)
% d = dist(X1, X2)
% returns geographical distance in kilometers
% input vectors must be the same size and shape!!
% Approximates earth as a sphere.
% Distance is calculated for

if nargin > 2
  nargin
end

% Quadratic mean radius from Wikipedia
R_avg = 6372.795477598;

% Convert to radians
q = pi / 180;
lon1 = coord1(1,:) * q;
lat1 = coord1(2,:) * q;
lon2 = coord2(1,:) * q;
lat2 = coord2(2,:) * q;

% Distance calculation (haversine law)
dlat = lat2 - lat1;
dlon = lon2 - lon1;
a = sin(dlat/2).^2 + cos(lat1).*cos(lat2).*(sin(dlon/2).^2);
c = 2 * atan2(sqrt(a), sqrt(1-a));
d = R_avg * c;

return

% $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ function [p,loglike] = gplearn(y, funcK, init_p)
% $$$ opts = optimset('GradObj', 'on');
% $$$ [p, negloglike] = fminunc(@(p) cost(y, funcK, p), init_p, opts);
% $$$ loglike = -negloglike;
% $$$
% $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ function [mu,Cov] = gppred(xh, x, y, funcK, noise)
% $$$ if nargin == 1
% $$$   x = [];
% $$$   y = [];
% $$$ end
% $$$
% $$$ Kxhx = funcK(xh,x);
% $$$ invKxx = inv(funcK(x,x)+noise^2*eye(length(x)));
% $$$ Kxhxh = funcK(xh,xh);
% $$$ if isempty(Kxhx)
% $$$   Kxhx = 0;
% $$$ end
% $$$ if isempty(invKxx)
% $$$   invKxx = 0;
% $$$ end
% $$$ mu = Kxhx*invKxx*y;
% $$$ if isempty(mu)
% $$$   mu = zeros(size(xh));
% $$$ end
% $$$ Cov = Kxhxh - Kxhx*invKxx*Kxhx';

% $$$ % $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ % $$$ function gpplot(x, mu, Cov)
% $$$ % $$$ e = sqrt(diag(Cov));
% $$$ % $$$ X = [x(1:end); x(end:-1:1)];
% $$$ % $$$ Y = [mu(1:end)+e; mu(end:-1:1)-e(end:-1:1)];
% $$$ % $$$ C = [.65,.65,.65];
% $$$ % $$$ fill(X,Y,C,'EdgeColor',C);
% $$$ % $$$ hold on
% $$$ % $$$ plot(x,mu,'k');
% $$$
% $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ function [K, dK] = gpK(varargin)
% $$$ K = 0;
% $$$ dK = [];
% $$$ for i=1:nargin
% $$$   if nargout == 1
% $$$     Knew = varargin{i}();
% $$$   else
% $$$     [Knew,dKnew] = varargin{i}();
% $$$     n = size(dKnew,3);
% $$$     if isempty(dK)
% $$$       dK = dKnew;
% $$$     else
% $$$       dK(:,:,end+(1:n)) = dKnew;
% $$$     end
% $$$   end
% $$$   K = K + Knew;
% $$$ end
% $$$
% $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ function [K, dK] = gpK_noise(n, p1);
% $$$ K = p1^2 * eye(n);
% $$$ if nargout >= 2
% $$$   dK = K * 2/p1;
% $$$ end
% $$$
% $$$ % $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ % $$$ function [K, dK] = gpK_sqexp(D, p1, p2)
% $$$ % $$$ K = p1^2*exp(-0.5*D.^2/(p2^2));
% $$$ % $$$
% $$$ % $$$ if nargout >= 2
% $$$ % $$$   dK = zeros([size(D), 2]);
% $$$ % $$$   dK(:,:,1) = K .* 2 / p1;
% $$$ % $$$   dK(:,:,2) = K .* (-0.5*D.^2) .* (-2*p2^(-3));
% $$$ % $$$ end
% $$$
% $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ function [K, dK] = gpK_decper(D, p1, p2, p3, p4)
% $$$ K = p1^2*exp(-0.5*(D.^2)/(p2^2)-2*sin(pi*D*p3).^2/(p4^2));
% $$$
% $$$ if nargout >= 2
% $$$   dK = zeros([size(D),4]);
% $$$   dK(:,:,1) = K .* 2 / p1;
% $$$   dK(:,:,2) = K .* (-0.5*D.^2) .* (-2*p2^(-3));
% $$$   dK(:,:,3) = K .* (-2*sin(pi*D*p3)*2) .* cos(pi*D*p3) .* (pi*D);
% $$$   dK(:,:,4) = K .* (-2*sin(pi*D*p3).^2) * (-2)*p4^(-3);
% $$$ end
% $$$
% $$$ % $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ % $$$ function [K,dK] = gpK_ratquad(D, p1, p2, p3)
% $$$ % $$$ %%% p3 = p3^2;
% $$$ % $$$ f = 1 + D.^2/(2*p3^2*p2^2);
% $$$ % $$$ K = p1^2*f.^(-p3^2);
% $$$ % $$$
% $$$ % $$$ if nargout >= 2
% $$$ % $$$   dK = zeros([size(D),3]);
% $$$ % $$$   dK(:,:,1) = K .* 2 / p1;
% $$$ % $$$   dK(:,:,2) = K .* (-p3^2).*f.^(-1) .* (-2).*D.^2/(2*p3^2)*p2^(-3);
% $$$ % $$$   dK(:,:,3) = K .* (-2*p3*log(f) + (-p3^2)./f.*D.^2/(2*p2^2) * (-2) * p3^(-3));
% $$$ % $$$ end
% $$$
% $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ function [f, df] = cost(y,funcK,p)
% $$$ [K,dK] = funcK(p);
% $$$ [f,df] = loglikelihood(y,K,dK);
% $$$ f = -f;
% $$$ df = -df;
% $$$
% $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ function [loglike, dloglike] = loglikelihood(y,K,dK)
% $$$ invK = inv(K);
% $$$ logdetK = log(det(K));
% $$$ if logdetK < -1e100;
% $$$   logdetK = -1e100;
% $$$ end
% $$$ n = length(y);
% $$$
% $$$ loglike = -0.5*y'*invK*y - 0.5*logdetK - 0.5*n*log(2*pi);
% $$$
% $$$ if nargout >= 2
% $$$   m = size(dK,3);
% $$$   dloglike = zeros(m,1);
% $$$   a = invK * y;
% $$$   W = a*a' - invK;
% $$$   for i = 1:m
% $$$     dloglike(i) = 0.5*sum(sum(W.*dK(:,:,i)));
% $$$   end
% $$$ end
% $$$
% $$$
% $$$ % $$$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% $$$ % $$$ function y = mnorm_rnd(mu, Cov, n)
% $$$ % $$$ m = max(length(mu), length(Cov));
% $$$ % $$$
% $$$ % $$$ opts.issym = true;
% $$$ % $$$ opts.isreal = true;
% $$$ % $$$ [V,D] = svd(Cov);
% $$$ % $$$ D(D<0) = 0;
% $$$ % $$$ D = sqrt(D);
% $$$ % $$$ A = V * D;
% $$$ % $$$
% $$$ % $$$ y = repmat(mu,1,n) + A*randn(m,n);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [output,extras] = urlread2(urlChar,method,body,headersIn,varargin)
%urlread2  Makes HTTP requests and processes response
%
%   [output,extras] = urlread2(urlChar, *method, *body, *headersIn, varargin)
%
%   * indicates optional inputs that must be entered in place
%
%   UNDOCUMENTED MATLAB VERSION
%
%   EXAMPLE CALLING FORMS
%   ... = urlread2(urlChar)
%   ... = urlread2(urlChar,'GET','',[],prop1,value1,prop2,value2,etc)
%   ... = urlread2(urlChar,'POST',body,headers)
%
%   FEATURES
%   =======================================================================
%   1) Allows specification of any HTTP method
%   2) Allows specification of any header. Very little is hard-coded
%      in for header handling.
%   3) Returns response status and headers
%   4) Should handle unicode properly ...
%
%   OUTPUTS
%   =======================================================================
%   output : body of the response, either text or binary depending upon
%            CAST_OUTPUT property
%   extras : (structure)
%       .allHeaders  - stucture, fields have cellstr values, HTTP headers may
%           may be repeated but will have a single field entry, with each
%           repeat's value another being another entry in the cellstr, for
%           example:
%               .Set_Cookie = {'first_value' 'second_value'}
%       .firstHeaders - (structure), variable fields, contains the first
%                      string entry for each field in allHeaders, this
%                      structure can be used to avoid dereferencing a cell
%                      for fields you expect not to be repeated ...
%                   EXAMPLE:
%                   .Response       : 'HTTP/1.1 200 OK'}
%                   .Server         : 'nginx'
%                   .Date           : 'Tue, 29 Nov 2011 02:23:16 GMT'
%                   .Content_Type   : 'text/html; charset=UTF-8'
%                   .Content_Length : '109155'
%                   .Connection     : 'keep-alive'
%                   .Vary           : 'Accept-Encoding, User-Agent'
%                   .Cache_Control  : 'max-age=60, private'
%                   .Set_Cookie     : 'first_value'
%       .status     - (structure)
%               .value  : numeric value of status, ex. 200
%               .msg    : message that goes along with status, ex. 'OK'
%       .url        - eventual url that led to output, this can change from
%                     the input with redirects, see FOLLOW_REDIRECTS
%       .isGood     - (logical) I believe this is an indicator of the presence of 400
%                      or 500 status codes (see status.value) but more
%                      testing is needed. In other words, true if status.value < 400.
%                      In code, set true if the response was obtainable without
%                      resorting to checking the error stream.
%
%   INPUTS
%   =======================================================================
%   urlChar   : The full url, must include scheme (http, https)
%   method    : examples: 'GET' 'POST' etc
%   body      : (vector)(char, uint8 or int8) body to write, generally used
%               with POST or PUT, use of uint8 or int8 ensures that the
%               body input is not manipulated before sending, char is sent
%               via unicode2native function with ENCODING input (see below)
%   headersIn : (structure array), use empty [] or '' if no headers are needed
%                   but varargin property/value pairs are, multiple headers
%                   may be passed in as a structure array
%       .name    - (string), name of the header, a name property is used
%                   instead of a field because the name must match a valid
%                   header
%       .value   - (string), value to use
%
%   OPTIONAL INPUTS (varargin, property/value pairs)
%   =======================================================================
%   CAST_OUTPUT      : (default true) output is uint8, useful if the body
%                       of the response is not text
%   ENCODING         : (default ''), ENCODING input to function unicode2native
%   FOLLOW_REDIRECTS : (default true), if false 3xx status codes will
%                       be returned and need to be handled by the user,
%                       note this does not handle javascript or meta tag
%                       redirects, just server based ones
%   READ_TIMEOUT     : (default 0), 0 means no timeout, value is in
%                       milliseconds
%
%   EXAMPLES
%   =======================================================================
%   GET:
%   --------------------------------------------
%   url    = 'http://www.mathworks.com/matlabcentral/fileexchange/';
%   query  = 'urlread2';
%   params = {'term' query};
%   queryString = http_paramsToString(params,1);
%   url = [url '?' queryString];
%   [output,extras] = urlread2(url);
%
%   POST:
%   --------------------------------------------
%   url    = 'http://posttestserver.com/post.php';
%   params = {'testChars' char([2500 30000]) 'new code' '?'};
%   [paramString,header] = http_paramsToString(params,1);
%   [output,extras] = urlread2(url,'POST',paramString,header);
%
%   From behind a firewall, use the Preferences to set your proxy server.
%
%   See Also:
%       http_paramsToString
%       unicode2native
%       native2unicode
%
%   Subfunctions:
%   fixHeaderCasing - small subfunction to fix case errors encountered in real
%       world, requires updating when casing doesn't match expected form, like
%       if someone sent the header content-Encoding instead of
%       Content-Encoding
%
%   Based on original urlread code by Matthew J. Simoneau
%
%   VERSION = 1.1

in.CAST_OUTPUT      = true;
in.FOLLOW_REDIRECTS = true;
in.READ_TIMEOUT     = 0;
in.ENCODING         = '';

%Input handling
%---------------------------------------
if ~isempty(varargin)
    for i = 1:2:numel(varargin)
        prop  = upper(varargin{i});
        value = varargin{i+1};
        if isfield(in,prop)
            in.(prop) = value;
        else
           error('Unrecognized input to function: %s',prop)
        end
    end
end

if ~exist('method','var') || isempty(method), method = 'GET'; end
if ~exist('body','var'), body = ''; end
if ~exist('headersIn','var'), headersIn = []; end

assert(usejava('jvm'),'Function requires Java')

import com.mathworks.mlwidgets.io.InterruptibleStreamCopier;
com.mathworks.mlwidgets.html.HTMLPrefs.setProxySettings %Proxy settings need to be set

%Create a urlConnection.
%-----------------------------------
urlConnection = getURLConnection(urlChar);
%For HTTP uses sun.net.www.protocol.http.HttpURLConnection
%Might use ice.net.HttpURLConnection but this has more overhead

%SETTING PROPERTIES
%-------------------------------------------------------
urlConnection.setRequestMethod(upper(method));
urlConnection.setFollowRedirects(in.FOLLOW_REDIRECTS);
urlConnection.setReadTimeout(in.READ_TIMEOUT);

for iHeader = 1:length(headersIn)
    curHeader = headersIn(iHeader);
    urlConnection.setRequestProperty(curHeader.name,curHeader.value);
end

if ~isempty(body)
    %Ensure vector?
    if size(body,1) > 1
        if size(body,2) > 1
            error('Input parameter to function: body, must be a vector')
        else
            body = body';
        end
    end

    if ischar(body)
        %NOTE: '' defaults to Matlab's default encoding scheme
        body = unicode2native(body,in.ENCODING);
    elseif ~(strcmp(class(body),'uint8') || strcmp(class(body),'int8'))
        error('Function input: body, should be of class char, uint8, or int8, detected: %s',class(body))
    end

    urlConnection.setRequestProperty('Content-Length',int2str(length(body)));
    urlConnection.setDoOutput(true);
    outputStream = urlConnection.getOutputStream;
    outputStream.write(body);
    outputStream.close;
else
    urlConnection.setRequestProperty('Content-Length','0');
end

%==========================================================================
%                   Read the data from the connection.
%==========================================================================
%This should be done first because it tells us if things are ok or not
%NOTE: If there is an error, functions below using urlConnection, notably
%getResponseCode, will fail as well
try
    inputStream = urlConnection.getInputStream;
    isGood = true;
catch ME
    isGood = false;
%NOTE: HTTP error codes will throw an error here, we'll allow those for now
%We might also get another error in which case the inputStream will be
%undefined, those we will throw here
    inputStream = urlConnection.getErrorStream;

    if isempty(inputStream)
        msg = ME.message;
        I = strfind(msg,char([13 10 9])); %see example by setting timeout to 1
        %Should remove the barf of the stack, at ... at ... at ... etc
        %Likely that this could be improved ... (generate link with full msg)
        if ~isempty(I)
            msg = msg(1:I(1)-1);
        end
        fprintf(2,'Response stream is undefined\n below is a Java Error dump (truncated):\n');
        error(msg)
    end
end

%POPULATING HEADERS
%--------------------------------------------------------------------------
allHeaders = struct;
allHeaders.Response = {char(urlConnection.getHeaderField(0))};
done = false;
headerIndex = 0;

while ~done
    headerIndex = headerIndex + 1;
    headerValue = char(urlConnection.getHeaderField(headerIndex));
    if ~isempty(headerValue)
        headerName = char(urlConnection.getHeaderFieldKey(headerIndex));
        headerName = fixHeaderCasing(headerName); %NOT YET FINISHED

        %Important, for name safety all hyphens are replace with underscores
        headerName(headerName == '-') = '_';
        if isfield(allHeaders,headerName)
            allHeaders.(headerName) = [allHeaders.(headerName) headerValue];
        else
            allHeaders.(headerName) = {headerValue};
        end
    else
        done = true;
    end
end

firstHeaders = struct;
fn = fieldnames(allHeaders);
for iHeader = 1:length(fn)
   curField = fn{iHeader};
   firstHeaders.(curField) = allHeaders.(curField){1};
end

status = struct(...
    'value',    urlConnection.getResponseCode(),...
    'msg',      char(urlConnection.getResponseMessage));

%PROCESSING OF OUTPUT
%----------------------------------------------------------
byteArrayOutputStream = java.io.ByteArrayOutputStream;
% This StreamCopier is unsupported and may change at any time. OH GREAT :/
isc = InterruptibleStreamCopier.getInterruptibleStreamCopier;
isc.copyStream(inputStream,byteArrayOutputStream);
inputStream.close;
byteArrayOutputStream.close;

if in.CAST_OUTPUT
    charset = '';

    %Extraction of character set from Content-Type header if possible
    if isfield(firstHeaders,'Content_Type')
        text = firstHeaders.Content_Type;
        %Always open to regexp improvements
        charset = regexp(text,'(?<=charset=)[^\s]*','match','once');
    end

    if ~isempty(charset)
        output = native2unicode(typecast(byteArrayOutputStream.toByteArray','uint8'),charset);
    else
        output = char(typecast(byteArrayOutputStream.toByteArray','uint8'));
    end
else
    %uint8 is more useful for later charecter conversions
    %uint8 or int8 is somewhat arbitary at this point
    output = typecast(byteArrayOutputStream.toByteArray','uint8');
end

extras              = struct;
extras.allHeaders   = allHeaders;
extras.firstHeaders = firstHeaders;
extras.status       = status;
%Gets eventual url even with redirection
extras.url          = char(urlConnection.getURL);
extras.isGood       = isGood;

end

function headerNameOut = fixHeaderCasing(headerName)
%fixHeaderCasing Forces standard casing of headers
%
%   headerNameOut = fixHeaderCasing(headerName)
%
%   This is important for field access in a structure which
%   is case sensitive
%
%   Not yet finished.
%   I've been adding to this function as problems come along

    switch lower(headerName)
        case 'location'
            headerNameOut = 'Location';
        case 'content_type'
            headerNameOut = 'Content_Type';
        otherwise
            headerNameOut = headerName;
    end
end

%==========================================================================
%==========================================================================
%==========================================================================

function urlConnection = getURLConnection(urlChar)
%getURLConnection
%
%   urlConnection = getURLConnection(urlChar)

% Determine the protocol (before the ":").
protocol = urlChar(1:find(urlChar==':',1)-1);

% Try to use the native handler, not the ice.* classes.
try
    switch protocol
        case 'http'
            %http://www.docjar.com/docs/api/sun/net/www/protocol/http/HttpURLConnection.html
            handler = sun.net.www.protocol.http.Handler;
        case 'https'
            handler = sun.net.www.protocol.https.Handler;
    end
catch ME
    handler = [];
end

% Create the URL object.
try
    if isempty(handler)
        url = java.net.URL(urlChar);
    else
        url = java.net.URL([],urlChar,handler);
    end
catch ME
    error('Failure to parse URL or protocol not supported for:\nURL: %s',urlChar);
end

% Get the proxy information using MathWorks facilities for unified proxy
% preference settings.
mwtcp = com.mathworks.net.transport.MWTransportClientPropertiesFactory.create();
proxy = mwtcp.getProxy();

% Open a connection to the URL.
if isempty(proxy)
    urlConnection = url.openConnection;
else
    urlConnection = url.openConnection(proxy);
end

end

author: plotly


text:
%% Copyright (C) 2014, 2016, 2018-2019, 2022 Colin B. Macdonald
%%
%% This file is part of OctSymPy.
%%
%% OctSymPy 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 software 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 software; see the file COPYING.
%% If not, see <http://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @documentencoding UTF-8
%% @defop  Method   @@sym mtimes {(@var{x}, @var{y})}
%% @defopx Operator @@sym {@var{x} * @var{y}} {}
%% Multiply symbolic matrices.
%%
%% Example:
%% @example
%% @group
%% syms x y
%% A = sym([1 2; 3 4])
%%   @result{} A = (sym 2×2 matrix)
%%       ⎡1  2⎤
%%       ⎢    ⎥
%%       ⎣3  4⎦
%% B = [x; y]
%%   @result{} B = (sym 2×1 matrix)
%%       ⎡x⎤
%%       ⎢ ⎥
%%       ⎣y⎦
%% A * B
%%   @result{} (sym 2×1 matrix)
%%       ⎡ x + 2⋅y ⎤
%%       ⎢         ⎥
%%       ⎣3⋅x + 4⋅y⎦
%% @end group
%% @end example
%% @end defop

function z = mtimes(x, y)

  cmd = { '(x,y) = _ins'
          'return x*y,' };

  z = pycall_sympy__ (cmd, sym(x), sym(y));

end

%!test
%! % scalar
%! syms x
%! assert (isa (x*2, 'sym'))
%! assert (isequal (2*sym(3), sym(6)))
%! assert (isequal (sym(2)*3, sym(6)))

%!test
%! % matrix-scalar
%! D = [0 1; 2 3];
%! A = sym(D);
%! assert (isa (2*A, 'sym'))
%! assert (isequal ( 2*A , 2*D  ))
%! assert (isequal ( A*2 , 2*D  ))

%!test
%! % matrix-matrix
%! D = [0 1; 2 3];
%! A = sym(D);
%! assert (isa (A*A, 'sym'))
%! assert (isequal ( A*A , D*D  ))

repo: octsympy


text:
function [im,header,timestamp,bb,mu] = mmf_read_frame(header,framei)

fp = header.fid;

% read from the current location if no frame input
if nargin < 2,
  loc = ftell(fp);
  framei = find(header.frame2file==loc,1);
  if isempty(framei),
    error('Could not find file location %d in frame index',loc);
  end
else
  fseek(fp,header.frame2file(framei),'bof');
end

[mu,header] = mmf_read_mean(header,'framei',framei);

im = mu;
im_idcode = fread(fp,1,'ulong'); %#ok<NASGU>
im_headersize = fread(fp,1,'int');
%im_depth = fread(fp,1,'int');
%im_nchannels = fread(fp,1,'int');
fseek(fp,8,'cof');
im_numblocks = fread(fp,1,'int');

% seek to the end of the header
fseek(fp,header.frame2file(framei)+im_headersize,'bof');

bb = nan(im_numblocks,4);
%for j = 1:im_numblocks,
  block_rect = fread(fp,4,'int');
  j=1;
while numel(block_rect) == 4 && block_rect(1)>=0
  %if block_rect(1)<0
      %disp(['skipping blocks ',num2str(framei),' ',mat2str(block_rect),sprintf(' block %d / %d',j,im_numblocks)])
      %break;
  %end
  bb(j,[2,1,4,3]) = block_rect;
  if block_rect(3)>header.height || block_rect(4) > header.width || any(block_rect<=0),
    break;
  end

  % then interlaced row ordered image data
  block_data = fread(fp,block_rect(3)*block_rect(4),'*uint8');
  block_data = reshape(block_data,[block_rect(3),block_rect(4)]);
  im(block_rect(1)+1:block_rect(1)+block_rect(3),block_rect(2)+1:block_rect(2)+block_rect(4)) = block_data;
  [block_rect,count] = fread(fp,4,'int');
  if count < 4,
    break;
  end
  j=j+1;
end
%disp([num2str(framei),' ',mat2str(block_rect),sprintf(' block %d / %d',j-1,im_numblocks)])
timestamp = nan;

author: kristinbranson

repo: JAABA


text:
function model = ivmEpUpdateM(model, index)

% IVMEPUPDATEM Update matrix M, L, varSigma and mu for EP.
% FORMAT
% DESC performs an EP update on the IVM model's represenations for
% a given point.
% ARG model : the mode for which the update will be done.
% ARG index : the index of the point for which the update will take
% place.
%
% SEEALSO : ivmEpUpdatePoint, ivmDowndateSites
%
% COPYRIGHT : Neil D. Lawrence, 2006

% IVM

d = length(model.I);
k = find(model.I == index);
if k == d % This point has just been included so EP update is irrelevant.
  return;
end
kvector = model.kern.Kstore(:, k);
for c = 1:length(model.Sigma)
  % Compute the kth row of  Sigma
  a = model.Sigma(c).M(:, index);
  s = kvector' - a'*model.Sigma(c).M;

  v = model.Sigma(c).Linv(k:end, k);
  sLambda_k = rocholhFactorise(v);

  % Update M
  model.Sigma(c).M(k:end, :) = ...
      rocholForeSub(sLambda_k, ...
		    model.Sigma(c).M(k:end, :));
  model.Sigma(c).M(k:end-1) = model.Sigma(c).M(k+1:end);

  % Update L
  t2 = model.Sigma(c).Linv(k, 1:k);
  slambda_11 = sLambda_k.u(1)*sLambda_k.v(1);
  slambda_21 = sLambda_k.v*sLambda_k.u(1);
  slambda_21(1) = [];
  sLambda_22 = sLambda_k;
  sLambda_22.s(1) = [];
  sLambda_22.u(1) = [];
  sLambda_22.v(1) = [];
  sLambda_22.n = sLambda_22.n-1;
  T3prime =  rocholForeSub(sLambda_22, model.Sigma(c).Linv(k+1:end, : ...
                                                    ));
  T3prime(:, 1:k) = T3prime(:, 1:k) ...
      - rocholForeSub(sLambda_22, slambda_21*1/slambda_11)*t2;
  v = T3prime(:, k+1:end)\T3prime(:, k);
  sVtilde = rocholFactorise(v);
  V = rocholTransMultiply(sVtilde, T3prime(:, k+1:end)')';
  T33inv = rocholTransMultiply(sLambda_22, model.Sigma(c).L(k+1:end, ...
                                                    k+1:end)')';
  invV = rocholForeSub(sVtilde, T33inv);
  model.Sigma(c).Linv(k:end-1, 1:end-1) = [T3prime(:, 1:k-1) V];
  model.Sigma(c).L(k:end-1, 1:end-1) = [-invV*T3prime(:, 1:k-1)* ...
                      model.Sigma(c).L(1:k-1, 1:k-1) invV];
  %/~
  oldVarSigma = model.varSigma;
  %~/
  model.varSigma(:, c) = model.varSigma(:, c) + ((model.nu(index, c)*s).*s)';
  %/~
  if any(model.varSigma(:, c)<0)
    warning('Variance less than zero')
  end
  %~/
  model.mu(:, c) = model.mu(:, c) + model.g(index, c)*s';

end
if length(model.Sigma)==1 & size(model.y, 2)>1
  for c = 2:size(model.y, 2)
    model.varSigma(:, c) = model.varSigma(:, c) ...
        - ((model.nu(index, c)*s).*s)';
    model.mu(:, c) = model.mu(:, c) + model.g(index, c)*s';
    %/~
    if any(model.varSigma(:, c)<0)
      warning('Variance less than zero')
    end
    %~/
  end
end

%model = ivmDowndateSites(model, index);
%model.J = [model.J index];

% Ensure nu and g now accurately reflect mu and varsigma.
model = ivmUpdateNuG(model, index);
% Update the site parameter.
model = ivmUpdateSites(model, index);
d = model.d;

% Compute the values of M for the new point
for c = 1:length(model.Sigma)
  a = model.Sigma(c).M(1:end-1, index);
  s = kvector' - a'*model.Sigma(c).M(1:end-1, ...
                                                    :);
  lValInv = sqrt(model.nu(index, c));
  %/~
  % If Nu is so low then the included data-point isn't really useful.
  if lValInv < 1e-16
    warning(['Square root of nu is ' num2str(lValInv)])
  end
  %~/
  model.Sigma(c).M(end, :) = lValInv*s;
  ainv = (-a*lValInv)'/tril(model.Sigma(c).L(1:end-1, 1:end-1));
  model.Sigma(c).L(end, :) = [a' 1/lValInv];
  model.Sigma(c).Linv(end, :) = [ainv lValInv];
  % make sure Matlab knows they are lower triangular.
  model.Sigma(c).L = tril(model.Sigma(c).L);
  model.Sigma(c).Linv = tril(model.Sigma(c).Linv);
  %/~
  oldVarSigma = model.varSigma(:, c);
  if any(model.varSigma(:, c)<0)
    warning('Variance less than zero')
  end
  %~/
  model.varSigma(:, c) = model.varSigma(:, c) - ((model.nu(index, c)*s).*s)';
  %/~
  if any(model.varSigma(:, c)<0)
    warning('Variance less than zero')
  end
  %~/
  model.mu(:, c) = model.mu(:, c) + model.g(index, c)*s';
end
if length(model.Sigma)==1 & size(model.y, 2)>1
  for c = 2:size(model.y, 2)
    model.varSigma(:, c) = model.varSigma(:, c) ...
        - ((model.nu(index, c)*s).*s)';
    model.mu(:, c) = model.mu(:, c) + model.g(index, c)*s';
    %/~
    if any(model.varSigma(:, c)<0)
      warning('Variance less than zero')
    end
    %~/
  end
end
% Swap the columns of the kernel matrix about.
model.kern.Kstore(:, k:end-1) = model.kern.Kstore(:, k+1:end);
model.kern.Kstore(:, end) = kvector;

% Move point to end of active set.
model.I(k:end-1) = model.I(k+1:end);
model.I(end) = index;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function result = polygon_integral_x ( n, v )

%*****************************************************************************80
%
%% POLYGON_INTEGRAL_X integrates the function X over a polygon.
%
%  Discussion:
%
%    The polygon is bounded by the points (X(1:N), Y(1:N)).
%
%    INTEGRAL = (1/6) * sum ( 1 <= I <= N )
%      ( X(I)^2 + X(I) * X(I-1) + X(I-1)^2 ) * ( Y(I) - Y(I-1) )
%
%    where X(0) and Y(0) should be replaced by X(N) and Y(N).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    S F Bockman,
%    Generalizing the Formula for Areas of Polygons to Moments,
%    American Mathematical Society Monthly,
%    1989, pages 131-132.
%
%  Parameters:
%
%    Input, integer N, the number of vertices of the polygon.
%    N should be at least 3 for a nonzero result.
%
%    Input, real V(2,N), the coordinates of the vertices
%    of the polygon.  These vertices should be given in counter-clockwise order.
%
%    Output, real RESULT, the value of the integral.
%
  result = 0.0;

  if ( n < 3 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'POLYGON_INTEGRAL_X - Warning!\n' );
    fprintf ( 1, '  The number of vertices must be at least 3.\n' );
    fprintf ( 1, '  The input value of N = %d\n', n );
    return
  end

  for i = 1 : n

    if ( i == 1 )
      im1 = n;
    else
      im1 = i - 1;
    end

    result = result + ( v(1,i).^2 + v(1,i) * v(1,im1) + v(1,im1).^2 ) ...
      * ( v(2,i) - v(2,im1) );

  end

  result = result / 6.0;

  return
end


text:
function x = sample_fs(xf, grid_sz)

% Samples the Fourier series

sz = [size(xf,1) size(xf,2)];

if nargin < 2
    x = prod(sz) * cifft2(xf);
else
    if any(grid_sz < sz)
        error('The grid size must be larger than or equal to the signal size')
    end
    tot_pad = grid_sz - sz;
    pad_sz = ceil(tot_pad/2);
    xf_pad = padarray(xf, pad_sz);
    if any(mod(tot_pad,2) == 1)
        % Handle odd padding
        xf_pad = xf_pad(1:end-mod(tot_pad(1),2), 1:end-mod(tot_pad(2),2), :, :);
    end
    x = prod(grid_sz) * cifft2(xf_pad);
end

repo: CFWCR


text:
function a = diagonal ( m, n, x )

%*****************************************************************************80
%
%% DIAGONAL returns the DIAGONAL matrix.
%
%  Formula:
%
%    if ( I = J )
%      A(I,J) = X(I)
%    else
%      A(I,J) = 0
%
%  Example:
%
%    M = 5, N = 5, X = ( 1, 2, 3, 4, 5 )
%
%    1 0 0 0 0
%    0 2 0 0 0
%    0 0 3 0 0
%    0 0 0 4 0
%    0 0 0 0 5
%
%  Square Properties:
%
%    A is banded, with bandwidth 1.
%
%    A is nonsingular if, and only if, each X(I) is nonzero.
%
%    The inverse of A is a diagonal matrix with diagonal values 1/X(I).
%
%    A is symmetric: A' = A.
%
%    Because A is symmetric, it is normal.
%
%    Because A is normal, it is diagonalizable.
%
%    LAMBDA(1:N) = X(1:N).
%
%    The matrix of eigenvectors of A is the identity matrix.
%
%    det ( A ) = product ( 1 <= I <= N ) X(I).
%
%    Because A is diagonal, it has property A (bipartite).
%
%    The family of matrices is nested as a function of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns of A.
%
%    Input, real X(min(M,N)), the diagonal entries of A.
%
%    Output, real A(M,N), the matrix.
%
  a = zeros ( m, n );

  for i = 1 : m
    for j = 1 : n

      if ( i == j )
        a(i,j) = x(i);
      else
        a(i,j) = 0.0;
      end

    end
  end

  return
end


text:
%% A Unit Test Class for mklJac
classdef mklJac_tests < matlab.unittest.TestCase

    properties
        absTol = 1e-8;
    end

    % Unit Tests
    methods (Test)

        %-- Valid Operation --%
        function scalarDiff(testCase)
            testCase.verifyEqual(1, mklJac(@(x) sin(x), 0), 'AbsTol', testCase.absTol);
            testCase.verifyEqual(1, mklJac(@(x) sin(x), 0, 1), 'AbsTol', testCase.absTol);
            testCase.verifyEqual(-1, mklJac(@(x) sin(x), pi), 'AbsTol', testCase.absTol);
            testCase.verifyEqual(0, mklJac(@(x) cos(x), pi), 'AbsTol', testCase.absTol);
            testCase.verifyEqual(0, mklJac(@(x) cos(x), -pi), 'AbsTol', testCase.absTol);
            testCase.verifyEqual(-1, mklJac(@(x) cos(x), pi/2), 'AbsTol', testCase.absTol);
            testCase.verifyEqual(1, mklJac(@(x) cos(x), -pi/2), 'AbsTol', testCase.absTol);
            testCase.verifyEqual(2, mklJac(@(x) 2*x, 0), 'AbsTol', testCase.absTol);
            testCase.verifyEqual(2, mklJac(@(x) 2*x, 1), 'AbsTol', testCase.absTol);
            testCase.verifyEqual(2, mklJac(@(x) 2*x, 2), 'AbsTol', testCase.absTol);
            testCase.verifyEqual(-2, mklJac(@(x) -2*x, 0), 'AbsTol', testCase.absTol);
            testCase.verifyEqual(-2, mklJac(@(x) -2*x, 1), 'AbsTol', testCase.absTol);
            testCase.verifyEqual(-2, mklJac(@(x) -2*x, 2), 'AbsTol', testCase.absTol);
        end

        function scalarDiffSweep(testCase)
            fun = @(x) 3*cos(x^2) + 0.5*sin(x/2);
            grad = @(x) cos(x/2)/4 - 6*x*sin(x^2);
            x = linspace(-pi,pi);
            for i = 1:length(x)
                testCase.verifyEqual(grad(x(i)), mklJac(fun, x(i)), 'AbsTol', testCase.absTol);
            end
        end

        function vectorDiff(testCase)
            testCase.verifyEqual(eye(3), mklJac(@(x) sin(x), zeros(3,1)), 'AbsTol', testCase.absTol);
            testCase.verifyEqual(eye(6), mklJac(@(x) sin(x), zeros(6,1)), 'AbsTol', testCase.absTol);
            testCase.verifyEqual(-eye(3), mklJac(@(x) sin(x), pi*ones(3,1)), 'AbsTol', testCase.absTol);
            testCase.verifyEqual(zeros(3), mklJac(@(x) sin(x), pi/2*ones(3,1)), 'AbsTol', testCase.absTol);
        end

        function vectorDiffSweep(testCase)
            fun = @(x) 3*cos(x(1)^2)*sin(x(2)) + 0.5*sin(x(2)/2);
            grad = @(x) [ -6*x(1)*sin(x(1)^2)*sin(x(2)), cos(x(2)/2)/4 + 3*cos(x(1)^2)*cos(x(2))];
            x1 = linspace(-pi,pi);
            x2 = linspace(pi,-pi);
            for i = 1:length(x1)
                testCase.verifyEqual(grad([x1(i);x2(i)]), mklJac(fun, [x1(i);x2(i)]), 'AbsTol', testCase.absTol);
            end
        end

        function vectorFunDiffSweep(testCase)
            fun = @(x) [100*(x(2)-x(1)^2); 1 - x(1)];
            grad = @(x) [-200*x(1) 100; -1 0];
            x = linspace(-1,1,10);
            for i = 1:length(x)
                for j = 1:length(x)
                    testCase.verifyEqual(grad([x(i);x(j)]), mklJac(fun, [x(i);x(j)]), 'AbsTol', testCase.absTol*10); %this is a hard one
                end
            end
        end

        function autoSizeIdentify(testCase)
            testCase.verifyEqual(1, numel(mklJac(@(x) sin(x), zeros(1,1))));
            testCase.verifyEqual(4, numel(mklJac(@(x) sin(x), zeros(2,1))));
            testCase.verifyEqual(9, numel(mklJac(@(x) sin(x), zeros(3,1))));
            testCase.verifyEqual(16, numel(mklJac(@(x) sin(x), zeros(4,1))));
        end

        %-- Input Args --%
        function inputArgs(testCase)
            testCase.verifyError(@() mklJac(@(x) sin(x)), 'OPTIMex:InputError');   % not enough args
            testCase.verifyError(@() mklJac(1, 1), 'OPTIMex:InputError');   % not fcn handle
            testCase.verifyError(@() mklJac(@(x) sin(x), int16(1)), 'OPTIMex:InputError'); % wrong input type
            testCase.verifyError(@() mklJac(@(x) sin(x), 1i), 'OPTIMex:InputError'); % wrong input type
            testCase.verifyError(@() mklJac(@(x) sin(x), [1 1; 1 1]), 'OPTIMex:InputError'); % wrong input type
            testCase.verifyError(@() mklJac(@(x) sin(x), 1, int16(1)), 'OPTIMex:InputError'); % wrong input type
            testCase.verifyError(@() mklJac(@(x) sin(x), 1, [1;1]), 'OPTIMex:InputError'); % wrong input type
            testCase.verifyError(@() mklJac(@(x) sin(x), 1, 0), 'OPTIMex:InputError'); % wrong val
            testCase.verifyError(@() mklJac(@(x) sin(x), 1, 1e9), 'OPTIMex:InputError'); % wrong val
            testCase.verifyError(@() mklJac(@(x) sin(x), 1, 1, int16(1)), 'OPTIMex:InputError'); % wrong input type
            testCase.verifyError(@() mklJac(@(x) sin(x), 1, 1, [1;1]), 'OPTIMex:InputError'); % wrong input type
            testCase.verifyError(@() mklJac(@(x) sin(x), 1, 1, 1e-18), 'OPTIMex:InputError'); % wrong val
            testCase.verifyError(@() mklJac(@(x) sin(x), 1, 1, 1.1), 'OPTIMex:InputError'); % wrong val
            testCase.verifyError(@() mklJac(@(x) sin(x), [1;1], 1), 'OPTIMex:DataError'); % wrong length
            testCase.verifyError(@() mklJac(@(x) sin(x), [1;1], 3), 'OPTIMex:DataError'); % wrong length
            testCase.verifyError(@() mklJac(@(x) sin(x), 1, 2), 'OPTIMex:DataError'); % wrong length (something odd about this ut)
        end
    end

end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
% Rounds a scalar, matrix or vector to a specified number of decimal places
% Format is roundoff(number,decimal_places)

function y = roundoff(number,decimal_places)

[INeg,JNeg] = find( number<0 ); % Negative numbers

if ~isempty(INeg)
   IndNeg = sub2ind(size(number),INeg,JNeg);
   Number = abs(number);
else
   Number = number;
end

decimals = 10.^decimal_places;
y1 = fix(decimals * Number + 0.5)./decimals;

if ~isempty(INeg)
   y1(IndNeg) = -y1(IndNeg);
end

y = y1;

repo: mcbench-benchmarks


text:
npts=1000;
ntaps = 5;
b =ones(1, ntaps) / ntaps; % create filter coefficients for 5- point moving average

x=(rand(npts,1)*2)-1; % raw data from -1 to +1
filtered_data=filter(b,1,x); % filter using 'b' coefficients

subplot(2,1,1); % 1st subplot
plot(x); % plot raw data
title('Raw Data');
subplot(2,1,2); % 2nd subplot
plot(filtered_data); %plot filtered data
title('Filtered Data');
xlabel('Time');

% Perform FFT on original and filtered signal
fftpts=npts; % number points in FFT
hpts=fftpts/2; % half number of FFT points
x_fft=abs(fft(x))/hpts; %scaled FFT of original signal
filtered_fft=abs(fft(filtered_data))/hpts; %scaled FFT of filtered signal

subplot(2,1,1) %1st subplot
plot(x_fft(1:hpts)); %plot first half of data points
title('Raw Data');
subplot(2,1,2) %2nd subplot
plot(filtered_fft(1:hpts));%plot first half data points
title('Filtered Data');
xlabel('Frequency');

author: yandld


text:
function g = sqrt(f)
%SQRT   Square root of a BALLFUN.
%   SQRT(F) is the square root of a BALLFUN F.

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

g = compose( f, @sqrt );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function tapas_hgf_ar1_binary_mab_plotTraj(r)
% Plots the estimated or generated trajectories for the binary HGF perceptual model for multi-armed
% bandit situations.
%
% Usage example:  est = tapas_fitModel(responses, inputs); tapas_hgf_binary_plotTraj(est);
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2013 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Optional plotting of standard deviations (true or false)
plotsd = true;

% Optional plotting of responses (true or false)
ploty = true;

% Set up display
scrsz = get(0,'screenSize');
outerpos = [0.2*scrsz(3),0.2*scrsz(4),0.8*scrsz(3),0.8*scrsz(4)];
figure(...
    'OuterPosition', outerpos,...
    'Name', 'HGF trajectories');

% Set up colors
colors = [1 0 0; 0.67 0 1; 0 0.67 1; 0.67 1 0];

% Number of bandits
b = r.c_prc.n_bandits;

% Number of trials
n = size(r.u,1);

% Time axis
t = ones(1,n);

ts = cumsum(t);
ts = [0, ts];

% Number of levels
try
    l = r.c_prc.n_levels;
catch
    l = (length(r.p_prc.p)+1)/5;
end

% Upper levels
for j = 1:l-2

    % Subplots
    subplot(l,1,j);

    if plotsd == true
        upperprior = r.p_prc.mu_0(l-j+1) +sqrt(r.p_prc.sa_0(l-j+1));
        lowerprior = r.p_prc.mu_0(l-j+1) -sqrt(r.p_prc.sa_0(l-j+1));
        upper = [upperprior; r.traj.mu(:,l-j+1,1)+sqrt(r.traj.sa(:,l-j+1,1))];
        lower = [lowerprior; r.traj.mu(:,l-j+1,1)-sqrt(r.traj.sa(:,l-j+1,1))];

        plot(0, upperprior, 'ob', 'LineWidth', 1);
        hold all;
        plot(0, lowerprior, 'ob', 'LineWidth', 1);
        fill([ts, fliplr(ts)], [(upper)', fliplr((lower)')], ...
             'b', 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
    end
    plot(ts, [r.p_prc.mu_0(l-j+1); r.traj.mu(:,l-j+1,1)], 'b', 'LineWidth', 2);
    hold all;
    plot(0, r.p_prc.mu_0(l-j+1), 'ob', 'LineWidth', 2); % prior
    xlim([0 ts(end)]);
    title(['Posterior expectation of x_' num2str(l-j+1)], 'FontWeight', 'bold');
    ylabel(['\mu_', num2str(l-j+1)]);
end

% Second level
subplot(l,1,l-1)

if plotsd == true
    for j=1:b
        upperprior = r.p_prc.mu_0(2) +sqrt(r.p_prc.sa_0(2));
        lowerprior = r.p_prc.mu_0(2) -sqrt(r.p_prc.sa_0(2));
        upper = [upperprior; r.traj.mu(:,2,j)+sqrt(r.traj.sa(:,2,j))];
        lower = [lowerprior; r.traj.mu(:,2,j)-sqrt(r.traj.sa(:,2,j))];

        plot(0, upperprior, 'o', 'Color', colors(j,:), 'LineWidth', 1);
        hold all;
        plot(0, lowerprior, 'o', 'Color', colors(j,:), 'LineWidth', 1);
        fill([ts, fliplr(ts)], [(upper)', fliplr((lower)')], ...
             colors(j,:), 'EdgeAlpha', 0, 'FaceAlpha', 0.15);
    end
end
for j=1:b
    plot(ts, [r.p_prc.mu_0(1); r.traj.mu(:,2,j)], 'Color', colors(j,:), 'LineWidth', 2);
    hold all;
    plot(0, r.p_prc.mu_0(1), 'o', 'Color', colors(j,:), 'LineWidth', 2); % prior
end
xlim([0 ts(end)]);
title('Posterior expectations of x_2', 'FontWeight', 'bold');
ylabel('\mu_2');

% Input level
subplot(l,1,l);

for j=1:b
    plot(ts, [tapas_sgm(r.p_prc.mu_0(2), 1); tapas_sgm(r.traj.mu(:,2,j), 1)], 'Color', colors(j,:), 'LineWidth', 2);
    hold all;
    plot(0, tapas_sgm(r.p_prc.mu_0(2), 1), 'o', 'Color', colors(j,:), 'LineWidth', 2); % prior
end
plot(ts(2:end), r.u(:,1), '.', 'Color', [0 0 0]); % inputs
plot(ts(2:end), r.traj.wt(:,1), 'k') % implied learning rate
if (ploty == true) && ~isempty(find(strcmp(fieldnames(r),'y'))) && ~isempty(r.y)
    y = r.y(:,1);
    if ~isempty(find(strcmp(fieldnames(r),'irr')))
        y(r.irr) = NaN; % weed out irregular responses
        plot(ts(r.irr)+1,  1.08.*ones([1 length(r.irr)]), 'x', 'Color', [1 0.7 0], 'Markersize', 11, 'LineWidth', 2); % irregular responses
    end
    for j=1:b
        plot(find(y==j), 1.08*ones([1 length(find(y==j))]), '.', 'Color', colors(j,:)); % responses
    end
    title(['Response y, input u (black dots), learning rate (fine black line), and posterior expectation of reward s(\mu_2) ', ...
           '(colour coded) for \phi=', num2str(r.p_prc.phi(2:end)), ', m=', num2str(r.p_prc.m(2:end)), ', \kappa=', ...
           num2str(r.p_prc.ka(2:end)), ', \omega=', num2str(r.p_prc.om(2:end))], ...
      'FontWeight', 'bold');
    ylabel('y, u, s(\mu_2)');
    axis([0 ts(end) -0.15 1.15]);
else
    title(['Input u (black dots), learning rate (fine black line), and posterior expectation of input s(\mu_2) ', ...
           '(red) for \rho=', num2str(r.p_prc.rho(2:end)), ', \kappa=', ...
           num2str(r.p_prc.ka(2:end)), ', \omega=', num2str(r.p_prc.om(2:end))], ...
      'FontWeight', 'bold');
    ylabel('u, s(\mu_2)');
    axis([0 ts(end) -0.1 1.1]);
end
plot(ts(2:end), 0.5, 'k');
xlabel('Trial number');
hold off;


text:
% StackExchange Signal Processing Q59325
% https://dsp.stackexchange.com/questions/59325
% Learning the Coefficients of Auto Regressive Model Using Least Mean
% Squares Filter for Signal Prediction
% References:
%   1.
% Remarks:
%   1.  B
% TODO:
% 	1.  C
% Release Notes Royi Avital RoyiAvital@yahoo.com
% - 1.0.000     14/04/2022
%   *   First release.

%% General Parameters

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Constants

%% Simulation Parameters

% Signal Generation
numSamples  = 2000;
numSignals = 30;

% Estimaiton
arModelOrder = 101;
numSamplesEst = round(0.75 * numSamples);

%% Generate / Load Data

vA = rand(numSignals, 1) + 0.35;
vF = 5 * rand(numSignals, 1) + 1;
vP = pi * rand(numSignals, 1);

vT = linspace(0, 5, numSamples).';
vX = zeros(numSamples, 1);

for ii = 1:length(vA)
    vX(:) = vX(:) + vA(ii) * sin(2 * pi  * vF(ii) * vT + vP(ii));
end

vW = zeros(arModelOrder, 1); %<! LMS initialization

%% LMS

vD = vX((arModelOrder + 1):(arModelOrder + 1 + numSamplesEst));
vY = vX(1:(numSamplesEst + 1));

vWW = LmsFilter(vW, vY, vD, arModelOrder, numSamplesEst, 5e-4, OFF);

vYY = filter(1, [1; vWW], vX); %<! Prediction by the AR Model as estimated by the LMS

%% Analysis

figureIdx = figureIdx + 1;

hFigure = figure('Position', figPosLarge);
hAxes   = axes(hFigure);
set(hAxes, 'NextPlot', 'add');
hLineObj = plot(vT(1:numSamplesEst), vX(1:numSamplesEst));
set(hLineObj, 'LineWidth', lineWidthNormal);
hLineObj = plot(vT((numSamplesEst + 1):end), vX((numSamplesEst + 1):end));
set(hLineObj, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', {['Linear Combination of ', num2str(numSignals), ' Sine Signals']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Time Index [Sec]']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['Value']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['Input Signal: Samples for LMS'], ['Input Signal: Samples to Predict']});

if(generateFigures == ON)
    % saveas(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
    print(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png'], '-dpng', '-r0'); %<! Saves as Screen Resolution
end

figureIdx = figureIdx + 1;

hFigure = figure('Position', figPosLarge);
hAxes   = axes(hFigure);
set(hAxes, 'NextPlot', 'add');
hLineObj = plot(vT, vX);
set(hLineObj, 'LineWidth', lineWidthNormal);
hLineObj = plot(vT, vYY);
set(hLineObj, 'LineWidth', lineWidthNormal, 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', {['LMS Filter Prediction'], ['The RMSE of the Estimation: ', num2str(sqrt(mean((vYY - vX) .^ 2)))]}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Time Index [Sec]']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['Value']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['Input Signal'], ['Predicted Signal (Model Order: ', num2str(arModelOrder), ')']});

if(generateFigures == ON)
    % saveas(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
    print(hFigure, ['Figure', num2str(figureIdx, figureCounterSpec), '.png'], '-dpng', '-r0'); %<! Saves as Screen Resolution
end

%% Restore Defaults

% set(0, 'DefaultFigureWindowStyle', 'normal');
% set(0, 'DefaultAxesLooseInset', defaultLoosInset);

author: RoyiAvital

repo: StackExchangeCodes


text:
function PlotAdaptiveContour2D(u, levels, tol)

% function PlotAdaptiveContour2D(u, levels, tol)
% Purpose: adaptively refine the mesh to approximately locate isocontours

Globals2D;

% build interpolation matrix (coarse->fine)
ri(:,1) = 0.5*(-(r+s)*( 0) + (1+r)*( 0) + (1+s)*(-1) );
si(:,1) = 0.5*(-(r+s)*(-1) + (1+r)*( 0) + (1+s)*( 0) );

ri(:,2) = 0.5*(-(r+s)*(-1) + (1+r)*( 0) + (1+s)*(-1) );
si(:,2) = 0.5*(-(r+s)*(-1) + (1+r)*(-1) + (1+s)*( 0) );

ri(:,3) = 0.5*(-(r+s)*( 1) + (1+r)*( 0) + (1+s)*( 0) );
si(:,3) = 0.5*(-(r+s)*(-1) + (1+r)*( 0) + (1+s)*(-1) );

ri(:,4) = 0.5*(-(r+s)*(-1) + (1+r)*(-1) + (1+s)*( 0) );
si(:,4) = 0.5*(-(r+s)*( 1) + (1+r)*( 0) + (1+s)*( 0) );

%interp = Vandermonde2D(N, ri(:), si(:))*invV;
interp = InterpMatrix2D(ri(:),si(:));

Nlevels = length(levels);

xref = x; yref = y; uref = u; Kref = K;

sk = 1;
F = spalloc(Np,Np,1);
for i=0:N % old ordering
  for j=0:N - i
    if(i+j<=1), F(sk,sk) = 1.; end;
    sk = sk+1;
  end
end

hold on
for nlev=1:Nlevels
  lev = levels(nlev);

  xref = x; yref = y; uref = u; Kref = K; Jref = J;

  err = 1;
  while(err > tol)

    umin = min(uref, [], 1);
    umax = max(uref, [], 1);

    refineflag = (umin <= lev & umax >= lev);

    toref = find( refineflag);
    Nref = length(toref);

    uref = reshape(interp*uref(:,toref), Np, 4*Nref);
    xref = reshape(interp*xref(:,toref), Np, 4*Nref);
    yref = reshape(interp*yref(:,toref), Np, 4*Nref);

    Kref = 4*Nref;

    ufilt = V*(F*(invV*uref));
    err = max(max(abs(ufilt-uref)));

  end

  ri = [-1;1;-1]; si = [-1;-1;1]; refNp = length(ri);
  interp1 = InterpMatrix2D(ri,si);
  xref = interp1*xref; yref = interp1*yref; uref = interp1*uref;

  ltri = delaunay(ri,si);
  Nltri = size(ltri,1);
  tri = zeros(Kref*Nltri, Nfaces);
  ks = (1:Nltri)';
  for k=1:Kref
    tri(ks,:) = ltri + (k-1)*refNp;
    ks = ks + Nltri;
  end

  PlotContour2D(tri, xref(:), yref(:), uref(:), levels)

end
hold off
return

repo: Aero-matlab


text:
function boxOut = bbox_scale2(boxIn,scale,szOut)
% Copyright (C) 2016 Hakan Bilen.
% 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).

if isempty(boxIn), boxOut = []; return; end

boxOut = scale * (boxIn-1) + 1;

boxOut = [max(1,round(boxOut(:,1))),...
  max(1,round(boxOut(:,2))),...
  min(szOut(1),round(boxOut(:,3))),...
  min(szOut(2),round(boxOut(:,4)))];

author: guosheng

repo: refinenet

save_path: github-repos/MATLAB/guosheng-refinenet


text:
function [x,rec] = bsearch(fun,x0,opt,varargin)
% BSEARCH Finds the minimum of a combinatorial function using backward search
%
%   X = BSEARCH(FUN, X0) attempts to find a combination of elements
%   of X0 which locally minimize the function FUN using backward
%   search (backward elimination). FUN accepts input X and returns
%   scalar function value F evaluated at X. X0 must be a vector of
%   initial indexes. Returned X contains indexes locally minimizing
%   the function FUN.
%
%   X = BSEARCH(FUN, X0, OPTIONS) allows use of optional
%   search parameters. See BSEARCH_OPT for details.
%
%   X = BSEARCH(FUN, X0, OPTIONS, P1, ..., Pn) P1,...,Pn are
%   additional parameters passed to the function FUN.
%
%   [X,REC] = BSEARCH(FUN, X0, ...) returns record of search as a
%   array of structs.
%
%   In backward search the elements are removed one at the time.
%   The element ehich removal minimizes the function value is
%   removed at each round. In a case of tie choose randomly.
%
% Copyright (c) Aki Vehtari (2007)

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

% 2007-02-23  Aki Vehtari  <Aki.Vehtari@hut.fi>
%             Use same parameter order as in fminunc etc.
%             Minimize instead of maximize
%             Rewrite.

% Columns
nx = size(x0,2);  % number of elements
nremoved = 0;     % number of removed elements
% Options
opt = bsearch_opt(opt);

% Base value
value = feval(fun,x0,varargin{:});
if (opt.display >= 1)
  fprintf(' Base value: %.4g\n Elements: %s\n', value, num2str(x0));
end
minvalue=value;
rec(nremoved+1).chosen = x0;
rec(nremoved+1).candidates = [];
rec(nremoved+1).values = value;

% The loop
while (nx > opt.nsel)
  values=zeros(nx,1);
  for i1=1:nx
    values(i1) = feval(fun,setdiff(x0,x0(i1)),varargin{:});
    if (opt.display >= 2)
      fprintf('%d %.5g\n', x0(i1),  values(i1));
    end
  end
  % Let's get the index of the best column
  % In a case of tie, choose randomly
  value=min(values);
  mini=randpick(find(values==value));
  if (opt.stop && value>=minvalue)
    if opt.display
      fprintf('Stopping because value not decreasing\n');
    end
    break
  end
  nremoved = nremoved + 1;
  rec(nremoved+1).candidates = x0;
  rec(nremoved+1).values = values;
  x0 = setdiff(x0,x0(mini));       % Remove the element from x0
  nx=nx-1;
  rec(nremoved+1).chosen = x0;
  if value<minvalue
    minvalue=value;
    minvaluei=nremoved+1;
  end

  if (opt.display >= 1)
    fprintf(' Value: %.4g\n Chosen: %s\n', value, num2str(x0));
  end

  if (minvalue<opt.stopvalue)
    if opt.display
      fprintf('Objective function value smaller than specified goal\n');
    end
    break
  end
end
x=rec(minvaluei).chosen;

author: gpstuff-dev

repo: gpstuff

save_path: github-repos/MATLAB/gpstuff-dev-gpstuff


text:
function p07_default ( dim_num )

%*****************************************************************************80
%
%% P07_DEFAULT sets default values for problem 07.
%
%  Discussion:
%
%    If a problem uses vector parameters, then the spatial dimension
%    DIM_NUM is needed as input, so that the vector parameters can be
%    properly dimensioned.
%
%    Every problem keeps a count of the number of function calls.  Calling
%    this routine causes that count to be reset to 0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 March 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the problem.
%
  i4 = 0;
  p07_i4 ( 'D', '*', i4 );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function y = beta_pdf(x,a,b)
%BETA_PDF    Beta probability density function (pdf).
%
%   Y = BETA_PDF(X,A,B) Returns the Beta pdf with
%   parameters A and B, at the values in X.
%
%   The size of Y is the common size of the input arguments. A
%   scalar input functions as a constant matrix of the same size as
%   the other inputs.
%
%   Default value for A and B is 1.

% Copyright (c) 2005 Aki Vehtari

% This software is distributed under the GNU General Public
% License (version 3 or later); please refer to the file
% License.txt, included with the software, for details.

if nargin < 3,
  a = 1;
end

if nargin < 2;
  b = 1;
end

if nargin < 1,
  error('Requires at least one input argument.');
end

y=exp((a-1).*log(x) +(b-1).*log(1-x) -betaln(a,b));

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [ y,dzdy ] = pad_data(  I,P,dzdy)
%PAD_DATA Summary of this function goes here
%   Detailed explanation goes here

    if(length(P)==1)
        P=ones(1,4)*P;
    end

    if isempty(dzdy)
        %%forward
        original_size=size(I);
        new_size=original_size;
        new_size(1)=new_size(1)+P(1)+P(2);
        new_size(2)=new_size(2)+P(3)+P(4);
        y=zeros(new_size,'like',I);
        y(1+P(1):P(1)+original_size(1),1+P(3):P(3)+original_size(2),:,:)=I;

    else
        %backward
        new_size=size(dzdy);
        original_size=new_size;
        original_size(1)=original_size(1)-P(1)-P(2);
        original_size(2)=original_size(2)-P(3)-P(4);
        y=I(1+P(1):P(1)+original_size(1),1+P(3):P(3)+original_size(2),:,:);
        dzdy=dzdy(1+P(1):P(1)+original_size(1),1+P(3):P(3)+original_size(2),:,:);

    end

end

author: yechengxi

repo: LightNet

save_path: github-repos/MATLAB/yechengxi-LightNet