text:
function result = HIDC_enterFunc(numIdentical,numAlike,numDistinct,numIter,inputPath)
    %% ======================================================================
    %%STEP 1: load the data
    fprintf('start load the data...\n');
%     fprintf(inputPath);
%     fprintf('\n');
    [TrainData, TestData, TrainLabel, TestLabel, numSource, numTrain, numTarget, numTest] = HIDC_loadData(inputPath);
      TrainData = double(TrainData);
      TestData = double(TestData);
      TrainLabel = double(TrainLabel);
      TestLabel = double(TestLabel);

      TrainData = sparse(TrainData);
      TestData = sparse(TestData);

    fprintf('training the model and testing...\n');
    [Results_TTL, Gt_TTL] = GenerativeTriTL(TrainData, TestData, TrainLabel, TestLabel, numIdentical, numAlike, numDistinct, numIter, numSource, numTrain, numTarget, numTest);

    result = Results_TTL(:,2)';
    result = result(end);
end

author: FuzhenZhuang

repo: Transfer-Learning-Toolkit


text:
function triangle_wandzura_rule_test01 ( )

%*****************************************************************************80
%
%% TRIANGLE_WANDZURA_RULE_TEST01 tests WANDZURA_RULE_NUM, WANDZURA_DEGREE, WANDZURA_ORDER_NUM.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 December 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TRIANGLE_WANDZURA_RULE_TEST01\n' );
  fprintf ( 1, '  WANDZURA_RULE_NUM returns the number of rules;\n' );
  fprintf ( 1, '  WANDZURA_DEGREE returns the degree of a rule;\n' );
  fprintf ( 1, '  WANDZURA_ORDER_NUM returns the order of a rule.\n' );

  rule_num = wandzura_rule_num ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of available rules = %d\n', rule_num );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      Rule    Degree     Order\n' );
  fprintf ( 1, '\n' );

  for rule = 1 : rule_num
    order_num = wandzura_order_num ( rule );
    degree = wandzura_degree ( rule );
    fprintf ( 1, '  %8d  %8d  %8d\n', rule, degree, order_num );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function x = singlehex2num(s)
%SINGLEHEX2NUM Convert single precision IEEE hexadecimal string to number.
%   SINGLEHEX2NUM(S), where S is a 8 character string containing
%   a hexadecimal number, returns a double type number
%   equal to the IEEE single precision
%   floating point number it represents.  Fewer than 8
%   characters are padded on the right with zeros.
%
%   If S is a character array, each row is interpreted as a single
%   precision number (and returned as a double).
%
%   NaNs, infinities and denorms are handled correctly.
%
%   Example:
%       hexsingle2num('40490fdb') returns Pi.
%       hexsingle2num('bf8') returns -1.
%
%   See also HEX2NUM.

% Based on Matlab's hex2num.
% Note: IEEE Standard 754 for floating point numbers
%
%  Floating point numbers are represented as:
%  x = +/- (1+f)*2^e
%
%  doubles: 64 bits
%           Bit 63       (1 bit)  = sign (0=positive, 1=negative)
%           Bit 62 to 52 (11 bits)= exponent biased by 1023
%           Bit 51 to 0  (52 bits)= fraction f of the number 1.f
%  singles: 32 bits
%           Bit 31       (1 bit)  = sign (0=positive, 1=negative)
%           Bit 30 to 23 (8 bits) = exponent biased by 127
%           Bit 22 to 0  (23 bits)= fraction f of the number 1.f
%
% Original file hexsingle2num from Mark Lubinski
% Changed on 19-may-05 by Matthias Noell: denormalized power set 2^-126

if iscellstr(s), s = char(s); end
if ~ischar(s)
    error('Input to hexsingle2num must be a string.')
end
if isempty(s), x = []; return, end

[row,col] = size(s);
blanks = find(s==' '); % Find the blanks at the end
if ~isempty(blanks), s(blanks) = '0'; end % Zero pad the shorter hex numbers.

% Convert characters to numeric digits.
% More than 8 characters are ignored
% For double: d = zeros(row,16);
d = zeros(row,8);
d(:,1:col) = abs(lower(s)) - '0';
d = d + ('0'+10-'a').*(d>9);
neg = d(:,1) > 7;
d(:,1) = d(:,1)-8*neg;

if any(d > 15) | any(d < 0)
    error('Input string to hexsingle2num should have just 0-9, a-f, or A-F.')
end

% Floating point exponent.
% For double: e = 16*(16*(d(:,1)-4) + d(:,2)) + d(:,3) + 1;
% For double: e = 256*d(:,1) + 16*d(:,2) + d(:,3) - 1023;
expBit = (d(:,3) > 7);
e = 32*d(:,1) + 2*d(:,2) + expBit - 127;
d(:,3) = d(:,3)-8*expBit;  % Remove most sig. bit of d(:,3) which belongs to exponent

% Floating point fraction.
% For double: sixteens = [16;256;4096;65536;1048576;16777216;268435456];
% For double: sixteens2 = 268435456*sixteens(1:6);
% For double: multiplier = 1./[sixteens;sixteens2];
% For double: f = d(:,4:16)*multiplier;
sixteens = [16;256;4096;65536;1048576;16777216];
multiplier = 2./[sixteens];
f = d(:,3:8)*multiplier;

x = zeros(row,1);
% Scale the fraction by 2 to the exponent.
% For double: overinf = find((e>1023) & (f==0));
overinf = find((e>127) & (f==0));
if ~isempty(overinf), x(overinf) = inf; end

% For double: overNaN = find((e>1023) & (f~=0));
overNaN = find((e>127) & (f~=0));
if ~isempty(overNaN), x(overNaN) = NaN; end

% For double: underflow = find(e<-1022);
underflow = find(e<-126);
if ~isempty(underflow), x(underflow) = pow2(f(underflow),-126); end

% For double: allothers = find((e<=1023) & (e>=-1022));
allothers = find((e<=127) & (e>=-126));
if ~isempty(allothers), x(allothers) = pow2(1+f(allothers),e(allothers)); end

negatives = find(neg);
if ~isempty(negatives), x(negatives) = -x(negatives); end

author: Sable

repo: mcbench-benchmarks


text:
function aistats2012_mohsst5_results()

%
% DATA
%

disp('--- Results for UNIFORM case ---')
show_results('data', 1);
disp('Kronecker GP')
show_results('gpkron', 1, 'results_aistats2012_mohsst5_gpkron_uniform_20110528');
disp('Kronecker GP (2 covfuncs)')
show_results('gpkron2', 1, 'results_aistats2012_mohsst5_gpkron-cov2_uniform_20111010');
disp('VB PCA')
show_results('pca', 1, 'results_aistats2012_mohsst5_vbpca_uniform_20111007');
disp('GPFA')
show_results('gpfa', 1, 'results_aistats2012_mohsst5_gpfa_uniform_20111006');

disp('--- Results for PATTERN case ---')
show_results('data', 2);
disp('Kronecker GP')
show_results('gpkron', 2, 'results_aistats2012_mohsst5_gpkron_pattern_20110528');
disp('Kronecker GP (2 covfuncs)')
show_results('gpkron2', 2, 'results_aistats2012_mohsst5_gpkron-cov2_pattern_20111010');
disp('VB PCA')
show_results('pca', 2, 'results_aistats2012_mohsst5_vbpca_pattern_20111007');
disp('GPFA')
show_results('gpfa', 2, 'results_aistats2012_mohsst5_gpfa_pattern_20111006');

return

% $$$ GP1 = load([dir 'results_aistats2012_mohsst5_gpkron-cov2_uniform_20110930']);
% $$$ % $$$ GP1 = load([dir 'results_aistats2012_mohsst5_gpkron-cov2_uniform_20110907']);
% $$$ % $$$ GP2 = load([dir 'results_aistats2012_mohsst5_gpkron-cov2_pattern_20110908']);
% $$$ PCA1 = load([dir 'results_aistats2012_mohsst5_vbpca_uniform_20110601']);
% $$$ PCA2 = load([dir 'results_aistats2012_mohsst5_vbpca_pattern_20110601']);
% $$$ GPFA1 = load([dir 'results_aistats2012_mohsst5_gpfa_uniform_20110921']);
% $$$ GPFA2 = load([dir 'results_aistats2012_mohsst5_gpfa_pattern_20111001']);
% $$$
% $$$ [gp1_rmse_train, gp1_rmse_test] = aistats2012_mohsst5_rmse(GP1.F, 1);
% $$$ fprintf('GP, uniform: train=%.4f and test=%.4f\n', gp1_rmse_train, gp1_rmse_test);
% $$$ [pca1_rmse_train, pca1_rmse_test] = aistats2012_mohsst5_rmse(PCA1.F, 1);
% $$$ fprintf('PCA, uniform: train=%.4f and test=%.4f\n', pca1_rmse_train, pca1_rmse_test);
% $$$ [pca2_rmse_train, pca2_rmse_test] = aistats2012_mohsst5_rmse(PCA2.F, 2);
% $$$ fprintf('PCA, pattern: train=%.4f and test=%.4f\n', pca2_rmse_train, pca2_rmse_test);
% $$$ [gpfa1_rmse_train, gpfa1_rmse_test] = aistats2012_mohsst5_rmse(GPFA1.F, 1);
% $$$ fprintf('GPFA, pattern: train=%.4f and test=%.4f\n', gpfa2_rmse_train, gpfa2_rmse_test);
% $$$ [gpfa2_rmse_train, gpfa2_rmse_test] = aistats2012_mohsst5_rmse(GPFA2.F, 2);
% $$$ fprintf('GPFA, pattern: train=%.4f and test=%.4f\n', gpfa2_rmse_train, gpfa2_rmse_test);
return

