text:
clear all, close all, clc
nu=0.001;   % Diffusion constant

% Define spatial domain
L = 20;             % Length of domain
N = 1000;           % Number of discretization points
dx = L/N;
x = -L/2:dx:L/2-dx; % Define x domain

% Define discrete wavenumbers
kappa = (2*pi/L)*[-N/2:N/2-1];
kappa = fftshift(kappa');    % Re-order fft wavenumbers

% Initial condition
u0 = sech(x);

% Simulate PDE in spatial domain
dt = 0.025;
t = 0:dt:100*dt;
[t,u] = ode45(@(t,u)rhsBurgers(t,u,kappa,nu),t,u0);

% Plot solution in time
subplot(1,2,1)
h=waterfall(real(u(1:10:end,:)));
set(h,'LineWidth',2,'FaceAlpha',.5);
colormap(jet/1.5)
view(5,55)

subplot(1,2,2)
imagesc(flipud(real(u)));
colormap jet

%% FIGURES (PRODUCTION)
figure
CC = colormap(jet(100));
dt = 0.1;
for k = 1:100
    if(mod(k-1,10)==0)
        plot(x,real(u(k,:)),'k','LineWidth',1.5)
%         plot(x,real(u(k,:)),'Color',CC(k,:),'LineWidth',1.5)
        hold on, grid on, drawnow
    end
end
% xlabel('Spatial variable, x')
% ylabel('Temperature, u(x,t)')
axis([-L/2 L/2 -.1 1.1])
set(gca,'LineWidth',1.2,'FontSize',12);
set(gcf,'Position',[100 100 550 220]);
set(gcf,'PaperPositionMode','auto')
print('-depsc2', '-loose', '../../figures/FFTBurgers1');

%
figure
subplot(1,2,1)
h=waterfall(real(u(1:10:end,:)));
set(h,'LineWidth',2,'FaceAlpha',.5);
colormap(jet/1.5)
view(5,55)
set(gca,'LineWidth',1.5)
set(gca,'XTick',[0 500 1000],'XTickLabels',{})
set(gca,'ZTick',[0 .5 1],'ZTickLabels',{})
set(gca,'YTick',[0 5 10],'YTickLabels',{})
set(gca,'YLim',[1 11])
set(gca,'ZLim',[-.1 1.1])

subplot(1,2,2)
imagesc(flipud(real(u)));
set(gca,'LineWidth',1.5)
set(gca,'XTick',[0 500 1000],'XTickLabels',{})
set(gca,'YTick',[0 50 100],'YTickLabels',{})

colormap jet
set(gcf,'Position',[100 100 600 250])
set(gcf,'PaperPositionMode','auto')
print('-depsc2', '-loose', '../../figures/FFTBurgers2');

author: dynamicslab


text:
function Q = metoffice_experiment_gpfa1(datanum, anomalies, validation)

%
% Load the data
%

if nargin < 3
  validation = true;
end

[data,dataset,folder,maskfile] = metoffice_getdata(datanum, anomalies, validation);

% Number of components
if anomalies
  comps = [0 20 80];
  %comps_spatial = [100 0 0]; % only independent components
  comps_spatial = [0 100 0];
else
  comps = [5 15 80];
  %comps_spatial = [100 0 0]; % only independent components
  comps_spatial = [5 95 0];
end
D = sum(comps);

maxiter = 2000;

%
% Model for temporal X
%

covfunc_x = cell(D,1);
theta_x = cell(D,1);
is_pseudos_x = false(D,1);

% Inputs (assume uniformly spaced time instances which is not exactly correct)
%in_x = data.time;
%pseudo_x = data.time(1:10:end);
in_x = linspace(data.time(1), data.time(end), length(data.time));
pseudo_x = in_x(1:10:end);

% Squared distances for covariance functions
D2_pp = sq_dist(pseudo_x);
D2_px = sq_dist(pseudo_x, in_x);
d2_x = zeros(size(in_x,2),1); %diag(D2_xx);
D_pp = sqrt(D2_pp);
D_px = sqrt(D2_px);
d_x = sqrt(d2_x);

% Update only the smooth components at the beginning
update_schedule = cell(maxiter,1);
update_schedule(:) = {1:D};
update_schedule(1:10) = {1:(comps(1)+comps(2))};
update_schedule_noise = 11;

% Distance matrices for covariance functions
%d_xx = sqrt(D2_xx(1,:));
d_xx = sqrt(sq_dist(in_x(:,1),in_x));

ind = 0;

% Periodic components (1 year period) with decay (rational quadratic)
if comps(1) > 0
  ind = ind(end) + (1:comps(1));
  fprintf('%d periodic components for X\n', length(ind));
  covfunc = @(D,D2) gp_cov_product(gp_cov_periodic(D, 'wavelength', 365.26), ...
                                   gp_cov_rq(D2));
  covfunc_x(ind) = {gp_cov_pseudo(gp_cov_jitter(covfunc(D_pp,D2_pp), 1e-4), ...
                                  covfunc(D_px,D2_px), ...
                                  covfunc(d_x,d2_x))};
  theta_x(ind) = columns_to_cells(...
      [linspace(1,0.1,length(ind))      % smoothness of the period
       365*linspace(500,50,length(ind)) % lengthscale of the decay (RQ)
       ones(1,length(ind))]);         % alpha (RQ)
  is_pseudos_x(ind) = true;
end

% Slow components (1-20 years): rational quadratic using pseudo inputs
if comps(2) > 0
  ind = ind(end) + (1:comps(2));
  fprintf('%d slow components for X\n', length(ind));
% $$$ covfunc = @(D2) gp_cov_se(D2); % SE
  covfunc = @(D2) gp_cov_rq(D2); % RQ
  covfunc_x(ind) = {gp_cov_pseudo(gp_cov_jitter(covfunc(D2_pp), 1e-3), ...
                                  covfunc(D2_px), ...
                                  covfunc(d2_x))};
% $$$ theta_x(ind) = columns_to_cells(...
% $$$     [365*linspace(20,1,length(ind))]); % lengthscale for SE
  theta_x(ind) = columns_to_cells(...
      [365*linspace(2,1,length(ind)) % lengthscale for RQ
       ones(1,length(ind))]);         % alpha for RQ
  is_pseudos_x(ind) = true;
end

% Fast components (4-18 months): piecewise polynomial in 1-D
% Take advantage of the Toeplitz structure of the covariance matrix
if comps(3) > 0
  ind = ind(end) + (1:comps(3));
  fprintf('%d fast components for X\n', length(ind));
  covfunc_x(ind) = {gp_cov_jitter(gp_cov_toeplitz(gp_cov_pp(d_xx,1)))};
  theta_x(ind) = columns_to_cells(...
      [30*linspace(6,3,length(ind))]); % lengthscale or cut-off
  is_pseudos_x(ind) = false;
end

%
% Model for spatial W
%

covfunc_w = cell(D,1);
theta_w = cell(D,1);
is_pseudos_w = false(D,1);

% Remove land area grid points
in_w = data.coordinates;
%in_w = mohsst5_remove_land(data.coordinates')';

%% Smooth components (using pseudo inputs)

% Pseudo inputs (uniformly with respect to area size)
pseudo_w = points_on_sphere(18); % uniform points by number of latitudes
% Remove pseudo inputs that are on land (the nearest grid point is land)
ind_pseudo_w = mohsst5_points_to_grid_index(pseudo_w);
mask = metoffice_get_mask(maskfile);
pseudo_w(:,~mask(ind_pseudo_w)) = [];
% $$$ % This code shows the pseudo inputs on the map
% $$$ figure
% $$$ map_projection('global-ellipse');
% $$$ map_plot(pseudo_w,'r+');
% $$$ map_coast()
% $$$ map_grid()
% $$$ return

% Transform inputs to 3-D Euclidean coordinates
in_w = geographic_to_euclidean(in_w);
pseudo_w = geographic_to_euclidean(pseudo_w);

% Squared distance matrices for the covariance functions
D2_ww = sq_dist(in_w);
D2_pp = sq_dist(pseudo_w);
D2_pw = sq_dist(pseudo_w, in_w);
d2_w = diag(D2_ww);

ind = 0;

if comps_spatial(1) > 0
  ind = ind(end) + (1:comps_spatial(1));
  fprintf('%d iid components for W\n', length(ind));

  covfunc_w(ind) = {gp_cov_scale(gp_cov_delta(size(in_w,2)))};
  theta_w(ind) = columns_to_cells(...
      [linspace(1,0.1,length(ind))]);       % magnitudes
end

if comps_spatial(2) > 0

  ind = ind(end) + (1:comps_spatial(2));
  fprintf('%d slow components for W (using %d pseudo inputs)\n', length(ind), ...
          size(pseudo_w,2));

  % Covariance function (scaled squared exponential) with pseudo inputs
  covfunc = @(D2) gp_cov_se(D2);
  covfunc_w(ind) = {gp_cov_pseudo(...
      gp_cov_scale(gp_cov_jitter(covfunc(D2_pp), 1e-3)), ...
      gp_cov_scale(covfunc(D2_pw)), ...
      gp_cov_scale(covfunc(d2_w)))};

  % Hyperparameters for the covariance functions
  theta_w(ind) = columns_to_cells(...
      [linspace(1,0.1,length(ind));       % magnitudes
       linspace(5000,1000,length(ind))]); % lengthscales
  is_pseudos_w(ind) = true;
end

if comps_spatial(3) > 0
% $$$ %% Short scale components: piecewise polynomial in 3-D
% $$$
% $$$ ind = (ind(end)+1):D;
% $$$ fprintf('%d fast components for W\n', length(ind));
% $$$
% $$$ % Use block-Toeplitz structure for the covariance function
% $$$ [lat,lon0] = meshgrid(data.latitude,data.longitude(1));
% $$$ in_w0 = geographic_to_euclidean([lon0(:)';lat(:)']);
% $$$ d_ww = sqrt(sq_dist(in_w0, ...
% $$$                     geographic_to_euclidean(data.coordinates)));
% $$$ covfunc = gp_cov_toeplitz_block(gp_cov_pp(d_ww,3));
% $$$ % Select only sea areas (discard land areas)
% $$$ sea = sum(~isnan(data.observations),2) > 0;
% $$$ covfunc = gp_cov_select(covfunc, sea);
% $$$ % Add scaling and jitter
% $$$ covfunc_w(ind) = {gp_cov_scale(gp_cov_jitter(covfunc))};
% $$$ % Hyperparameters
% $$$ theta_w(ind) = columns_to_cells(...
% $$$     [linspace(0.5,0.1,length(ind));     % magnitude
% $$$      linspace(3000,2000,length(ind))]); % lengthscale
% $$$ is_pseudos_w(ind) = false;
end

%
% Process data
%

% Form the data matrix
Y = data.data;
%Y = mohsst5_remove_land(data.data);
[M,N] = size(Y);
Obs = ~isnan(Y);

%
% GPFA inference
%

% Filename for saving the results
folder = [folder '/gpfa'];
mkdir(folder);
filename = sprintf('%s/results_rectest_%s_gpfa1_D=%d_anomalies=%d_remval=%d_%s', ...
                   folder, ...
                   dataset, ...
                   D, ...
                   anomalies, ...
                   validation, ...
                   datestr(now,'yyyymmdd'));

% Component-wise factorization for X
X_module = factor_module_gp_factorized(N, covfunc_x, theta_x, ...
                                       'update_hyperparameters', [5 10:10:100 100:25:2000], ...
                                       'maxiter_hyperparameters', 5, ...
                                       'is_pseudo', is_pseudos_x, ...
                                       'init', zeros(D,N), ...
                                       'update_schedule', {update_schedule});

% Component-wise factorization for W
W_module = factor_module_gp_factorized(M, covfunc_w, theta_w, ...
                                       'update_hyperparameters', [5 10:10:100 100:25:2000], ...
                                       'maxiter_hyperparameters', 5, ...
                                       'is_pseudo', is_pseudos_w, ...
                                       'update_schedule', {update_schedule});

% Isotropic noise (precisions weighted proportionally to grid size)
weights = repmat(data.gridsize, [1, N]);
noise_module = noise_module_isotropic(M, N, 1e-3, 1e-3, ...
                                      'init', 4, ...
                                      'weights', weights, ...
                                      'update_schedule', update_schedule_noise);

% Run GPFA
Q = gpfa(D, Y, W_module, X_module, noise_module, ...
         'maxiter', maxiter, ...
         'rotate', 1:50, ... %[1:50 60:10:2000], ...
         'autosavefile', filename, ...
         'autosave', [10:100:2000]);

% Reconstruct
Yh = Q.W'*Q.X;

% Some performance measures
fprintf('Weighted training RMSE of the reconstruction: %f\n',  ...
        rmsew(Y(Obs)-Yh(Obs),weights(Obs)));

% Save the results
save(filename, '-struct', 'Q');

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function pass = test_feval(pref)
% Test feval

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

tol = 100*pref.cheb2Prefs.chebfun2eps;

f = chebfun2(@(x,y) x);
pass(1) = (abs(f(pi/6,pi/12)-pi/6) < tol);

f = chebfun2(@(x,y) x, [-1 2 -pi/2 pi]);
pass(2) = (abs(f(0,0)) < 1e-14);
pass(3) = (abs(f(pi/6,pi/12)-pi/6) < tol);

f = chebfun2(@(x,y) y, [-1 2 -pi/2 pi]);
pass(4) = (abs(f(0,0)) < 1e-14);
pass(5) = (abs(f(pi/6,pi/12)-pi/12) < tol);

% some harder tests.
f = @(x,y) cos(x) + sin(x.*y);
g = chebfun2(f);

r = 0.126986816293506; s = 0.632359246225410; % two fixed random number in domain.
pass(6) = (abs(f(r,s) - g(r,s))<tol);

% Are we evaluating on arrays correctly
seedRNG(0)
r = rand(10,1);
s = rand(10,1);
[rr, ss]=meshgrid(r,s);
pass(7) = (norm((f(r,s) - g(r,s))) < tol);
pass(8) = (norm((f(rr,ss) - g(rr,ss))) < tol);  % on arrays as well.

% Does this work off [-1,1]^2
g = chebfun2(f,[-pi/6 pi/2 -pi/12 sqrt(3)]); % strange domain.
r = 0.126986816293506; s = 0.632359246225410; % two fixed random number in domain.
pass(9) = (abs(f(r,s) - g(r,s))<tol);

% Are we evaluating on arrays correctly
r = rand(10,1); s = rand(10,1); [rr, ss]=meshgrid(r,s);
pass(10) = (norm((f(r,s) - g(r,s)))<tol);
pass(11) = (norm((f(rr,ss) - g(rr,ss)))<2*tol); % on arrays as well.

% Evaluation at complex arguments.
f = chebfun2(@(x,y) x+1i*y);
pass(12) = ( norm( f(-1,-1) - (-1-1i) )  < tol );
pass(13) = ( norm( f(r,s) - (r + 1i*s) )  < tol );
pass(14) = (norm((f(rr,ss) - (rr + 1i*ss)))<tol); % on arrays as well.

% Evaluation at complex arguments, different syntax.
f = chebfun2(@(z) z);
pass(15) = ( norm( f(-1-1i) - (-1-1i) )  < tol );
pass(16) = ( norm( f(r+1i*s) - (r + 1i*s) )  < tol );
pass(17) = (norm((feval(f,rr,ss) - (rr + 1i*ss)))<tol); % on arrays as well.

% Evaluation at transposed meshgrid
op = @(x,y) cos(pi*(x+y))+y.*sin(pi*x);
[yy,xx] = meshgrid(linspace(-1,1,1001));
f = chebfun2( op );
pass(18) = (norm( feval(f,xx,yy) - op(xx,yy) ) < 100*tol);

% Check sizes:
n = 10;
f = chebfun2(@(x,y) cos(x.*y));
x = ones(1,n);
[xx, yy] = meshgrid(x);
pass(19) = ( all( size( feval(f, 1, 1) ) == [1 1] ) );
pass(20) = ( all( size( feval(f, [1 1], [1 1]) ) == [1 2] ) );
pass(21) = ( all( size( feval(f, [1;1], [1;1]) ) == [2 1] ) );
pass(22) = ( all( size( feval(f, [1 1;1 1], [1 1; 1 1]) ) == [2 2] ) );
pass(23) = ( all( size( feval(f, x, x) ) == size(x) ) );
pass(24) = ( all( size( feval(f, xx, yy)  ) == [n n] ) );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
%kigeowarp 'Perform Direct Bilinear Geometric Warping'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros igeowarp.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'Input data object'
% Toggle: planes 'Process by planes', default: 0: 'Process data by full planes instead of small prisms'
% OutputFile: o 'Output', required: 'Output data object'
% InputFile: wcoeffs 'W warp coeffs object', optional: 'W warp function coeffs object'
% InputFile: hcoeffs 'H warp coeffs object', optional: 'H warp function coeffs object'
%
% Example: o = kigeowarp({i, wcoeffs, hcoeffs}, {'i','';'planes',0;'o','';'wcoeffs','';'hcoeffs',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% igeowarp - Perform Direct Bilinear Geometric Warping
%
%  DESCRIPTION
% .I igeowarp
% is used to apply a direct bilinear warping to an image. In this mode, the
% bilinear mapping equations give the mapping of output data coordinates to
% input data coordinates. A given output data coordinate is operated on by the
% mapping functions to give the input data coordinate. The value of the input
% data coordinate is interpolated using bilinear interpolation, and the resulting
% value is assigned to the output data coordinate.
%
% All warping is assumed to occur only in the WxH plane.
% If warping in a different plane is desired, it is necessary to use the kaxis
% program to reorient the data such that the desired plane is the WxH plane before
% applying igeowarp. This limitation is due to computational complexity and
% performance issues, particularly when dealing with large data sets.
%
% The mapping functions used by igeowarp are:
%
%
%             f(w,h) = output image
%
%             g(w',h') = input image coordinate
%
%             w' = a00 + a01*w + a10*h +a11*w*h
%             h' = b00 + b01*w + b10*h +b11*w*h
%
%
% Note that in general, w' and h' are not integers for a given integer pair
% w and h. The
% data value at g(w',h') is obtained by bilinear interpolation of the four
% nearest neighbors to coordinate (w',h').
%
% The mapping is called direct because it directly gives the location in the
% input data that corresponds to a given location in the output data. The
% opposite case gives the location in the output data that corresponds to a
% given input data location; this requires much trickier and computationally
% expensive interpolation methods to obtain good results.
%
% The input parameters a00, a01, a10, a11, b00, b01, b10, b11 are exactly as
% specified in the mapping functions above. These parameters can also be supplied
% in the form of object value data where the data is assumed to be organized with
% the parameters a00, a01, a10, and a11 present in that storage order for the
% W coefficients and likewise for the H coefficients.
%
% If the input object has a mask, a new mask will be computed for the output
% object indicating which data points contain reliable data.
%
% If the input object has a map, the data is pulled through the map prior to
% rotation, and the output object will have no map.
%
% The -planes flag is used to change the way the data is accessed for processing.
% If the data set is small enough that individual planes of data (including the
% E data) parallel to the WxH plane will fit in memory, then use of the -planes
% flag will permit processing to happen in a plane-by-plane basis; this can
% occur with efficiency and speed. Otherwise,
% processing will occur in prisms down the D axis, which is much slower but will
% work on data sets of any size. If you have lots of memory, you may be able to
% get away with using -planes even for quite large data sets, say around 2Kx2K
% or more points per plane.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
% irotate(1)
%
%  RESTRICTIONS
%
% For this initial release, the E dimension data is "not\fR handled
% correctly. Data in WxHxD "is\fR handled correctly. This deficiency
% will be corrected in the next version of
% .I igeowarp.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kigeowarp(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kigeowarp(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'planes', 0;'o', '__output';'wcoeffs', '__input';'hcoeffs', '__input'};
maxval={0,0,0,1,1};
minval={0,0,0,1,1};
istoggle=[0,1,0,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','Toggle','OutputFile','InputFile','InputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'igeowarp"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function [textures] = getGLSZMtextures_STS(GLSZM)
% -------------------------------------------------------------------------
% function [textures] = getGLSZMtextures_STS(GLSZM)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes texture features from an input Gray-Level Size
% Zone Matrix (GLSZM). (STS study)
% -------------------------------------------------------------------------
% REFERENCES:
% [1] Galloway, M. M. (1975). Texture analysis using gray level run lengths.
%     Computer Graphics and Image Processing, 4(2), 172–179.
% [2] Chu, A., Sehgal, C. M., & Greenleaf, J. F. (1990). Use of gray value
%     distribution of run lengths for texture analysis. Pattern Recognition
%     Letters, 11(6), 415-419.
% [3] Dasarathy, B. V., & Holder, E. B. (1991). Image characterizations
%     based on joint gray level-run length distributions. Pattern
%     Recognition Letters, 12(8), 497-502.
% [4] Thibault, G., Fertil, B., Navarro, C., Pereira, S., Cau, P., Levy,
%     N., Mari, J.-L. (2009). Texture Indexes and Gray Level Size Zone
%     Matrix. Application to Cell Nuclei Classification. In Pattern
%     Recognition and Information Processing (PRIP) (pp. 140–145).
% -------------------------------------------------------------------------
% INPUTS:
% - GLSZM: Gray-Level Size Zone Matrix.
%
% ** 'GLSZM' should be the output from 'getGLSZM.m' **
% -------------------------------------------------------------------------
% OUTPUTS:
% - textures: Struture specifying the values of different GLSZM texture
%             features as defined below.
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: January 2013
% - Revision: May 2015
% -------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

% USEFUL MATRICES, VECTORS AND QUANTITIES
sz = size(GLSZM); % Size of GLSZM
nRuns = sum(GLSZM(:));
cVect = 1:sz(2); rVect = 1:sz(1);% Row and column vectors
[cMat,rMat] = meshgrid(cVect,rVect); % Column and row indicators for each entry of the GLSZM
pg = sum(GLSZM,2)'; % Gray-Level Run-Number Vector
pr = sum(GLSZM); % Run-Length Run-Number Vector

% COMPUTATION OF TEXTURE FEATURES
% 1. Small Zone Emphasis (SZE), Ref.[1,4]
textures.SZE = (pr*(cVect.^(-2))')/nRuns;

% 2. Large Zone Emphasis (LZE), Ref.[1,4]
textures.LZE = (pr*(cVect.^2)')/nRuns;

% 3. Gray-Level Nonuniformity (GLN), adapted from Ref.[1,4]
textures.GLN = sum(pg.^2)/nRuns;

% 4. Zone-Size Nonuniformity (ZSN), adapted from Ref.[1,4]
textures.ZSN = sum(pr.^2)/nRuns;

% 5. Zone Percentage (ZP), adapted from Ref.[1,4]
textures.ZP = nRuns/(pr*cVect');

% 6. Low Gray-Level Zone Emphasis (LGZE), Ref.[2,4]
textures.LGZE = (pg*(rVect.^(-2))')/nRuns;

% 7. High Gray-Level Zone Emphasis (HGZE), Ref.[2,4]
textures.HGZE = (pg*(rVect.^2)')/nRuns;

% 8. Small Zone Low Gray-Level Emphasis (SZLGE), Ref.[3,4]
textures.SZLGE = sum(sum(GLSZM.*(rMat.^(-2)).*(cMat.^(-2))))/nRuns;

% 9. Small Zone High Gray-Level Emphasis (SZHGE), Ref.[3,4]
textures.SZHGE = sum(sum(GLSZM.*(rMat.^2).*(cMat.^(-2))))/nRuns;

% 10. Large Zone Low Gray-Level Emphasis (LZLGE), Ref.[3,4]
textures.LZLGE = sum(sum(GLSZM.*(rMat.^(-2)).*(cMat.^2)))/nRuns;

% 11. Large Zone High Gray-Level Emphasis (LZHGE), Ref.[3,4]
textures.LZHGE = sum(sum(GLSZM.*(rMat.^2).*(cMat.^2)))/nRuns;

% New features according to Ref.[4]
GLSZM = GLSZM./nRuns; % In the future, this operation will be applied at the beginning of the function
pg=sum(GLSZM,2)'; pr=sum(GLSZM);
ug = (pg*rVect')/(sz(1)*sz(2));
ur = (pr*cVect')/(sz(1)*sz(2));

% 12. Gray-Level Variance (GLV), adapted from Ref.[4]
GLV = 0;
for g = 1:sz(1)
    for r = 1:sz(2)
        GLV = GLV + (GLSZM(g,r)*g-ug)^2;
    end
end
textures.GLV = GLV/(sz(1)*sz(2));

% 13. Zone-Size Variance (ZSV), adapted from Ref.[4]
ZSV = 0;
for g = 1:sz(1)
    for r = 1:sz(2)
        ZSV = ZSV + (GLSZM(g,r)*r-ur)^2;
    end
end
textures.ZSV = ZSV/(sz(1)*sz(2));

end

author: mvallieres

repo: radiomics


text:
function [y]=funcrs(tt,fun,eps,y,nswp,varargin)
%Cross approximation of a function of a TT-tensor, Method 1
%   [Y]=FUNCRS(TT,FUN,EPS,Y,NSWP)
%   Computes approximation to the function FUN(TT) with accuracy EPS
%   Auxiliary parameters:  Y (initial approximation), NSWP
%   (number of sweeps in the method
%   Much faster then usual cross by vectorized computation of subtensors
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

%PARAMETERS SECTION
rmin=1;
dropsweeps = 4;
kick_rank=10;
verb=true;
if (~isempty(y))
   yold=y;
end
%For this procedure we need only local (!) indices, since it
%is more or less equivalent to orthogonalization;
d=tt.d;
ps=tt.ps;
core=tt.core;
n=tt.n;
r=tt.r;

ry=y.r;
psy=y.ps;
cry=y.core;
phi=cell(d+1,1);
phi{d+1}=1;
phi{1}=1;
%Warmup procedure: orthogonalize from right-to-left & maxvol
pos1=psy(d);
%Cores i
for i=d-1:-1:1
   cr=cry(pos1:pos1+ry(i+1)*n(i+1)*ry(i+2)-1);
   cr2=cry(pos1-ry(i)*n(i)*ry(i+1):pos1-1);
   cr2=reshape(cr2,[ry(i)*n(i),ry(i+1)]);
   cr=reshape(cr,[ry(i+1),n(i+1)*ry(i+2)]);
   cr=cr.';
   [cr,rm]=qr(cr,0);
   ry(i+1)=size(cr,2);
   ind=maxvol2(cr);
   r1=cr(ind,:);
   cr=cr/r1;
   cr=cr.';
   cry(pos1:pos1+ry(i+1)*n(i+1)*ry(i+2)-1)=cr(:);
   pos1=pos1-ry(i)*n(i)*ry(i+1);
   cr2=cr2*(r1*rm).';
   cry(pos1:pos1+ry(i)*n(i)*ry(i+1)-1)=cr2(:);
   %Take phi matrix; convolve from right with current cors of V
   cr0=core(ps(i+1):ps(i+2)-1);
   cr0=reshape(cr0,[r(i+1)*n(i+1),r(i+2)]);
   cr0=cr0*phi{i+2}; %cr0 is now r(i)*n(i)*ry(i+1);
   cr0=reshape(cr0,[r(i+1),n(i+1)*ry(i+2)]);
   phi{i+1}=cr0(:,ind);
 %  psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
%y.core=cry;
%y.r=ry;
%y.ps=psy;
%keyboard

   %Orthogonalization & maxvol is "local operation" (touches two
   %cores)
   %Computation of elements requries storage of phi matrices
end
%Truncate cry
%pos1=pos1-n(1)*r(2);
cry=cry(pos1:numel(cry));
y.core=cry;
y.r=ry;
y.ps=psy;
%keyboard
swp=1;

yold=[];
not_converged = true;
pos1=1;
last_sweep = false;
while ( swp < nswp && not_converged )
    max_er=0;
cry_old=cry;
psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
pos1=1;
 for i=1:d-1
     %fprintf('i=%d \n',i);
     %We care for two cores, with number i & number i+1, and use
     %psi(i) and psi(i+2) as a basis; also we will need to recompute
     %psi(i+1)
     ps1=phi{i}; ps2=phi{i+2};
     %Take current core; convolve it with
     %core=core(ps
     %Compute (!) superblock and function (!) of it
     cr1=core(ps(i):ps(i+1)-1);
     cr2=core(ps(i+1):ps(i+2)-1);
     cr1=reshape(cr1,[r(i),n(i)*r(i+1)]);
     cr1=ps1*cr1;
     cr2=reshape(cr2,[r(i+1)*n(i+1),r(i+2)]);
     cr2=cr2*ps2;
     cr1=reshape(cr1,[ry(i)*n(i),r(i+1)]);
     cr2=reshape(cr2,[r(i+1),n(i+1)*ry(i+2)]);
     cr=cr1*cr2;
     cr=reshape(cr,[ry(i)*n(i),n(i+1)*ry(i+2)]);
     cr=fun(cr); %Elements are evaluated here!
     cr=reshape(cr,[ry(i)*n(i),n(i+1)*ry(i+2)]);
     %Check for local approximation of cr for the error
     cry1=cry(pos1:pos1+ry(i)*n(i)*ry(i+1)-1);

     cry2=cry_old(psy(i+1):psy(i+2)-1);

     cry1=reshape(cry1,[ry(i)*n(i),ry(i+1)]);
     cry2=reshape(cry2,[ry(i+1),n(i+1)*ry(i+2)]);
     appr=cry1*cry2;
     er=norm(appr-cr,'fro')/norm(cr,'fro');
     max_er=max(er,max_er);
     %Compute SVD of cr

     if (mod(swp,dropsweeps)~=0)&&(swp>1)&&(~last_sweep)
         [u,s,v]=svd(cr-appr,'econ');
     else
         [u,s,v]=svd(cr,'econ');
     end;
     s=diag(s);
     r2=my_chop2(s,eps*norm(cr,'fro')/sqrt(d-1));
     s=s(1:r2); u=u(:,1:r2); v=v(:,1:r2);

     v=conj(v)*diag(s);

     if (mod(swp,dropsweeps)~=0)&&(swp>1)&&(~last_sweep)
         u = [cry1, u];
         v = [cry2.', v];
         [u,rv]=qr(u,0);
         v = v*(rv.');
     else
         if (~last_sweep)
             %Kick rank of u
             ur=randn(size(u,1),kick_rank);
             %Orthogonalize ur to u by Golub-Kahan reorth
             u=reort(u,ur);
             radd=size(u,2)-r2;
             if ( radd > 0 )
                 vr=zeros(size(v,1),radd);
                 v=[v,vr];
             end
         end;
     end;

     %vr=randn(size(v,1),kick_rank);
     %v=reort(v,vr);
     %radd=size(v,2)-rnew;
     %if ( radd > 0 )
     %    ur=zeros(size(u,1),radd);
     %    u=[u,ur];
     %end
     %rnew=rnew+radd;

     %u=[u,uadd]; v=[v,vadd];
     %[u,ru]=qr(u,0);
     %v=v*(ru.');
     r2=size(u,2);

     ind=maxvol2(u);
     r1=u(ind,:);
     u=u/r1; v=v*r1'; v=v.';
     %Compute new phi
     ry(i+1)=r2;
     cry(pos1:pos1+ry(i)*n(i)*ry(i+1)-1)=u(:);
     pos1=pos1+ry(i)*n(i)*ry(i+1);
     cry(pos1:pos1+ry(i+1)*n(i+1)*ry(i+2)-1)=v(:);

     %Now we have to: cr1 with phi from the left
     phi{i+1}=cr1(ind,:); %phi{i+1}=phi{i+1};

 end
  %Truncate local memory
  cry=cry(1:pos1+ry(d)*n(d)*ry(d+1)-1);
  psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);
y.core=cry;
y.r=ry;
y.ps=psy;
%keyboard;
cry_old=cry;
%m=numel(cry);
%cry=[zeros(size(cry)),cry];
%pos1=pos1+m;
%cry2=cry_old(psy(i+1):psy(i+2)-1);
cry=cry(psy(d):psy(d+1)-1); %Start--only two cores
 for i=d-1:-1:1
     %fprintf('i=%d \n',i);
     %We care for two cores, with number i & number i+1, and use
     %psi(i) and psi(i+2) as a basis; also we will need to recompute
     %psi(i+1)
     ps1=phi{i}; ps2=phi{i+2};
     %Take current core; convolve it with
     %core=core(ps
     %Compute (!) superblock and function (!) of it
     cr1=core(ps(i):ps(i+1)-1);
     cr2=core(ps(i+1):ps(i+2)-1);
     cr1=reshape(cr1,[r(i),n(i)*r(i+1)]);
     cr1=ps1*cr1;
     cr2=reshape(cr2,[r(i+1)*n(i+1),r(i+2)]);
     cr2=cr2*ps2;
     cr1=reshape(cr1,[ry(i)*n(i),r(i+1)]);
     cr2=reshape(cr2,[r(i+1),n(i+1)*ry(i+2)]);
     cr=cr1*cr2;
     cr=reshape(cr,[ry(i)*n(i),n(i+1)*ry(i+2)]);
     cr=fun(cr); %Elements are evaluated here!
     cr=reshape(cr,[ry(i)*n(i),n(i+1)*ry(i+2)]);
     %Check for local approximation of cr for the error
     %cry1=cry(pos1-ry(i)*n(i)*ry(i+1):pos1-1);
     cry1=cry_old(psy(i):psy(i+1)-1);
     %cry2=cry(ry(:ry(i+1)*n(i+1)*ry(i+2));
     %cry1=cry(1:ry(i)*n(i)*ry(i+1));
     %cry2=cry(ry(i)*n(i)*ry(i+1)+1:ry(i)*n(i)*ry(i+1)+ry(i+1)*n(i+1)*ry(i+2));
     cry2=cry(1:ry(i+1)*n(i+1)*ry(i+2));
     cry(1:ry(i+1)*n(i+1)*ry(i+2))=[];
     %cry2=cry_old(psy(i+1):psy(i+2)-1);
     %cry(1:(ry(i+1)*n(i+1)*ry(i+2)))=[]; %Delete both of first cores
     %cry(1:ry(i)*n(i)*ry(i+1)+ry(i+1)*n(i+1)*ry(i+2))=[];
     cry1=reshape(cry1,[ry(i)*n(i),ry(i+1)]);
     cry2=reshape(cry2,[ry(i+1),n(i+1)*ry(i+2)]);
     appr=cry1*cry2;
     er=norm(appr-cr,'fro')/norm(cr,'fro');
     %er
     max_er=max(er,max_er);

     if (mod(swp,dropsweeps)~=0)&&(swp>1)&&(~last_sweep)
         [u,s,v]=svd(cr-appr,'econ');
     else
         %Compute SVD of cr
         [u,s,v]=svd(cr,'econ');
     end;
     s=diag(s);
     r2=my_chop2(s,eps*norm(cr,'fro')/sqrt(d-1));
     s=s(1:r2); u=u(:,1:r2); v=conj(v(:,1:r2));
     %Make it standard
     u=u*diag(s);

     if (mod(swp,dropsweeps)~=0)&&(swp>1)&&(~last_sweep)
         u = [cry1, u];
         v = [cry2.', v];
         [v,rv]=qr(v,0);
         u = u*(rv.');
     else
         if (~last_sweep)
             %Kick rank
             vr=randn(size(v,1),kick_rank);
             v=reort(v,vr);
             radd=size(v,2)-r2;
             if ( radd > 0 )
                 ur=zeros(size(u,1),radd);
                 u=[u,ur];
             end
         end;
     end;

     r2=size(v,2);
     ind=maxvol2(v);
     r1=v(ind,:);
     v=v/r1; v=v.';
     u=u*r1';
     %Compute new phi;
     ry(i+1)=r2;
     u=u(:); u=u'; v=v(:); v=v';
     cry=[u,v,cry];
     %keyboard;
     %We need new memory for
     %cry(pos1:pos1+ry(i+1)*n(i+1)*ry(i+2)-1)=v(:); %Here memory has to be (?) enlarged
     %if ( pos1 <= ry(i)*n(i)*ry(i+1) ) %Have to enlarge memory
     %  cry=cry(pos1:numel(cry));
     %  cry=[u(:)',cry];
     %  pos1=ry(i)*n(i)*ry(i+1)+1;
     %else
     %   pos1=pos1-ry(i)*n(i)*ry(i+1);
      %  cry(pos1:pos1+ry(i)*n(i)*ry(i+1)-1)=u(:);
     %end
     %Now we have to: cr1 with phi from the left
     phi{i+1}=cr2(:,ind); phi{i+1}=phi{i+1};

 end
 %Truncate local memory
 %cry=cry(pos1:numel(cry));

% keyboard;
  psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);

y.core=cry;
y.r=ry;
y.ps=psy;
%keyboard;
if ( isempty(yold) )
 yold=y;
 er_nrm=1;
else

   er_nrm=norm(yold-y)/norm(y);
   yold=y;
end
if ( verb )
 fprintf('sweep=%d, er=%3.2e er_nrm=%3.2e \n',swp,max_er,er_nrm);
end
if (last_sweep)
    break;
end;
if (er_nrm<eps)
    last_sweep=true;
end;
swp=swp+1;
end
  psy=cumsum([1;n.*ry(1:d).*ry(2:d+1)]);

y.core=cry;
y.r=ry;
y.ps=psy;

end

author: oseledets

repo: TT-Toolbox


text:
function y = msmoothboxlogpdf(x,a,b,sigma)
%MSMOOTHBOXLOGPDF Multivariate smooth-box log probability density function.
%   Y = MSMOOTHBOXLOGPDF(X,A,B,SIGMA) returns the logarithm of the pdf of
%   the multivariate smooth-box distribution with pivots A and B and scale
%   SIGMA, evaluated at the values in X. The multivariate smooth-box pdf is
%   the product of univariate smooth-box pdfs in each dimension.
%
%   For each dimension i, the univariate smooth-box pdf is defined as a
%   uniform distribution between pivots A(i), B(i) and Gaussian tails that
%   fall starting from p(A(i)) to the left (resp., p(B(i)) to the right)
%   with standard deviation SIGMA(i).
%
%   X can be a matrix, where each row is a separate point and each column
%   is a different dimension. Similarly, A, B, and SIGMA can also be
%   matrices of the same size as X.
%
%   The log pdf is typically preferred in numerical computations involving
%   probabilities, as it is more stable.
%
%   See also MSMOOTHBOXPDF, MSMOOTHBOXRND.

% Luigi Acerbi 2022

[N,D] = size(x);

if any(sigma(:) <= 0)
    error('msmoothboxpdf:NonPositiveSigma', ...
        'All elements of SIGMA should be positive.');
end

if D > 1
    if isscalar(a); a = a*ones(1,D); end
    if isscalar(b); b = b*ones(1,D); end
    if isscalar(sigma); sigma = sigma*ones(1,D); end
end

if size(a,2) ~= D || size(b,2) ~= D || size(sigma,2) ~= D
    error('msmoothboxpdf:SizeError', ...
        'A, B, SIGMA should be scalars or have the same number of columns as X.');
end

if size(a,1) == 1; a = repmat(a,[N,1]); end
if size(b,1) == 1; b = repmat(b,[N,1]); end
if size(sigma,1) == 1; sigma = repmat(sigma,[N,1]); end

if any(a(:) >= b(:))
    error('msmoothboxpdf:OrderError', ...
        'For all elements of A and B, the order A < B should hold.');
end

y = -inf(size(x));
lnf = log(1/sqrt(2*pi)./sigma) - log1p(1/sqrt(2*pi)./sigma.*(b - a));

for ii = 1:D
    idx = x(:,ii) < a(:,ii);
    y(idx,ii) = lnf(idx,ii) - 0.5*((x(idx,ii) - a(idx,ii))./sigma(idx,ii)).^2;

    idx = x(:,ii) >= a(:,ii) & x(:,ii) <= b(:,ii);
    y(idx,ii) = lnf(idx,ii);

    idx = x(:,ii) > b(:,ii);
    y(idx,ii) = lnf(idx,ii) - 0.5*((x(idx,ii) - b(idx,ii))./sigma(idx,ii)).^2;
end

y = sum(y,2);

author: acerbilab

repo: vbmc

save_path: github-repos/MATLAB/acerbilab-vbmc


text:
clear all; close all; clc;

S = spherefactory(5, 2);
M = powermanifold(S, 7);
checkmanifold(M);
checkretraction(M);

author: NicolasBoumal

repo: manopt


text:
function [ mat ] = tsgReadMatrix( filename )
%
% [ mat ] = tsgReadMatrix( filename )
%
% reads a matrix from a file format
%
% 3 4
% 1 2 3 4
% 5 6 7 8
% 9 10 11 12
%
% results in the matrix [ 1 2 3 4; 5 6 7 8; 9 10 11 12; ]
%

fid = fopen ( filename, 'rt' );

[ s ] = fscanf( fid, ' %d ', [1, 2 ] ); % load the number of points

Ni = s(1);
Nj = s(2);

if ( (Ni>0) && (Nj>0) )
    mat = zeros( Ni, Nj );
else
    mat = [];
end

for i = 1:Ni

    [ s ] = fscanf( fid, ' %f ', [1, Nj ] );

    mat(i,:) = s;

end;

fclose( fid );

author: johannesgerer

repo: jburkardt-m


text:
function [ y, m, d, f ] = jed_to_ymdf_julian2 ( jed )

%*****************************************************************************80
%
%% JED_TO_YMDF_JULIAN2 converts a JED to a Julian YMDF date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real JED, the Julian Ephemeris Date.
%
%    Output, integer Y, M, D, real F, the YMDF date.
%

%
%  Check the input.
%
  ierror = jed_check ( jed );

  if ( ierror ~= 0 )
    y = -1;
    m = -1;
    d = -1;
    f = -1.0;
    return
  end

  j = floor ( jed + 0.5 );
  f = ( jed + 0.5 ) - j;

  jd = floor ( ( j + 1524 - 122.1 ) / 365.25 );
  je = floor ( 365.25 * jd );
  jg = floor ( ( j + 1524 - je ) / 30.6001 );
%
%  Now compute D, M and Y.
%
  d = j + 1524 - je - floor ( 30.6001 * jg );

  if ( jg <= 13 )
    m = jg - 1;
  else
    m = jg - 13;
  end

  if ( 2 < m )
    y = jd - 4716;
  else
    y = jd - 4715;
  end
%
%  Any year before 1 AD must be moved one year further back, since
%  this calendar does not include a year 0.
%
  y = y_astronomical_to_common ( y );

  return
end


text:
function [Xt,Yt] = elec_3d_2d(X,Y,Z,Zrad)

% elec_3d_2d - Project Cartesian 3D coordinates to a 2D plane.
%
% [Xt,Yt] = elec_3d_2d(X,Y,Z,Zrad)
%
% Project 3D electrode positions onto a 2D plane, using
% gnomonic projection, which draws a line from a point
% halfway between equator and south pole, through electrode,
% to a plane tangential to north pole.
%
% Given:  Set of 3D Cartesian coordinates (X,Y,Z with Z > 0)
%         X,Y midpoint at 0
%         Zrad is radius in Z
%
% See also elec_2d_3d.m
%

% $Revision: 1.1 $ $Date: 2009-04-28 22:13:54 $

% Licence:  GNU GPL, no implied or express warranties
% History:  10/1999, Chris Harvey
%           07/2001, Darren.Weber_at_radiology.ucsf.edu
%                    - using matrix algebra rather than indexed looping
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rad = Zrad / 2;

t = (Z + rad) * ( 1/rad );

Xt = X .* t;
Yt = Y .* t;

return

author: PatternRecognition

repo: OpenBMI


text:
% Two-input NAND gate sizing (GP)
% Boyd, Kim, Patil, and Horowitz, "Digital circuit optimization
% via geometric programming"
% Written for CVX by Almir Mutapcic 02/08/06
% (a figure is generated)
%
% This is an example taken directly from the paper:
%
%   Digital circuit optimization via geometrical programming
%   by Boyd, Kim, Patil, and Horowitz
%   Operations Research 53(6): 899-932, 2005.
%
% Solves the problem of choosing device widths w_i for the given
% NAND2 gate in order to achive minimum Elmore delay for different
% gate transitions, subject to limits on the device widths,
% gate area, power, and so on. The problem is a GP:
%
%   minimize   D = max( D_1, ..., D_k )  for k transitions
%       s.t.   w_min <= w <= w_max
%              A <= Amax, etc.
%
% where variables are widths w.
%
% This code is specific to the NAND2 gate shown in figure 19
% (page 926) of the paper. All the constraints and the objective
% are hard-coded for this particular circuit.

%********************************************************************
% problem data and hard-coded GP specs (evaluate all transitions)
%********************************************************************
N = 4;       % number of devices
Cload = 12;  % load capacitance
Vdd = 1.5;   % voltage

% device specs
NMOS = struct('R',0.4831, 'Cdb',0.6, 'Csb',0.6, 'Cgb',1, 'Cgs',1);
PMOS = struct('R',2*0.4831, 'Cdb',0.6, 'Csb',0.6, 'Cgb',1, 'Cgs',1);

% maximum area and power specification
wmin = 1;

% varying parameters for the tradeoff curve
Npoints = 25;
Amax = linspace(5,45,Npoints);
Dopt = [];

disp('Generating the optimal tradeoff curve...')
for k = 1:Npoints
    fprintf(1,'  Amax = %5.2f:', Amax(k));
    cvx_begin gp quiet
        % device width variables
        variable w(N)

        % device specs
        device(1:2) = PMOS; device(3:4) = NMOS;

        for num = 1:N
            device(num).R   = device(num).R/w(num); %#ok
            device(num).Cdb = device(num).Cdb*w(num); %#ok
            device(num).Csb = device(num).Csb*w(num); %#ok
            device(num).Cgb = device(num).Cgb*w(num); %#ok
            device(num).Cgs = device(num).Cgs*w(num); %#ok
        end

        % capacitances
        C1 = sum([device(1:3).Cdb]) + Cload;
        C2 = device(3).Csb + device(4).Cdb;

        % input capacitances
        Cin_A = sum([ device([2 3]).Cgb ]) + sum([ device([2 3]).Cgs ]);
        Cin_B = sum([ device([1 4]).Cgb ]) + sum([ device([1 4]).Cgs ]);

        % resistances
        R = [device.R]';

        % area definition
        area = sum(w);

        % delays and dissipated energies for all six possible transitions
        % transition 1 is A: 1->1, B: 1->0, Z: 0->1
        D1 = R(1)*(C1 + C2);
        E1 = (C1 + C2)*Vdd^2/2;
        % transition 2 is A: 1->0, B: 1->1, Z: 0->1
        D2 = R(2)*C1;
        E2 = C1*Vdd^2/2;
        % transition 3 is A: 1->0, B: 1->0, Z: 0->1
        % D3 = C1*R(1)*R(2)/(R(1) + R(2)); % not a posynomial
        E3 = C1*Vdd^2/2;
        % transition 4 is A: 1->1, B: 0->1, Z: 1->0
        D4 = C1*R(3) + R(4)*(C1 + C2);
        E4 = (C1 + C2)*Vdd^2/2;
        % transition 5 is A: 0->1, B: 1->1, Z: 1->0
        D5 = C1*(R(3) + R(4));
        E5 = (C1 + C2)*Vdd^2/2;
        % transition 6 is A: 0->1, B: 0->1, Z: 1->0
        D6 = C1*R(3) + R(4)*(C1 + C2);
        E6 = (C1 + C2)*Vdd^2/2;

        % objective is the worst-case delay
        minimize( max( [D1 D2 D4] ) )
        subject to
            area <= Amax(k); %#ok
            w >= wmin; %#ok
    cvx_end
    % display and store computed values
    fprintf(1,' delay = %3.2f\n',cvx_optval);
    Dopt = [Dopt cvx_optval]; %#ok
end

% plot the tradeoff curve
plot(Dopt,Amax);
xlabel('Dmin'); ylabel('Amax');
disp('Optimal tradeoff curve plotted.')

author: yu-jiang

repo: radpbook

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


text:
function varargout = pnopt_curvtrack( x, d, t, f_old, dg_x, smoothF, nonsmoothF, ...
  desc_param, xtol, maxIter )
% pnopt_curvtrack : Curve search for step that satisfies the Armijo condition
%
%   $Revision: 0.8.0 $  $Date: 2012/12/01 $
%
% ------------ Initialize ------------
  % Set line search parameters
  beta = 0.5;

  % Set termination flags
  FLAG_SUFFDESC = 1;
  FLAG_TOLX     = 2;
  FLAG_MAXFUNEV = 3;

  iter = 0;

  % ------------ Main Loop ------------
  while 1
    iter = iter + 1;

    % Evaluate trial point and function value.
    [ h_y, y ]  = nonsmoothF( x + t * d, t );
    if nargout > 6
      [ g_y, Dg_y, D2g_y ] = smoothF( y );
    else
      [ g_y, Dg_y ] = smoothF( y );
    end
    f_y = g_y + h_y;

    % Check termination criteria
    desc = 0.5 * norm( y - x ) ^2;
    if f_y < max( f_old ) + desc_param * t * desc    % Sufficient descent condition satisfied
      flag = FLAG_SUFFDESC;
      break
    elseif t <= xtol            % Step length too small
      flag = FLAG_TOLX;
      break
    elseif iter >= maxIter      % Too many line search iterations
      flag = FLAG_MAXFUNEV;
      break
    end

    % Backtrack if objective value not well-defined of function seems linear
    if isnan( f_y ) || isinf( f_y ) || abs( f_y - f_old(end) - t * dg_x ) <= 1e-9
      t = beta * t;
    % Safeguard quadratic interpolation
    else
      t_interp = - ( dg_x * t ^2) / ( 2 * ( f_y - f_old(end) - t * dg_x ) );
      if 0.1 <= t_interp || t_interp <= 0.9*t
        t = t_interp;
      else
        t = beta * t;
      end
    end
  end

  if nargout > 6
    varargout = { y, f_y, Dg_y, D2g_y, t, flag ,iter };
  else
    varargout = { y, f_y, Dg_y, t, flag ,iter };
  end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
function TessMat = in_tess_off(TessFile)
% IN_TESS_OFF: Read Geomview .off mesh files
%
% USAGE:  TessMat = in_tess_off(TessFile);
%
% INPUT:
%     - TessFile : full path to a tesselation file
% OUTPUT:
%     - TessMat:  Brainstorm tesselation structure
% FORMAT:
%     ASCII file with four blocks:
%     - Header line 1: "OFF"
%     - Header line 2: "nVertices nFaces ?"
%     - Vertices     : nVertices lines ("x y z"), values in miliimeters
%     - Faces        : nFaces lines ("nvert vertex1 vertex2 vertex3")
%
% SEE ALSO: in_tess

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2012

% Open tesselation file
fid = fopen(TessFile, 'r');
if fid < 0
    error(['Cannot open file ', TessFile])
end

% ===== READ HEADER =====
% Read format name
h1 = fgetl(fid);
if ~strcmpi(h1, 'OFF')
    error('Not a valid Geomview .off file.');
end
% Read number of vertices and faces
h2 = fgetl(fid);
dim = sscanf(h2, '%d', 2);
nVertices = dim(1);
nFaces    = dim(2);

% ===== READ MESH =====
% Read vertices
Vertices = double(fscanf(fid, '%f', [3 nVertices]));
% Go to next line
fgetl(fid);
% Read faces
Faces = double(fscanf(fid, '%f',[4 nFaces]) + 1);
Faces = Faces(2:4,:);
% Close file
fclose(fid);

%% ===== CONVERT IN BRAINSTORM FORMAT =====
TessMat.Vertices = Vertices';
TessMat.Faces    = Faces';

author: brainstorm-tools


text:
function sz = tsize(a,idx)
%TSIZE Tensor size of sptenmat.
%
%   D = TSIZE(X) returns the size of the tensor being stored as a
%   matrix.
%
%   M = TSIZE(X,DIM) returns the length of the dimension(s) specified
%   by DIM.  For example, SIZE(X,1) returns the size of the first
%   dimension of the tensor.
%
%   See also SPTENMAT, SPTENMAT/SIZE.
%
%MATLAB Tensor Toolbox.
%Copyright 2015, Sandia Corporation.

% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
% http://www.sandia.gov/~tgkolda/TensorToolbox.
% Copyright (2015) Sandia Corporation. Under the terms of Contract
% DE-AC04-94AL85000, there is a non-exclusive license for use of this
% work by or on behalf of the U.S. Government. Export of this data may
% require a license from the United States Government.
% The full license terms can be found in the file LICENSE.txt

if isempty(a.tsize)
    sz = [];
    return;
end

if exist('idx', 'var')
    sz = a.tsize(idx);
else
    sz = a.tsize;
end

return;

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function vw = drawROIs(vw, whichRois, method)
%
% vw = drawROIs(vw, <whichRois=get from view>, <method=get from vw>)
%
% Draw the specified ROIs on the view's current display.
%
% whichRois: vector of ROIs in the vw.ROIs field to draw. If omitted,
% parses the view's ROI settings (in vw.ui.showROIs)
% to decide the ROI order:
%   -2: draw all ROIs perimeter
%   -1: draw selected ROIs perimeter
%    0: don't draw any ROIs
%    1: draw selected ROI, boxes around each pixel
%    2: draw all ROIs, boxes around each pixel.
%
% whichRois can also be a cell array of ROI names
%
% method: can be one of the following strings or flags:
%   1 or 'perimeter', outline pixel perimeters
%   2 or 'boxes': boxes around each pixel
%   3 or 'patches': semi-transparent patches. (Matlab 7+ only)
%  If omitted, this is taken from the view struct, or (for hidden views)
%  defaults to 1.
%
%
% ras 05/06: updated to be much more efficient. Added patches option.
% TODO: re-implement the 'Filled Perimeter' option.
if ~exist('whichRois', 'var') | isempty(whichRois)
    % get from view
    if checkfields(vw, 'ui', 'showROIs')
        if length(vw.ui.showROIs) <= 1
            switch vw.ui.showROIs
                case -2, whichRois = 1:length(vw.ROIs);
                case -1, whichRois = vw.selectedROI;
                case 0, whichRois = [];
                otherwise, whichRois = vw.ui.showROIs;  % manually specify
            end
        else
            whichRois = vw.ui.showROIs;
        end
    else
        whichRois = 1:length(vw.ROIs);
    end
end

if ~exist('method', 'var') | isempty(method)
    % get from view
    if checkfields(vw, 'ui', 'roiDrawMethod')
        method = vw.ui.roiDrawMethod;

    else
        method = 'perimeter';

    end
end

% ML 7 check
if isequal(method, 'patches') | isequal(method, 3)
    ver = version;
    if str2num(ver(1)) < 7
        warning('Need MATLAB 7+ to draw patches ... drawing boxes instead')
        method = 'boxes';
    end
end

if iscell(whichRois)    % allow names of ROIs to be passed in
    roiNames = {vw.ROIs.name};
    whichRois = ismember(roiNames, whichRois);
    drawROIs(vw, whichRois, method);
    return
end

%%%%%ensure a field is set for the ROI object handles
N = length(vw.ROIs);
% if ~checkfields(vw, 'ui', 'roiHandles') | length(vw.ui.roiHandles) < N
%     vw.ui.roiHandles = cell(1, N);
% elseif N > length(vw.ui.roiHandles)
%     vw.ui.roiHandles = vw.ui.roiHandles(1:N);
% end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Draw the ROIs -- updated algorithms %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch vw.viewType
    case 'Inplane'
        vw = inplaneDrawROIs(vw, whichRois, method); % attached below

    case {'Volume' 'Gray'}
        % test for 3-axes view vs. single-orientation view:
        if checkfields(vw, 'ui', 'axiAxesHandle')
            vw = volumeDrawROIs(vw, 1, 1, whichRois, method);
            vw = volumeDrawROIs(vw, 2, 1, whichRois, method);
            vw = volumeDrawROIs(vw, 3, 1, whichRois, method);
        else
            ori = getCurSliceOri(vw);
			slice = viewGet(vw, 'Current Slice');
            vw = volumeDrawROIs(vw, ori, slice, whichRois, method);
        end

    case 'Flat'
        vw = flatDrawROIs(vw, whichRois, method); % attached below

end

return
% /------------------------------------------------------------------/ %

% /------------------------------------------------------------------/ %
function vw = inplaneDrawROIs(vw, whichRois, method);
% Draw specified ROIs according to the specified method, for inplanes
delete( findobj('Parent', vw.ui.mainAxisHandle, 'Type', 'line') )
if isempty(vw.ROIs) | whichRois==0, return; end
prefs.method = method;
for r = whichRois
    if isempty(vw.ROIs(r).coords), continue; end

    % get color
    color = vw.ROIs(r).color;
    if r==vw.selectedROI, color = viewGet(vw, 'selRoiColor'); end

    % delete old ROI hanldes
    try,   delete(vw.ui.roiHandles{r}); end
    % delete(findobj('Parent', vw.ui.mainAxisHandle, 'Type', 'line'))

    % find coords of ROI within the displayed slice
    pixels = coords2Montage(vw, vw.ROIs(r).coords);

    prefs.color = color;
    vw.ui.roiHandles{r} = outline(pixels, prefs);
end
return
% /------------------------------------------------------------------/ %

% /------------------------------------------------------------------/ %
function vw = volumeDrawROIs(vw, ori, slice, whichRois, method);
% Draw specified ROIs in an image of the specified orientation,
% according to the specified method, for volume/gray views
delete( findobj('Parent', vw.ui.mainAxisHandle, 'Type', 'line') )
if isempty(vw.ROIs) | whichRois==0, return; end
prefs.method = method;
if isempty(whichRois), return; end

for r = whichRois
    % get color
    color = vw.ROIs(r).color;
    if r==vw.selectedROI, color = viewGet(vw, 'selRoiColor'); end

    % delete old ROI hanldes
    try, delete(vw.ui.roiHandles{r}); end
    delete(findobj('Parent', vw.ui.mainAxisHandle, 'Type', 'line'))

    % find coords of ROI within the displayed slice
    coords = canOri2CurOri(vw, vw.ROIs(r).coords, ori);
	inSlice = find(coords(3,:)==slice);
	if isempty(inSlice)
		vw.ui.roiHandles{r} = [];
		continue
	end

    prefs.color = color;
    vw.ui.roiHandles{r} = outline(coords(1:2,inSlice), prefs);
end
% /------------------------------------------------------------------/ %

% /------------------------------------------------------------------/ %
function vw = flatDrawROIs(vw, whichRois, method);
% Draw specified ROIs according to the specified method, for flat views
prefs.method = method;
delete(findobj('Parent', vw.ui.mainAxisHandle, 'Type', 'line'))
if isempty(vw.ROIs) | whichRois==0, return; end
for r = whichRois
    % get color
    color = vw.ROIs(r).color;
    if r==vw.selectedROI, color = viewGet(vw, 'selRoiColor'); end

    % delete old ROI handles
    if checkfields(vw, 'ui', 'roiHandles') & length(vw.ui.roiHandles) >= r
        delete(vw.ui.roiHandles{r});
    end
    coords = vw.ROIs(r).coords;

	if isempty(coords), continue; end

    % The FLAT view has a 'rotateImageDegrees' field that
    % specifies a rotation angle for each slice (L or R).
    % If this is set, then we have to transform the ROIs by this amount as
    % well to make them register with the anatomy and functional data
    coords = rotateCoords(vw, coords);

    % restrict to coords in the slice, zoom
    slice = viewGet(vw, 'Current Slice');  % should really make this work for hidden views
    xZoom = vw.ui.zoom(1:2);
    yZoom = vw.ui.zoom(3:4);
    ok = find(coords(1,:)>=yZoom(1) & coords(1,:)<=yZoom(2) & ...
              coords(2,:)>=xZoom(1) & coords(2,:)<=xZoom(2) & ...
              coords(3,:)==slice);
    coords = coords(:,ok);

    prefs.color = color;
    prefs.method = method;
    vw.ui.roiHandles{r} = outline(coords, prefs);
end

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% solve_chol - solve linear equations from the Cholesky factorization.
% Solve A*X = B for X, where A is square, symmetric, positive definite. The
% input to the function is R the Cholesky decomposition of A and the matrix B.
% Example: X = solve_chol(chol(A),B);
%
% NOTE: The program code is written in the C language for efficiency and is
% contained in the file solve_chol.c, and should be compiled using matlabs mex
% facility. However, this file also contains a (less efficient) matlab
% implementation, supplied only as a help to people unfamiliar with mex. If
% the C code has been properly compiled and is available, it automatically
% takes precendence over the matlab code in this file.
%
% Copyright (c) by Carl Edward Rasmussen and Hannes Nickisch 2010-09-18.

function X = solve_chol(L, B)

if nargin ~= 2 || nargout > 1
  error('Wrong number of arguments.');
end

if size(L,1) ~= size(L,2) || size(L,1) ~= size(B,1)
  error('Wrong sizes of matrix arguments.');
end

X = L\(L'\B);

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function phi=pt2MaxEntropyCoords(x,v,algorithm,maxIter,epsScal)
%%PT2MAXENTROPYCOORDS Convert a point in a polygon to maximum entropy
%       coordinates. This is a type of barycentric coordinate system in 2D
%       that applies to points in polygons and is described in Section 7 of
%       [1]. These coordinates can be useful for Barycentric mapping and
%       for interpolation.
%
%INPUTS: x The 2XnumPts set of points to convert into maximum entropy
%          coordinates.
%        v The 2XnumBases set of vertices of the polygon defining the
%          coordinate system. The vertices can be given in clockwise or
%          counterclockwise order. It is assumed that the first vertex is
%          not repeated at the end.
% algorithm An optional parameter specifying the edge weight function to
%          use. Possible values are:
%          0 (The default if omitted or an empty matrix is passed) Use the
%            edge weight of Hormann and Sukumar, which is the first one
%            given in [1].
%          1 Use the second edge weight given in [1].
%  maxIter The maximum number of iterations of Newton's method to use. The
%          default if omitted or an empty matrix is passed is 100 and
%          generally, far fewer iterations are actually needed.
%  epsScal Convergence is determined when
%          abs(stepSize)<=epsScal*eps(lambda), where lambda is the 2X1
%          vector parameter being estimated from which the weights phi are
%          derived in [1]. The default if omitted or an empty matrix is
%          passed is 1.
%
%OUTPUTS: phi A numBasesXnumPts set of the maximum entropy coordinates of
%             the points given in x. Note that sum(phi(:,k))=1 for any k.
%             Points on the edge and outside of the polygon will typically
%             return a vector of NaNs.
%
%EXAMPLE 1:
%In this example, the coordinates of two points are found and the points
%are recreated from the phi values. The residual error of the reverse
%conversion is seen to be on the order of finite precision errors.
% numVertices=5;
% v=zeros(2,numVertices);
% v(:,1)=[-2;3];
% v(:,2)=[11;0];
% v(:,3)=[9;10];
% v(:,4)=[7;11];
% v(:,5)=[0;10];
% x=[[6;5],[8;2]];
% phi=pt2MaxEntropyCoords(x,v);
% xBack=barycentricCoords2Pt(phi,v);
% ResidualErr=max(max(abs(xBack-x)))
%
%EXAMPLE 2:
%This draws a convex polygon and a set of horizontal and vertical lines in
%the polygon. Then, the points of the lines are converted to maximum
%entropy coordinates and the vertices of the polygon are moved. After
%conversion back to Cartesian coordinates, one can see how the original
%grid in the polygon has been warped.
% numVertices=5;
% v=zeros(2,numVertices);
% v(:,1)=[-2;3];
% v(:,2)=[11;0];
% v(:,3)=[9;10];
% v(:,4)=[7;11];
% v(:,5)=[0;10];
%
% vertLines=-1:1:10;
% numVertLines=length(vertLines);
% horizLines=1:10;
% numHorizLines=length(horizLines);
% numLinPts=99;
% xVert=zeros(2,numLinPts,numVertLines);
% xHoriz=zeros(2,numLinPts,numHorizLines);
% for k=1:numVertLines
%     xVert(1,:,k)=vertLines(k);
%     xVert(2,:,k)=linspace(0,11,numLinPts);
%     sel=pointIsInPolygon(v,xVert(:,:,k),false)==0;
%     xVert(:,sel,k)=NaN;
% end
% for k=1:numHorizLines
%     xHoriz(1,:,k)=linspace(-2,11,numLinPts);
%     xHoriz(2,:,k)=horizLines(k);
%     sel=pointIsInPolygon(v,xHoriz(:,:,k),false)==0;
%     xHoriz(:,sel,k)=NaN;
% end
%
% figure(1)
% clf
% hold on
% plot([v(1,:),v(1,1)],[v(2,:),v(2,1)],'-c','linewidth',2)
% axis([-2, 12, 0, 12])
% for k=1:numVertLines
%     plot(xVert(1,:,k),xVert(2,:,k),'-r')
% end
% for k=1:numHorizLines
%     plot(xHoriz(1,:,k),xHoriz(2,:,k),'-b')
% end
% phiVert=pt2MaxEntropyCoords(xVert(:,:),v);
% phiHoriz=pt2MaxEntropyCoords(xHoriz(:,:),v);
%
% %Move the vertices, but keep the shape convex and the vertices still in
% %counterclockwise order.
% v(:,1)=[3;2];
% v(:,2)=[8;0];
% v(:,3)=[10;10];
% v(:,4)=[9;11];
% v(:,5)=[3;10];
%
% figure(2)
% clf
% hold on
% plot([v(1,:),v(1,1)],[v(2,:),v(2,1)],'-c','linewidth',2)
% axis([-2, 12, 0, 12])
% %Synthesize the corresponding points in the transformed coordinate
% %system.
% xTransVert=reshape(barycentricCoords2Pt(phiVert,v),[2,numLinPts,numVertLines]);
% xTransHoriz=reshape(barycentricCoords2Pt(phiHoriz,v),[2,numLinPts,numHorizLines]);
% for k=1:numVertLines
%     plot(xTransVert(1,:,k),xTransVert(2,:,k),'-r')
% end
% for k=1:numHorizLines
%     plot(xTransHoriz(1,:,k),xTransHoriz(2,:,k),'-b')
% end
%
%REFERENCES:
%[1] M. S. Floater, "Generalized barycentric coordinates and applications,"
%    Acta Numerica, vol. 24, pp. 161-214, 1 May 2015.
%
%September 2021 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

if(nargin<5||isempty(epsScal))
    epsScal=1;
end

if(nargin<4||isempty(maxIter))
    maxIter=100;
end

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

numPts=size(x,2);
numBases=size(v,2);
phi=zeros(numBases,numPts);
for curPt=1:numPts
    if(any(isnan(x(:,curPt))))
       phi(:,curPt)=NaN;
       continue;
    end

    diffs=bsxfun(@minus,x(:,curPt),v);
    mags=sqrt(sum(diffs.*diffs,1));

    rhoVals=zeros(numBases,1);
    if(algorithm==0)
        for k=1:(numBases-1)
            rhoVals(k)=mags(k)+mags(k+1)-norm(v(:,k)-v(:,k+1));
        end
        rhoVals(numBases)=mags(numBases)+mags(1)-norm(v(:,numBases)-v(:,1));
    elseif(algorithm==1)
        for k=1:(numBases-1)
            rhoVals(k)=mags(k)*mags(k+1)+dot(diffs(:,k),diffs(:,k+1));
        end
        rhoVals(numBases)=mags(numBases)*mags(1)+dot(diffs(:,numBases),diffs(:,1));
    else
       error('Unknown algorithm specified.')
    end

    piVals=zeros(1,numBases);
    piVals(1)=prod(rhoVals(2:(numBases-1)));
    for k=2:(numBases-1)
        piVals(k)=prod(rhoVals([1:(k-2),(k+1):numBases]));
    end
    piVals(numBases)=prod(rhoVals(1:(numBases-2)));

    mi=piVals/sum(piVals);
    d=-diffs;

    %lambda can start uniform.
    lambda=[1/2;1/2];
    %Use Newton's method.
    for curIter=1:maxIter
        [F,H]=getFAndH(lambda,d,mi);

        if(any(~isfinite(H(:))))
           lambda=[NaN;NaN];
           break;
        end

        stepVal=pinv(H)*F;

        if(all(abs(stepVal)<=epsScal*eps(lambda)))
            %Convergence to some multiple of finite precision limits.
            break;
        end
        lambda=lambda-stepVal;
    end

    wi=mi.*exp(sum(bsxfun(@times,lambda,d),1));
    phi(:,curPt)=wi/sum(wi);
end
end

function [gradF,H]=getFAndH(lambda,d,mi)
%%GETFANDH Get the the gradient and the Hessian of the cost function.

mExpdTerms=bsxfun(@times,mi.*exp(sum(bsxfun(@times,lambda,d),1)),d);

%The gradient.
gradF=sum(mExpdTerms,2);
H=zeros(2,2);
H(1,1)=sum(mExpdTerms(1,:).*d(1,:));
H(2,1)=sum(mExpdTerms(1,:).*d(2,:));
H(1,2)=H(2,1);
H(2,2)=sum(mExpdTerms(2,:).*d(2,:));

end

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

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
% USERDATAGRAPH  User data for SLAMTB_GRAPH.
%   Edit this script to enter the information you need for SLAM. Variable
%   names and comments should make this file easy to understand. Follow
%   these guidelines:
%
%   * Specify site and estimation details for the current run in ExpOpt.
%   * Specify sampling time and start and end frames in Tim.
%   * Use as many robots and sensors as you wish with Robot{} and Sensor{}.
%   * Assign sensors to robots via Sensor{i}.robot.
%   * Use field Sensor{i}.distortion for radial distortion parameters if
%   desired.
%   * Use the field Opt.map.numLmk and .lmkSize to specify the maximum
%   number of landmarks that the SLAM map must support.
%   * Use Opt.init.initType to select the type of landmarks to use. Try
%   with one in this list:
%       'eucPnt'. -- only this one is supported so far.
%   * Use World.points to create artificial worlds of points.
%   Check function THICKCLOISTER.
%
%   See further comments within the file for more detailed information.
%
%   NOTE: You can have multiple copies of this file with different names to
%   store different simulation conditions. Just modify the call in
%   SLAMTB_GRAPH to point to the particular 'USERDATAGRAPH' file you want.
%
%   See also SLAMTB_GRAPH, EULERANGLES, THICKCLOISTER.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.
%   Copyright 2015-     Joan Sola @ IRI-CSIC-UPC.

% Time variables
%   - sampling time, first and last frames
Time = struct(...
  'dt',                   .1,...          % sampling time, seconds
  'firstFrame',           1,...           % first frame #
  'lastFrame',            800);           % last frame #

% Simulated world
%   - Simulation landmark sets, playground dimensions
World = struct(...
  'points',           thickCloister(-6,6,-6,6,1,7),... % 3d point landmarks - see THICKCLOISTER.
  'segments',         []);  % 3D segments - see HOUSE.

% Robot things with their controls
%   - each robot's type and initial configuration, and controls.
%   - motion models (add new model strings if you need more):
%       'constVel'    6D Constant velocity model
%       'odometry'    6D Odometry model
%       'inertial'    6D IMU-based model
%   - See EULERANGLES for orientations specifications.
Robot{1} = struct(...                  % ODOMETRY EXAMPLE
  'id',                 1,...           % robot identifier
  'name',               'Dala',...      % robot name
  'type',               'atrv',...      % type of robot
  'motion',             'odometry',...  % motion model
  'position',           [0;-5;0],...     % robot position in map
  'orientationDegrees', [0;0;0],...     % orientation, in degrees, [roll; pitch; yaw].
  'positionStd',        [0;0;0],...     % position error, std
  'orientationStd',     [0;0;0],...     % orient. error, std, in degrees
  'dx',                 [.08;0;0],...     % position increment
  'daDegrees',          [0;0;.9],...     % angle increment, degrees
  'dxStd',              0.03*[1;1;1],...  % odo linear error std
  'daStd',              0.3*[1;1;1]);      % odo ang error std, degrees

% Robot{2} = struct(...                  % ODOMETRY EXAMPLE
%   'id',                 2,...           % robot identifier
%   'name',               'Dala',...      % robot name
%   'type',               'atrv',...      % type of robot
%   'motion',             'odometry',...  % motion model
%   'position',           [0;5;0],...     % robot position in map
%   'orientationDegrees', [0;0;180],...     % orientation, in degrees, [roll; pitch; yaw].
%   'positionStd',        [0;0;0],...     % position error, std
%   'orientationStd',     [0;0;0],...     % orient. error, std, in degrees
%   'dx',                 [.08;0;0],...     % position increment
%   'daDegrees',          [0;0;.9],...     % angle increment, degrees
%   'dxStd',              0.005*[1;1;1],...  % odo linear error std
%   'daStd',              0.05*[1;1;1]);      % odo ang error std, degrees

% Sensor things
%   - each sensor's type and parameters, noise, non-measurable prior.
%   - Sensor types (add new type strings if you need more):
%       'pinHole'       Pin-hole camera
%       'pinHoleDepth'  Pin-hole camera with depth measurement (RGBD)
%   - See EULERANGLES for orientations specifications.
% Sensor{1} = struct(...
%   'id',                 1,...           % sensor identifier
%   'name',               'Micropix',...  % sensor name
%   'type',               'pinHole',...   % type of sensor
%   'robot',              1,...           % robot where it is mounted
%   'position',           [0;0;.6],...    % position in robot
%   'orientationDegrees', [-90;0;-90],... % orientation in robot, [roll; pitch; yaw]
%   'positionStd',        [0;0;0],...     % position error std
%   'orientationStd',     [0;0;0],...     % orient. error std
%   'imageSize',          [640;480],...   % image size
%   'pixErrorStd',        1.0,...         % pixel error std
%   'intrinsic',          [320;240;320;320],... % intrinsic params [u0 v0 au av]
%   'distortion',         [-0.3;0.1],...          % distortion params
%   'frameInMap',         false,...       % add sensor frame in slam map?
%   'imGrid',               struct(...      % grid for Active Search
%     'numCells',         [8;6],...         % number of H and V grid cells
%     'skipOuter',        true));           % skip outer cells for initialization?

Sensor{1} = struct(...
  'id',                 1,...           % sensor identifier
  'name',               'Micropix',...  % sensor name
  'type',               'pinHoleDepth',...   % type of sensor
  'robot',              1,...           % robot where it is mounted
  'position',           [0;0;.6],...    % position in robot
  'orientationDegrees', [-90;0;-90],... % orientation in robot, [roll; pitch; yaw]
  'positionStd',        [0;0;0],...     % position error std
  'orientationStd',     [0;0;0],...     % orient. error std
  'imageSize',          [640;480],...   % image size
  'pixErrorStd',        1.0,...         % pixel error std [pixels]
  'depthErrorStd',      0.5,...         % depth error std [m]
  'intrinsic',          [320;240;320;320],... % intrinsic params [u0 v0 au av]
  'distortion',         [],...          % distortion params, e.g. [-0.3;0.1]
  'frameInMap',         false,...       % add sensor frame in slam map?
  'imGrid',             struct(...      % grid for Active Search
    'numCells',         [8;6],...         % number of H and V grid cells
    'skipOuter',        true));           % skip outer cells for initialization?

% Estimation options
Opt = struct(...
  'map',              struct(...    % options for the map
    'type',           'graph',...    % type of map {'ekf','graph'}
    'numLmks',        73,...         % number of 3d landmarks
    'lmkSize',        3,...          % Size of landmark state
    'lmkDSize',       3,...          % Size of lmk error state
    'numFrames',      25,...         % number of frames in graph
    'kfrmPeriod',     20),...        % period between keyframes
  'solver',           struct(...    % graph solver
    'decomposition',  'QR',...        % Matrix decomposition {'Cholesky', 'QR', 'Schur'}
    'niterations',    1,...           % number of iterations
    'target_dres',    1e-2,...        % exit criterion for residual variation
    'target_res',     1e-3),...       % exit criterion for current residual
  'correct',          struct(...    % options for lmk correction
    'reprojectLmks',  false,...       % reproject lmks after active search?
    'reparametrize',  true,...       % reparametrize lmk?
    'nUpdates',       10,...          % max simultaneus updates
    'MD2th',          9,...          % Threshold on Mahalanobis distance squared
    'linTestIdp',     0.1,...        % threshold on IDP linearity test
    'lines',          struct(...     % options for line corrections
      'innType',      'ortDst',...    % innovation type for lines
      'extPolicy',    false,...       % line extending policy ?
      'extSwitch',    10)),...        % extension policy switch point in pixels
  'init',             struct(...    % Options for initialization
    'nbrInits',       [5 3],...      % number of inits [firstFrame, otherFrames]
    'initType',       'eucPnt',...   % Type of lmk to use for init
    'idpPnt',         struct(...     % options for lmk initialization
      'nonObsMean',   .1,...         % mean of non obs
      'nonObsStd',    .5),...         % std of non obs
    'plkLin',         struct(...     % opt. for Plucker and anchored Plucker lines init
      'nonObsMean',   [.1;0],...      % mean of non obs
      'nonObsStd',    [.25;1])),...   % std of non obs
  'obs',              struct(...    % Observation options
    'lines',          struct(...     % lines options
      'minLength',    20)));          % minimum segment length

% Simulation options
%   - random
SimOpt = struct(...
  'random',           struct(...      % random generator options
    'newSeed',        false,...         % select new random seed?
    'fixedSeed',      209065,...       % random seed for non-random runs
    'seed',           []),...          % current seed
  'obs',              Opt.obs);       % Observation options

% Figure options
%   - view, projection, video, ellipses.
%   - figure projections - mapProj:
%       'persp'     Perspective
%       'ortho'     Orthographic
%   - 3D figure views - mapView - see MAPOBSERVER.
%       [a,e,f]     Custom azimuth/elevation/FOV vector. Distance automatic
%       [a,e,f,d]   custom az/el/fov/distance vector.
%   - 3D figure predefined views (edit mapObserver.m to create/edit views):
%       'top'       Top view
%       'side'      Side view
%       'view'      Generic view
%       'normal'    Normal view
%   - objects colors - two options for color specification:
%       'rgbcmykw'  1-char predifined Matlab colors
%       [r g b]     RGB color vector. [0 0 0] is black, [1 1 1] is white.
FigOpt = struct(...
  'renderer',       'zbuffer',...    % renderer
  'rendPeriod',     10,...           % frames to skip for faster processing
  'createVideo',    false,...       % create video sequences?
  'map',            struct(...      % map figure options
    'size',         [320 240],...   % map figure size
    'lims',         struct(...      % playground limits
      'xMin',            -10,...
      'xMax',             10,...
      'yMin',            -10,...
      'yMax',             10,...
      'zMin',            -10,...
      'zMax',             10),...
    'proj',         'persp',...     % projection of the 3d figure {'persp', 'ortho'}
    'view',         'view',...      % viewpoint of the 3d figure [30 45 40 20]
    'orbit',        [0 0],...       % Azimuth and Elevation orbit angle increments - use to animate figure
    'showSimLmk',   false,...        % show simulated landmarks?
    'showLmkId',    true,...        % Show landmark id?
    'showEllip',    true,...        % show ellipsoids?
    'showMotFac',   true,...        % show motion factors?
    'showMeaFac',   false,...       % show measurement factors?
    'colors',       struct(...      % map figure colors
      'border',     [1 1 1],...      %   [r g b]
      'axes',       [0 0 0],...      % with:
      'bckgnd',     [1 1 1],...      %   [0 0 0] black
      'simLmk',     .3*[1 1 1],...   %   [1 1 1] white
      'defPnt',     struct(...       % euclidean point colors
        'mean',     'b',...           % mean dot
        'ellip',    [.7 .7 1]),...    % ellipsoid
      'othPnt',     struct(...       % other point colors
        'mean',     'r',...           % mean dot
        'ellip',    [1 .7 .7]),...    % ellipsoid
      'defLin',     struct(...       % Plucker line colors
        'mean',     [0 .8 0],...      % mean line
        'ellip',    [.6 1 .6]),...    % ellipsoid
      'othLin',     struct(...       % Plucker line colors
        'mean',     [.8 0 0],...      % mean line
        'ellip',    [1 .6 .6]),...    % ellipsoid
      'graph',      struct(...       % graph
        'frame',    'b',...           % frame
        'motion',   'b',...           % motion factors
        'meas',     [1 .8 .6]),...          % measurement factors
      'simu',       'b',...          %   or 'r', 'b', etc.
      'est',        'g',...          % estimated robots and sensors
      'ground',     [.8 .8 .8],...   % simulated robots and sensors
      'label',      [.0 .5 0])),...  % landmark ID labels
  'sensor',         struct(...      % sensor figures options
    'size',         [320 240],...    % sensor figure size
    'showEllip',    false,...        % show ellipses?
    'colors',       struct(...       % Sensor figure colors:
      'border',     .8*[1 1 1],...    %
      'axes',       [0 0 0],...       %
      'bckgnd',     [1 1 1],...       %
      'raw',        .3*[1 1 1],...    %
      'defPnt',     struct(...       % Default point colors
        'updated',  'c',...           % updated
        'predicted','b'),...          % predicted
      'othPnt',     struct(...       % other point colors
        'updated',  'r',...           % updated
        'predicted','m'),...          % predicted
      'defLin',     struct(...       % Default line colors
        'meas',     'b',...           % measurement
        'mean',     'g',...           % mean line
        'ellip',    'y'),...          % ellipsoid
      'othLin',     struct(...       % other line colors
        'meas',     'b',...           % measurement
        'mean',     'm',...           % mean line
        'ellip',    'r'),...          % ellipsoid
      'label',      [.5 .5 .5])));    %

% Experiment options
%   - site name, series gathered, estimation run number
ExpOpt = struct(...
  'root',         '~/SLAM/',...         % root directory
  'site',         'simu',...            % Name of the site
  'dataRun',      1,...                 % Run # on this site
  'estimateRun',  1,...                 % slam run for data and site
  'lmkTypes',     Opt.init.initType,... % types of landmarks used
  'sensingType',  'mono',...            % sensing mode
  'mappingType',  'single');            % mapping mode

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB 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.
%
%   SLAMTB 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 SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
% APM Retrieve an Option Value (Tag)
%
% response = apm_tag(server,app,name)
%
% Function apm_tag retrieves an option specification from
%   the APM server with the following arguments:
%
%   server = address of server
%      app = application name
%     name = option name
% response = option value
%
% Either global options (NLC.option_name) or parameter/variable
%   options (NAME.option_name) are able to be retrieved with this
%   function. A full list of available options with their default
%   values are available here:
%
% Global Options:
%    http://apmonitor.com/wiki/index.php/Main/DbsGlobal
%
% Variable Options:
%    http://apmonitor.com/wiki/index.php/Main/DbsVariable
%
function response = apm_tag(server,app,name)
   app = lower(deblank(app));

    % Web-server URL base
    url_base = [deblank(server) '/online/get_tag.php'];

    % Send request to web-server
    params = ['?p=' urlencode(app) '&n=' urlencode(name)];
    url = [url_base params];
    % Send request to web-server
    response = str2num(urlread_apm(url));

author: APMonitor

repo: arduino

save_path: github-repos/MATLAB/APMonitor-arduino


text:
function idx = GetParameterIndex(modelname, parametername)
%
% function idx = GetParameterIndex(modelname, parametername)
%
% Given an input modelname and parametername, this function returns the index
% of the parameter within the model.
%
% author: Gary Hui Zhang (gary.zhang@ucl.ac.uk)
%

strings = GetParameterStrings(modelname);

for i=1:length(strings)
    if (strcmp(strings(i), parametername))
        idx = i;
        return;
    end
end

idx = -1;

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function [ a_lu, info ] = r8cbb_fa ( n1, n2, ml, mu, a )

%*****************************************************************************80
%
%% R8CBB_FA factors a R8CBB matrix.
%
%  Discussion:
%
%    Note that in C++ and FORTRAN, we can look at A as an abstract
%    vector, but then look at parts of A as storing a two dimensional
%    array.  MATLAB assigns an inherent dimensionality to a data object,
%    and gets very unhappy when you try to manipulate the data yourself.
%    This means that the MATLAB implementation of this routine requires
%    the use of temporary 2D arrays.
%
%    The R8CBB storage format is for a compressed border banded matrix.
%    Such a matrix has the logical form:
%
%      A1 | A2
%      ---+---
%      A3 | A4
%
%    with A1 a (usually large) N1 by N1 banded matrix, while A2, A3 and A4
%    are dense rectangular matrices of orders N1 by N2, N2 by N1, and N2 by N2,
%    respectively.
%
%    The R8CBB format is the same as the DBB format, except that the banded
%    matrix A1 is stored in compressed band form rather than standard
%    banded form.  In other words, we do not include the extra room
%    set aside for fill in during pivoting.
%
%    A should be defined as a vector.  The user must then store
%    the entries of the four blocks of the matrix into the vector A.
%    Each block is stored by columns.
%
%    A1, the banded portion of the matrix, is stored in
%    the first (ML+MU+1)*N1 entries of A, using the obvious variant
%    of the LINPACK general band format.
%
%    The following formulas should be used to determine how to store
%    the entry corresponding to row I and column J in the original matrix:
%
%    Entries of A1:
%
%      1 <= I <= N1, 1 <= J <= N1, (J-I) <= MU and (I-J) <= ML.
%
%      Store the I, J entry into location
%      (I-J+MU+1)+(J-1)*(ML+MU+1).
%
%    Entries of A2:
%
%      1 <= I <= N1, N1+1 <= J <= N1+N2.
%
%      Store the I, J entry into location
%      (ML+MU+1)*N1+(J-N1-1)*N1+I.
%
%    Entries of A3:
%
%      N1+1 <= I <= N1+N2, 1 <= J <= N1.
%
%      Store the I, J entry into location
%      (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
%
%    Entries of A4:
%
%      N1+1 <= I <= N1+N2, N1+1 <= J <= N1+N2
%
%      Store the I, J entry into location
%      (ML+MU+1)*N1+N1*N2+(J-1)*N2+(I-N1).
%      (same formula used for A3).
%
%
%    Once the matrix has been factored by R8CBB_FA, R8CBB_SL may be called
%    to solve linear systems involving the matrix.
%
%    R8CBB_FA uses special non-pivoting versions of LINPACK routines to
%    carry out the factorization.  The special version of the banded
%    LINPACK solver also results in a space saving, since no entries
%    need be set aside for fill in due to pivoting.
%
%    The linear system must be border banded, of the form:
%
%      ( A1 A2 ) (X1) = (B1)
%      ( A3 A4 ) (X2)   (B2)
%
%    where A1 is a (usually big) banded square matrix, A2 and A3 are
%    column and row strips which may be nonzero, and A4 is a dense
%    square matrix.
%
%    The algorithm rewrites the system as:
%
%         X1 + inverse(A1) A2 X2 = inverse(A1) B1
%
%      A3 X1 +             A4 X2 = B2
%
%    and then rewrites the second equation as
%
%      ( A4 - A3 inverse(A1) A2 ) X2 = B2 - A3 inverse(A1) B1
%
%    The algorithm will certainly fail if the matrix A1 is singular,
%    or requires pivoting.  The algorithm will also fail if the A4 matrix,
%    as modified during the process, is singular, or requires pivoting.
%    All these possibilities are in addition to the failure that will
%    if the total matrix A is singular.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    24 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N1, N2, the order of the banded and dense blocks.
%    N1 and N2 must be nonnegative, and at least one must be positive.
%
%    Input, integer ML, MU, the lower and upper bandwidths.
%    ML and MU must be nonnegative, and no greater than N1-1.
%
%    Input, real A( (ML+MU+1)*N1 + 2*N1*N2 + N2*N2), the compact
%    border-banded coefficient matrix.
%
%    Output, real A_LU( (ML+MU+1)*N1 + 2*N1*N2 + N2*N2).
%    information describing a partial factorization
%    of the original coefficient matrix.
%
%    Output, integer INFO, singularity flag.
%    0, no singularity detected.
%    nonzero, the factorization failed on the INFO-th step.
%
  nband = (ml+mu+1)*n1;

  a_lu(1:nband+2*n1*n2+n2*n2) = a(1:nband+2*n1*n2+n2*n2);
%
%  Factor the A1 band matrix, overwriting A1 by its factors.
%
  if ( 0 < n1 )

    a1(1:ml+mu+1,1:n1) = r8vec_to_r8cb ( n1, n1, ml, mu, a(1:nband) );

    [ a1_lu, info ] = r8cb_np_fa ( n1, ml, mu, a1 );

    if ( info ~= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8CBB_FA - Fatal error!\n' );
      fprintf ( 1, '  R8CB_NP_FA returned INFO = %d\n', info );
      fprintf ( 1, '  Factoring failed for column INFO.\n' );
      fprintf ( 1, '  The band matrix A1 is singular.\n' );
      fprintf ( 1, '  This algorithm cannot continue!\n' );
      error ( 'R8CBB_FA - Fatal error!' );
    end

    a_lu(1:nband) = r8cb_to_r8vec ( n1, n1, ml, mu, a1_lu );

  end

  if ( 0 < n1 & 0 < n2 )
%
%  Set A2 := -inverse(A1) * A2.
%
    job = 0;

    for j = 1 : n2
      b2(1:n1) = -a(nband+(j-1)*n1+1:nband+(j-1)*n1+n1);
      x2 = r8cb_np_sl ( n1, ml, mu, a1_lu, b2, job );
      a_lu(nband+(j-1)*n1+1:nband+(j-1)*n1+n1) = x2(1:n1);
    end
%
%  Set A4 := A4 + A3*A2
%
    for i = 1 : n2
      for j = 1 : n1
        ij = nband + n1*n2 + (j-1)*n2 + i;
        for k = 1 : n2
          ik = nband + 2*n1*n2 + (k-1)*n2 + i;
          jk = nband + (k-1)*n1 + j;
          a_lu(ik) = a_lu(ik) + a_lu(ij) * a_lu(jk);
        end
      end
    end

  end
%
%  Factor A4.
%
  if ( 0 < n2 )

    a4 = r8vec_to_r8ge ( n2, n2, a_lu(nband+2*n1*n2+1:nband+2*n1*n2+n2*n2) );

    [ a4_lu, info ] = r8ge_np_fa ( n2, a4 );

    if ( info ~= 0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'R8CBB_FA - Fatal error!\n' );
      fprintf ( 1, '  R8GE_NP_FA returned INFO = %d\n',info );
      fprintf ( 1, '  This indicates singularity in column INFO\n' );
      info = n1 + info;
      fprintf ( 1, '  of the A4 submatrix, which is column %d\n', info );
      fprintf ( 1, '  of the full matrix.\n' );
      fprintf ( 1, '\n' );
      fprintf ( 1, '  It is possible that the full matrix is \n' );
      fprintf ( 1, '  nonsingular, but the algorithm R8CBB_FA may\n' );
      fprintf ( 1, '  not be used for this matrix.\n' );
      return
    end

    a_lu(nband+2*n1*n2+1:nband+2*n1*n2+n2*n2) = r8ge_to_r8vec ( n2, n2, a4_lu );

  end

  return
end


text:
function [matching,optimal_cost]=hungarian(A, options)
%HUNGARIAN Solve the Assignment problem using the Hungarian method.
%
%[matching,optimal_cost]=hungarian(A, options)
%A - a square cost matrix.
%matching - the optimal assignment (from column to row)
%optimal_cost - the cost of the optimal assignment.
if nargin < 2
    options = struct;
end
if ~isfield(options, 'verbose')
    options.verbose = 0;
end
[matching,optimal_cost] = mex_hungarian(A, options.verbose);
end

author: indigits

repo: sparse-plex


text:
%tstoolbox/mex/range_search
%   Syntax:
%
%     * [count, neighbors] = range_search(pointset, atria, query_points,
%       r)
%     * [count, neighbors] = range_search(pointset, atria, query_indices,
%       r, exclude)
%
%   Input arguments:
%
%     * pointset - a N by D double matrix containing the coordinates of
%       the point set, organized as N points of dimension D
%     * atria - output of (cf. Section )nn_prepare for pointset
%     * query_points - a R by D double matrix containing the coordinates
%       of the query points, organized as R points of dimension D
%     * query_indices - query points are taken out of the pointset,
%       query_indices is a vector of length R which contains the indices
%       of the query points
%     * r - range or search radius (r > 0)
%     * exclude - in case the query points are taken out of the pointset,
%       exclude specifies a range of indices which are omitted from
%       search. For example if the index of the query point is 124 and
%       exclude is set to 3, points with indices 121 to 127 are omitted
%       from search. Using exclude = 0 means: exclude self-matches
%
%   Output arguments:
%
%     * count - a vector of length R contains the number of points within
%       distance r to the corresponding query point
%     * neighbors - a Matlab cell structure of size R by 2 which contains
%       vectors of indices and vectors of distances to the neighbors for
%       each given query point. This output argument can not be stored in
%       a standard Matlab matrix because the number of neighbors within
%       distance r is not the same for all query points. The vectors if
%       indices and distances for one query point have exactly the length
%       that is given in count. The values in the distances vectors are
%       not sorted..
%
% Copyright 1997-2001 DPI Goettingen, License http://www.physik3.gwdg.de/tstool/gpl.txt

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa