text:
if ~exist('n_ima'),
   fprintf(1,'No data to process.\n');
   return;
end;

if n_ima == 0,
    fprintf(1,'No image data available\n');
    return;
end;

if ~exist('active_images'),
	active_images = ones(1,n_ima);
end;
n_act = length(active_images);
if n_act < n_ima,
   active_images = [active_images ones(1,n_ima-n_act)];
else
   if n_act > n_ima,
      active_images = active_images(1:n_ima);
   end;
end;

ind_active = find(active_images);

% I did not call check_active_images, because I want to prevent a break
%check_active_images;

fprintf(1,'\nThis function is useful to select a subset of images to calibrate\n');

   fprintf(1,'\nThere are currently %d active images selected for calibration (out of %d):\n',length(ind_active),n_ima);

   if ~isempty(ind_active),

      if length(ind_active) > 2,

   		for ii = 1:length(ind_active)-2,

         	fprintf(1,'%d, ',ind_active(ii));

      	end;

      	fprintf(1,'%d and %d.',ind_active(end-1),ind_active(end));

      else

         if length(ind_active) == 2,

            fprintf(1,'%d and %d.',ind_active(end-1),ind_active(end));

         else

            fprintf(1,'%d.',ind_active(end));

         end;

      end;

   end;

   fprintf(1,'\n');

   if length(ind_active)==0,
      fprintf(1,'\nYou probably want to add images\n');
      choice = 1;
   else
      if length(ind_active)==n_ima,
         fprintf(1,'\nYou probably want to suppress images\n');
         choice = 0;
      else
         choice = 2;
      end;
   end;

   if (choice~=0) & (choice ~=1),
   	fprintf(1,'\nDo you want to suppress or add images from that list?\n');
   end;

while (choice~=0)&(choice~=1),
   choice = input('For suppressing images enter 0, for adding images enter 1 ([]=no change): ');
   if isempty(choice),
      fprintf(1,'No change applied to the list of active images.\n');
      return;
   end;
   if (choice~=0)&(choice~=1),
      disp('Bad entry. Try again.');
   end;
end;

if choice,

	ima_numbers = input('Number(s) of image(s) to add ([] = all images) = ');

if isempty(ima_numbers),
	   fprintf(1,'All %d images are now active\n',n_ima);
   	ima_proc = 1:n_ima;
	else
   	ima_proc = ima_numbers;
	end;

else

	ima_numbers = input('Number(s) of image(s) to suppress ([] = no image) = ');

	if isempty(ima_numbers),
      fprintf(1,'No image has been suppressed. No modication of the list of active images.\n',n_ima);
   	ima_proc = [];
	else
   	ima_proc = ima_numbers;
	end;

end;

if ~isempty(ima_proc),

   active_images(ima_proc) = choice * ones(1,length(ima_proc));

end;

   check_active_images;

   fprintf(1,'\nThere is now a total of %d active images for calibration:\n',length(ind_active));

   if ~isempty(ind_active),

      if length(ind_active) > 2,

   		for ii = 1:length(ind_active)-2,

         	fprintf(1,'%d, ',ind_active(ii));

      	end;

      	fprintf(1,'%d and %d.',ind_active(end-1),ind_active(end));

      else

         if length(ind_active) == 2,

            fprintf(1,'%d and %d.',ind_active(end-1),ind_active(end));

         else

            fprintf(1,'%d.',ind_active(end));

         end;

      end;

   end;

   fprintf(1,'\n\nYou may now run ''Calibration'' to recalibrate based on this new set of images.\n');

repo: ekfmonoslam


text:
% Change to your downloaded location
clear
addpath('C:\liblinear\matlab')
addpath('../training_code/');
addpath('../utilities/');
addpath('../../data extraction/');
%% load shared definitions and AU data
shared_defs;

% Set up the hyperparameters to be validated
hyperparams.c = 10.^(-7:0.5:1);
hyperparams.e = 10.^(-3);

hyperparams.validate_params = {'c', 'e'};

% Set the training function
svm_train = @svm_train_linear;

% Set the test function (the first output will be used for validation)
svm_test = @svm_test_linear;

pca_loc = '../../pca_generation/generic_face_rigid.mat';

hog_data_dir_BP4D = hog_data_dir;

aus = [1, 2, 4, 6, 7, 10, 12, 14, 15, 17, 23];
%%
for a=1:numel(aus)

    au = aus(a);

    rest_aus = setdiff(all_aus, au);

    % load the training and testing data for the current fold
    [train_samples, train_labels, valid_samples, valid_labels, ~, PC, means, scaling] = Prepare_HOG_AU_data_generic_dynamic(train_recs, devel_recs, au, BP4D_dir, hog_data_dir_BP4D);

    train_samples = sparse(train_samples);
    valid_samples = sparse(valid_samples);

    %% Cross-validate here
    [ best_params, ~ ] = validate_grid_search_no_par(svm_train, svm_test, false, train_samples, train_labels, valid_samples, valid_labels, hyperparams);
    model = svm_train(train_labels, train_samples, best_params);

    [~, predictions_all] = svm_test(valid_labels, valid_samples, model);

    name = sprintf('results_BP4D_devel/AU_%d_dynamic.mat', au);

    [ accuracies, F1s, corrs, ccc, rms, classes ] = evaluate_regression_results( predictions_all, valid_labels );

    save(name, 'model', 'F1s', 'accuracies', 'predictions_all', 'valid_labels');

    % Write out the model
    name = sprintf('models/AU_%d_dynamic.dat', au);

    pos_lbl = model.Label(1);
    neg_lbl = model.Label(2);

    w = model.w(1:end-1)';
    b = model.w(end);

    svs = bsxfun(@times, PC, 1./scaling') * w;

    write_lin_dyn_svm(name, means, svs, b, pos_lbl, neg_lbl);
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function objectiveAbbr = checkObjective(model)
warning('checkObjective will be depreciated because the function name is a misnomer, please use printObjective instead')
% Prints out the Stoichiometric Coefficients for each
% Metabolite, with the name of the objective
%
% USAGE:
%
%    objectiveAbbr = checkObjective(model)
%
% INPUT:
%    model:            COBRA model structure
%
% OUTPUT:
%    objectiveAbbr:    Objective reaction abbreviation
%
% .. Authors:
%       - Ronan Fleming 22/10/2008
%       - Thomas Pfau 15/12/2015 - Made the function compatible with sparse S matrices
%       - Laurent Heirendt March 2017 - Compatibility with large models and conversion to table

objRxnInd = find(model.c ~= 0);
objectiveAbbr = model.rxns(objRxnInd);
T = cell(length(objRxnInd), 1);

if isempty(objRxnInd)
    warning('There is no objective!');
else
    objMetVect = {};
    objRxnVect = {};
    objCoeffVect = {};
    for k = 1:length(objRxnInd)
        objMetInd = find(model.S(:, objRxnInd(k)));
        objMetVect{k} = model.mets(objMetInd);
        rxnName = model.rxns(objRxnInd(k));
        objCoeffVect{k} = full(model.S(objMetInd, objRxnInd(k)));

        objRxnVect = {};
        objRxnIDVect = {};
        % fill the list with the reaction name and ID
        for m = 1:length(objMetInd)
            objRxnVect{m} = char(rxnName);
            objRxnIDVect{m} = objRxnInd(k);
        end

        % save the table for reaction k
        T{k} = table(objCoeffVect{k}, categorical(objMetVect{k}), objMetInd, categorical(objRxnVect'), cell2mat(objRxnIDVect'), ...
                     'VariableNames', {'Coefficient', 'Metabolite', 'metID', 'Reaction', 'RxnID'});

        % concatenate the tables
        if k == 1
            summaryT = T{1};
        else
            summaryT = vertcat(summaryT, T{k});
        end
    end

    % display a summary
    display(summaryT);
end

author: opencobra

repo: cobratoolbox


text:
function grid = reshapeGrid(~, grid)
%RESHAPEGRID   Add the repeated endpoint to a 1D periodic grid.

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

% Note: We use periodic discretizations in 1D/2D/3D. These discretizatioins
% do not include the repeated endpoints. Before plotting data with
% SPINOPERATOR/INITIALIZEMOVIE, we add these repeated endpoints to the grid.

% Get the points:
xx = grid{1};

% Add the endpoint:
xx = [xx; 2*xx(end) - xx(end-1)];

% Output the new grid:
grid{1} = xx;

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
clear;
clc;
close all;

RP = [ -2175464.65976786          4387261.15978363          4072912.71678943]';
SP = [  15987741.24878          3203328.88432069          21122447.2672636]';

this_TOW = 225119;

[lat,  lon, h] = ch_ECEF2LLA(RP);
[az, el] = satellite_az_el(SP, RP);

el = linspace(0.1, pi*0.9, 100);

 dtrop = tropo_correction(el, 1*1000);

plot(dtrop)
title("对流层误差");
ylabel("m");

author: yandld


text:
function tetrahedron_arbq_rule_test02 ( degree, n, header )

%*****************************************************************************80
%
%% TETRAHEDRON_ARBQ_RULE_TEST02 gets a rule and writes it to a file.
%
%  Licensing:
%
%    This code is distributed under the GNU GPL license.
%
%  Modified:
%
%    08 July 2014
%
%  Author:
%
%    Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Hong Xiao, Zydrunas Gimbutas,
%    A numerical algorithm for the construction of efficient quadrature
%    rules in two and higher dimensions,
%    Computers and Mathematics with Applications,
%    Volume 59, 2010, pages 663-676.
%
%  Parameters:
%
%    Input, integer DEGREE, the desired total polynomial degree exactness
%    of the quadrature rule.  0 <= DEGREE <= 15.
%
%    Input, integer N, the number of nodes to be used by the rule.
%
%    Input, string HEADER, an identifier for the filenames.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TETRAHEDRON_ARBQ_RULE_TEST02\n' );
  fprintf ( 1, '  Get a quadrature rule for the tetrahedron.\n' );
  fprintf ( 1, '  Then write it to a file.\n' );
  fprintf ( 1, '  Polynomial exactness degree DEGREE = %d\n', degree );
%
%  Retrieve a symmetric quadrature rule.
%
  [ x, w ] = tetrahedron_arbq ( degree, n );
%
%  Write the points and weights to a file.
%
  rule_filename = strcat ( header, '.txt' );

  rule_unit = fopen ( rule_filename, 'wt' );
  for i = 1 : n
    fprintf ( rule_unit, '%g  %g  %g  %g\n', x(1,i), x(2,i), x(3,i), w(i) );
  end
  fclose ( rule_unit );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Quadrature rule written to file "%s"\n', rule_filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [jj,tx,j,a,v,p,tt]=profile3(t,j,acc,plt)

% function [jj,tx,j,a,v,p,tt]=profile3(t,j,acc)
%
% Calculate symmetrical third order profiles from times:
%
%  Inputs:
%
%      t(1) = constant jerk phase duration
%      t(2) = constant acceleration phase duration (default 0)
%      t(3) = constant velocity phase duration (default 0)
%
%      j    = bound on jerk
%      acc  = continuous time: accuracy for profiles: t(1)*acc = minimal timestep
%             discrete time:   sample time
%
%  Outputs:
%
%      jj  = derivative of jerk profile suitable for simulink
%
%      tx  = time sequence for plotting profiles
%      j   = jerk profile
%      a   = acceleration profile
%      v   = velocity profile
%      p   = position profile
%
%      tt  = 8 switching times for profile:
%
%       0 1              6 7
%       .-.              .-.
%       | |              | |
%       | |  2 3    4 5  | |
%       '-'--.-.----.-.--' '--
%            | |    | |
%            | |    | |
%            '-'    '-'
%
%  Note: coinciding switching times are not removed

%
% Copyright 2004, Paul Lambrechts, The MathWorks, Inc.
%

if nargin < 3 || nargin > 4
    help profile3
    return
end
if nargin==3
    plt=1;
end

if length(t)==1  % min distance with max jerk
    tt=   [0 1 1 2 2 3 3 4 ]*t;

elseif length(t)==2 % constant acceleration phase
    tt=   [0 1 1 2 2 3 3 4 ]*t(1) ...
        + [0 0 1 1 1 1 2 2 ]*t(2);

elseif length(t)==3 % constant velocity phase
    tt=   [0 1 1 2 2 3 3 4 ]*t(1) ...
        + [0 0 1 1 1 1 2 2 ]*t(2) ...
        + [0 0 0 0 1 1 1 1 ]*t(3) ;
else
    return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate Simulink look-up table
jt=[];
for i=1:8
    jt =  [jt   [1 1] * tt(i) ]  ;
end
jt  = [jt 1.5*tt(8)];

jj  = [jt ; [ 0 j j 0 0 -j -j 0 0 -j -j 0 0 j j 0 0 ] ];

if plt==0  % no plot required
    tx=[];j=[];a=[];v=[];p=[];    % dummy outputs
    return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate profiles for plotting

% Determine continuous or discrete
if max(abs( round(t/acc)-t/acc )) > 1e-12 % continuous

   disp('Calculating continuous time profiles')
   step = t(1)*acc;
   tx=0:step:1.2*tt(8);
   x=[];
   for i=0:step:1.2*tt(8)
       j=find(i<=jj(1,:));
       x=[x ; jj(2,j(1))];
   end
   j=x;
   a=cumsum(j)*step;
   v=cumsum(a)*step;
   p=cumsum(v)*step;

else % discrete

   disp('Calculating discrete time profiles')
   Ts=acc;
   ttest=[tt 1.5*tt(8)];
   len = round(1.2*tt(8)/Ts + 1); % length of profiles
   xj = zeros(len,1);
   xa = xj;
   xv = xj;
   xp = xj;
   xj(1) = j;
   tx=0:Ts:1.2*tt(8)+Ts/2;
   for time=Ts:Ts:(1.2*tt(8)+Ts/2)
      i = find( (time + Ts/2) <= ttest ); i = i(1)-1;
      k = round(time/Ts);
      if i==1 || i==7
          xj(k+1) =  j;
      elseif i==3 || i==5
          xj(k+1) = -j;
      else
          xj(k+1) =  0;
      end
      xa(k+1) = xa(k) + xj(k)*Ts;
      xv(k+1) = xv(k) + xa(k)*Ts;
      xp(k+1) = xp(k) + xv(k)*Ts;
   end
   j=xj;a=xa;v=xv;p=xp;

end

%close all
%figure
subplot(411);plot(tx,j,'k','LineWidth',1.5);hold on;plot([0 0],[-1 1]*max(j),'k--',[1 1]*max(tt),[-1 1]*max(j),'k--','LineWidth',1.5);grid on; axis([ [-0.01 1]*max(tx) [-1.1 1.1]*max(j)]);
title('Third order trajectory profiles');ylabel('j [m/s3]');
subplot(412);plot(tx,a,'k','LineWidth',1.5);hold on;plot([0 0],[-1 1]*max(a),'k--',[1 1]*max(tt),[-1 1]*max(a),'k--','LineWidth',1.5);grid on; axis([ [-0.01 1]*max(tx) [-1.1 1.1]*max(a)]);
ylabel('a [m/s2]');
subplot(413);plot(tx,v,'k','LineWidth',1.5);hold on;plot([0 0],[0 1]*max(v),'k--',[1 1]*max(tt),[0 1]*max(v),'k--','LineWidth',1.5);grid on; axis([ [-0.01 1]*max(tx) [-0.1 1.1]*max(v)]);
ylabel('v [m/s]');
subplot(414);plot(tx,p,'k','LineWidth',1.5);hold on;plot([0 0],[0 1]*max(p),'k--',[1 1]*max(tt),[0 1]*max(p),'k--','LineWidth',1.5);grid on; axis([ [-0.01 1]*max(tx) [-0.1 1.1]*max(p)]);
xlabel('time [s]');ylabel('x [m]');
set(1,'position',[700 200 500 680])
set(1,'paperposition',[0 0 5 6.8])

subplot(411);
text(tt(1)+tt(8)/200,max(j)/5,'t_0');
text(tt(2)+tt(8)/200,max(j)/5,'t_1');
text(tt(3)          ,max(j)/5,'t_2');
text(tt(4)          ,max(j)/5,'t_3');
text(tt(5)          ,max(j)/5,'t_4');
text(tt(6)          ,max(j)/5,'t_5');
text(tt(7)+tt(8)/200,max(j)/5,'t_6');
text(tt(8)+tt(8)/200,max(j)/5,'t_7');

author: Sable

repo: mcbench-benchmarks


text:
function res = computeDirectionWeights3d13(varargin)
%COMPUTEDIRECTIONWEIGHTS3D13 Direction weights for 13 directions in 3D
%
%   C = computeDirectionWeights3d13
%   Returns an array of 13-by-1 values, corresponding to directions:
%   C(1)  = [+1  0  0]
%   C(2)  = [ 0 +1  0]
%   C(3)  = [ 0  0 +1]
%   C(4)  = [+1 +1  0]
%   C(5)  = [-1 +1  0]
%   C(6)  = [+1  0 +1]
%   C(7)  = [-1  0 +1]
%   C(8)  = [ 0 +1 +1]
%   C(9)  = [ 0 -1 +1]
%   C(10) = [+1 +1 +1]
%   C(11) = [-1 +1 +1]
%   C(12) = [+1 -1 +1]
%   C(13) = [-1 -1 +1]
%   The sum of the weights in C equals 1.
%   Some values are equal whatever the resolution:
%   C(4)==C(5);
%   C(6)==C(7);
%   C(8)==C(9);
%   C(10)==C(11)==C(12)==C(13);
%
%   C = computeDirectionWeights3d13(DELTA)
%   With DELTA = [DX DY DZ], specifies the resolution of the grid.
%
%   Example
%   c = computeDirectionWeights3d13;
%   sum(c)
%   ans =
%       1.0000
%
%   c = computeDirectionWeights3d13([2.5 2.5 7.5]);
%   sum(c)
%   ans =
%       1.0000
%
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-10-18,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

%% Initializations

% grid resolution
delta = [1 1 1];
if ~isempty(varargin)
    delta = varargin{1};
end

% If resolution is [1 1 1], return the pre-computed set of weights
if all(delta == [1 1 1])
    area1 = 0.04577789120476 * 2;
    area2 = 0.03698062787608 * 2;
    area3 = 0.03519563978232 * 2;
    res = [...
        area1; area1; area1; ...
        area2; area2; area2; area2; area2; area2;...
        area3; area3; area3; area3 ];
    return;
end

% Define points of interest in the 26 discrete directions
% format is pt[Xpos][Ypos][Zpos], with [X], [Y] or [Z] being one of
% 'N' (for negative), 'P' (for Positive) or 'Z' (for Zero)

% points below the OXY plane
ptPNN = normalizeVector3d([+1 -1 -1].*delta);
ptPZN = normalizeVector3d([+1  0 -1].*delta);
ptNPN = normalizeVector3d([-1 +1 -1].*delta);
ptZPN = normalizeVector3d([ 0 +1 -1].*delta);
ptPPN = normalizeVector3d([+1 +1 -1].*delta);

% points belonging to the OXY plane
ptPNZ = normalizeVector3d([+1 -1  0].*delta);
ptPZZ = normalizeVector3d([+1  0  0].*delta);
ptNPZ = normalizeVector3d([-1 +1  0].*delta);
ptZPZ = normalizeVector3d([ 0 +1  0].*delta);
ptPPZ = normalizeVector3d([+1 +1  0].*delta);

% points above the OXY plane
ptNNP = normalizeVector3d([-1 -1 +1].*delta);
ptZNP = normalizeVector3d([ 0 -1 +1].*delta);
ptPNP = normalizeVector3d([+1 -1 +1].*delta);
ptNZP = normalizeVector3d([-1  0 +1].*delta);
ptZZP = normalizeVector3d([ 0  0 +1].*delta);
ptPZP = normalizeVector3d([+1  0 +1].*delta);
ptNPP = normalizeVector3d([-1 +1 +1].*delta);
ptZPP = normalizeVector3d([ 0 +1 +1].*delta);
ptPPP = normalizeVector3d([+1 +1 +1].*delta);

%% Spherical cap type 1, direction [1 0 0]

% Compute area of voronoi cell for a point on the Ox axis, i.e. a point
% in the 6-neighborhood of the center.
refPoint = ptPZZ;

% neighbours of chosen point, sorted by CCW angle
neighbors = [ptPNN; ptPNZ; ptPNP; ptPZP; ptPPP; ptPPZ; ptPPN; ptPZN];

% compute area of spherical polygon
area1 = sphericalVoronoiDomainArea(refPoint, neighbors);

%% Spherical cap type 1, direction [0 1 0]

% Compute area of voronoi cell for a point on the Oy axis, i.e. a point
% in the 6-neighborhood of the center.
refPoint    = ptZPZ;

% neighbours of chosen point, sorted by angle
neighbors   = [ptPPZ; ptPPP; ptZPP; ptNPP; ptNPZ; ptNPN; ptZPN; ptPPN];

% compute area of spherical polygon
area2 = sphericalVoronoiDomainArea(refPoint, neighbors);

%% Spherical cap type 1, direction [0 0 1]

% Compute area of voronoi cell for a point on the Oz axis, i.e. a point
% in the 6-neighborhood of the center.
refPoint = ptZZP;

% neighbours of chosen point, sorted by angle
neighbors = [ptPZP; ptPPP; ptZPP; ptNPP; ptNZP; ptNNP; ptZNP; ptPNP];

% compute area of spherical polygon
area3 = sphericalVoronoiDomainArea(refPoint, neighbors);

%% Spherical cap type 2, direction [1 1 0]

% Compute area of voronoi cell for a point on the Oxy plane, i.e. a point
% in the 18-neighborhood
refPoint = ptPPZ;

% neighbours of chosen point, sorted by angle
neighbors = [ptPZZ; ptPPP; ptZPZ; ptPPN];

% compute area of spherical polygon
area4 = sphericalVoronoiDomainArea(refPoint, neighbors);

%% Spherical cap type 2, direction [1 0 1]

% Compute area of voronoi cell for a point on the Oxz plane, i.e. a point
% in the 18-neighborhood
refPoint = ptPZP;
% neighbours of chosen point, sorted by angle
neighbors = [ptPZZ; ptPPP; ptZZP; ptPNP];

% compute area of spherical polygon
area5 = sphericalVoronoiDomainArea(refPoint, neighbors);

%% Spherical cap type 2, direction [0 1 1]

% Compute area of voronoi cell for a point on the Oxy plane, i.e. a point
% in the 18-neighborhood
refPoint = ptZPP;
% neighbours of chosen point, sorted by angle
neighbors = [ptZPZ; ptNPP; ptZZP; ptPPP];

% compute area of spherical polygon
area6 = sphericalVoronoiDomainArea(refPoint, neighbors);

%% Spherical cap type 3 (all cubic diagonals)

% Compute area of voronoi cell for a point on the Oxyz diagonal, i.e. a
% point in the 26 neighborhood only
refPoint = ptPPP;
% neighbours of chosen point, sorted by angle
neighbors = [ptPZP; ptZZP; ptZPP; ptZPZ; ptPPZ; ptPZZ];

% compute area of spherical polygon
area7 = sphericalVoronoiDomainArea(refPoint, neighbors);

%% Concatenate results

% return computed areas, formatted as fraction of sphere surface
res = [...
    area1 area2 area3 ...
    area4 area4 area5 area5 area6 area6...
    area7 area7 area7 area7...
    ]/(2*pi);

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function hole_point = p00_hole_point ( test, hole_index, m )

%*****************************************************************************80
%
%% P00_HOLE_POINT returns a point inside a given hole.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer TEST, the index of the test problem
%
%    Input, integer HOLE_INDEX, the index of the hole.
%
%    Input, integer M, the spatial dimension.
%
%    Output, real HOLE_POINT(M), a point in the hole
%
  if ( test == 1 )
    hole_point = p01_hole_point ( hole_index, m );
  elseif ( test == 2 )
    hole_point = p02_hole_point ( hole_index, m );
  elseif ( test == 3 )
    hole_point = p03_hole_point ( hole_index, m );
  elseif ( test == 4 )
    hole_point = p04_hole_point ( hole_index, m );
  elseif ( test == 5 )
    hole_point = p05_hole_point ( hole_index, m );
  elseif ( test == 6 )
    hole_point = p06_hole_point ( hole_index, m );
  elseif ( test == 7 )
    hole_point = p07_hole_point ( hole_index, m );
  elseif ( test == 8 )
    hole_point = p08_hole_point ( hole_index, m );
  elseif ( test == 9 )
    hole_point = p09_hole_point ( hole_index, m );
  elseif ( test == 10 )
    hole_point = p10_hole_point ( hole_index, m );
  elseif ( test == 11 )
    hole_point = p11_hole_point ( hole_index, m );
  elseif ( test == 12 )
    hole_point = p12_hole_point ( hole_index, m );
  elseif ( test == 13 )
    hole_point = p13_hole_point ( hole_index, m );
  elseif ( test == 14 )
    hole_point = p14_hole_point ( hole_index, m );
  elseif ( test == 15 )
    hole_point = p15_hole_point ( hole_index, m );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_HOLE_POINT - Fatal error!\n' );
    fprintf ( 1, '  Input value of TEST is out of range.\n' );
    error ( 'P00_HOLE_POINT - Fatal error!' );
  end

  return
end


text:
%--- help for parameter_parsing_tool ---
%
%  INTERNAL FUNCTION
%
%

author: jmaih


text:
%% Copyright (C) 2014-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
%% @defun  vpasolve (@var{e})
%% @defunx vpasolve (@var{e}, @var{x})
%% @defunx vpasolve (@var{e}, @var{x}, @var{x0})
%% Numerical solution of a symbolic equation.
%%
%% Variable-precision numerical solution of the equation @var{e}
%% for variable @var{x} using initial guess of @var{x0}.
%%
%% Example:
%% @example
%% @group
%% syms x
%% eqn = exp(x) == x + 2;
%% vpasolve(eqn, x, 0.1)
%%   @result{} (sym) 1.1461932206205825852370610285214
%% @end group
%% @end example
%%
%% Systems of equations are supported:
%% @example
%% @group
%% syms x y
%% eqns = [y*exp(x) == 16; x^5 == x + x*y^2]
%%   @result{} eqns = (sym 2×1 matrix)
%%
%%       ⎡     x       ⎤
%%       ⎢  y⋅ℯ  = 16  ⎥
%%       ⎢             ⎥
%%       ⎢ 5      2    ⎥
%%       ⎣x  = x⋅y  + x⎦
%% @end group
%%
%% @group
%% vpasolve(eqns, [x; y], [1; 1])
%%   @result{} (sym 2×1 matrix)
%%
%%       ⎡1.7324062670465659633407456995303⎤
%%       ⎢                                 ⎥
%%       ⎣2.8297332667835974266598942031498⎦
%% @end group
%% @end example
%%
%% Complex roots can be found but you must provide a complex initial
%% guess:
%% @example
%% @group
%% vpasolve(x^2 + 2 == 0, x, 1i)
%%   @result{} (sym) 1.4142135623730950488016887242097⋅ⅈ
%% @end group
%% @end example
%%
%% @seealso{vpa}
%% @end defun

function r = vpasolve(e, x, x0)

  if (nargin < 1 || nargin > 3)
    print_usage ();
  end

  if (nargin < 3)
    x0 = sym(0);
  end
  if (nargin < 2)
    x = symvar(e, 1);
  end

  n = digits();

  % everything to column vector
  e = e(:);
  x = x(:);
  x0 = x0(:);

  if (isscalar (x0)) && (numel (x) >= 2)
    x0 = x0 * ones (size (x));
  end

  % IPC cannot pass double matrices yet
  if (isnumeric (x0) && (numel (x0) > 1))
    x0 = num2cell (x0);
  end

  cmd = {
      '(e, x, x0, n) = _ins'
      'import mpmath'
      'mpmath.mp.dps = n'
      'r = nsolve(e, x, x0)'
      'return r' };
  r = pycall_sympy__ (cmd, sym(e), x, x0, n);

end

%!test
%! syms x
%! vpi = vpa(sym(pi), 64);
%! e = tan(x/4) == 1;
%! q = vpasolve(e, x, 3.0);
%! w = q - vpi ;
%! assert (double(w) < 1e-30)

%!test
%! syms x
%! vpi = vpa(sym(pi), 64);
%! e = tan(x/4) == 1;
%! q = vpasolve(e, x);
%! w = q - vpi;
%! assert (double(w) < 1e-30)
%! q = vpasolve(e);
%! w = q - vpi;
%! assert (double(w) < 1e-30)

%!test
%! % very accurate pi
%! syms x
%! e = tan(x/4) == 1;
%! m = digits(256);
%! q = vpasolve(e, x, 3);
%! assert (double(abs(sin(q))) < 1e-256)
%! digits(m);

%!test
%! % very accurate sqrt 2
%! syms x
%! e = x*x == 2;
%! m = digits(256);
%! q = vpasolve(e, x, 1.5);
%! assert (double(abs(q*q - 2)) < 1e-256)
%! digits(m);

%!test
%! % very accurate sqrt pi
%! % (used to fail https://github.com/sympy/sympy/issues/8564)
%! syms x
%! e = x*x == sym(pi);
%! m = digits(256);
%! q = vpasolve(e, x, 3);
%! assert (double(abs(sin(q*q))) < 1e-256)
%! digits(m);

%!test
%! syms x
%! r = vpasolve(x^2 + 2 == 0, x, 1i);
%! assert (double (imag(r)^2 - 2), 0, 1e-32)
%! assert (double (real(r)^2), 0, 1e-32)
%! r = vpasolve(x^2 + 2 == 0, x, -3i + 5);
%! assert (double (imag(r)^2 - 2), 0, 1e-32)
%! assert (double (real(r)^2), 0, 1e-32)

%!test
%! % system
%! syms x y
%! f = 3*x^2 - 2*y^2 - 1;
%! g = x^2 - 2*x + y^2 + 2*y - 8;
%! r = vpasolve([f; g], [x; y], sym([-1; 1]));
%! assert (isa (r, 'sym'))
%! assert (numel (r) == 2)

%!test
%! % system, double guess
%! syms x y
%! f = 3*x^2 - 2*y^2 - 1;
%! g = x^2 - 2*x + y^2 + 2*y - 8;
%! r = vpasolve([f; g], [x; y], [-1.1 1.2]);

%!test
%! % system, double guess
%! syms x y
%! f = 3*x^2 - 2*y^2 - 1;
%! g = x^2 - 2*x + y^2 + 2*y - 8;
%! r1 = vpasolve([f; g], [x; y], [-1.1]);
%! r2 = vpasolve([f; g], [x; y], [-1.1 -1.1]);
%! assert (isequal (r1, r2))

%!test
%! % system, more eqns than unknowns
%! syms x y
%! eqns = [x^3 - x - y == 0; y*exp(x) == 16; log(y) + x == 4*log(sym(2))];
%! r = vpasolve (eqns, [x; y], [1; 1]);
%! A = subs (lhs (eqns), [x; y], r);
%! err = A - [0; 16; 4*log(sym(2))];
%! assert (double (err), zeros (size (err)), 1e-31)

repo: octsympy


text:
% Deterministic error model.
% Demonstrates how the MSE various as the perturbation level increases.
% Produces figures similar to Fig. 4, 5.

%% Configure
clear();

wavelength = 1; % normalized
d_0 = wavelength / 2;
design_set = 'SameAperture'; % can also be 'SameAperture'
designs = get_design_set(design_set, d_0);

n_designs = length(designs);
if strcmpi(design_set, 'SameNumSensor')
    n_doas = 11;
else
    n_doas = 6;
end
doas = linspace(-pi/3, pi/3, n_doas);

n_snapshots = 1000;
source_power = 1;
noise_power = 1;

perturb_type = 'gaussian';
n_params_perturb = 20;
max_perturb = 0.08*d_0;
params_perturb = linspace(0, max_perturb, n_params_perturb);
n_repeats = 1000;

%% Run
mse_em = zeros(n_designs, n_params_perturb, n_repeats);
mse_an = zeros(n_designs, n_params_perturb);

fprintf('Fixed snapshot count, varying pos err std parameters:\n');
fprintf('    DOAs: [%s]\n', num2str(doas, '%f '));
fprintf('    Number of snapshots: %d\n', n_snapshots);
fprintf('    SNR: %.1f\n', 10*log10(source_power / noise_power));
fprintf('    Perturbation range: [0 %f]d_0\n', max_perturb/d_0);
fprintf('    Repeat: %d\n', n_repeats);

for dd = 1:n_designs
    cur_design = designs{dd};
    fprintf('Running simulations for %s:\n', cur_design.name);
    progressbar('reset', n_params_perturb);
    for kk = 1:n_params_perturb
        % collect empirical results
        cur_mses = zeros(n_repeats, 1);
        perturb_std = params_perturb(kk);
        for rr = 1:n_repeats
            pos_err = gen_pos_err(perturb_std, cur_design.element_count, perturb_type);
            perturbed_design = cur_design;
            perturbed_design.position_errors = pos_err;
            [~, R] = snapshot_gen_sto(perturbed_design, doas, wavelength, n_snapshots, noise_power, source_power);
            [Rv, ~, ~] = virtual_ula_cov_1d(cur_design, R, 'SS');
            sp = rmusic_1d(Rv, n_doas, 2*pi * cur_design.element_spacing / wavelength);
            mse_em(dd,kk,rr) = sum((sp.x_est - doas).^2) / n_doas;
        end
        % compute analytical results
        [~, perturb_cov] = gen_pos_err(params_perturb(kk), cur_design.element_count, perturb_type);
        error_stat = struct;
        error_stat.PositionErrorCov = perturb_cov;
        mse_an(dd, kk) = sum(ecov_perturbed_coarray_music_1d(cur_design, wavelength, ...
            doas, source_power, noise_power, n_snapshots, error_stat, 'DiagonalsOnly')) / n_doas;
        progressbar('advance');
    end
    progressbar('end');
end

%% Plot
markers = {'x', 'o', 's', 'd', '*', '^', '>', '<', 'h'};
figure;
for dd = 1:n_designs
    hp = plot(params_perturb/d_0, rad2deg(sqrt(mean(squeeze(mse_em(dd,:,:)), 2))), ...
        ['' markers{mod(dd, length(markers))+1}], ...
        'DisplayName', [designs{dd}.name ' empirical']);
    hold on;
    plot(params_perturb/d_0, rad2deg(sqrt(mse_an(dd,:))), '--', 'Color', get(hp, 'Color'), ...
        'DisplayName', [designs{dd}.name ' analytical']);
end
hold off;
xlabel('\delta_p/d_0');
ylabel('RMSE/deg');
grid on;
legend('show');
legend('Location', 'northwest');
title('RMSE vs Perturbation Level - Deterministic Error Model');

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
% APM Index 2+ DAE Integrator
%
% y = apm_solve(app)
%
% Function apm_solve uploads the model file (apm) and optionally
%   a data file (csv) with the same name to the web-server and performs
%   a forward-time stepping integration of ODE or DAE equations
%   with the following arguments:
%
%  Input:      app = model (apm) and data file (csv) name
%
% Output: y.names  = names of all variables
%         y.values = tables of values corresponding to y.names
%         y.nvar   = number of variables
%         y.x      = combined variables and values but variable
%                      names may be modified to make them valid MATLAB
%                      characters (e.g. replace '[' with '_')
%
function y = apm_solve(app)
    % server and application file names
    server = 'http://xps.apmonitor.com';
    app_model = [app '.apm'];
    app_data = [app '.csv'];

    % clear previous application
    apm(server,app,'clear all');

    % check that model file exists (required)
    if (~exist(app_model,'file')),
        disp(['Error: file ' app_model ' does not exist']);
        y = [];
        return
    else
        % load model file
        apm_load(server,app,app_model);
    end

    % check if data file exists (optional)
    if (exist(app_data,'file')),
        % load data file
        csv_load(server,app,app_data);
    end

    % default options
    % use or don't use web viewer
    web = false;
    if web,
        apm_option(server,app,'nlc.web',2);
    else
        apm_option(server,app,'nlc.web',0);
    end
    % internal nodes in the collocation (between 2 and 6)
    apm_option(server,app,'nlc.nodes',3);
    % sensitivity analysis (default: 0 - off)
    apm_option(server,app,'nlc.sensitivity',0);
    % simulation mode (1=ss,  2=mpu, 3=rto)
    %                 (4=sim, 5=est, 6=nlc, 7=sqs)
    apm_option(server,app,'nlc.imode',7);

    % attempt solution
    solver_output = apm(server,app,'solve');

    % check for successful solution
    status = apm_tag(server,app,'nlc.appstatus');

    if status==1,
        % open web viewer if selected
        if web,
            apm_web(server,app);
        end
        % retrieve solution and solution.csv
        y = apm_sol(server,app);
        return
    else
        apm_web_root(server,app);
        disp(solver_output);
        disp('Error: Did not converge to a solution');
        y = solver_output;
        return
    end

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function [Fao,Fso] = blockframepairaccel(Fa, Fs, Lb, varargin)
%BLOCKFRAMEPAIRACCEL Precompute structures for block processing
%   Usage: F = blockframepairaccel(Fa,Fs,Lb);
%
%   `[Fao,Fso]=blockframepairaccel(Fa,Fs,Lb)` works similar to
%   |blockframeaccel| with a pair of frames. The only difference from
%   calling |blockframeaccel| separatelly for each frame is correct
%   default choice of the slicing windows. Frame objects `Fa,Fs` will be
%   accelerated for length `2*Lb`.
%
%   The following optional arguments are recognized:
%
%      `'anasliwin',anasliwin`  : Analysis slicing window. `sliwin` have to
%                                 be a window of length *2Lb* or a string
%                                 accepted by the |firwin| function. It is
%                                 used only in the slicing window approach.
%                                 The default is `'hann'`.
%
%      `'synsliwin',synsliwin`  : Synthesis slicing window. The same as the
%                                 previous one holds. The default is `'rect'`.
%
%      `'zpad',zpad`   : Number of zero samples the block will be padded
%                        after it is windowed by a slicing window. Note the
%                        frames will be accelerated for length
%                        `2*Lb+2*kv.zpad`.
%

complainif_notenoughargs(nargin,3,'BLOCKFRAMEPAIRACCEL');
complainif_notvalidframeobj(Fa,'BLOCKFRAMEPAIRACCEL');
complainif_notvalidframeobj(Fs,'BLOCKFRAMEPAIRACCEL');

definput.flags.blockalg = {'naive','sliced','segola'};
definput.keyvals.anasliwin = [];
definput.keyvals.synsliwin = [];
definput.keyvals.zpad = 0;
[flags,kv]=ltfatarghelper({},definput,varargin);

isSliProp = ~isempty(kv.anasliwin) || ~isempty(kv.synsliwin) || kv.zpad~=0;
assert(~(~flags.do_sliced && isSliProp),...
   sprintf(['%s: Definig slicing window properties without setting the',...
            ' ''sliced'' flag.'], mfilename));

if flags.do_sliced
   if isempty(kv.anasliwin)
      kv.anasliwin = 'hann';
   end

   if isempty(kv.synsliwin)
      kv.synsliwin = 'hann';
   end

   Fao = blockframeaccel(Fa,Lb,'sliced','sliwin',kv.anasliwin,...
                         'zpad',kv.zpad);
   Fso = blockframeaccel(Fs,Lb,'sliced','sliwin',kv.synsliwin,...
                         'zpad',kv.zpad);
else
   Fao = blockframeaccel(Fa,Lb,flags.blockalg);
   Fso = blockframeaccel(Fs,Lb,flags.blockalg);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function d = gsp_hop_distanz(G,i,j)
%GSP_HOP_DISTANZ Compute the hop distance between two node
%   Usage:  d = gsp_hop_distanz(G,i,j);
%
%   Input parameters:
%       G   : Graph
%       i   : node
%       j   : node
%   Output parameters:
%       d   : hop distanz
%
%   This code computes the hop distance between node i and node j. It uses
%   a naive greedy algorithm and has to be improved.
%

% Author: Nathanael Perraudin
% Date  : 15 septembre 2015
% Testing: test_gsp_hope_distanz

M = double(logical(G.W));
s = zeros(G.N,1);
s(i) = 1;
s = logical(s);
d = 0;
while s(j)==0 && d<=G.N+1
    d = d+1;
    s = logical(M*double(s)) +s;
end

if d == G.N+1
    d = inf;
end

end

repo: gspbox


text:
%  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., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
%function tfrrppat
%TFRRPPAT Unit test for the function TFRRPPAG.

%	O. Lemoine - April 1996.

N=128;

% Reality of the TFR
sig=noisecg(N);
[tfr,rtfr]=tfrrppag(sig);
if sum(any(abs(imag(rtfr))>sqrt(eps)))~=0,
 error('tfrrppag test 1 failed');
end

% Energy conservation
sig=fmlin(N);
[tfr,rtfr]=tfrrppag(sig);
Es=norm(sig)^2;
Etfr=sum(mean(rtfr));
if abs(Es-Etfr)>sqrt(eps),
 error('tfrrppag test 2 failed');
end

% time localization
t0=30; sig=((1:N)'==t0);
[tfr,rtfr]=tfrrppag(sig);
[ik,jk]=find(abs(rtfr)>sqrt(eps));
if any(jk~=t0)|any(ik'-(1:N)),
 error('tfrrppag test 3 failed');
end;

% frequency localization
f0=30;
sig=fmconst(N+6,f0/N);
[tfr rtfr]=tfrrppag(sig,N/2+2,N,tftb_window(N+1,'rect'));
if any(find(rtfr>2*max(rtfr)/N)~=f0+1)|(abs(mean(rtfr)-1.0)>sqrt(eps)),
 error('tfrrppag test 4 failed');
end;

N=131;

% Reality of the TFR
sig=noisecg(N);
[tfr,rtfr]=tfrrppag(sig);
if sum(any(abs(imag(rtfr))>sqrt(eps)))~=0,
 error('tfrrppag test 5 failed');
end

% Energy conservation
sig=fmlin(N);
[tfr,rtfr]=tfrrppag(sig);
Es=norm(sig)^2;
Etfr=sum(mean(rtfr));
if abs(Es-Etfr)>sqrt(eps),
 error('tfrrppag test 6 failed');
end

% time localization
t0=30; sig=((1:N)'==t0);
[tfr,rtfr]=tfrrppag(sig);
[ik,jk]=find(abs(rtfr)>sqrt(eps));
if any(jk~=t0)|any(ik'-(1:N)),
 error('tfrrppag test 7 failed');
end;

% frequency localization
f0=30;
sig=fmconst(N+6,f0/N);
[tfr rtfr]=tfrrppag(sig,round(N/2)+2,N,tftb_window(N,'rect'));
if any(find(rtfr>2*max(rtfr)/N)~=f0+1)|(abs(mean(rtfr)-1.0)>sqrt(eps)),
 error('tfrrppag test 8 failed');
end;

author: HeLiangHIT


text:
function [C,h,Ph,F,Fa,Fc,Eh,Ch,hE,hC,Q] = spm_reml_sc(YY,X,Q,N,hE,hC,V)
% ReML estimation of covariance components from y*y' - proper components
% FORMAT [C,h,Ph,F,Fa,Fc,Eh,Ch,hE,hC,Q] = spm_reml_sc(YY,X,Q,N,[hE,hC,V])
%
% YY  - (m x m) sample covariance matrix Y*Y'  {Y = (m x N) data matrix}
% X   - (m x p) design matrix
% Q   - {1 x q} covariance components
% N   - number of samples
%
% hE  - hyperprior expectation in log-space [default = -32]
% hC  - hyperprior covariance  in log-space [default = 256]
% V   - fixed covariance component
%
% C   - (m x m) estimated errors = h(1)*Q{1} + h(2)*Q{2} + ...
% h   - (q x 1) ReML hyperparameters h
% Ph  - (q x q) conditional precision of log(h)
%
% hE  - prior expectation of log scale parameters
% hC  - prior covariances of log scale parameters
% Eh  - posterior expectation of log scale parameters
% Ch  - posterior covariances of log scale parameters
%
% Q   - scaled covariance components
%
% F   - [-ve] free energy F = log evidence = p(Y|X,Q) = ReML objective
%
% Fa  - accuracy
% Fc  - complexity (F = Fa - Fc)
%
% Performs a Fisher-Scoring ascent on F to find MAP variance parameter
% estimates.  NB: uses weakly informative log-normal hyperpriors.
% See also spm_reml for an unconstrained version that allows for negative
% hyperparameters.
%
%__________________________________________________________________________
%
% SPM ReML routines:
%
%      spm_reml:    no positivity constraints on covariance parameters
%      spm_reml_sc: positivity constraints on covariance parameters
%      spm_sp_reml: for sparse patterns (c.f., ARD)
%
%__________________________________________________________________________
% Copyright (C) 2007-2017 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_reml_sc.m 7305 2018-05-07 13:35:06Z karl $

% assume a single sample if not specified
%--------------------------------------------------------------------------
try, N; catch, N  = 1;  end
try, V; catch, V  = 0;  end

% initialise h
%--------------------------------------------------------------------------
n     = length(Q{1});
m     = length(Q);
h     = zeros(m,1);
dFdh  = zeros(m,1);
dFdhh = zeros(m,m);
Inn   = speye(n,n);

[PQ{1:m}] = deal(zeros(n,n));

% ortho-normalise X
%--------------------------------------------------------------------------
if isempty(X)
    X = sparse(n,0);
    R = Inn;
else
    X = spm_svd(X,0);
    R = Inn - X*X';
end

% check fixed component
%--------------------------------------------------------------------------
if length(V) == 1
    V = V*Inn;
end

% initialise and specify hyperpriors
%==========================================================================

% scale Q and YY
%--------------------------------------------------------------------------
sY = spm_trace(R,YY)/(N*n);
YY = YY/sY;
V  = V/sY;
for i = 1:m
    sh(i,1) = spm_trace(R,Q{i})/n;
    Q{i}    = Q{i}/sh(i);
end

% hyperpriors
%--------------------------------------------------------------------------
try, hE = hE(:);        catch, hE = -32;   end
try, hP = spm_inv(hC);  catch, hP = 1/256; end

% check sise
%--------------------------------------------------------------------------
if length(hE) < m, hE = hE(1)*ones(m,1);   end
if length(hP) < m, hP = hP(1)*speye(m,m);  end

% intialise h: so that sum(exp(h)) = 1
%--------------------------------------------------------------------------
if any(diag(hP) > exp(16))
    h = hE;
end

% ReML (EM/VB)
%--------------------------------------------------------------------------
dF    = Inf;
as    = 1:m;
t     = 4;
for k = 1:32

    % compute current estimate of covariance
    %----------------------------------------------------------------------
    C     = V;
    for i = as
        C = C + Q{i}*exp(h(i));
    end
    iC    = spm_inv(C);

    % E-step: conditional covariance cov(B|y) {Cq}
    %======================================================================
    iCX    = iC*X;
    if ~isempty(X)
        Cq = inv(X'*iCX);
    else
        Cq = sparse(0);
    end

    % M-step: ReML estimate of hyperparameters
    %======================================================================

    % Gradient dF/dh (first derivatives)
    %----------------------------------------------------------------------
    P     = iC  - iCX*Cq*iCX';
    U     = Inn - P*YY/N;
    for i = as

        % dF/dh = -trace(dF/diC*iC*Q{i}*iC)
        %------------------------------------------------------------------
        PQ{i}   = P*Q{i};
        dFdh(i) = -spm_trace(PQ{i},U)*N/2;

    end

    % Expected curvature E{dF/dhh} (second derivatives)
    %----------------------------------------------------------------------
    for i = as
        for j = as

            % dF/dhh = -trace{P*Q{i}*P*Q{j}}
            %--------------------------------------------------------------
            dFdhh(i,j) = -spm_trace(PQ{i},PQ{j})*N/2;
            dFdhh(j,i) =  dFdhh(i,j);

        end
    end

    % modulate
    %----------------------------------------------------------------------
    dFdh  = dFdh.*exp(h);
    dFdhh = dFdhh.*(exp(h)*exp(h)');

    % add hyperpriors
    %----------------------------------------------------------------------
    e     = h     - hE;
    dFdh  = dFdh  - hP*e;
    dFdhh = dFdhh - hP;

    % Fisher scoring: update dh = -inv(ddF/dhh)*dF/dh
    %----------------------------------------------------------------------
    dh    = spm_dx(dFdhh(as,as),dFdh(as),{t});
    h(as) = h(as) + dh;

    % predicted change in F - increase regularisation if increasing
    %----------------------------------------------------------------------
    pF    = dFdh(as)'*dh;
    if pF > dF
        t = t - 1;
    else
        t = t + 1/8;
    end
    dF    = pF;

    % convergence
    %----------------------------------------------------------------------
    fprintf('%s %-23d: %10s%e [%+3.2f]\n','  ReML Iteration',k,'...',full(dF),t);
    if dF < 1e-2
        break
    else
        % eliminate redundant components (automatic selection)
        %------------------------------------------------------------------
        as  = find(h > hE);
        as  = as(:)';
    end
end

% log evidence = ln p(y|X,Q) = ReML objective = F = trace(R'*iC*R*YY)/2 ...
%--------------------------------------------------------------------------
Ph    = -dFdhh;
if nargout > 3

    % tr(hP*inv(Ph)) - nh (complexity KL cost of parameters = 0)
    %----------------------------------------------------------------------
    Ft = trace(hP/Ph) - length(Ph);

    % complexity - KL(Ph,hP)
    %----------------------------------------------------------------------
    Fc = Ft/2 + e'*hP*e/2 + spm_logdet(Ph/hP)/2;

    % Accuracy - ln p(Y|h)
    %----------------------------------------------------------------------
    Fa = Ft/2 - spm_trace(C*P,YY*P)/2 - N*n*log(2*pi)/2  - N*spm_logdet(C)/2;

    % Free-energy
    %----------------------------------------------------------------------
    F  = Fa - Fc - N*n*log(sY)/2;

end

% priors and posteriors of log parameters (with scaling)
%--------------------------------------------------------------------------
if nargout > 7

    hE = hE + log(sY) - log(sh);
    hC = spm_inv(hP);
    Eh = h  + log(sY) - log(sh);
    Ch = spm_inv(Ph);

end

% return exp(h) hyperpriors and rescale
%--------------------------------------------------------------------------
h  = sY*exp(h)./sh;
C  = sY*C;

author: spm


text:
function sydney_add_labels_to_list(source_list_file,dest_list_file,rotation_start_index,all_singlerots_to_zero)

list = textread(source_list_file,'%s');
parfor i = 1 : numel(list)
    [class_name,rotation] = sydney_extract_classname_and_rot_from_filename(list{i});
    [class_label,pose_label] = sydney_generate_class_and_pose_labels(class_name,rotation-rotation_start_index,all_singlerots_to_zero);
    list{i} = sprintf('%s %d %d %d' ,list{i},1,class_label,pose_label);
end

%--- write the modified list to the destination file
fp = fopen(dest_list_file,'wt');
fprintf(fp,'%s\n',list{:});
fclose(fp);

author: lmb-freiburg

repo: orion

save_path: github-repos/MATLAB/lmb-freiburg-orion


text:
function [g, kern] = pathKernDiagGradient(kern, x, covDiag)

% PATHKERNDIAGGRADIENT Compute the gradient of the PATH kernel's diagonal wrt parameters.
% FORMAT
% DESC computes the gradient of functions of the diagonal of the
% path kernel matrix with respect to the parameters of the kernel. The
% parameters' gradients are returned in the order given by the
% pathKernExtractParam command.
% ARG kern : the kernel structure for which the gradients are
% computed.
% ARG x : the input data for which the gradient is being computed.
% ARG factors : partial derivatives of the function of interest with
% respect to the diagonal elements of the kernel.
% RETURN g : gradients of the relevant function with respect to each
% of the parameters. Ordering should match the ordering given in
% pathKernExtractParam.
% RETURN kern : the updated kernel structure
%
% SEEALSO : pathKernParamInit, pathDiagGradient, pathKernExtractParam, pathKernGradient
%
% COPYRIGHT : Andrea Baisero, Carl Henrik Ek, 2013

% SHEFFIELDML

maxl=max(cellfun(@(x)size(x,1),x));
kern=pathKernUpdateWMat(kern,maxl);

num=length(x);

gwd=zeros(1,num);
gwhv=zeros(1,num);
ggk=zeros(num,kern.gkern.nParams);
for i=1:num
    li=size(x{i},1);
    w=kern.wmat(1:li,1:li);
    w=.5*(w+w(end:-1:1,end:-1:1));
    dwd=kern.dwdmat(1:li,1:li);
    dwd=.5*(dwd+dwd(end:-1:1,end:-1:1));
    dwhv=kern.dwhvmat(1:li,1:li);
    dwhv=.5*(dwhv+dwhv(end:-1:1,end:-1:1));

    gk=kernCompute(kern.gkern,x{i});
    gwd(i)=sum(sum(dwd.*gk));
    gwhv(i)=sum(sum(dwhv.*gk));
    ggk(i,:)=kernGradient(kern.gkern,x{i},w);
end

g(1,1)=sum(covDiag.*gwd);
g(1,2)=sum(covDiag.*gwhv);
for i=1:kern.gkern.nParams
    g(1,i+2)=sum(covDiag.*ggk(:,i));
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
classdef MimStorageFormat
    % MimStorageFormat. Describes how a cached file will be stored
    %
    %   When storing data in the cache, this enumeration can be used to
    %   choose the file format
    %
    %    Mat           - Store as a Matlab .mat file
    %    Xml           - Convert to an XML file using a structured format
    %                    that can be loaded back into the original data
    %                    using CoreLoadXml(). Use this format when you wat
    %    SimplifiedXml - Save using a more compact XML format that uses
    %                    less disk space and is easier for people to parse,
    %                    but loading the data from the simplified XML file
    %                    will not necessarily produce exactly the same
    %                    structure that was saved
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    enumeration
        Mat           % MATLAB .MAT file
        Xml           % XML file that can be automatically reloaded
        XmlSimplified % XML file that may require conversion when loading
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
% this is a demo low-performance EPI sequence;
% it doesn't use ramp-samping and is only good for educational purposes.
% in addition, it demonstrated how LABEL extension can be used to set data
% header values, which can be used either in combination with integrated
% image reconstruction or to guide the off-line reconstruction tools
%
seq=mr.Sequence();              % Create a new sequence object
fov=220e-3; Nx=96; Ny=Nx;       % Define FOV and resolution
thickness=3e-3;                 % slice thinckness
Nslices=7;
Nreps = 4;
Navigator = 3;

% Set system limits
lims = mr.opts('MaxGrad',32,'GradUnit','mT/m',...
               'MaxSlew',130,'SlewUnit','T/m/s', ...
               'rfRingdownTime', 30e-6, 'rfDeadTime', 100e-6);

% Create 90 degree slice selection pulse and gradient
[rf, gz] = mr.makeSincPulse(pi/2,'system',lims,'Duration',3e-3,...
    'SliceThickness',thickness,'apodization',0.5,'timeBwProduct',4);

% define the trigger to play out
trig=mr.makeTrigger('physio1','duration', 2000e-6); % duration after

% Define other gradients and ADC events
deltak=1/fov;
kWidth = Nx*deltak;
dwellTime = 4e-6; % I want it to be divisible by 2
readoutTime = Nx*dwellTime;
flatTime=ceil(readoutTime*1e5)*1e-5; % round-up to the gradient raster
gx = mr.makeTrapezoid('x',lims,'Amplitude',kWidth/readoutTime,'FlatTime',flatTime);
adc = mr.makeAdc(Nx,'Duration',readoutTime,'Delay',gx.riseTime+flatTime/2-(readoutTime-dwellTime)/2);

% Pre-phasing gradients
preTime=8e-4;
gxPre = mr.makeTrapezoid('x',lims,'Area',-gx.area/2,'Duration',preTime); % removed -deltak/2 to aligh the echo between the samples
gzReph = mr.makeTrapezoid('z',lims,'Area',-gz.area/2,'Duration',preTime);
gyPre = mr.makeTrapezoid('y',lims,'Area',+Ny/2*deltak,'Duration',preTime); % phase area should be Kmax for clin=0 and -Kmax for clin=Ny... strange

% Phase blip in shortest possible time
dur = ceil(2*sqrt(deltak/lims.maxSlew)/10e-6)*10e-6;
gy = mr.makeTrapezoid('y',lims,'Area',-deltak,'Duration',dur); % phase area should be Kmax for clin=0 and -Kmax for clin=Ny... strange

gz_spoil=mr.makeTrapezoid('z',lims,'Area',deltak*Nx*4);

% Define sequence blocks
for r=1:Nreps
    seq.addBlock(trig, mr.makeLabel('SET','SLC', 0));
    for s=1:Nslices
        rf.freqOffset=gz.amplitude*thickness*(s-1-(Nslices-1)/2);
        rf.phaseOffset=-rf.freqOffset*mr.calcRfCenter(rf); % compensate for the slice-offset induced phase
        seq.addBlock(rf,gz);
        seq.addBlock(gxPre,gzReph, ...
                     mr.makeLabel('SET','NAV',1),...
                     mr.makeLabel('SET','LIN', round(Ny/2)));
        for n=1:Navigator
            seq.addBlock(gx,adc, ...
                         mr.makeLabel('SET','REV', gx.amplitude<0), ...
                         mr.makeLabel('SET','SEG', gx.amplitude<0), ...
                         mr.makeLabel('SET','AVG',n==3));
            if (n~=Navigator)
                seq.addBlock(mr.makeDelay(mr.calcDuration(gy))); % we need this dummy blip pulse to maintain ientical RO gradient timing anf the correspnding eddy currents
            end
            gx.amplitude = -gx.amplitude;   % Reverse polarity of read gradient
        end
        %seq.addBlock(gxPre,gyPre,gzReph);
        seq.addBlock(gyPre, ...
                     mr.makeLabel('SET','LIN', 0), ...
                     mr.makeLabel('SET','NAV', 0), ...
                     mr.makeLabel('SET','AVG', 0) );% lin/nav/avg reset
        for i=1:Ny
            seq.addBlock(mr.makeLabel('SET','REV', gx.amplitude<0), ...
                         mr.makeLabel('SET','SEG', gx.amplitude<0));
            seq.addBlock(gx,adc);   % Read one line of k-space
            seq.addBlock(gy,mr.makeLabel('INC','LIN', 1));  % Phase blip
            gx.amplitude = -gx.amplitude;   % Reverse polarity of read gradient
        end
        seq.addBlock(gz_spoil,mr.makeDelay(0.1),mr.makeLabel('INC','SLC', 1));
        if rem(Navigator+Ny,2)~=0
            gx.amplitude = -gx.amplitude;
        end
    end
    seq.addBlock(mr.makeLabel('INC','REP', 1)); % zero-duration block, but it's OK
end

%% check whether the timing of the sequence is correct
[ok, error_report]=seq.checkTiming;

if (ok)
    fprintf('Timing check passed successfully\n');
else
    fprintf('Timing check failed! Error listing follows:\n');
    fprintf([error_report{:}]);
    fprintf('\n');
end

%% prepare sequence export
seq.setDefinition('FOV', [fov fov thickness*Nslices]);
seq.setDefinition('Name', 'epi_lbl');

seq.write('epi_label.seq');   % Output sequence for scanner

return

%% plots, etc.
seq.plot('TimeRange',[0 0.1], 'TimeDisp', 'ms', 'Label', 'SEG,LIN,SLC');
% trajectory calculation
[ktraj_adc, t_adc, ktraj, t_ktraj, t_excitation, t_refocusing] = seq.calculateKspacePP();

% plot k-spaces
figure; plot(t_ktraj, ktraj'); % plot the entire k-space trajectory
hold; plot(t_adc,ktraj_adc(1,:),'.'); % and sampling points on the kx-axis
figure; plot(ktraj(1,:),ktraj(2,:),'b'); % a 2D plot
axis('equal'); % enforce aspect ratio for the correct trajectory display
hold; plot(ktraj_adc(1,:),ktraj_adc(2,:),'r.');

% seq.sound(); % simulate the seq's tone
% test read-write
%seq.read('epi_label.seq');
%seq.write('epi_label2.seq');

author: pulseq

repo: pulseq

save_path: github-repos/MATLAB/pulseq-pulseq


text:
clear all; close all;
I=imread('circles.png');
J=bwulterode(I);
figure;
subplot(121);  imshow(I);
subplot(122);  imshow(J);

author: HuangCongQing


text:
function [winner, roundnr] = connect4(player, nmvs, posval)
%CONNECT4([player, nmvs, posval])   Connect-4 game played on a 7x7 board
%
% Optional input arguments:
%  player .. e.g. 'HC' human vs. computer (default); other options: 'CH', 'CC', 'HH'
%  nmvs ..  minimum number of moves to be evaluated
%  posval .. weights for 1,2, or 3 same-color pieces in an otherwise empty arrays of 4
%
% The program uses a minimax algorithm with brute-force search to a flexible depth of the search tree (however usually not more than 2-3).
% Evaluation of positions follows the tactic of maximizing the weighted amount of same-color pieces in empty arrays of 4
% It takes a brute-force minimax approach in evaluating all possible move-combinations down to a certain depth
% The actual depth of evaluation corresponds to the depth reached after computing a predefined minimum number of available moves.
% Despite the usage of only one simple tactic the program appears to show some adequate or even intelligent behavior.
% At each move, the figure title gives information on the round-nr, computed depth, nr of computed moves, assessment of current position, and
% evaluation of the best-possible position at computed depth. Positive/negative values indicate an assumed advantage for player 1/2.

% By Mathias Benedek 09/2011

global INDEX4 EVAL_INDEX CNTR NMVS iDPTH movelist

if nargin < 1
    player = 'HC';  %stabdard setting: human vs. computer
else
    if ~any(strcmpi(player(1),{'C','H'})) && any(strcmpi(player(2),{'C','H'}))
        error('Player-argument must be e.g. for human vs computer ''HC''!')
    end
end
if nargin < 2
    NMVS = [2000 2000];     % minimum number of moves to compute; val < 50 corresponds to depth==1, val < 2800 to depth==2, val> 2800 to depth==3
else
    NMVS = nmvs;
end
if nargin < 3
    v = [1 3 10];   % evaluation v for 1x1+3x0, 2x1+2x0, 3x1+1x0, 4x1 (=4 wins)
    w = [1 3 10];
else
    v = posval{1};
    w = posval{2};
end
if strcmp(player,'CC')
    gamov_announced = 1;    %do not preannounce soon gameover in C-C game
else
    gamov_announced = 0;
end
messbxs = 1;

board = zeros(7,7);

%INDEX4: list all possible sets of 4 adjacent positions in board
row1 = [1:7:22; 8:7:29; 15:7:36; 22:7:43];
col1 = [1:4; 2:5; 3:6; 4:7];
p6 = [0 6 12 18];
p8 = [0 8 16 24];
INDEX4 = [row1; row1+1; row1+2; row1+3; row1+4; row1+5; row1+6;...
    col1; col1+7; col1+14; col1+21; col1+28; col1+35; col1+42;...
    4+p6; 5+p6; 11+p6; 6+p6; 12+p6; 18+p6; 7+p6; 13+p6; 19+p6; 25+p6; 14+p6; 20+p6; 26+p6; 21+p6; 27+p6; 28+p6;...
    4+p8; 3+p8; 11+p8; 2+p8; 10+p8; 18+p8; 1+p8; 9+p8; 17+p8; 25+p8; 8+p8; 16+p8; 24+p8; 15+p8; 23+p8; 22+p8];

% Evaluation indices of player 1 and 2
EVAL_INDEX = {};
mxv = 10^4; %maximal value
EVAL_INDEX{1} = [0 v(1) -v(1) v(1) v(2) 0 -v(1) 0 -v(2) v(1) v(2) 0 v(2) v(3) 0 0 0 0 -v(1) 0 -v(2) 0 0 0 -v(2) 0 -v(3) ...
    v(1) v(2) 0 v(2) v(3) 0 0 0 0 v(2) v(3) 0 v(3) mxv zeros(1,13) ...
    -v(1) 0 -v(2) 0 0 0 -v(2) 0 -v(3) zeros(1,9) -v(2) 0 -v(3) 0 0 0 -v(3) 0 -mxv];
% Evaluation of sets [0 0 0 0], [0 0 0 1], [0 0 0 2], [0 0 1 0], ...
EVAL_INDEX{2} = [0 w(1) -w(1) w(1) w(2) 0 -w(1) 0 -w(2) w(1) w(2) 0 w(2) w(3) 0 0 0 0 -w(1) 0 -w(2) 0 0 0 -w(2) 0 -w(3) ...
    w(1) w(2) 0 w(2) w(3) 0 0 0 0 w(2) w(3) 0 w(3) mxv zeros(1,13) ...
    -w(1) 0 -w(2) 0 0 0 -w(2) 0 -w(3) zeros(1,9) -w(2) 0 -w(3) 0 0 0 -w(3) 0 -mxv];

%Plot game
figure('Units','normalized','Position',[.2 .2 .6 .6], 'NumberTitle','off', 'MenuBar','none', 'Name','4-Wins','Color',[.5 .5 .5]); %,'WindowKeyPressFcn',@keypfcn);
axes('Units','normalized','Position',[.1 .1 .8 .8],'Color',[0 0 .8],'Box','on');
set(gca, 'XTick',[],'YTick',[],'XLim',[0,80],'YLim',[0,80])
hold on;

for row = 1:7
    for col = 1:7
        field(row,col) = plot(10+(col-1)*10, 10+(row-1)*10, 'wo','MarkerSize',35,'LineWidth',2,'MarkerFaceColor',[1 1 1]','MarkerEdgeColor', [0 0 0],'ButtonDownFcn',@click_col);
    end
end
col = {[1 0 0], [1 1 0]};
cnames = {'Red','Yellow'};
drawnow;

%Play game
gameover = 0;
move_rc = [1,1];
roundnr = 0;    %Round number
turn = 1;       %player red => turn==1, player yellow => turn==2
movelist = [];
while ~gameover

    if turn == 1, roundnr = roundnr + 1; end
    CNTR = 0;
    move_rc_last = move_rc;
    [board, move_rc, valD] = next_move(board, turn, player(turn));
    movelist(roundnr, turn) = move_rc(2);

    %Show move
    set(field(move_rc(1), move_rc(2)),'MarkerFaceColor',col{turn},'MarkerEdgeColor',[1 1 1])
    set(field(move_rc_last(1), move_rc_last(2)),'MarkerEdgeColor',[0 0 0])
    v = eval_board(board, turn);
    set(gcf,'Name',['Connect-4   (R: ',num2str(roundnr),', D: ',num2str(iDPTH),', #M: ',num2str(CNTR),', C: ',num2str(v),', P: ',num2str(valD),')'])
    drawnow;

    %Announce soon game-over
    cwins = (turn==2 && valD<-10^3) || (turn==1 && valD>10^3);
    if cwins && strcmp(player(turn), 'C')  && ~gamov_announced && messbxs && iDPTH > 1
        msgbox(['Game over in ',num2str(iDPTH-1),' move(s)'])
        gamov_announced = 1;
    end

    %Check if game is over
    if abs(v) > 10^3
        gameover = 1;
        if messbxs, msgbox([cnames{turn},' wins!']); end
        winner = turn;
    elseif ~any(any(board==0))
        gameover = 1;
        if messbxs, msgbox('Patt!'); end
        winner = 0;
    end

    turn = 3 - turn;
end

function [board, move_rc, valD] = next_move(board, turn, player_act)

if strcmp(player_act, 'H')      %Human move
    [board, move_rc] = get_move(board, turn);
    valD = nan;

elseif strcmp(player_act, 'C')  %Computer-Move
    [board, move_rc, valD] = make_move(board, turn);

end

function [board, move_rc] = get_move(board, turn)
global mymove

mymove = -1;

valid_moves =  find(board(1,:) == 0);
isvalidmove = 0;
while ~isvalidmove
    isvalidmove = any(mymove == valid_moves);
    pause(.01)
end

row = find(board(:,mymove) == 0, 1, 'last');
board(row, mymove) = turn;
move_rc = [8-row, mymove];

function [board, move_rc, valD] = make_move(board, turn)
global CNTR NMVS iDPTH

for iDPTH = 1:10 	% usually depth is no more than 3, but in cases of restricted choices it might go deeper until at least minimum number of moves are used

    v = MTree(board, turn, turn, 1, iDPTH*2);   %a depth of 1 always includes 2 steps/moves
    if turn == 2
        [valD, move] = min(v);
        nopts = length(find(v < 10^3));   %how many reasonable options for moves
    else
        [valD, move] = max(v);
        nopts = length(find(v > 10^3));   %how many reasonable options for moves
    end

    %compute again for increased depth if ..
    if CNTR > NMVS(turn) || abs(valD) > 10^3 || nopts == 1    %.. not already more than k comps, not already lost/won, and more than one option
        break;
    end

end

row = find(board(:,move) == 0, 1, 'last');
board(row, move) = turn;
move_rc = [8-row, move];

function v = MTree(board, turn, turnD, curD, maxD)   %actual turn in game, and turn in current depth, current depth, maximum depth
global CNTR

finalD = curD == maxD;
nfreecell = length(find(board == 0));

v = nan(1,7);   %standard value nan represents full column
for move = 1:7

    board0 = board;
    row = find(board0(:,move) == 0, 1, 'last');  %find lowest free row in selected column

    if ~isempty(row) %column not full
        board0(row, move) = turnD;            % place piece of current color
        v(move) = eval_board(board0, turn);   % evaluate position
        CNTR = CNTR + 1;

        %Go one step deeper..
        if  ~finalD && nfreecell > 1 && ~(abs(v(move)) > 10^3)  % .. if not already maximum depth, free cells available, position is not won

            new_turnD = 3-turnD;
            if turnD == 1
                v(move) = min(MTree(board0, turn, new_turnD, curD+1, maxD));
            else
                v(move) = max(MTree(board0, turn, new_turnD, curD+1, maxD));
            end
        end

    end
end

function val = eval_board(board, i)
global INDEX4 EVAL_INDEX

p = board(INDEX4)*[27; 9; 3; 1]+1;  % Each of the 88 arrays of 4 fields is related to a number of 1-81 indicating the actual combination of pieces

val = sum(EVAL_INDEX{i}(p));        % A weighted sum of all arrays yields the evaluation of the given game

function click_col(scr, event)    % Record click on column
global mymove

p = get(gca,'CurrentPoint');
mymove = round(p(1,1)/10);


text:
clear
close all

lambda = 1;
iter = 4;
p = 0.8;
eps = 0.0001;

Img = im2double((imread('flower.png')));
Smoothed = ILS_LNorm(Img, lambda, p, eps, iter);
% Smoothed = ILS_LNorm_GPU(Img, lambda, p, eps, iter);

Diff = Img - Smoothed;
ImgE = Img + 3 * Diff;

figure; imshow(Smoothed)
figure; imshow(ImgE)

%%
Img = im2double(imread('clip_art.jpg'));
lambda = 30;
gamma = 10/255;
iter = 10;

Smoothed = ILS_Welsch(Img, lambda, gamma, iter);
% Smoothed = ILS_Welsch_GPU(Img, lambda, gamma, iter);

figure; imshow(Img)
figure; imshow(Smoothed)

author: wliusjtu


text:
function model = gdetect_dp(pyra, model)
% Compute dynamic programming tables used for finding detections.
%   model = gdetect_dp(pyra, model)
%
%   This function implements the dynamic programming algorithm for
%   computing high-scoring derivations using an Object Detection Grammar.
%   It is assumed that the detection grammar is an Isolated Deformation
%   Grammar and therefore contains only structural schemas and deformation
%   schemas.
%
% Return value
%   model   Object model augmented to store the dynamic programming tables
%
% Arguments
%   pyra    Feature pyramid returned by featpyramid.m
%   model   Object model

% AUTORIGHTS
% -------------------------------------------------------
% Copyright (C) 2009-2012 Ross Girshick
%
% This file is part of the voc-releaseX code
% (http://people.cs.uchicago.edu/~rbg/latent/)
% and is available under the terms of an MIT-like license
% provided in COPYING. Please retain this notice and
% COPYING if you use this file (or a portion of it) in
% your project.
% -------------------------------------------------------

% cache filter response
model = filter_responses(model, pyra);

% compute detection scores
L = model_sort(model);
for s = L
  for r = model.rules{s}
    model = apply_rule(model, r, pyra.pady, pyra.padx);
  end
  model = symbol_score(model, s, pyra);
end

% done
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute score pyramid for symbol s
function model = symbol_score(model, s, pyra)
% model  object model
% s      grammar symbol

% take pointwise max over scores for each rule with s as the lhs
rules = model.rules{s};
score = rules(1).score;

for r = rules(2:end)
  for i = 1:length(r.score)
    score{i} = max(score{i}, r.score{i});
  end
end
model.symbols(s).score = score;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute score pyramid for rule r
function model = apply_rule(model, r, pady, padx)
% model  object model
% r      structural|deformation rule
% pady   number of rows of feature map padding
% padx   number of cols of feature map padding

if r.type == 'S'
  model = apply_structural_rule(model, r, pady, padx);
else
  model = apply_deformation_rule(model, r);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute score pyramid for structural rule r
function model = apply_structural_rule(model, r, pady, padx)
% model  object model
% r      structural rule
% pady   number of rows of feature map padding
% padx   number of cols of feature map padding

% structural rule -> shift and sum scores from rhs symbols
% prepare score for this rule
score      = model.scoretpt;
bias       = model_get_block(model, r.offset) * model.features.bias;
loc_w      = model_get_block(model, r.loc);
loc_f      = loc_feat(model, length(score));
loc_scores = loc_w * loc_f;
for i = 1:length(score)
  score{i}(:) = bias + loc_scores(i);
end

% sum scores from rhs (with appropriate shift and down sample)
for j = 1:length(r.rhs)
  ax = r.anchor{j}(1);
  ay = r.anchor{j}(2);
  ds = r.anchor{j}(3);
  % step size for down sampling
  step = 2^ds;
  % amount of (virtual) padding to halucinate
  virtpady = (step-1)*pady;
  virtpadx = (step-1)*padx;
  % starting points (simulates additional padding at finer scales)
  starty = 1+ay-virtpady;
  startx = 1+ax-virtpadx;
  % score table to shift and down sample
  s = model.symbols(r.rhs(j)).score;
  for i = 1:length(s)
    level = i - model.interval*ds;
    if level >= 1
      % ending points
      endy = min(size(s{level},1), starty+step*(size(score{i},1)-1));
      endx = min(size(s{level},2), startx+step*(size(score{i},2)-1));
      % y sample points
      iy = starty:step:endy;
      oy = sum(iy < 1);
      iy = iy(iy >= 1);
      % x sample points
      ix = startx:step:endx;
      ox = sum(ix < 1);
      ix = ix(ix >= 1);
      % sample scores
      sp = s{level}(iy, ix);
      sz = size(sp);
      % sum with correct offset
      stmp = -inf(size(score{i}));
      stmp(oy+1:oy+sz(1), ox+1:ox+sz(2)) = sp;
      score{i} = score{i} + stmp;
    else
      score{i}(:) = -inf;
    end
  end
end
model.rules{r.lhs}(r.i).score = score;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute score pyramid for deformation rule r
function model = apply_deformation_rule(model, r)
% model  object model
% r      deformation rule

% deformation rule -> apply distance transform
def_w      = model_get_block(model, r.def);
score      = model.symbols(r.rhs(1)).score;
bias       = model_get_block(model, r.offset) * model.features.bias;
loc_w      = model_get_block(model, r.loc);
loc_f      = loc_feat(model, length(score));
loc_scores = loc_w * loc_f;
for i = 1:length(score)
  score{i} = score{i} + bias + loc_scores(i);
  % Bounded distance transform with +/- 4 HOG cells (9x9 window)
  [score{i}, Ix{i}, Iy{i}] = fast_bounded_dt(score{i}, def_w(1), def_w(2), ...
                                             def_w(3), def_w(4), 4);
  % Unbounded distance transform
  %[score{i}, Ix{i}, Iy{i}] = dt(score{i}, def_w(1), def_w(2), ...
  %                              def_w(3), def_w(4));
end
model.rules{r.lhs}(r.i).score = score;
model.rules{r.lhs}(r.i).Ix    = Ix;
model.rules{r.lhs}(r.i).Iy    = Iy;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute all filter responses (filter score pyramids)
function model = filter_responses(model, pyra)
% model    object model
% pyra     feature pyramid

% gather filters for computing match quality responses
filters = cell(model.numfilters, 1);
for i = 1:model.numfilters
  filters{i} = single(model_get_block(model, model.filters(i)));
end

for level = 1:length(pyra.feat)
  if ~pyra.valid_levels(level)
    % not processing this level, so set default values
    model.scoretpt{level} = 0;
    for i = 1:model.numfilters
      model.symbols(model.filters(i).symbol).score{level} = -inf;
    end
    continue;
  end

  % Determine the convolution routine to use
  % Try to use an unrolled SSE routine if available
  % 4 floats are processed at a time, so round the feature dim
  % up to the next largest multiple of 4
  fconv_num = 4*ceil(size(pyra.feat{level},3)/4);
  fconv_name = ['fconv_' num2str(fconv_num)];
  if exist(fconv_name) == 3  % 3 ==> MEX function
    fconv_fun = str2func(fconv_name);
  else
    % Fall back to the slower version that works with any dimension
    fconv_fun = @fconv_var_dim;
  end
  % compute filter response for all filters at this level
  r = fconv_fun(pyra.feat{level}, filters, 1, length(filters));

  % find max response array size for this level
  s = cellfun(@size, r, 'UniformOutput', false);
  s = max(cat(1, s{:}), [], 1);

  % set filter response as the score for each filter terminal
  for i = 1:length(r)
    % normalize response array size so all responses at this
    % level have the same dimension
    spady = s(1) - size(r{i},1);
    spadx = s(2) - size(r{i},2);
    if spady > 0 || spadx > 0
      r{i} = post_pad(r{i}, spady, spadx, -inf);
      %r{i} = padarray(r{i}, [spady spady], -inf, 'both');
    end
    fsym = model.filters(i).symbol;
    model.symbols(fsym).score{level} = r{i};
  end
  model.scoretpt{level} = zeros(s);
end

author: rbgirshick

repo: voc-dpm

save_path: github-repos/MATLAB/rbgirshick-voc-dpm


text:
function variance = semicircular_variance ( a, b )

%*****************************************************************************80
%
%% SEMICIRCULAR_VARIANCE returns the variance of the Semicircular PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < B.
%
%    Output, real VARIANCE, the variance of the PDF.
%
  variance = b * b / 4.0;

  return
end


text:
function i4vec_direct_product_test ( )

%*****************************************************************************80
%
%% I4VEC_DIRECT_PRODUCT_TEST tests I4VEC_DIRECT_PRODUCT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 April 2009
%
%  Author:
%
%    John Burkardt
%
  factor_num = 3;
  point_num = 24;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'I4VEC_DIRECT_PRODUCT_TEST\n' );
  fprintf ( 1, '  I4VEC_DIRECT_PRODUCT forms the entries of a\n' );
  fprintf ( 1, '  direct product of a given number of I4VEC factors.\n' );

  x(1:factor_num,1:point_num) = 0;

  for factor_index = 1 : factor_num

    if ( factor_index == 1 )
      factor_order = 4;
      factor_value = [ 1, 2, 3, 4 ];
    elseif ( factor_index == 2 )
      factor_order = 3;
      factor_value = [ 50, 60, 70 ];
    elseif ( factor_index == 3 )
      factor_order = 2;
      factor_value = [ 800, 900 ];
    end

    x = i4vec_direct_product ( factor_index, factor_order, factor_value,  ...
      factor_num, point_num, x );

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '     J     X(1)  X(2)  X(3)\n' );
  fprintf ( 1, '\n' );

  for j = 1 : point_num
    fprintf ( 1, '  %4d    %4d  %4d  %4d\n', j, x(1:factor_num,j) );
  end

  return
end


text:
function [ xmin, ymin ] = minquad ( x1, y1, x2, y2, x3, y3 )

%*****************************************************************************80
%
%% MINQUAD finds a local minimum of F(X) = A * X * X + B * X + C.
%
%  Discussion:
%
%    MINQUAD is primarily intended as a utility routine.
%    The square of the distance function between a point
%    and a line segment has the form of F(X).  Hence, we can seek
%    the line on the second segment which minimizes the square of
%    the distance to the other line segment.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X1, Y1, X2, Y2, X3, Y3, three sets of data
%    of the form ( X, F(X) ).  The three X values must be distinct.
%
%    Output, real XMIN, YMIN.  XMIN is a point within the interval
%    spanned by X1, X2 and X3, at which F takes its local minimum value YMIN.
%

%
%  Refuse to deal with coincident data.
%
  if ( x1 == x2 | x2 == x3 | x3 == x1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MINQUAD - Fatal error!\n' );
    fprintf ( 1, '  X values are equal.\n' );
    error ( 'MINQUAD - Fatal error!' );
  end
%
%  Find the interval endpoints.
%
  xleft = min ( x1, min ( x2, x3 ) );
  xrite = max ( x1, max ( x2, x3 ) );
%
%  Find the minimizer and its function value, over the three input points.
%
  if ( y1 <= y2 & y1 <= y3 )
    xmin = x1;
    ymin = y1;
  elseif ( y2 <= y1 & y2 <= y3 )
    xmin = x2;
    ymin = y2;
  else
    xmin = x3;
    ymin = y3;
  end
%
%  Find the minimizer and its function value over the real line.
%
  [ x, y, ierror ] = parabola_ex ( x1, y1, x2, y2, x3, y3 );
%
%  If F is linear, then take the already computed min.
%
  if ( ierror == 2 )
%
%  If F has a maximum, then take the already computed min.
%
  elseif ( ymin < y )
%
%  If the minimizer is to the left, take the already computed min.
%
  elseif ( x < xleft )
%
%  If the minimizer is to the right, take the already computed min.
%
  elseif ( xrite < x )

  else

    xmin = x;
    ymin = y;

  end

  return
end


text:
function r8vec_concatenate_test ( )

%*****************************************************************************80
%
%% R8VEC_CONCATENATE_TEST tests R8VEC_CONCATENATE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 October 2014
%
%  Author:
%
%    John Burkardt
%
  n1 = 5;
  n2 = 3;
  n3 = n1 + n2;

  a1 = [ 91.1, 31.2, 71.3, 51.4, 31.5 ];
  a2 = [ 42.6, 22.7, 12.8 ];

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_CONCATENATE_TEST\n' );
  fprintf ( 1, '  R8VEC_CONCATENATE concatenates two R8VECs\n' );

  r8vec_print ( n1, a1, '  Array 1:' );
  r8vec_print ( n2, a2, '  Array 2:' );
  a3 = r8vec_concatenate ( n1, a1, n2, a2 );
  r8vec_print ( n3, a3, '  Array 3 = Array 1 + Array 2:' );

  return
end


text:
function value = alnfac ( n )

%*****************************************************************************80
%
%% ALNFAC computes the logarithm of the factorial of N.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the argument of the factorial.
%
%    Output, real VALUE, the logarithm of the factorial of N.
%
  value = alngam ( n + 1 );

  return
end


text:
function test_validity_of_datasource(ds)

% A helper function that tests whether a given data source returns any of the same data
%  in the training and test data.  Obviously if some of the same data is in the training
%  and test sets then there are major problems with the datasource!
%
% The input of this function should be a datasource that has been initialized with
%  real data.  This function replaces the real data with unique numbers and tests
%  whether any of the same numbers are in the training and test sets.  The function
%  prints a message that says whether the datasource has problems or seems ok.

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

%     This code is part of the Neural Decoding Toolbox.
%     Copyright (C) 2011 by Ethan Meyers (emeyers@mit.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 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

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

% replace the original data with unique numbers for data point

the_original_data = ds.the_data;

start_ind = 1;

for iSite = 1:length(the_original_data)

    end_ind = start_ind - 1 + prod(size(the_original_data{iSite}));

    the_new_data{iSite} = reshape(start_ind:end_ind, [size(the_original_data{iSite}, 1) size(the_original_data{iSite}, 2)]);

    start_ind = end_ind + 1;

end

clear the_original_data

ds.the_data = the_new_data;

% get the training and test data with the new data that has unique labels at each time point, and
%   make sure that no unique label is in both the training and test sets

[XTr_all_time_cv YTr_all_cv XTe_all_time_cv YTe_all_cv] = get_data(ds);

num_probs_with_ds = 0;

for iTime = 1:size(XTr_all_time_cv, 2)
    for iCV = 1:ds.num_cv_splits

        if length(intersect(XTr_all_time_cv{iTime}{iCV}(:), XTe_all_time_cv{iTime}{iCV}(:)))
            num_probs_with_ds = num_probs_with_ds + 1;
        end

    end
end

if num_probs_with_ds > 1
    'This datasource has problems - some of the same data is in the training and test sets!'
else
    'This datasource seems ok'
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
function [ distMoment ] = CalcNormalDistributionMoments( paramMu, paramSigmaSquared, momentOrder )
% See https://en.wikipedia.org/wiki/Normal_distribution#Moments.

switch(momentOrder)
    case(1)
        distMoment = (paramMu ^ momentOrder);
    case(2)
        distMoment = (paramMu ^ momentOrder) + paramSigmaSquared;
    case(3)
        distMoment = (paramMu ^ momentOrder) + (3 * paramMu * paramSigmaSquared);
    case(4)
        distMoment = (paramMu ^ momentOrder) + (6 * (paramMu ^ 2) * paramSigmaSquared) + (3 * (paramSigmaSquared ^ 2));
    case(5)
        distMoment = (paramMu ^ momentOrder) + (10 * (paramMu ^ 3) * paramSigmaSquared) + (15 * paramMu * (paramSigmaSquared ^ 2));
    case(6)
        distMoment = (paramMu ^ momentOrder) + (15 * (paramMu ^ 4) * paramSigmaSquared) + (45 * (paramMu ^ 2) * (paramSigmaSquared ^ 2)) + (15 * (paramSigmaSquared ^ 3));
    case(7)
        distMoment = (paramMu ^ momentOrder) + (21 * (paramMu ^ 5) * paramSigmaSquared) + (105 * (paramMu ^ 3) * (paramSigmaSquared ^ 2)) + (105 * paramMu * (paramSigmaSquared ^ 3));
    case(8)
        distMoment = (paramMu ^ momentOrder) + (28 * (paramMu ^ 6) * paramSigmaSquared) + (210 * (paramMu ^ 4) * (paramSigmaSquared ^ 2)) + (420 * (paramMu ^ 2) * (paramSigmaSquared ^ 3)) + (105 * (paramSigmaSquared ^ 4));
    otherwise
        gridRadius      = 10 * sqrt(paramSigmaSquared);
        gridNaumSamles  = 1e6;
        vX = linspace(paramMu - gridRadius, paramMu + gridRadius, gridNaumSamles);
        dX = mean(diff(vX));
        vNormalPdf = (1 / sqrt(2 * pi * paramSigmaSquared)) * exp(-( (vX - paramMu) .^ 2 ) ./ (2 * paramSigmaSquared));
        distMoment = sum( (vX .^ momentOrder) .* vNormalPdf * dX);
end

end

author: RoyiAvital

repo: StackExchangeCodes


text:
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

clear workspace

visual_odo_initial( ...
    'ODO_IMG_TRANS_Y_RANGE', 1:270, ...
    'ODO_IMG_TRANS_X_RANGE', 1:480, ...
    'ODO_IMG_HEIGHT_V_Y_RANGE', 1:270, ...
    'ODO_IMG_HEIGHT_V_X_RANGE', 1:480, ...
    'ODO_IMG_YAW_ROT_Y_RANGE', 1:270, ...
    'ODO_IMG_YAW_ROT_X_RANGE', 1:480, ...
    'ODO_IMG_TRANS_RESIZE_RANGE', [130, 240], ...
    'ODO_IMG_YAW_ROT_RESIZE_RANGE', [130, 240], ...
    'ODO_IMG_HEIGHT_V_RESIZE_RANGE', [130, 240], ...
    'ODO_TRANS_V_SCALE', 30, ...
    'ODO_YAW_ROT_V_SCALE', 1, ...
    'ODO_HEIGHT_V_SCALE', 5, ...
    'MAX_TRANS_V_THRESHOLD', 0.4, ...
    'MAX_YAW_ROT_V_THRESHOLD', 4.2, ...
    'MAX_HEIGHT_V_THRESHOLD', 0.4, ...
    'ODO_SHIFT_MATCH_HORI', 30, ...
    'ODO_SHIFT_MATCH_VERT', 30, ...
    'FOV_HORI_DEGREE', 75, ...
    'FOV_VERT_DEGREE', 20, ...
    'KEY_POINT_SET', [3750, 4700, 8193, 9210], ...
    'ODO_STEP', 5);

% vt initial
vt_image_initial('*.png', ...
    'VT_MATCH_THRESHOLD', 3.4, ...
    'VT_IMG_CROP_Y_RANGE', 1:270, ...
    'VT_IMG_CROP_X_RANGE', 1:480, ...
    'VT_IMG_RESIZE_X_RANGE', 48, ...
    'VT_IMG_RESIZE_Y_RANGE', 27, ...
    'VT_IMG_X_SHIFT', 5, ...
    'VT_IMG_Y_SHIFT', 3, ...
    'VT_GLOBAL_DECAY', 0.5, ...
    'VT_ACTIVE_DECAY', 0.5, ...
    'PATCH_SIZE_Y_K', 5, ...
    'PATCH_SIZE_X_K', 5, ...
    'VT_PANORAMIC', 0, ...
    'VT_STEP', 1);

% yaw_height_hdc initial
yaw_height_hdc_initial( ...
    'YAW_HEIGHT_HDC_Y_DIM', 36, ...
    'YAW_HEIGHT_HDC_H_DIM', 36, ...
    'YAW_HEIGHT_HDC_EXCIT_Y_DIM', 8, ...
    'YAW_HEIGHT_HDC_EXCIT_H_DIM', 8, ...
    'YAW_HEIGHT_HDC_INHIB_Y_DIM', 5, ...
    'YAW_HEIGHT_HDC_INHIB_H_DIM', 5, ...
    'YAW_HEIGHT_HDC_EXCIT_Y_VAR', 1.9, ...
    'YAW_HEIGHT_HDC_EXCIT_H_VAR', 1.9, ...
    'YAW_HEIGHT_HDC_INHIB_Y_VAR', 3.1, ...
    'YAW_HEIGHT_HDC_INHIB_H_VAR', 3.1, ...
    'YAW_HEIGHT_HDC_GLOBAL_INHIB', 0.0002, ...
    'YAW_HEIGHT_HDC_VT_INJECT_ENERGY', 0.1, ...
    'YAW_ROT_V_SCALE', 1, ...
    'HEIGHT_V_SCALE', 1, ...
    'YAW_HEIGHT_HDC_PACKET_SIZE', 5);

% 3d gc initial
gc_initial( ...
    'GC_X_DIM', 36, ...
    'GC_Y_DIM', 36, ...
    'GC_Z_DIM', 36, ...
    'GC_EXCIT_X_DIM', 7, ...
    'GC_EXCIT_Y_DIM', 7, ...
    'GC_EXCIT_Z_DIM', 7, ...
    'GC_INHIB_X_DIM', 5, ...
    'GC_INHIB_Y_DIM', 5, ...
    'GC_INHIB_Z_DIM', 5, ...
    'GC_EXCIT_X_VAR', 1.5, ...
    'GC_EXCIT_Y_VAR', 1.5, ...
    'GC_EXCIT_Z_VAR', 1.5, ...
    'GC_INHIB_X_VAR', 2, ...
    'GC_INHIB_Y_VAR', 2, ...
    'GC_INHIB_Z_VAR', 2, ...
    'GC_GLOBAL_INHIB', 0.0002, ...
    'GC_VT_INJECT_ENERGY', 0.1, ...
    'GC_HORI_TRANS_V_SCALE', 1, ...
    'GC_VERT_TRANS_V_SCALE', 1, ...
    'GC_PACKET_SIZE', 4);

% 3d em initial 1/2/5/
exp_initial( ...
    'DELTA_EXP_GC_HDC_THRESHOLD', 30, ...
    'EXP_LOOPS', 5, ...
    'EXP_CORRECTION', 0.5 );

% main process
% file path: visualDataFile, expMapHistoryFile, odoMapHistoryFile, groundTruthFile,vtHistoryFile, emHistoryFile, varargin

imageFolderPath = 'C:\NeuroSLAM_Datasets\01_NeuroSLAM_Datasets\03_QUTCarparkData';
SynPerDataFile = fullfile (imageFolderPath, '');

gtDataFolderPath = '';
gtDataFile = fullfile(gtDataFolderPath, '');

experiDataPath = 'C:\NeuroSLAM_Datasets\03_NeuroSLAM_Experiments_Results\QUTCarparData';
expMapFile  = fullfile(experiDataPath, '01_exp_map_ml.txt');
odoMapFile  = fullfile(experiDataPath, '02_odo_map_ml.txt');
emHisFile   = fullfile(experiDataPath, '03_em_history_ml.txt');
vtHisFile   = fullfile(experiDataPath, '04_vt_history_ml.txt');
gcTrajFile  = fullfile(experiDataPath, '05_gc_trajectory.txt');
hdcTrajFile = fullfile(experiDataPath, '06_hdc_trajectory.txt');

main(SynPerDataFile, ...
    gtDataFile, ...
    expMapFile, ...
    odoMapFile, ...
    emHisFile, ...
    vtHisFile, ...
    gcTrajFile, ...
    hdcTrajFile, ...
    'BLOCK_READ', 50, ...
    'RENDER_RATE', 2, ...
    'GT_ODO_X_SCALING', 20, ...
    'GT_ODO_Y_SCALING', 20, ...
    'GT_ODO_Z_SCALING', 20, ...
    'GT_EXP_X_SCALING', 20, ...
    'GT_EXP_Y_SCALING', 20, ...
    'GT_EXP_Z_SCALING', 20, ...
    'ODO_MAP_X_SCALING', -0.8, ...
    'ODO_MAP_Y_SCALING', 0.8, ...
    'ODO_MAP_Z_SCALING', 0.1, ...
    'EXP_MAP_X_SCALING', 0.8, ...
    'EXP_MAP_Y_SCALING', 0.8, ...
    'EXP_MAP_Z_SCALING', 0.1);

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


text:
function polarPlot(Z,arg1)
% polarPlot(Z,params)
%   creates a polar plot using the parameters
%   determined by the structure 'params'.
%
% PolarPlot(Z,'PropertyName1',PropertyValue1,'PropertyName2',... );
%   creates a polar plot using the property values and names
%
%   default fields of params or properties:
%      params.grid            = 'on';
%      params.line            = 'on';
%      params.gridColor       = [0.6,0.6,0.6];
%      params.gridLineWidth   = 1;
%      params.fontSize        = 12;
%      params.symbol          = 'os^dp';
%      params.size            = 20;
%      params.color           = 'k';
%      params.fillColor       = [1,1,1;0.6,0.6,0.6;0,0,0];
%      params.lineWidth       = 2;
%      params.backgroundColor = 'w';
%      params.maxAmp          = ceil(max(abs(Z(:)))*10)/10;
%      params.ringTicks       = [0:0.1:defParams.maxAmp];
%	   params.sigFigs		  = 1;
%
%  fillColor loops through the COLUMNS of Z
%  symbols loop through the ROWS of Z
%
%  Example:
%  [x,y] = meshgrid(linspace(0,2*pi,11),linspace(0,2*pi,11));
%  Z = x.*exp(sqrt(-1)*(y+x/5));
%  polarPlot(Z,'fillColor',hsv(11),'grid','off','line','off','size',30)

% 4/9/98 gmb wrote it.
% 3/20/09 ras updated many years later; more flexibility about specifying
% the maximum amplitude of the data, small cleanup using newer matlab
% features.
% 9/17/2016 RL minor clean up to unused input variables, add option to not
% print tick labels (when making figures for paper, often cleaner to add own labels)

%set up default parameters
defParams.grid            = 'on';
defParams.line            = 'on';
defParams.gridColor       = [0.6,0.6,0.6];
defParams.gridLineWidth   = 1;
defParams.fontSize        = 12;
defParams.symbol          = 'os^dp';
defParams.size            = 20;
defParams.color           = 'k';
defParams.fillColor       = [1,1,1;0.6,0.6,0.6;0,0,0];
defParams.lineWidth       = 2;
defParams.backgroundColor = 'w';
defParams.maxAmp          = ceil(max(abs(Z(:)))*10)/10;
defParams.ringTicks       = [0:0.1:defParams.maxAmp];
defParams.sigFigs		  = 1;
defParams.tickLabel       = true;

% also allow the user to override the default params with whatever fields
% are provided
if ~exist('arg1', 'var')
	params = struct;  % empty structure
else
	params = arg1;
end

% If a string is passed as second argument it must be a property name
% so build the params file from the arguments
if strcmp(class(params),'char')
	for arg=1:2:nargin-1
		estr =['propertyName = arg',int2str(arg),';'];   eval(estr)
		estr = ['params.',propertyName,' = arg',int2str(arg+1),';'];
		eval(estr)
	end
end

% for any parameters not specified by the user, use the default value.
for f = fieldnames(defParams)'
	if ~isfield(params, f{1}) | isempty(params.(f{1}))
		params.(f{1}) = defParams.(f{1});
	end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Ready to plot

set(gcf,'Color',params.backgroundColor);

hold on
if strcmp(params.grid,'on')
	%% radial lines
	maxRad = max(params.ringTicks);
	for i=linspace(0,pi,5);
		line([-maxRad*cos(i),maxRad*cos(i)],[-maxRad*sin(i),maxRad*sin(i)],...
			'Color',params.gridColor,'lineWidth',params.gridLineWidth);
	end

	%% concentric circles.
	arg = linspace(0,pi*2,64);
	ringTicks = setdiff(params.ringTicks, 0); % we don't want a ring at r=0
	for i=ringTicks
		x = i*cos(arg);
		y = i*sin(arg);
		plot(x,y,'Color',params.gridColor,'lineWidth',params.gridLineWidth);
	end

	%% labels
	% offset from grid

    % check that we want labels to begin with
    if params.tickLabel
        if isempty(params.maxAmp), dx = 0.075;
        else,	dx = 0.075 * params.maxAmp;
        end

        % make labels
        for i=params.ringTicks(2:length(params.ringTicks));
            pattern = sprintf('%%3.%if', params.sigFigs);
            text(dx, i+dx, sprintf(pattern,i), 'Color', params.gridColor, ...
                'FontSize', params.fontSize);
        end
    end
end

% %hack:  place four white points in the corners to fix image size
%
% plot(params.maxAmp*exp(sqrt(-1)*[45,135,225,315]*pi/180),'w.')

%loop through the columns of Z

%lines first
if strcmp(params.line,'on')
	for i=1:size(Z,2)
		color = params.color(mod(i-1,length(params.color))+1);
		for j=1:size(Z,1)
			plot([0,real(Z(j,i))],[0,imag(Z(j,i))],color,...
				'LineWidth',params.lineWidth)
		end
	end
end

%then symbols
for i=1:size(Z,1)
	symbol = params.symbol(mod(i-1,length(params.symbol))+1);
	color = params.color(mod(i-1,length(params.color))+1);
	for j=1:size(Z,2)
		fillColor = params.fillColor(mod(j-1,size(params.fillColor,1))+1,:);

		plot(real(Z(i,j)),imag(Z(i,j)),[symbol,color],...
			'LineWidth',params.lineWidth,...
			'MarkerFaceColor',fillColor,...
			'MarkerSize',params.size)
	end
end

axis off
axis equal

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function value = r4_chu ( a, b, x )

%*****************************************************************************80
%
%% R4_CHU evaluates the confluent hypergeometric function of R4 arguments.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 October 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real A, B, the parameters.
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the function value.
%
  persistent eps

  if ( isempty ( eps ) )
    eps = r4_mach ( 3 );
  end

  if ( x < 0.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_CHU - Fatal error!\n' );
    fprintf ( 1, '  X < 0.\n' );
    error ( 'R4_CHU - Fatal error!' )
  end

  if ( x == 0.0 )
    if ( 1.0 <= b )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R4_CHU - Fatal error!\n' );
      fprintf ( 1, '  X = 0 and 1 <= B.\n' );
      error ( 'R4_CHU - Fatal error!' )
    end
    value = r4_gamma ( 1.0 - b ) / r4_gamma ( 1.0 + a - b );
    return
  end

  if ( max ( abs ( a ), 1.0 ) * max ( abs ( 1.0 + a - b ), 1.0 ) < 0.99 * abs ( x ) )
    value = x^( - a ) * r4_chu_scaled ( a, b, x );
    return
  end
%
%  The ascending series will be used, because the descending rational
%  approximation (which is based on the asymptotic series) is unstable.
%
  if ( b < 0.0 )
    aintb = r4_aint ( b - 0.5 );
  else
    aintb = r4_aint ( b + 0.5 );
  end
  beps = b - aintb;
  n = aintb;

  alnx = log ( x );
  xtoeps = exp ( - beps * alnx );
%
%  Evaluate the finite sum.
%
%  Consider the case b < 1.0 first.
%
  if ( n < 1 )

    sum = 1.0;
    t = 1.0;
    m = - n;
    for i = 1 : m
      xi1 = i - 1;
      t = t * ( a + xi1 ) * x / ( ( b + xi1 ) * ( xi1 + 1.0 ) );
      sum = sum + t;
    end

    sum = r4_poch ( 1.0 + a - b, - a ) * sum;
%
%  Now consider the case b .ge. 1.0.
%
  else

    sum = 0.0;
    m = n - 2;

    if ( 0 <= m )

      t = 1.0;
      sum = 1.0;

      for i = 1 : m
        xi = i;
        t = t * ( a - b + xi ) * x / ( ( 1.0 - b + xi ) * xi );
        sum = sum + t;
      end

      sum = r4_gamma ( b - 1.0 ) * r4_gamr ( a ) ...
        * x^( 1 - n ) * xtoeps * sum;

    end

  end
%
%  Now evaluate the infinite sum.
%
  if ( n < 1 )
    istrt = 1 - n;
  else
    istrt = 0;
  end

  xi = istrt;

  factor = r4_mop ( n ) * r4_gamr ( 1.0 + a - b ) * x^istrt;

  if ( beps ~= 0.0 )
    factor = factor * beps * pi / sin ( beps * pi );
  end

  pochai = r4_poch ( a, xi );
  gamri1 = r4_gamr ( xi + 1.0 );
  gamrni = r4_gamr ( aintb + xi );
  b0 = factor * r4_poch ( a, xi - beps ) * gamrni ...
    * r4_gamr ( xi + 1.0 - beps );
%
%  x^(-beps) is close to 1.0, so we must be careful in evaluating
%  the differences.
%
  if ( abs ( xtoeps - 1.0 ) <= 0.5 )

    pch1ai = r4_poch1 ( a + xi, - beps );
    pch1i = r4_poch1 ( xi + 1.0 - beps, beps );
    c0 = factor * pochai * gamrni * gamri1 * ( ...
      - r4_poch1 ( b + xi, -beps ) + pch1ai ...
      - pch1i + beps * pch1ai * pch1i );
%
%  xeps1 = (1.0 - x^(-beps)) / beps = (x^(-beps) - 1.0)/(-beps)
%
    xeps1 = alnx * r4_exprel ( - beps * alnx );
    value = sum + c0 + xeps1 * b0;
    xn = n;

    for i = 1 : 1000
      xi = istrt + i;
      xi1 = istrt + i - 1;
      b0 = ( a + xi1 - beps ) * b0 * x ...
        / ( ( xn + xi1 ) * ( xi - beps ) );
      c0 = ( a + xi1 ) * c0 * x / ( ( b + xi1 ) * xi ) ...
        - ( ( a - 1.0 ) * ( xn + 2.0 * xi - 1.0 )...
        + xi * ( xi - beps ) ) * b0 ...
        / ( xi * ( b + xi1 ) * ( a + xi1 - beps ) );
      t = c0 + xeps1 * b0;
      value = value + t;
      if ( abs ( t ) < eps * abs ( value ) )
        return
      end
    end

    fprintf ( 1, '\n' );
    fprintf ( 1, 'R4_CHU - Fatal error!\n' );
    fprintf ( 1, '  No convergence in 1000 terms.\n' );
    error ( 'R4_CHU - Fatal error!' )

  end
%
%  x^(-beps) is very different from 1.0, so the straightforward
%  formulation is stable.
%
  a0 = factor * pochai * r4_gamr ( b + xi ) * gamri1 / beps;
  b0 = xtoeps * b0 / beps;

  value = sum + a0 - b0;

  for i = 1 : 1000
    xi = istrt + i;
    xi1 = istrt + i - 1;
    a0 = ( a + xi1 ) * a0 * x / ( ( b + xi1 ) * xi );
    b0 = ( a + xi1 - beps ) * b0 * x ...
      / ( ( aintb + xi1 ) * ( xi - beps ) );
    t = a0 - b0;
    value = value + t;
    if ( abs ( t ) < eps * abs ( value ) )
      return
    end
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R4_CHU - Fatal error!\n' );
  fprintf ( 1, '  No convergence in 1000 terms.'\n' );

  error ( 'R4_CHU - Fatal error!' )
end

author: johannesgerer

repo: jburkardt-m


text:
function value = p23_autonomous ( )

%*****************************************************************************80
%
%% P23_AUTONOMOUS reports whether problem P23 is autonomous.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer VALUE, is 1 if the equation is autonomous.
%
  value = 1;

  return
end