% $$$ Y = GP1.Ytest;
% $$$ I = ~isnan(GP1.Ytrain);
% $$$ Y(I) = GP1.Ytrain(I);
% $$$
% $$$ data = mohsst5_loaddata();
% $$$ [LON,LAT] = meshgrid(data.longitude, data.latitude);
% $$$ weights = 1./cosd(LAT(:));
% $$$ weights = mohsst5_remove_land(weights);
% $$$
% $$$ %
% $$$ % Means (uniform)
% $$$ %
% $$$
% $$$ t = 800;
% $$$ % t = 1200;
% $$$
% $$$ cl = max(max(abs(Y(:,t))), ...
% $$$          max(abs(GP1.F(:,t))));
% $$$ % $$$ cl = max(max(abs(Y(:,t))), ...
% $$$ % $$$          max(max(abs(GP1.F(:,t))), ...
% $$$ % $$$              max(abs(PCA1.F(:,t)))));
% $$$
% $$$ dir = '/home/jluttine/papers/gpkron/figures_aistats2012/';
% $$$
% $$$ fig();
% $$$ mohsst5_mapplot(GP1.F(:,t));
% $$$ colormap_redblue();
% $$$ polish(gca);
% $$$ set(gca, 'clim', [-cl cl]);
% $$$ print('-depsc2', [dir 'fig_sst_gp_mean']);
% $$$
% $$$ % $$$ figure
% $$$ % $$$ mohsst5_mapplot(PCA1.F(:,t))
% $$$ % $$$ colormap_redblue();
% $$$ % $$$ polish(gca);
% $$$ % $$$ set(gca, 'clim', [-cl cl]);
% $$$ % $$$ print('-depsc2', [dir 'fig_sst_pca_mean']);
% $$$
% $$$ %
% $$$ % DATA
% $$$ %
% $$$
% $$$ fig()
% $$$ mohsst5_mapplot(Y(:,t))
% $$$ colormap_redblue();
% $$$ polish(gca);
% $$$ set(gca, 'clim', [-cl cl]);
% $$$ print('-depsc2', [dir 'fig_sst_data']);
% $$$
% $$$ %
% $$$ % STD (uniform)
% $$$ %
% $$$
% $$$ burnin = 3;
% $$$ GPstd = sqrt(GP1.stdF(:,t).^2 + weights*mean(GP1.results.theta(end,burnin:end).^2));
% $$$ % $$$ PCAstd = sqrt(PCA1.stdF(:,t).^2 + 1./PCA1.results.Tau(:,t));
% $$$
% $$$ % $$$ fig()
% $$$ % $$$ plot(GPstd(:))
% $$$ % $$$ return
% $$$
% $$$ % $$$ fig()
% $$$ % $$$ %plot(GPstd(:))
% $$$ % $$$ plot(GP1.stdF(:))
% $$$ % $$$ min(GP1.stdF(:))
% $$$ % $$$ size(GP1.results.theta)
% $$$ % $$$ min(weights(:)*mean(GP1.results.theta(end,1000:end).^2))
% $$$
% $$$
% $$$ cl = max(GPstd(:));
% $$$ % $$$ cl = max( max(GPstd), max(PCAstd) );
% $$$
% $$$ fig()
% $$$ mohsst5_mapplot(GPstd);
% $$$ colormap_scale();
% $$$ polish(gca);
% $$$ set(gca, 'clim', [0 cl]);
% $$$ print('-depsc2', [dir 'fig_sst_gp_std']);
% $$$
% $$$ % $$$ figure
% $$$ % $$$ mohsst5_mapplot(PCAstd);
% $$$ % $$$ colormap_scale();
% $$$ % $$$ polish(gca);
% $$$ % $$$ set(gca, 'clim', [0 cl]);
% $$$ % $$$ print('-depsc2', [dir 'fig_sst_pca_std']);
% $$$
% $$$ disp(datestr(data.time(t)));
% $$$
% $$$ % return
% $$$
% $$$ %return
% $$$
% $$$ %
% $$$ % Means (pattern)
% $$$ %
% $$$
% $$$ t = 1429;
% $$$
% $$$ cl = max(abs(GP2.F(:,t)));
% $$$ % $$$ cl = max( max(abs(GP2.F(:,t))), max(abs(PCA2.F(:,t))) );
% $$$
% $$$ % GP
% $$$
% $$$ fig()
% $$$ mohsst5_mapplot(GP2.F(:,t));
% $$$ colormap_redblue();
% $$$ polish(gca);
% $$$ set(gca, 'clim', [-cl cl]);
% $$$ print('-depsc2', [dir 'fig_sst_pattern_gp_mean']);
% $$$
% $$$ % $$$ % PCA
% $$$ % $$$
% $$$ % $$$ figure
% $$$ % $$$ mohsst5_mapplot(PCA2.F(:,t))
% $$$ % $$$ colormap_redblue();
% $$$ % $$$ polish(gca);
% $$$ % $$$ set(gca, 'clim', [-cl cl]);
% $$$ % $$$ print('-depsc2', [dir 'fig_sst_pattern_pca_mean']);
% $$$
% $$$ % Data
% $$$
% $$$ fig()
% $$$ mohsst5_mapplot(Y(:,t))
% $$$ colormap_redblue();
% $$$ polish(gca);
% $$$ set(gca, 'clim', [-cl cl]);
% $$$ print('-depsc2', [dir 'fig_sst_pattern_data']);
% $$$
% $$$ % $$$ figure
% $$$ % $$$ mohsst5_mapplot(PCA2.Ytrain(:,t))
% $$$ % $$$ colormap_redblue();
% $$$ % $$$ polish(gca);
% $$$ % $$$ set(gca, 'clim', [-cl cl]);
% $$$ % $$$ print('-depsc2', [dir 'fig_sst_pattern_train']);
% $$$
% $$$ disp(datestr(data.time(t)));
% $$$
% $$$ % $$$ figure
% $$$ % $$$ EGP = abs(Y-GP2.F);
% $$$ % $$$ EPCA = abs(Y-PCA2.F);
% $$$ % $$$ mohsst5_mapplot(EGP(:,t)-EPCA(:,t))
% $$$ % $$$ colormap_redblue();
% $$$ % $$$ polish(gca);
% $$$ % $$$ %set(gca, 'clim', [-cl cl]);
% $$$
% $$$ % $$$ figure
% $$$ % $$$ mohsst5_mapplot(E2(:,t))
% $$$ % $$$ colormap_scale();
% $$$ % $$$ polish(gca);
% $$$ %set(gca, 'clim', [-cl cl]);
% $$$
% $$$ figure()
% $$$ burnin = ceil(size(GP2.results.theta,2)/2);
% $$$ plot_scatterhist(GP2.results.theta(:,burnin:end)');
% $$$
% $$$ figure()
% $$$ subplot(1,2,1)
% $$$ semilogy(GP2.results.theta');
% $$$ subplot(1,2,2)
% $$$ plot(acorr(GP2.results.theta(:,burnin:end)'));
% $$$
% $$$ return
% $$$
% $$$ %
% $$$ % STD (pattern)
% $$$ %
% $$$
% $$$ GPstd = sqrt(GP2.stdF(:,t).^2 + weights*mean(GP2.results.theta(4,1000:end).^2));
% $$$ PCAstd = sqrt(PCA2.stdF(:,t).^2 + 1./PCA2.results.Tau(:,t));
% $$$
% $$$ cl = max( max(GPstd), max(PCAstd) );
% $$$ %cl = max( max(GP2.stdF(:,t)), max(PCA2.stdF(:,t)) );
% $$$
% $$$ figure
% $$$ mohsst5_mapplot(GPstd);
% $$$ colormap_scale();
% $$$ polish(gca);
% $$$ set(gca, 'clim', [0 cl]);
% $$$ print('-depsc2', [dir 'fig_sst_pattern_gp_std']);
% $$$
% $$$ figure
% $$$ mohsst5_mapplot(PCAstd);
% $$$ colormap_scale();
% $$$ polish(gca);
% $$$ set(gca, 'clim', [0 cl]);
% $$$ print('-depsc2', [dir 'fig_sst_pattern_pca_std']);
% $$$
% $$$ return

function show_results(model, testset, filename)

dir = '/home/jluttine/papers/gpkron/figures/';
printfile = [dir 'fig_sst'];
if testset == 1
  teststring = 'uniform';
  t = 800;
else
  teststring = 'pattern';
  t = 1429;
end

if nargin >= 3
  % Load results
  dir = '/home/jluttine/matlab/publications/aistats2012/';
  Q = load([dir filename], 'F');
  % RMSE
  [rmse_train, rmse_test] = aistats2012_mohsst5_rmse(Q.F, testset);
  fprintf('train=%.4f and test=%.4f\n', rmse_train, rmse_test);
  F = Q.F(:,t);
  plot_map(F);
  printfile = sprintf('%s_%s_%s', printfile, teststring, model);
  print('-depsc2', printfile);
  title([teststring ' ' model]);
else
  [Ytrain, Ytest, Itrain, Itest] = aistats2012_mohsst5_sets(testset);
  F = Ytrain;
  F(Itest) = Ytest(Itest);
  plot_map(F(:,t));
  print('-depsc2', sprintf('%s_%s_data', printfile, teststring));
  title([teststring ' data']);
  plot_map(Ytrain(:,t));
  print('-depsc2', sprintf('%s_%s_train', printfile, teststring));
  title([teststring ' train']);
end

function plot_map(F)
% Plot reconstruction
cl = 3.3;
f = fig();
mohsst5_mapplot(F);
colormap_redblue();
polish(gca);
set(gca, 'clim', [-cl cl]);
% NaNs to grey
whitebg(f, [0.7 0.7 0.7]);

function polish(hax)
hcb = colorbar('SouthOutside');
set(hcb, 'Units', 'normalized');
set_figure_size(10,8);
set_colorbar_position(hcb, hax, 0.03, 0.07);

function f = fig()
f = figure('Position', [1 300 400 300]);

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
% Creates input parameter structure array needed to run "treeqsm" function
% and "filtering" function.
% NOTE: use this code to define all the parameters but the PatchDiam and
% BallRad parameters can be conveniently defined by "define_input"
% function.
%
% Last update 11 May 2022

clear inputs

%% QSM reconstruction parameters
%%% THE THREE INPUT PARAMETERS TO BE OPTIMIZED.
% These CAN BE VARIED AND SHOULD BE OPTIMIZED
% One possibility to define these is to use "define_input" code
% (These can have multiple values given as vectors, e.g. [0.01 0.02]).
% Patch size of the first uniform-size cover:
inputs.PatchDiam1 = [0.08 0.12];
% Minimum patch size of the cover sets in the second cover:
inputs.PatchDiam2Min = [0.02 0.03];
% Maximum cover set size in the stem's base in the second cover:
inputs.PatchDiam2Max = [0.07 0.1];

%%% ADDITIONAL PATCH GENERATION PARAMETERS.
% The following parameters CAN BE VARIED BUT CAN BE USUALLY KEPT AS SHOWN
% (i.e. little bigger than PatchDiam parameters).
% One possibility to define these is to use "define_input" code
% Ball radius in the first uniform-size cover generation:
inputs.BallRad1 = inputs.PatchDiam1+0.015;
% Maximum ball radius in the second cover generation:
inputs.BallRad2 = inputs.PatchDiam2Max+0.01;

% The following parameters CAN BE USUALLY KEPT FIXED as shown.
% Minimum number of points in BallRad1-balls, generally good value is 3:
inputs.nmin1 = 3;
% Minimum number of points in BallRad2-balls, generally good value is 1:
inputs.nmin2 = 1;
% Does the point cloud contain points only from the tree (if 1, then yes):
inputs.OnlyTree = 1;
% Produce a triangulation of the stem's bottom part up to the first main
% branch (if 1, then yes):
inputs.Tria = 0;
% Compute the point-model distances (if 1, then yes):
inputs.Dist = 1;

%%% RADIUS CORRECTION OPTIONS FOR MODIFYING TOO LARGE AND TOO SMALL CYLINDERS.
% These parameters CAN BE USUALLY KEPT FIXED as shown.
% Traditional TreeQSM choices:
% Minimum cylinder radius, used particularly in the taper corrections:
inputs.MinCylRad = 0.0025;
% Radius correction based on radius of the parent. If 1, radii in a branch
% are always smaller than the radius of the parent in the parent branch:
inputs.ParentCor = 1;
% Parabola taper correction of radii inside branches. If 1, use the
% correction:
inputs.TaperCor = 1;

% Growth volume correction approach introduced by Jan Hackenberg,
% allometry: Radius = a*GrowthVol^b+c
inputs.GrowthVolCor = 0; % If 1, use growth volume (GV) correction
% fac-parameter of the GV-approach, defines upper and lower bound. When
% using GV-approach, consider setting TaperCorr = 0, ParentCorr = 0,
% MinCylinderRadius = 0.
inputs.GrowthVolFac = 1.5; % Defines the allowed radius:
%  1/fac*predicted_radius <= radius <= fac*predicted_radius
%  However, the radii of the branch tip cylinders are not increased.

%% Filtering parameters
% NOTE: These are all optional, but needed to run the "filtering" function.
% Statistical k-nearest neighbor distance outlier filtering, applied if
% filter.k > 0. The value filter.k is the number of nearest neighbors.
inputs.filter.k = 10;
% Statistical point density outlier filtering, applied if filter.radius > 0.
% The value filter.radius is the radius of the ball neighborhood. This is
% usually meant as alternative to the above knn-filtering.
inputs.filter.radius = 0.00;
% The value filter.nsigma is the multiplier of the standard deviation of
% the kth-nearest neighbor distance/point density and points whose
% kth-nearest neighbor distance/point density is larger/lower than the
% average +/- filter.nsigma * std are removed:
inputs.filter.nsigma = 1.5;
% Small component filtering is applied if filter.ncomp > 0. This filter is
% based on cover whose patches are defined by filter.PatchDiam1 and
% filter.BallRad1. The points which are included in components that have
% less than filter.ncomp patches are removed:
inputs.filter.PatchDiam1 = 0.05;
inputs.filter.BallRad1 = 0.075;
inputs.filter.ncomp = 2;
% Cubical downsampling is applied if filter.EdgeLength > 0.
% The value filter.EdgeLength is the length of the cube edges:
inputs.filter.EdgeLength = 0.004;
% Plot the filtering results automatically after the filtering if
% filter.plot > 0
inputs.filter.plot = 1;

%% Other inputs
% These parameters don't affect the QSM-reconstruction but define what is
% saved, plotted, and displayed and how the models are named/indexed
% Name string for saving output files and naming models:
inputs.name = 'tree';
% Tree index. If modelling multiple trees, then they can be indexed uniquely:
inputs.tree = 1;
% Model index, can separate models if multiple models with the same inputs:
inputs.model = 1;
% Save the output struct QSM as a matlab-file into \result folder.
% If name = 'pine', tree = 2, model = 5, the name of the saved file is
% 'QSM_pine_t2_m5.mat':
inputs.savemat = 1;
% Save the models in .txt-files (check "save_model_text.m"):
inputs.savetxt = 1;
% What are plotted during reconstruction process:
% 2 = plots the QSM, the segmentated point cloud and distributions,
% 1 = plots the QSM and the segmentated point cloud
% 0 = plots nothing
inputs.plot = 2;
% What are displayed during the reconstruction: 2 = display all;
% 1 = display name, parameters and distances; 0 = display only the name:
inputs.disp = 2;

author: InverseTampere

repo: TreeQSM


text:
function RHOTPH2O = h2o_rhotp(P,X)
% H2O_RHOTP  Two-phase density of H2O in kg/m^3
% H2O_RHOTP(P,X) Returns the two-phase density
% of H2O at a given pressure and quality.
%  Called function: h2o_rhof(P), h2o_rhog(P)
%  Required Inputs are: P - pressure (kPa), X - quality (fraction)
% ---------------------------------------------------------------
%  The MATLAB function was created by Tibor Balint, December 1998
% TBoreal Research Corporation, Toronto, Ont. Canada
% (tibor@netcom.ca) and also, University of Warwick, UK
% ---------------------------------------------------------------

format long g;  % set the format of the calculations

RF=h2o_rhof(P);
RG=h2o_rhog(P);

INVRHOTP=(X/RG)+((1-X)/RF);
RHOTPH2O=1/INVRHOTP;
return


text:
% MCMCSUMM - Summary Statistics
% Copyright (c) 1998, Harvard University. Full copyright in the file Copyright
%
%   [S] = mcmcsumm(A)
%
% A = r x c x s array of s samples of an r x c matrix of parameters
%
% S = structure containing returned values (mean, median, etc.)
%   select components with S.mean, S.median, etc.
%
% Note: all summary statistics are marginal, there are no multivariate
%   summaries at this time.
%
% These routines use the last dimension of an array as the
% sample index.  So an array with dimension (nr,nc,ns)
% will be a collection of ns samples of an nr by nc matrix of
% parameters.  An array with dimension (nr,nc) will
% be nc samples of an nr-vector of parameters.
% When the summary statistics are calculated, the last dimension
% is dropped.
%
% See also: MCMCTRACE, MCMCLT
%

function [S] = mcmcsumm(A)

if isnan(A),
 S.mean = NaN;
 S.min = NaN;
 S.max = NaN;
 S.std = NaN;
 S.sorted = NaN;
 S.median = NaN;
 S.meanvec = NaN;
 S.cov = NaN;
 S.acf= NaN;
 S.acf10max = NaN;
 S.acf10med = NaN;
 S.gr2 = NaN ;
 S.gr2max = NaN ;
else

dd = size(A) ;
ll = length(dd) ;
if (ll==2),
  aa = reshape(A, [dd(1),1,dd(2)]) ;
else
  aa = A ;
end

[nr,nc,ns] = size(aa) ;

maxlag = min(100,ns-1) ;

Z = zeros(nr,nc,ns) ;
S = struct('mean',Z) ;

S.mean = mean(aa,3) ;
S.min = min(aa,[],3) ;
S.max = max(aa,[],3) ;
S.std = std(aa,0,3) ;
S.sorted = NaN*zeros(nr,nc,ns) ;
S.median = NaN*zeros(nr,nc) ;

tmpvec = reshape(S.mean, nr*nc, 1) ;
sel = ~isnan(tmpvec) ;
S.meanvec = tmpvec(sel,:) ;

aavec = reshape(aa, nr*nc, ns) ;
aavec = aavec(sel,:) ;

if nr>0 & nc>0 & ns>0,
% then there's something to work with

S.cov = cov(aavec') ;

for ir = 1:nr,
for ic = 1:nc,
  xx = reshape(aa(ir,ic,:),1,ns) ;
  S.sorted(ir,ic,:) = sort(xx) ;
  S.median(ir,ic) = median(xx) ;
  xx0 = xx - mean(xx) ;
  if S.max(ir,ic)-S.min(ir,ic) < .0000000001,
    xc = NaN * zeros(1,2*maxlag+1) ;
  else
    xc = xcorr(xx0,xx0,maxlag,'coeff');
  end
  S.acf(ir,ic,:) = [xc(maxlag+(1:(maxlag+1)))] ;
  S.acf1 = S.acf(:,:,2) ;
end
end

if ns>10,
  S.acf10 = S.acf(:,:,11) ;
  tmpacf = reshape(S.acf10,1,nr*nc) ;
  if any(~isnan(tmpacf)),
    tmpacf = tmpacf(~isnan(tmpacf)) ;
    S.acf10max = max(max( tmpacf )) ;
    S.acf10med = median(median( tmpacf )) ;
  else
    S.acf10max = NaN ;
    S.acf10med = NaN ;
  end
else
  S.acf10 = NaN*S.acf(:,:,1) ;
  S.acf10max = NaN ;
  S.acf10med = NaN ;
end

else
  % one of nr nc ns == 0
  S.sorted(:,:,:) = A ;
  S.median = A(:,:,1) ;
  S.acf10max = NaN ;
  S.acf10med = NaN ;
  S.acf = NaN * zeros(nr,nc,maxlag+1) ;
  S.acf1 = S.acf(:,:,2) ;
  S.acf10 = S.acf(:,:,11) ;
  S.cov = A(:,:,1) ;
end

S.gr2 = mcmcgr(A,2) ;
S.gr2max = max(max(S.gr2)) ;

end
% end of NaN branch


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Max manipulability index ALONG A LINE.
% Use stomp like to optimize along a surface/line
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function experiment1B_K10_N30
close all;
global parameters

%STOMP PARAMETERS
%number of particles
K = 10;

%repeat experiment number of times
parameters.n_repeat = 3;
parameters.experiment_name = 'experiment1B_K10_N30_time.mat';
parameters.experiment_number = 1;
parameters.animate = 0;

%repeat the experiment E times
random_manips=[];
Gout = [];
tic;
for i=1:parameters.n_repeat
    close all
    [pk, final_manip] = path_planning_SCO_4DOF(K);
    Gout{i}=pk;
    random_manips = [random_manips; final_manip];
    save(parameters.experiment_name)
end
elapsed = toc;
fprintf('TIC TOC: %g\n', elapsed/parameters.n_repeat);

'ended'
parameters.experiment_name

repo: ARTE


text:
% fig2: anat averages from paired regressors: sensory/motor correlations
% plots a pair of regressors in a graded double colormap, saves as tiff
% stacks, and also saves anat average plot as tiff file.
% (compare to fig2A_correlation_1reg_with_hist for more details)

clear all; close all; clc

outputDir = GetOutputDataDir;

%% Init load
hfig = figure;
InitializeAppData(hfig);
ResetDisplayParams(hfig);

%% set params
thres_reg_const = 0.5;
% M_reg_thres = {0.5,0.5,0.5};

for setflag = 7
    % reset
    ResetDisplayParams(hfig);

    switch setflag
        case 1 % Phototaxis
%             stimrange = 1;
M_stimrange = GetStimRange('P');
            M_regtypeflag = [1,1,2]; % 0 for VAR, 1 for stim and 2 for motor
            M_reg_name = {'PT','HalfField','PT_SwimLR'};
            M_reg_range = {[3,2],[5,6],[1,3]};
            M_fishrange = {[1:18],[1:7],[1:3,5:18]};
            n_reg = length(M_reg_name);

        case 2 % OMR
%             stimrange = 2;
M_stimrange = GetStimRange('O');
            M_regtypeflag = [1,1,2]; % 1 for stim and 0 for motor
            M_reg_name = {'OMR_FwBw','OMR_LR','OMR_SwimLR'};
            M_reg_range =  {[7,6],[9,8],[1,3]};
            M_fishrange = {[8:18],[8:18],[8:18]};
            n_reg = length(M_reg_name);

        case 3 % Looming
%             stimrange = 5;
M_stimrange = GetStimRange('L');
            M_regtypeflag = [1,2]; % 1 for stim and 0 for motor
            M_reg_name = {'Loom_LR','Loom_SwimLR'};
            M_reg_range =  {[11,12],[1,3]};
            M_fishrange = {[9:15,17:18],[9:15,17:18]};
            n_reg = length(M_reg_name);

        case 4 % Dark Flash (black-white)
%             stimrange = 3;
M_stimrange = GetStimRange('D');
            M_regtypeflag = [1,2]; % 1 for stim and 0 for motor
            M_reg_name = {'DF_BW','DF_SwimLR'};
            M_reg_range =  {[1,4],[1,3]};
            M_fishrange = {[12:15,17:18],[12:15,17:18]};
            n_reg = length(M_reg_name);

%             % from VAR:

        case 5 % ABN
            M_regtypeflag = [0];
            M_stimrange = GetStimRange();%('2');
            M_reg_name = {'ABN_reg0.5_defstimrange_test'};
            M_clus_range = {[12,1]};
            M_fishrange = {[1:12,14:18]};
            n_reg = length(M_reg_name);
        case 6 % Fw
            M_regtypeflag = [0];
            M_stimrange = GetStimRange();%('2');
            M_reg_name = {'Fw_reg0.5_defS'};
            M_clus_range = {[11,4]};
            M_fishrange = {[1:18]};
            n_reg = length(M_reg_name);

        case 7 % left HBO
            M_regtypeflag = [0,0];
            M_stimrange = GetStimRange();%('2');
            M_reg_name = {'HBO-L2_reg0.5_defS','HBO-R2_reg0.5_defS'};
            M_clus_range = {[10,2],[10,3]};
            M_fishrange = {[1:3,5:18],[1:3,5:18]};
            n_reg = length(M_reg_name);
    end
    M_fishrange_im = M_fishrange;

    %% Load fish
    range = 1:18;
    IM_full = cell(n_reg,max(range));
    %%

%     [M_fishrange_im,fishrange_load] = CheckIfLoadFish(M_fishrange,M_ClusterIDs,M_stimrange);
% [cIX_load,gIX_load,M,stim,behavior,M_0] = LoadSingleFishDefault(i_fish,hfig,M_ClusterIDs{i_set},M_stimrange{i_set});
    for i_fish = range
        if ismember(i_fish,cell2mat(M_fishrange))
            ClusterIDs = [1,1];%GetClusterIDs('all');
            stimrange = M_stimrange{i_fish};
            % load fish data
            if isempty(stimrange)
                [cIX_load,gIX_load,M,stim,behavior,M_0] = LoadSingleFishDefault(i_fish,hfig,ClusterIDs);
            else
                [cIX_load,gIX_load,M,stim,behavior,M_0] = LoadSingleFishDefault(i_fish,hfig,ClusterIDs,stimrange);
            end
        end

        %% Load stim/motor
        for i_set = 1:n_reg

            if ~ismember(i_fish,M_fishrange{i_set})
                continue;
            end

            reg_thres = thres_reg_const; %M_reg_thres{i_set};

            % get stim/motor regressors
            if M_regtypeflag(i_set)==0
                ClusterIDs = M_clus_range{i_set};%[12,1];% GetClusterIDs('all');
                [cIX,gIX] = LoadCluster_Direct(i_fish,ClusterIDs(1),ClusterIDs(2));
                M = UpdateIndices_Manual(hfig,cIX,gIX);
                Reg = FindClustermeans(gIX,M);

            elseif M_regtypeflag(i_set)==1
                fishset = getappdata(hfig,'fishset');
                [~,names,regressors] = GetStimRegressor(stim,fishset,i_fish);

                reg_range = M_reg_range{i_set}; % left/right pair
                Reg = regressors(reg_range,:);

            elseif M_regtypeflag(i_set)==2
                isMotorseed = 0;
                setappdata(hfig,'isMotorseed',isMotorseed);
                [~,~,behavior] = UpdateTimeIndex(hfig);
                [~,names,regressors] = GetMotorRegressor(behavior,i_fish);

                reg_range = M_reg_range{i_set}; % left/right pair
                Reg = regressors(reg_range,:);
            end

            % code adapted from 'best regressor regression' code 'AllRegsRegression'

            Corr = corr(Reg',M_0');
            [corr_max,IX_regtype] = max(Corr,[],1);
            cIX = find(corr_max>reg_thres)';
            gIX_offset = IX_regtype(cIX)';
            clrIX = MapXto1Dcolormap(corr_max(cIX),[reg_thres,1],64);

            if isempty(cIX)
                M_fishrange_im{i_set} = setdiff(M_fishrange_im{i_set} ,i_fish);
                continue;
            end

            gIX = clrIX+(gIX_offset-1)*64;
            numK = length(unique(gIX));

            %% make double colormap - ??
            clr1 = [1,0,0];
            clr1_ = [0.7,0.5,0.5];
            clr2 = [0,1,1];
            clr2_ = [0.5,0.7,0.7];
            numC = 64;
            clrmap1 = Make1DColormap([clr1_;clr1],numC);
            clrmap2 = Make1DColormap([clr2_;clr2],numC);
            clrmap = [clrmap1;clrmap2];

            %% make figure
            I = LoadCurrentFishForAnatPlot(hfig,cIX,gIX,clrmap);
            [h,im_full] = DrawCellsOnAnat(I);

            %         % add 2 colorbars
            % %         AddColorbarToAnat(clrmap,cmin,cmax)
            % %         colormap(clrmap1);
            % %         caxis([reg_thres,1])
            % %         colorbar('Location','manual','Position',[0.8,0.7,0.05,0.15],'Units','normalized')
            %         ax = axes('Position',[0.75,0.8,0.05,0.15],'Units','normalized');
            %         DrawCustomColorbar(clrmap1,[reg_thres,1],2,ax);
            %
            %         ax = axes('Position',[0.9,0.8,0.05,0.15],'Units','normalized');
            %         DrawCustomColorbar(clrmap2,[reg_thres,1],2,ax);
            %
            %% save figure
            close(h);
            IM_full{i_set,i_fish} = im_full;

        end
    end

    %% save as tiff stack
    for i_set = 1:n_reg
        range_im = M_fishrange_im{i_set};
        tiffdir = fullfile(outputDir,[M_reg_name{i_set},'_allfish.tiff']);
        IM = IM_full(i_set,range_im);

        SaveImToTiffStack(IM,tiffdir);
    end

    %% [for later] plot from tiff stack
    isPlotfromtiffstack = 0;

    if isPlotfromtiffstack
        IM_full = cell(n_reg,18);
        for i_set = 1:n_reg
            %% get tiff-stack path
            tiffdir = fullfile(outputDir,[M_reg_name{i_set},'_allfish.tiff']);

            %% load

            for i_fish = 1:18
                im = double(imread(tiffdir,i_fish))./255;
                IM_full{i_set,i_fish} = im(317:1236,1:621,:);
            end
        end
    end

    %% make average plots
    % M_k_scale = {1,1.5,1};
    % M_k_contrast = {1.2,1.5,1.2};

    for i_set = 1:n_reg;
        range_im = M_fishrange_im{i_set};%[1:3,5:7];%[1:3,5:18];
        cellarray = IM_full(i_set,range_im);

        % adjust params for visualization
        k_scale = 0.4;%1/1.5;%M_k_scale{i_set};
        k_contrast = 1;%M_k_contrast{i_set};

        [h_anat,im_avr] = AverageAnatPlot(cellarray,k_contrast,k_scale);

        tiffdir = fullfile(outputDir,[M_reg_name{i_set},'_avr.tiff']);
        imwrite(im_avr, tiffdir, 'compression','none','writemode','overwrite');
    end

end

author: xiuyechen

repo: FishExplorer


text:
function imdb = getImdbDCFNet(varargin)
opts = [];
opts.dataDir = fullfile('..','data');

opts.visualization = false;
opts.output_size = 125;
opts.padding = 1.5;
opts = vl_argparse(opts, varargin);

bbox_mode = 'axis_aligned';
if strcmp(bbox_mode,'axis_aligned')
    get_bbox = @get_axis_aligned_BB;
elseif strcmp(bbox_mode,'minmax')
    get_bbox = @get_minmax_BB;
end
% -------------------------------------------------------------------------
%   full dataset:
%           vid2015: 464873
% -------------------------------------------------------------------------
set_name = {'vid2015'};
num_all_frame = 464873;
%% Be careful!!! It takes a HUGE RAM for fast speed!
imdb.images.set = int8(ones(1, num_all_frame));
imdb.images.set(randperm(num_all_frame, 1000)) = int8(2);
imdb.images.images = zeros(opts.output_size, opts.output_size, 3, num_all_frame, 'uint8');
imdb.images.up_index = zeros(1, num_all_frame, 'double'); % The farthest frame can it touch
now_index = 0;

% -------------------------------------------------------------------------
%                                                                   VID2015
% -------------------------------------------------------------------------
if any(strcmpi(set_name, 'vid2015'))
    disp('VID2015 Data:');
    if exist('vid_2015_seg.mat', 'file')
        load('vid_2015_seg.mat');%% use dataPreprocessing;
    else
        error('You should generate <vid_2015_seg.mat> according ''dataPreprocessing'' at first.')
    end
    videos = seg;
    n_videos = numel(videos);
    for  v = 1:n_videos
        video = videos{v};fprintf('%3d / %3d \n', v, n_videos);

        img_files = video.path;

        ground_truth = cell2mat((video.rect)');
        ground_truth_4xy = [ground_truth(:,1),ground_truth(:,2),...
            ground_truth(:,1),ground_truth(:,2)+ground_truth(:,4),...
            ground_truth(:,1)+ground_truth(:,3),ground_truth(:,2)+ground_truth(:,4),...
            ground_truth(:,1)+ground_truth(:,3),ground_truth(:,2)];
        invaildindex = any(isnan(ground_truth), 2);
        img_files(invaildindex) = [];
        ground_truth_4xy(invaildindex,:) = [];

        im_frist = vl_imreadjpeg(img_files(1));
        [H, W, ~] = size(im_frist{1});
        im_bank = vl_imreadjpeg(img_files, 'Pack', 'Resize', [H, W], 'numThreads', 32);
        bbox_gt = get_bbox(ground_truth_4xy);
        n_frames = size(bbox_gt,1);
        imdb.images.images(:,:,:,now_index+(1:n_frames)) = uint8(...
            imcrop_pad(im_bank{1}, bbox_gt, opts.padding, opts.output_size([1,1])));
        imdb.images.up_index(now_index+(1:n_frames)) = (n_frames:-1:1)-1;
        imdb.images.set(now_index+n_frames) = 4; %should not be selected as x.
        now_index = now_index + n_frames;
    end %%end v
end %%end VID2017

% dataMean = single(mean(imdb.images.images,4));
dataMean(1,1,1:3) = single([123,117,104]);
imdb.images.data_mean(1, 1, 1:3) = dataMean;
imdb.meta.sets = {'train', 'val'} ;
end %%end function

author: foolwood

repo: DCFNet

save_path: github-repos/MATLAB/foolwood-DCFNet


text:
function c=pfilt(f,g,varargin)
%PFILT  Apply filter with periodic boundary conditions
%   Usage:  h=pfilt(f,g);
%           h=pfilt(f,g,a,dim);
%
%   `pfilt(f,g)` applies the filter *g* to the input *f*. If *f* is a
%   matrix, the filter is applied along each column.
%
%   `pfilt(f,g,a)` does the same, but downsamples the output keeping only
%   every a'th sample (starting with the first one).
%
%   `pfilt(f,g,a,dim)` filters along dimension dim. The default value of
%   [] means to filter along the first non-singleton dimension.
%
%   The filter *g* can be a vector, in which case the vector is treated
%   as a zero-delay FIR filter.
%
%   The filter *g* can be a cell array. The following options are
%   possible:
%
%     * If the first element of the cell array is the name of one of the
%       windows from |firwin|, the whole cell array is passed onto
%       |firfilter|.
%
%     * If the first element of the cell array is `'bl'`, the rest of the
%       cell array is passed onto |blfilter|.
%
%     * If the first element of the cell array is `'pgauss'`, `'psech'`,
%       the rest of the parameters is passed onto the respective
%       function. Note that you do not need to specify the length *L*.
%
%   The coefficients obtained from filtering a signal *f* by a filter *g* are
%   defined by
%
%   ..          L-1
%      c(n+1) = sum f(l+1) * g(an-l+1)
%               l=0
%
%   .. math:: c\left(n+1\right)=\sum_{l=0}^{L-1}f\left(l+1\right)g\left(an-l+1\right)
%
%   where $an-l$ is computed modulo $L$.
%
%   See also: pconv

% Assert correct input.
if nargin<2
  error('%s: Too few input parameters.',upper(mfilename));
end;

definput.import={'pfilt'};
definput.keyvals.a=1;
definput.keyvals.dim=[];
[flags,kv,a,dim]=ltfatarghelper({'a','dim'},definput,varargin);

[f,L,Ls,W,dim,permutedsize,order]=assert_sigreshape_pre(f,[],dim,upper(mfilename));

[g,info] = comp_fourierwindow(g,L,upper(mfilename));

outIsReal = isreal(f) &&...
              (isfield(g,'h') && isreal(g.h) && ~(isfield(g,'fc') && g.fc~=0) || isfield(g,'H') && g.realonly);

asan = comp_filterbank_a(a,1);
g = comp_filterbank_pre({g},a,L,kv.crossover);

c = comp_filterbank(f,g,a);
c = c{1};

permutedsize(1)=size(c,1);

c=assert_sigreshape_post(c,dim,permutedsize,order);

if outIsReal
   c = real(c);
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function  [cct, ta1, tetael] = srcint1(la, nn1, ct4, pp1, pp2, ctc0, ctad, nprimtilt, lmax, K, stfname)
% srcint1.m integrates the source including polarisation over conformal time using the approach of Seljak :
% integration over the photon past light cone, see M Zaldarriaga et al., ApJ nr 494, 491 (1998) for K > 0
% the function needs an l-range (la), a k-range (nn1) a ctime-range (ct4), a series of splines of the source
% in conformal time for the temperature (pp1) and for the polarisation (pp2),the present conformal time (ctc0),
% decoupling time (ctad), the index of the primary power spectrum (nprimtilt),
% the maximum la of the anistropy spectrum (lmax) and the curvature energy content from the Friedman equation (K),
% the startfile name of the ultra-spherical function.
%
% the result is the square of the temperature spectrum (1) resp the E-polarisation spectrum (2)
% resp the cross-correlation of temperature and polarisation (3).
%
% D Vangheluwe 13 june 2005
% remark 1a: attention :nn1 is an integer array :wavenumber array k1 is made from it by taking k1 = nn1 * sqrt(K)
% remark 1: we use ultra-spherical bessel values which have two parameters and must be found by integration:
%   for the calculation of the ultra-spherical bessel functions see: cmb/usphint.m
% remark 2: for the ode equation of the u-function, see equation (36) of Zaldarriaga & Seljak.
% remark 3: take attention the value of nn1 should lie within the range of kb1 values : there is no check!!
%   see spline of start  values x10 and pbd0.
% remark 4: how do we get the ultra-spherical bessel function values? see development and test routine usphpar.m
%  and usphpar1.m
% remark 5: this function is a fast version of srcpint.m (the source need not to be in the loop
% if there is enough memory, probably 80Mbytes will be used)
% remark6 on 13 june: this routine has been obtained by modifying srcpintf.m according to usphpar1.m

if (K < 0),
   error('the space curvature constant K should be > 0 for this routine');
   return;
end
kk0 = sign(K);
lla = size(la, 2);
k1 = nn1 * sqrt(K);
lk1 = size(k1, 2);
lct4 = size(ct4, 2);
k1step = k1(2) - k1(1);
ct4step = ct4(2) - ct4(1);
k1max = k1(lk1);
k1min = k1(1);
ct4sym = 0.5 * pi/sqrt(K);

% if necessary adapt the source to the period of the ultra-spherical function for K >0
ct0 = ct4(1);
ict1 = 1 : lct4;
if ct4sym < ct4(end)
   ctc4 = ctc0 - ct4;
   ict1 = find(ctc4 < ct4sym);
   lct1 = size(ict1,2);
%   ct0 = ctc0 - ct4sym;
   ct0 = ct4(ict1(1));
% map the region ctc4 >= ct4sym onto ict1 by mirroring it wrt the point: ctc4 = ct4sym
   ict2 = find(ctc4 < ct4sym  &  ctc4 > 2*ct4sym - ctc0 + ct4step);
% values of ct4 to be used to calculate (by interpolation) the mirrored source values:
   ct44 = 2*(ctc0 - ct4sym) - ct4(ict2);
% define a (small) region (index icts) where the source changes rapidely and should be splined for reconstruction:
% ct4(icts(1..end)) is projected onto ct44s(end..1) :notice that the index is inverted, the phaseshift between
% the two amounts to : 2*rem(ctc0 - ct4sym - ct4(1), ct4step).
% The same phaseshift and index inversion exists also between ct4(ict1(1..end) and ct44(end..1)
% icts and ict3 have the same size and take attention that they can be empty for ctc0 - ctad +100 < ct4sym!!
   ctcsmin = 2*ct4sym - ctc0 + (ctad - 100);
   ctcsmax = 2*ct4sym - ctc0 + (ctad + 100);
   icts = find((ct4 < ctad + 100)  &  (ct4 >= ctad - 100)  &  (ctc0 - ct4 >= ct4sym));
   ict3 = find(ctc4(ict2) < ctcsmax  &  ctc4(ict2) >= ctcsmin);
   if ~isempty(ict3),    ct44s = ct44(ict3); end
end

% load the table and find the start values for the integration of the ode for u-functions
% the startfile can be obtained by running usphst.m
%  stv1 = load('usphst.dat');
stv1 = load(stfname);
lmax0 = stv1.ust.la(end);

nrsteps = 10000;   % default 10000
%lmax = 1500;
xmax = k1(lk1) * ctc0;   % =default 3000
xstep = 2*lmax0/nrsteps;
kx = 1e-12 : xstep : xmax;

% allocate the data in order to speed up the routine
kctc0 = zeros(1, lct4);
kctc = zeros(1, lct4);
jl = zeros(lct4, 1);
htable = zeros(lct4, 1);
cmbtable = zeros(1, lk1);

% calculate the curvature length times wavenumber : a parameter for the ultra_spherical bessel function
kb = sqrt(abs(K)) ./k1;

%#########
%kb(1)
%kb(end)
if all(kb == 0), kbzero = 1; x10 = 1e-12 * ones(1, lk1);
else  % not all kb are zero

   kbzero = 0;
   ikb = find(kb > 10);
   if ~isempty(ikb),
      kb(ikb) = 10 * ones(1, size(ikb, 2));
      message('boundary kb > 10 reached')
   end

% split the structure stv1 from the startfile : x1: start values where pbi_beta = 1e-6,
% pbd : dphi_beta/dx at the start values, all data exact within 1e-6.
% the la-values resp kb-vector should be the same as in the program usphst.m
   la1 = stv1.ust.la;
   lla1 = size(la1, 2);
   kb2 = stv1.ust.kb;
   x1v1 = stv1.ust.x1;
   pbdv1 = stv1.ust.pbd;

end %all(kb == 0)

%############################# start
for il = 1:lla

    l = la(il)
    laa = sqrt(l * (l + 1));

% make a range of wavenumber and kb starting at 1/(l+1)
    clear('k2')
%    ik2 = (l+1) : lk1;
    ik2 = find(nn1 >= l+1);
    k2 = k1(ik2);
    lk2 = size(k2, 2);
    nn2 = nn1(ik2);
    cmbtable = zeros(1,lk2);

% the range of indices in src for k2 will be ik2 : take care, do not neglect!!:
    k2min = k2(1);
    k2max = k2(lk2);
    tetatl = zeros(1, lk2);
    tetael = zeros(1, lk2);
    kb = sqrt(abs(K)) ./k2;
    odd_la = xor(rem(nn2,2), rem(l,2));

% the kb-range is made more progressive towards 1/la
    nrkbsteps = 100;
    kbmax = 1/laa;
    ekbmax = -7;
    ekbmin = log(kbmax - 1e-5)/log(10);
    ekbstep = (ekbmax - ekbmin)/nrkbsteps;
    ekb1 = ekbmin : ekbstep : ekbmax;
    kb1 = kbmax - 10 .^ekb1;
    lkb1 = size(kb1, 2);
    if l == la1(1)
        if any(kb1 ~= kb2)  error('kb vector is unequal to the start vector'); return; end
    end

% find the start values (hs) for the ode integration step as a function of la and kb (is a surface):
% the long sought magic formula : la1 is the la-vector from usphst1.m!!
    hs = 0.5 * x1v1 .* repmat((la1 .^-0.9)', 1, lkb1);
    if ~kbzero
       il1 = find(la1 == l);
       if isempty(il1)
           error('la not in the range of la1')
           break;
       else
           x10 = spline(kb1, x1v1(il1,:), kb);
           pbd0 = spline(kb1, pbdv1(il1,:), kb);
           hstart = spline(kb1, hs(il1,:), kb);
       end

% calculate the start values for the u-function : u0 = [u, du/dx] with u = r(x) * phi_beta(x) and
%  du/dx = phi_beta(x) * dr(x)/dx + r(x) * dphi_beta(x)/dx :
       r10 = sin(kb .* x10) ./kb;
       rd10 = cos(kb .* x10);
       u0 = [r10 *1e-6; r10 .* pbd0 + 1e-6 * rd10];

% define the step and set 'odeint' to a constant number of steps for all cases (kb)
       toi = 1e-4;
       hmax = 0.3;   %we take hmax = 0.3 as the default value
       maxstp = 150;   % default maxstp = 150
       x1e = ones(1, lk2) * xmax;
% solve the u-function instead of the phi_beta function : xs2 does not have a constant step
       [xs2, u2, nok, nbad, nfev] = odeintp(u0, x10, x1e, toi, hstart, 0, hmax, maxstp, @uspheq1, l, kb, kk0);

% prepare a spline of the solution for a constant step in xs (the values are found with ppval1):
       y2der = splinep(xs2, u2);
% calculate the number of constant steps, xstep in the solution xs2: ilast is the last one
%       ilast = ceil((xs2(end,:) - x10)/xstep);
% set the number of overlapping steps for the matching to the asymptotic solution (default= 40):
       noverflow = 85;
       ncsteps = 0;

    end  % kbzero

% calculate also the spherical bessel function as we may need it for kb < 1e-5
    table_bv = sphbes(l, kx);
% set some constants needed for the polarisation formula
    gl = sqrt((l + 2) * (l + 1) * l * (l - 1));

% interpolate and integrate (sum with the Simpson rule) over conformal time :
    for i = 1:lk2
%    for i = 26:26

       clear('htable', 'src4', 'src5')
       src4 = ppval(pp1, k2(i));
       src5 = ppval(pp2, k2(i));

       if kb(i) > 1e-5 * (1500/l)

% define a vector of argument values for the ultra-spherical bessel function (x1): reverse ct4
%          clear('xs', 'x1', 'x2', 'yspl2','iover')
          x1 = x10(i) : xstep : xmax;
% define some parameters for the calculation ; x1=x1sym up to the symmetry point, x1* kb = pi/2 needed for K>0 :
          x1sym = min(0.5*pi/kb(i), xmax);
% in the next 10 lines make a table over x1 (ul) of ultra-spherical bessel values:

          if  xs2(end,i) >=  x1sym

% asymptotic expansion is not needed here as we found already the complete solution: include one point xs2 > x1sym+xstep
% adding xs2step + xstep to x1sym in the find logic statement of xs2 makes sure that xs(end) > x1(ilast) > x1sym
             xs2step = max(diff(xs2(:,i)));
             ix2 = find(xs2(:,i) <= (x1sym + xs2step + xstep));
% include one point xs2 > x1sym in the interpolation of the ode solution:
%             ix2 = [ix2', (ix2(end) + 1)];
             [xs, yspl2] = ppval1(xs2(ix2,i), u2(ix2,i), y2der(ix2,i), xstep);
%             ilast = floor((x1sym - x10(i))/xstep) + 1;
             ilast = floor((x1sym - x10(i))/xstep) + 2;

          else  % all steps in odeint are used (full interpolation) : an asymptotic solution is needed

             [xs, yspl2] = ppval1(xs2(:,i), u2(:,i), y2der(:,i), xstep);
             ilast = floor((xs(end) - x10(i))/xstep) + 1 - ncsteps;
% define the overflow index (iover): default 85 (or 40) steps should at least include one zero and one maximum of u2
             iover = (ilast - noverflow) : ilast;

% find the asymptotic values and the phase correction in the overflow region:
             dphi = phdif(x1(iover), usphas(l, kb(i), x1(iover), 0, 0, kk0), yspl2(iover));

          end

% calculate the ultra-spherical bessel values
          kctc0 = ctc0 * k2(i) * ones(1, lct4);
          kctc = k2(i) * ct4;
          xctc = kctc0 - kctc;
% xctc > x1(1) is necessary as x1(1)=x10 is not the origin of ct, but the point where odeint started
          if   xs2(end,i) >=  x1sym
% debug action 12-6-2005 : interpolation over xs limits our range to xs(end) at most if it happens xs(end) < x1sym
%             ixode = find(xctc > x1(1)  &  xctc <= x1sym);
             ixode = find(xctc > x1(1)  &  xctc <= min(x1sym, xs(end)));
          else
             ixode = find(xctc > x1(1)  &  xctc <= x1(ilast));
          end
          ul = zeros(1, lct4);
          if  ~isempty(ixode)
             r1 = sin(kb(i) * xs)/kb(i);

% interpolate the ultra-spherical bessel function when we are in the k-range of the ode solution
             ul(ixode) = interpl(xs, yspl2 ./ r1, xctc(ixode));
% calculate the necessary values of the asymptotic expansion of the ultra-spherical bessel function
             if  (xs2(end,i) <  x1sym)   % asymptotic expansion needed for K>0
                ixas = find(xctc > x1(ilast)  &  xctc <= x1sym);
                ul(ixas) = usphas1(l, kb(i), xctc(ixas), dphi, 1, kk0);
             end

% if necessary adapt the source to the period of the ultra-spherical function for K >0 by mirroring:
% we should have: find(k2(i)* (ctc0 - ct4(ict1)) < x1sym) == find(xctc < x1sym) == ict1
             if  ct4sym < ct4(end)
% reconstruct the source by linear interpolation : not so accurate but sufficient where changes are slow
                src42 = interpl(ct4, src4', ct44)';
                src52 = interpl(ct4, src5', ct44)';
% for efficieny reason reconstruct the source by splines only in the region (ict3) where it changes fast
                if ~isempty(ict3)
                  src42(ict3) = spline(ct4(icts), src4(icts), ct44s)';
                  src52(ict3) = spline(ct4(icts), src5(icts), ct44s)';
                end
                if odd_la(i) == 1
                   src4(ict2) = src4(ict2) + src42;
                   src5(ict2) = src5(ict2) + src52;
                else
                   src4(ict2) = src4(ict2) - src42;
                   src5(ict2) = src5(ict2) - src52;
                end
             end

debug = 0;
if debug == 1
ct4step
nn2(i)
ctc0-ct4sym
ul(ict1(1))
src4(ict1(1))
ct0
odd_la(i)
x10(i)
xs2(end,i)
xs(1)
xs(end)
x1(1)
x1(ilast)
x1sym
xstep

figure(1)
%plot(ct4(ict1), src4(ict1), ct4, ul, 'k--', ct4, src4(:,ik2(i)), 'b--', ct4(ict2(ict3)), src42(ict3),'r')
plot(ct4(ict1), src4(ict1), ct4, ul, 'k--', ct4, src4, 'b--')
%axis([6500, 7100, -0.05, 0.05])
end

% integrate the source function over ct4, following formula (40) of Zaldarriaga et al.(1998)
             htable = src4(ict1) .* ul(ict1)';
             tetatl(i) = simpsint(ct0, ct4(end), htable) + (ct0 - ctc0 + ct4sym) * src4(ict1(1)) * ul(ict1(1));
             htable = src5(ict1) .* ul(ict1)';
             tetael(i) = simpsint(ct0, ct4(end), htable) + (ct0 - ctc0 + ct4sym) * src5(ict1(1)) * ul(ict1(1));
          else
             tetatl(i) = 0;
             tetael(i) = 0;
          end

       else  %if kb(i) <= 1e-5 *(1500/l)

% interpolate the spherical bessel function
          kctc0 = ctc0 * k2(i) * ones(1, lct4);
          kctc = k2(i) * ct4;
          xctc = kctc0 - kctc;
          jl = interpl(kx, table_bv, xctc)';
% integrate the source function following (12) and (13) of Seljak resp (18) of Zaldarriaga and Seljak
% the integration interval is limited to ct4 where the source <> 0
          htable = src4 .* jl;
          tetatl(i) = simpsint(ct4(1), ct4(end), htable);
          htable = src5 .* jl;
          tetael(i) = simpsint(ct4(1), ct4(end), htable);

       end  %if kb(i)

   end  % forloop k2

ta1 = tetatl;
f2 = 1;
%if (K ~= 0),  f2 = coth(pi*k2/sqrt(abs(K))); end
factor = sqrt((k2 .^2 - 4*K) ./ (k2 .^2 - K));
%factor = 1;
% perform the integration over k2min-k2max, following (9) of Seljak resp (19) of Zaldarriaga and Seljak
   cmbtable = f2 * factor .* (tetatl .^2) .* (k2 ./ (k2 .^2 - K)) .^ (2 - nprimtilt);
%ta1 = cmbtable;
%    cct.ctt(il) = l*(l + 1) * sqrt(K) * sum(cmbtable');
    cct.ctt(il) = l*(l + 1) * k1step * sum(cmbtable');
%   cct.ctt(il) = l*(l + 1) * simpsint(k2min, k2max, cmbtable')
   cmbtable = factor .* (tetael .^2) .* (k2 ./ (k2 .^2 - K)) .^ (2-nprimtilt);
   cct.cee(il) = gl^2 * l*(l + 1) * k1step * sum(cmbtable');
   cmbtable = factor .* (tetael .* tetatl) .* (k2 ./ (k2 .^2 - K)) .^ (2-nprimtilt);
   cct.cte(il) = l*(l + 1) * gl * k1step * sum(cmbtable');
end  % forloop la


text:
function imdb_cnn_train(imdb, opts, varargin)
% Train a CNN model on a dataset supplied by imdb

opts.lite = false ;
opts.numFetchThreads = 0 ;
opts.train.batchSize = opts.batchSize ;
opts.train.numEpochs = 25 ;
opts.train.continue = true ;
opts.train.useGpu = true ;
opts.train.prefetch = false ;
opts.train.learningRate = [0.001*ones(1, 10) 0.0001*ones(1, 10) 0.00001*ones(1,10)] ;
opts.train.expDir = opts.expDir ;
opts = vl_argparse(opts, varargin) ;

% -------------------------------------------------------------------------
%                                                    Network initialization
% -------------------------------------------------------------------------

net = initializeNetwork(imdb, opts) ;

% Initialize average image
if isempty(net.meta.normalization.averageImage),
    % compute the average image
    averageImagePath = fullfile(opts.expDir, 'average.mat') ;
    if exist(averageImagePath, 'file')
      load(averageImagePath, 'averageImage') ;
    else
      train = find(imdb.images.set == 1) ;
      bs = 256 ;
      fn = getBatchWrapper(net.meta.normalization, opts.numFetchThreads) ;
      for t=1:bs:numel(train)
        batch_time = tic ;
        batch = train(t:min(t+bs-1, numel(train))) ;
        fprintf('computing average image: processing batch starting with image %d ...', batch(1)) ;
        temp = fn(imdb, batch) ;
        im{t} = mean(temp, 4) ;
        batch_time = toc(batch_time) ;
        fprintf(' %.2f s (%.1f images/s)\n', batch_time, numel(batch)/ batch_time) ;
      end
      averageImage = mean(cat(4, im{:}),4) ;
      save(averageImagePath, 'averageImage') ;
    end

    net.meta.normalization.averageImage = averageImage ;
    clear averageImage im temp ;
end

% -------------------------------------------------------------------------
%                                               Stochastic gradient descent
% -------------------------------------------------------------------------
fn = getBatchWrapper(net.meta.normalization, opts.numFetchThreads) ;
[net,info] = cnn_train(net, imdb, fn, opts.train, 'conserveMemory', true) ;

% Save model
net = vl_simplenn_move(net, 'cpu');
saveNetwork(fullfile(opts.expDir, 'final-model.mat'), net);

% -------------------------------------------------------------------------
function saveNetwork(fileName, net)
% -------------------------------------------------------------------------
layers = net.layers;

% Replace the last layer with softmax
layers{end}.type = 'softmax';
layers{end}.name = 'prob';

% Remove fields corresponding to training parameters
ignoreFields = {'filtersMomentum', ...
                'biasesMomentum',...
                'filtersLearningRate',...
                'biasesLearningRate',...
                'filtersWeightDecay',...
                'biasesWeightDecay',...
                'class'};
for i = 1:length(layers),
    layers{i} = rmfield(layers{i}, ignoreFields(isfield(layers{i}, ignoreFields)));
end
classes = net.meta.classes;
normalization = net.meta.normalization;
save(fileName, 'layers', 'classes', 'normalization');

% -------------------------------------------------------------------------
function fn = getBatchWrapper(opts, numThreads)
% -------------------------------------------------------------------------
fn = @(imdb,batch) getBatch(imdb,batch,opts,numThreads) ;

% -------------------------------------------------------------------------
function [im,labels] = getBatch(imdb, batch, opts, numThreads)
% -------------------------------------------------------------------------
images = strcat([imdb.imageDir '/'], imdb.images.name(batch)) ;
im = imdb_get_batch(images, opts, ...
                            'numThreads', numThreads, ...
                            'prefetch', nargout == 0);
labels = imdb.images.label(batch) ;

% -------------------------------------------------------------------------
function net = initializeNetwork(imdb, opts)
% -------------------------------------------------------------------------
scal = 1 ;
init_bias = 0.1;
numClass = length(imdb.classes.name);
if ~isempty(opts.model)
    net = load(fullfile('data/models', opts.model)); % Load model if specified
    net.meta.normalization.keepAspect = opts.keepAspect;
    fprintf('Initializing from model: %s\n', opts.model);

    % Replace the last but one layer with random weights
    net.layers{end-1} = struct('type', 'conv', ...
                           'filters', 0.01/scal * randn(1,1,4096,numClass,'single'), ...
                           'biases', zeros(1, numClass, 'single'), ...
                           'stride', 1, ...
                           'pad', 0, ...
                           'filtersLearningRate', 10, ...
                           'biasesLearningRate', 20, ...
                           'filtersWeightDecay', 1, ...
                           'biasesWeightDecay', 0);

    % Last layer is softmaxloss (switch to softmax for prediction)
    net.layers{end} = struct('type', 'softmaxloss') ;

    % Rename classes
    net.meta.classes.name = imdb.classes.name;
    net.meta.classes.description = imdb.classes.name;

    % TODO: add dropout layers if initializing from previous models
    return;
end

% Else initial model randomly
net.layers = {} ;

% Block 1
net.layers{end+1} = struct('type', 'conv', ...
                           'filters', 0.01/scal * randn(11, 11, 3, 96, 'single'), ...
                           'biases', zeros(1, 96, 'single'), ...
                           'stride', 4, ...
                           'pad', 0, ...
                           'filtersLearningRate', 1, ...
                           'biasesLearningRate', 2, ...
                           'filtersWeightDecay', 1, ...
                           'biasesWeightDecay', 0) ;
net.layers{end+1} = struct('type', 'relu') ;
net.layers{end+1} = struct('type', 'pool', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'normalize', ...
                           'param', [5 1 0.0001/5 0.75]) ;

% Block 2
net.layers{end+1} = struct('type', 'conv', ...
                           'filters', 0.01/scal * randn(5, 5, 48, 256, 'single'), ...
                           'biases', init_bias*ones(1, 256, 'single'), ...
                           'stride', 1, ...
                           'pad', 2, ...
                           'filtersLearningRate', 1, ...
                           'biasesLearningRate', 2, ...
                           'filtersWeightDecay', 1, ...
                           'biasesWeightDecay', 0) ;
net.layers{end+1} = struct('type', 'relu') ;
net.layers{end+1} = struct('type', 'pool', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'normalize', ...
                           'param', [5 1 0.0001/5 0.75]) ;

% Block 3
net.layers{end+1} = struct('type', 'conv', ...
                           'filters', 0.01/scal * randn(3,3,256,384,'single'), ...
                           'biases', init_bias*ones(1,384,'single'), ...
                           'stride', 1, ...
                           'pad', 1, ...
                           'filtersLearningRate', 1, ...
                           'biasesLearningRate', 2, ...
                           'filtersWeightDecay', 1, ...
                           'biasesWeightDecay', 0) ;
net.layers{end+1} = struct('type', 'relu') ;

% Block 4
net.layers{end+1} = struct('type', 'conv', ...
                           'filters', 0.01/scal * randn(3,3,192,384,'single'), ...
                           'biases', init_bias*ones(1,384,'single'), ...
                           'stride', 1, ...
                           'pad', 1, ...
                           'filtersLearningRate', 1, ...
                           'biasesLearningRate', 2, ...
                           'filtersWeightDecay', 1, ...
                           'biasesWeightDecay', 0) ;
net.layers{end+1} = struct('type', 'relu') ;

% Block 5
net.layers{end+1} = struct('type', 'conv', ...
                           'filters', 0.01/scal * randn(3,3,192,256,'single'), ...
                           'biases', init_bias*ones(1,256,'single'), ...
                           'stride', 1, ...
                           'pad', 1, ...
                           'filtersLearningRate', 1, ...
                           'biasesLearningRate', 2, ...
                           'filtersWeightDecay', 1, ...
                           'biasesWeightDecay', 0) ;
net.layers{end+1} = struct('type', 'relu') ;
net.layers{end+1} = struct('type', 'pool', ...
                           'method', 'max', ...
                           'pool', [3 3], ...
                           'stride', 2, ...
                           'pad', 0) ;

% Block 6
net.layers{end+1} = struct('type', 'conv', ...
                           'filters', 0.01/scal * randn(6,6,256,4096,'single'),...
                           'biases', init_bias*ones(1,4096,'single'), ...
                           'stride', 1, ...
                           'pad', 0, ...
                           'filtersLearningRate', 1, ...
                           'biasesLearningRate', 2, ...
                           'filtersWeightDecay', 1, ...
                           'biasesWeightDecay', 0) ;
net.layers{end+1} = struct('type', 'relu') ;
net.layers{end+1} = struct('type', 'dropout', ...
                           'rate', 0.5) ;

% Block 7
net.layers{end+1} = struct('type', 'conv', ...
                           'filters', 0.01/scal * randn(1,1,4096,4096,'single'),...
                           'biases', init_bias*ones(1,4096,'single'), ...
                           'stride', 1, ...
                           'pad', 0, ...
                           'filtersLearningRate', 1, ...
                           'biasesLearningRate', 2, ...
                           'filtersWeightDecay', 1, ...
                           'biasesWeightDecay', 0) ;
net.layers{end+1} = struct('type', 'relu') ;
net.layers{end+1} = struct('type', 'dropout', ...
                           'rate', 0.5) ;

% Block 8
net.layers{end+1} = struct('type', 'conv', ...
                           'filters', 0.01/scal * randn(1,1,4096,numClass,'single'), ...
                           'biases', zeros(1, numClass, 'single'), ...
                           'stride', 1, ...
                           'pad', 0, ...
                           'filtersLearningRate', 1, ...
                           'biasesLearningRate', 2, ...
                           'filtersWeightDecay', 1, ...
                           'biasesWeightDecay', 0) ;

% Block 9
net.layers{end+1} = struct('type', 'softmaxloss') ;

% Other details
net.meta.normalization.imageSize = [227, 227, 3] ;
net.meta.normalization.interpolation = 'bicubic' ;
net.meta.normalization.border = 256 - net.meta.normalization.imageSize(1:2) ;
net.meta.normalization.averageImage = [] ;
net.meta.normalization.keepAspect = true ;

repo: DBCNN


text:
function vbrfa2011_testbed_measure(model, Q)

if nargin < 2 || isempty(Q)
  % Load default results
  disp('Loading default results..')
  files = {'testbed_results_gaussian_D=30_20110617', ...
           'testbed_results_multi-t_D=30_20110617', ...
           'testbed_results_ind-t_D=30_20110617', ...
           'testbed_results_laplace_D=30_20110617'};
  Q = load(files{model}, 'W_struct', 'W', 'CovW', 'X', 'CovX', 'Tau_struct');
end

% Load test data
disp('Loading test data..')
data_test = load(['/share/climate/jluttine/testbed/' ...
                  'testbed_vbrfa2011_testdata']);

% Compute predictive log-likelihood
disp('Computing predictive log-likelihood..')
samples = 10;
[M,N] = size(data_test.observations);
Ytest = data_test.observations;
Itest = ~isnan(Ytest);
loglike = nan(samples,1);
for i=1:samples

  % Sample a noise level
  tau = gamrnd(Q.W_struct.rho_struct.a_tau, 1/Q.W_struct.rho_struct.b_tau);

  % Sample a reconstruction
  for m=1:M
    W(:,m) = gaussian_rand(Q.W(:,m), 1/tau * Q.CovW(:,:,m));
  end
  for n=1:N
    X(:,n) = gaussian_rand(Q.X(:,n), Q.CovX(:,:,n));
  end
% $$$   W = mvnrnd(Q.W', 1/tau * Q.CovW)';
% $$$   X = mvnrnd(Q.X', Q.CovX)';
  F = W'*X;

  % Log-likelihood
  switch model
   case 1 % GAUSSIAN
    Z = Ytest(Itest) - F(Itest);
    loglike(i) = gaussian_logpdf(tau*(Z'*Z), ...
                                 0, ...
                                 0, ...
                                 -numel(Z)*log(tau), ...
                                 numel(Z));

   case 2 % MULTIVARIATE T
    alpha = Q.Tau_struct.a_u;
    beta = Q.Tau_struct.b_u;
    tau = tau * (alpha./beta);
    nu = 2*alpha;
    Z = Ytest - F;
    Z(~Itest) = 0;
    Z2 = tau .* dot(Z,Z,1);
    Nmv = sum(Itest,1)>0;
    loglike(i) = sum(t_logpdf(Z2(Nmv), ...
                              -log(tau(Nmv)), ...
                              nu(Nmv), ...
                              sum(Itest(:,Nmv),1)));

   case 3 % INDEPENDENT T
    [I,J] = find(Itest);
    Z2 = tau * (Ytest(Itest) - F(Itest)).^2;
    loglike(i) = sum(t_logpdf(Z2, ...
                              -log(tau), ...
                              Q.Tau_struct.nu(I), ...
                              1));

   case 4 % LAPLACE
    loglike(i) = sum(laplace_logpdf(Ytest(Itest), ...
                                    F(Itest), ...
                                    sqrt(tau)));
  end

end
%loglike = loglike / samples;
fprintf('Mean predictive log-density: %.4e\n', mean(loglike));
bias = max(loglike);
fprintf('Log predictive density: %.4e\n', log(mean(exp(loglike-bias))) + bias);


text:
%
% read_smooth_eccen.m
%
% Original Author: Bruce Fischl
% CVS Revision Info:
%    $Author: nicks $
%    $Date: 2011/03/02 00:04:12 $
%    $Revision: 1.3 $
%
% Copyright © 2011 The General Hospital Corporation (Boston, MA) "MGH"
%
% Terms and conditions for use, reproduction, distribution and contribution
% are found in the 'FreeSurfer Software License Agreement' contained
% in the file 'LICENSE' found in the FreeSurfer distribution, and here:
%
% https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
%
% Reporting: freesurfer@nmr.mgh.harvard.edu
%

if (exist('nfig') == 0)
   nfig = 0;
end;

% Should run 'left or right' scripts to set up fname_* variable,s
% and set nfig so that plotting starts in the right location.
fclose('all');
fid=fopen(fname_patch,'rt');
s = fgetl(fid);
s = fgetl(fid);
[poly,count] = sscanf(s,'%d');
numvert = poly(1);
numquad = poly(2);
vertex_coordinates = zeros(3,1);
face_index = zeros(4,1);
vertx_list = zeros(numvert,3);
face_list = zeros(numquad,5);
tic;
for vert = 1:1:(numvert),
    s = fgetl(fid);
    vertx = sscanf(s,'%d');
    s = fgetl(fid);
    vertx_coordinates = sscanf(s,'%f');
    vertx_list(vert,:) = [vertx_coordinates(1:2)' vertx];
end;
toc
tic;
for face = 1:1:(numquad),
    s = fgetl(fid);
    facenum = sscanf(s,'%d');
    s = fgetl(fid);
    face_vertx = sscanf(s,'%d');
    face_list(face,:) = [face_vertx' facenum];
end;
toc
fclose(fid);
full_vertx=zeros(max(max(vertx_list))+1,3);
full_vertx(vertx_list(:,3)+1,1)=vertx_list(:,1);
full_vertx(vertx_list(:,3)+1,2)=vertx_list(:,2);
mesh_real=File2Var(fname_real);
vertx_values=zeros(max(max(vertx_list))+1,1);
vertx_values(vertx_list(:,3)+1)=mesh_real(:,5);
mesh_imag=File2Var(fname_imag);
vertx_complex=zeros(max(max(vertx_list))+1,1);
vertx_complex(vertx_list(:,3)+1)=mesh_imag(:,5);
v_complex = vertx_values(:,1) + i*vertx_complex;
%v_phase=angle(v_complex');

strpwd=pwd;
subplot(2,1,nfig+1);
title([strpwd((length(strpwd)-25):length(strpwd)) ' ' fname_real]);
p_eccen_handle=patch('Vertices',full_vertx,'Faces',face_list(:,1:4)+1,'FaceVertexCData',angle(v_complex),'FaceColor','interp','EdgeColor','none');
%colormap(rgb(256));
colorbar;

short_complex = mesh_real(:,5) + (i*mesh_imag(:,5));
new_val = avg_vertx(vertx_list(:,1),vertx_list(:,2),short_complex,4);
new_idx = insert_vals(new_val,vertx_list(:,3));

subplot(2,1,nfig+2);
p_eccen_handle=patch('Vertices',full_vertx,'Faces',face_list(:,1:4)+1,'FaceVertexCData',angle(new_idx),'FaceColor','interp','EdgeColor','none');
title(['Smoothed: 'strpwd((length(strpwd)-25):length(strpwd)) ' ' fname_real]);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [history, stopReason] = lars(yin, xin, XTX, type, stopCriterion, regularization, trace, quiet)
% %{
% This program implements "LARS" algorithm and its lasso modification
% introduced by Efron et. al. 2003.  Read the paper to understand codes
% of this function.  Each line of this file has corresponding equation
% number in Efron et. al. 2003 for reader's convenience.
%
%
% *** CAUTION
% history(1).mu_OLS contains original 'yin' to provide convinience in
% writing user defined stop criterion function. Actual mu_OLS of the first
% step should be just mean of yin which is a simple array of copy of
% history(1).mu. So, if history(1).mu_OLS contains that information, it is
% redundant. In fact, this is also contained in history(1).b, which is a
% bias of the output. Therefore, to provide more information to user who
% want to write his/her own stop criterion function, history(1).mu_OLS
% contains yin.
%
%
% Example 1: moderate size x.
%     stopCrioterion = {};
%     stopCrioterion{1,1} = 'maxKernels';
%     stopCrioterion{1,2} = 100;
%
%     XTX = lars_getXTX(x_original);      % this takes long time.
%     sol = lars(y, x, 'lasso', XTX, stopCriterion);
%
% Example 2: very small size x, or a really really big size x
%     stopCrioterion = {};
%     stopCrioterion{1,1} = 'maxKernels';
%     stopCrioterion{1,2} = 100;
%
%     sol = lars(y, x, 'lasso', XTX, stopCriterion);
%
% Note:
%     Users can add any kind of stop criterion by editing
%     the corresponding portion of this file.  See the code
%     for existing examples.
%
% Note 2:
%     This m-file does not implement routine for missing data.
% %}
%

global USING_CLUSTER;
global RESOLUTION_OF_LARS;
global REGULARIZATION_FACTOR;
lars_init();

regularization_factor = REGULARIZATION_FACTOR; % Tikhonov regularization factor (or the ridge regression factor)
                              % -> This should be small enough in this case to get
                              %    a reasonable pseudoinverse.

stopReason = {};

%%% Check parameters
if length(yin)==0 | length(xin)==0
    warning('\nInput or Output has zero length.\n');
    history.active_set = [];
    stopReason{1} = 'Parameter error';
    stopReason{2} = 0;
    return;
end
if size(yin,1) ~= size(xin,1)
    warning('\nSize of y does not match to that of x.\n');
    history.active_set = [];
    stopReason{1} = 'Parameter error';
    stopReason{2} = 0;
    return;
end
if ~strcmp(type, 'lasso') & ~strcmp(type, 'lars') & ~strcmp(type, 'forward_stepwise')
    warning('\nUnknown type of regression.\n');
    history.active_set = [];
    stopReason{1} = 'Parameter error';
    stopReason{2} = 0;
    return;
end
if strcmp(type, 'forward_stepwise')
    warning('\nForward_stepwise is not implemented.\n');
    history.active_set = [];
    stopReason{1} = 'Parameter error';
    stopReason{2} = 0;
    return;
end

if exist('regularization','var') & ~isempty(regularization)
    regularization = 10;
else
    regularization  = 0;
end

if ~exist('trace','var') | isempty(trace)
    trace=0;
end

if ~exist('quiet','var') | isempty(quiet)
    quiet=0;
elseif quiet==1
    trace=0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Data preparation

% Program automatically centers and standardizes predictors.
if ~exist('XTX','var')
    XTX=[];
end
no_xtx = 0;
if ~isempty(XTX)
    if ~quiet & trace >=0
        fprintf('\nLars is using the provided xtx.\n');
    end
elseif size(xin,2)^2 > 10^6
    if ~quiet & trace >=0
        fprintf('Too large matrix (size(x,2)^2 > 10^6). lars will not pre-calculate xtx.\n');
    end
    no_xtx      = 1;
    XTX         = lars_getXTX(xin,no_xtx);
else
%     fprintf('\nCalculating xtx.\n');
    XTX         = lars_getXTX(xin);
end

x               = XTX.x;            % normalized xin
mx              = XTX.mx;           % mean xin
sx              = XTX.sx;           % length of each column of xin
ignores         = XTX.ignores;      % indices for constant terms
all_candidate   = XTX.all_candidate;% indices for all possible columns
if ~no_xtx
    xtx         = XTX.xtx;          % xtx matrix
    dup_columns = XTX.dup_columns;  % duplicated columns which will be automatically ignored.
end

my              = mean(yin);
y               = yin-my;

n               = size(x,1);        % # of samples
m               = size(x,2);        % # of predictors

% Now, we can determine the maximum number of kernels.
%maxKernels = min(maxKernels, min(size(xin,1)-1, length(all_candidate)));
%maxKernels = min(maxKernels, min(rank(xin), length(all_candidate)));
existMaxKernels = 0;
existMSE        = 0;
for is = 1:size(stopCriterion,1)
    if strcmp(stopCriterion{is,1},'maxKernels')
        existMaxKernels = 1;
%        stopCriterion{is,2} = min(stopCriterion{is,2}, min(size(xin,1)-1, length(all_candidate)));
        stopCriterion{is,2} = min(stopCriterion{is,2}, min(rank(xin), length(all_candidate)));
        if stopCriterion{is,2}<1
            warning('Max Kernel is less than 1. It must be larger than 0.\n');
            stopCriterion{is,2} = 1;
        end
    end
    if strcmp(stopCriterion{is,1},'MSE')
        existMSE        = 1;
        if stopCriterion{is,2}<1.0e-10
            warning('Maximum MSE is too small. Automatically set to 1.0e-10\n');
            stopCriterion{is,2} = 1.0e-10;
        end
    end
end
if ~existMaxKernels
    is = size(stopCriterion,1);
    stopCriterion{is+1,1} = 'maxKernels';
%    stopCriterion{is+1,2} = min(size(xin,1)-1, length(all_candidate));  % Stop when size of active set is data.maxKernels.
    stopCriterion{is+1,2} = min(rank(xin), length(all_candidate));  % Stop when size of active set is data.maxKernels.
end
if ~existMSE
    is = size(stopCriterion,1);
    stopCriterion{is+1,1} = 'MSE';
    stopCriterion{is+1,2} = 1.0e-10;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialization

active      = [];               % active set
inactive    = all_candidate;    % inactive set

mu_a        = zeros(n,1);       % current estimate (eq. 2.8)
mu_a_plus   = 0;                % next estimate (eq. 2.12)
mu_a_OLS    = 0;                % OLS estimate  (eq. 2.19)

beta        = zeros(1,size(x,2));
beta_new    = beta;
beta_OLS    = beta;

history.active_set          = [];
history.add                 = [];
history.drop                = [];
history.beta_norm           = [];
history.beta                = [];
history.b                   = my;
history.mu                  = my;
history.beta_OLS_norm       = [];
history.beta_OLS            = [];
history.b_OLS               = my;
history.mu_OLS              = my*ones(size(yin));
history.MSE                 = sum(y.^2)/length(y);
history.R_square            = 0;
history.resolution_warning  = [];

if var(yin)==0
    stopReason{1} = 'zeroVarY';
    stopReason{2} = var(yin);
    return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Main loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c               = 0;                % correlation vector
C_max           = max(abs(c));
C_max_ind       = [];
C_max_ind_pl    = [];
drop            = [];               % used for 'lasso'
k               = 1;                % iteration index
if ~quiet & trace >= 0
    fprintf('Active predictors / total  :::::  Current iteration\n                            ');
end
while 1

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%           Exit Criterions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if exist('stopCriterion','var')     % If there is any stop criterion...
        % Any of these is satisfied, algorithm stops.
        %
        % Other criterions: Make your own cost function.

        for is = 1:size(stopCriterion,1)  % there can be a number of stop criterions.

            % Default Criterions.
            % Maximum number of consecutive drops or maximum number of drops within a window.
            if strcmp(stopCriterion{is,1},'maxDrops')
                drop_window = stopCriterion{is,2}(1);   % number of backward history to be considerd. 0 means all history.
                if drop_window==0
                    drop_window=k;
                end
                drop_n      = min(drop_window,stopCriterion{is,2}(2));   % number of maximum drops within the window.
                drop_vector = [];
                for z = max(k-drop_window+1,1):k
                    drop_vector=[drop_vector, history(z).drop];
                end
                if length(drop_vector)>=drop_n
                    stopReason{1} = 'maxDrops';
                    stopReason{2} = drop_n;
                    break;
                end
            end
            % Maximum number of kernels.
            if strcmp(stopCriterion{is,1},'maxKernels')
                if length(active) >= min(stopCriterion{is,2}, min(size(xin,1)-1, length(all_candidate)))
                    stopReason{1} = 'maxKernels';
                    stopReason{2} = length(active);
                    break;
                end
            end
            % Maximum number of iterations.
            if strcmp(stopCriterion{is,1},'maxIterations')
                if k >= stopCriterion{is,2}
                    stopReason{1} = 'maxIterations';
                    stopReason{2} = k;
                    break;
                end
            end
            % MSE.
            if strcmp(stopCriterion{is,1},'MSE')
                if history(k).MSE <= stopCriterion{is,2}
                    stopReason{1} = 'MSE';
                    stopReason{2} = history(k).MSE;
                    break;
                end
            end

            % User defined stop criterion.
            if strcmp(stopCriterion{is,1},'userDefinedCriterion')
                fhandle = stopCriterion{is,2}.fhandle;
                r_fhandle = fhandle(history, stopCriterion{is,2}.data);
                if r_fhandle.stop
                    stopReason{1} = 'userDefinedCriterion';
                    stopReason{2} = r_fhandle;
                    break;
                end
            end

        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    end % end of stop criterion checking :  if exist('stopCriterion','var')
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if length(stopReason)>0 % if there is any reason of stopping the algorithm, exit loop.
        break;
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%           LARS Algorithm
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %if ~USING_CLUSTER
    %fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%5d/%5d    :::::   %5d',length(active),length(all_candidate),k);
    %end
    % start of algorithm
    c                   = x'*(y-mu_a);                      % eq 2.8
    [C_max,C_max_ind]   = max(abs(c(inactive)));            % eq 2.9
    C_max_ind           = inactive(C_max_ind);
    % active          = sort(union(active,C_max_ind));      % If there is no machine limit, this can be used.
    % But because of machine limit, there can be multiple new predictors.
    % This dramatically improves the overall precision of the result,
    % and speeds up the whole process.
    C_max_ind_pl        = abs(c(inactive))>C_max-RESOLUTION_OF_LARS;
    C_max_ind_pl        = inactive(C_max_ind_pl);
    active          = sort(union(active,C_max_ind_pl));
    inactive        = setdiff(all_candidate, active);   % eq 2.9
    if strcmp(type,'lasso')
        if ~isempty(drop) & length(find(drop==C_max_ind))==0           % If there is a drop, that must have the maximum correlation in inactive set.
            if ~quiet & trace >=0
                fprintf('\n');
                warning('Dropped item and index of maximum correlation is not the same. But it is being ignored here...');
                fprintf('\n                            ');
            end
            active(find(active==C_max_ind))=[];
        end
        if ~isempty(drop)
            C_max_ind   = [];
            C_max_ind_pl= [];
        end
        active          = setdiff(active,drop);             % eq 3.6
        inactive        = sort(union(inactive,drop));       % eq 3.6
    end

    s       = sign(c(active));              % eq 2.10
    xa      = x(:,active).*repmat(s',n,1);  % eq 2.4
    if ~no_xtx
        ga  = xtx(active,active).*(s*s');   % eq 2.5
    else
        ga  = xa'*xa;                       % eq 2.5
    end
    if regularization > 2
        ga  = ga+eye(length(ga))*regularization_factor; % This routine will make the test below
    end
    invga   = ga\eye(size(ga,1));           % eq 2.5
    aa      = sum(sum(invga))^(-1/2);       % eq 2.5
    wa      = aa*sum(invga,2);              % eq 2.6
    ua      = xa*wa;                        % eq 2.6

    % test using eq 2.7
    test_1      = xa'*ua;
    test_2      = aa*ones(size(test_1));
    test_1_2    = sum(sum(abs(test_1-test_2)));
    test_3      = norm(ua) - 1;

    history(k+1).resolution_warning=0;
    if test_1_2 > RESOLUTION_OF_LARS*100 | abs(test_3 ) > RESOLUTION_OF_LARS*100
        if regularization <=2
            if ~quiet & trace>0
                fprintf('\n');
                warning('Eq 2.7 test failure.');
                fprintf('\n                            ');
            end
            regularization = regularization + 1;
            if regularization > 2
                if ~quiet & trace>0
                    fprintf('\n');
                    warning('Lots of Eq 2.7 test failure. Regularization will be applied from now on.');
                    fprintf('\n                            ');
                end
            end
        end
        history(k+1).resolution_warning=1;
    end

    a               = x'*ua;                % eq 2.11
    tmp_1           = (C_max - c(inactive))./(aa - a(inactive));
    tmp_2           = (C_max + c(inactive))./(aa + a(inactive));
    tmp_3           = [tmp_1, tmp_2];
    tmp             = tmp_3(find(tmp_3>0));
    gamma = min(tmp);                       % eq 2.13
    if length(gamma)==0 % if this is the last step (i.e. length(active)==maxKernels)
        gamma = C_max/aa;                   % eq 2.19, eq 2.21 and 5 lines below eq 2.22
    end

    d               = zeros(1,m);
    d(active)       = s.*wa;

    if length(find(d(active)==0))
        fprintf('\n');
        warning('Something wrong with vector d: eq 3.4.');
        fprintf('\n                            ');
    end

    tmp             = zeros(1,m);
    tmp(active)     = -1*beta(active)./d(active);   % eq 3.4
    tmp2            = tmp(find(tmp>0));

    drop            = [];
    gamma_tilde = inf;                              % eq 3.5
    if ~isempty(tmp2) & gamma >= min(tmp2)
        gamma_tilde = min(tmp2);                    % eq 3.5
        drop        = find(tmp==gamma_tilde);       % eq 3.6
    end

    if strcmp(type, 'lars')
        mu_a_plus   = mu_a + gamma*ua;              % eq 2.12
        beta_new    = beta + gamma*d;               % eq 3.3
        drop        = [];
    elseif strcmp(type, 'lasso')
        mu_a_plus   = mu_a + min(gamma, gamma_tilde)*ua;    % eq 3.6
        beta_new    = beta + min(gamma, gamma_tilde)*d;     % eq 3.3
        active      = setdiff(active,drop);                 % eq 3.6
        inactive    = setdiff(all_candidate,active);
        beta_new(drop) = 0;
    elseif strcmp(type, 'forward_stepwise')
        drop        = [];
        error('forward.stepwise has not been implemented yet.');
        return;
    end

    mu_a_OLS    = mu_a + C_max/aa*ua;          % eq 2.19, 2.21
    beta_OLS    = beta + C_max/aa*d;           % eq 2.19, 2.21
    MSE         = sum((y - mu_a_OLS).^2)/length(y);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % update and save
    mu_a = mu_a_plus;
    beta = beta_new;

    % history with scale correction
    k = k+1;
    history(k).active_set   = active;
    history(k).drop         = drop;
    history(k).add          = C_max_ind_pl;
    history(k).beta_norm    = beta(active);
    history(k).beta         = beta(active)./sx(active);
    history(k).b            = my - sum(mx./sx.*beta);
    history(k).mu           = xin * (beta./sx)' + history(k).b;
    history(k).beta_OLS_norm= beta_OLS(active);
    history(k).beta_OLS     = beta_OLS(active)./sx(active);
    history(k).b_OLS        = my - sum(mx./sx.*beta_OLS);
    history(k).mu_OLS       = xin * (beta_OLS./sx)' + history(k).b_OLS;
    history(k).MSE          = MSE;
    history(k).R_square     = 1 - var(yin  - history(k).mu_OLS)/var(yin);

    % exit if exact mathing is achieved.
    if abs(C_max/aa - min(gamma,gamma_tilde)) < RESOLUTION_OF_LARS
        stopReason{1}   = 'ExactMatching';
        stopReason{2}   = 0;
        break;
    end

end % end of while loop
if ~quiet & trace >=0
    fprintf('\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%5d/%5d    :::::   %5d\n',length(active),length(all_candidate),k);
end

return;