thermalblock.m

Go to the documentation of this file.

function res = thermalblock(step, params)
% thermalblock example
%
%  step 1: detailed simulation with different parameters
%  step 2: reduced basis simulation
%  step 3: demo_rb_gui
%  step 4: time measurements comparing OO model with old model.

% B. Haasdonk 27.2.2011

if nargin <= 1
  params.numintervals_per_block = 5;
end

if nargin<1
  step = 1;
end;

switch step
 case 1 % different detailed simulations

  disp('detailed simulation:');
  descr      = thermalblock_model(params);
  dmodel     = gen_detailed_model(descr);
  %  dmodel = thermalblock_model;

  model_data = gen_model_data(dmodel);
  dmodel = set_mu(dmodel,[1,1,1,1,1,1,1,1,1]);

  sim_data = detailed_simulation(dmodel,model_data);
  plot_params.title = ['detailed solution , \mu_i=1, output s=',...
                       num2str(sim_data.s)];
  figure, plot_sim_data(dmodel,model_data,sim_data,plot_params);
  dmodel = set_mu(dmodel,[0.1,0.1,1,1,0.1,0.1,1,0.1,1]);

  sim_data = detailed_simulation(dmodel,model_data);
  plot_params.title = ['detailed solution , \mu_i=0.1/1, output s=',...
                       num2str(sim_data.s)];
  figure, plot_sim_data(dmodel,model_data,sim_data,plot_params);

 case 2 % rb treatment and comparison to detailed

  % for non-snapshot parameter nice error behaviour is observed
  % for snapshot parameter error is numerically zero, but
  % estimators due to numerical noise not zero.

  disp('detailed simulation:');
  descr            = thermalblock_model(params);
  [dmodel, rmodel] = gen_models(descr);
%  dmodel            = thermalblock_model;
%  rmodel            = dmodel;
  model_data       = gen_model_data(dmodel);
  % non-snapshot parameter, error <= Delta
  dmodel = set_mu(dmodel,[0.5,1,1,1,1,1,1,1,1]);
  % non-snapshot parameter with excellent Delta==error!!!
  %model = model.set_mu(model,[1,1,1,1,1,1,1,1,1]);
  % snapshot parameter with error, Delta approximately 0.
  %model = model.set_mu(model,[0.1,1,1,1,1,1,1,1,1]);
  sim_data = detailed_simulation(dmodel,model_data);
  sim_data.mu = get_mu(dmodel);
  plot_params.title = 'detailed solution';
  disp(['output s(mu) = ',num2str(sim_data.s)]);
  figure, plot_sim_data(dmodel,model_data,sim_data,plot_params);
  disp('offline computations I (gen_detailed_data: basis):');
  detailed_data = gen_detailed_data(rmodel,model_data);
  disp('offline computations II (gen_reduced_data: operator components):');
  reduced_data = gen_reduced_data(rmodel,detailed_data);
%  model.N = 1;
%  reduced_data = model.reduced_data_subset(model,reduced_data);

  disp('online reduced simulation:');
  rmodel = set_mu(rmodel, sim_data.mu);
  rb_sim_data = rb_simulation(rmodel,reduced_data);
  rb_sim_data = rb_reconstruction(rmodel, detailed_data, rb_sim_data);
  plot_params.title = 'reduced solution';
  figure, plot_sim_data(rmodel,model_data,rb_sim_data,plot_params);
  disp(['output sN(mu) = ',num2str(rb_sim_data.s)]);

  disp('-------------------------------------------');
  eh = sim_data.uh - rb_sim_data.uh;
  l2_err = fem_l2_norm(eh);
  h10_err = fem_h10_norm(eh);
  disp(['L2-error: ',num2str(l2_err)]);
  disp(['H10-error: ',num2str(h10_err)]);
  disp(['error s(mu)-sN(mu) = ',num2str(sim_data.s-rb_sim_data.s)]);
  disp(['error estimator Delta_u(mu) = ',num2str(rb_sim_data.Delta)]);
  disp(['error estimator Delta_s(mu) = ',num2str(rb_sim_data.Delta_s)]);

  % save detailed_data for demo_rb_gui:
  save(fullfile(rbmatlabresult,'thermal_block_offline_data.mat'),...
       'dmodel','detailed_data');

 case 3 % demo_rb_gui

  fn = fullfile(rbmatlabresult,'thermal_block_offline_data.mat');
  if ~exist(fn) % call step 2 first
    disp('calling thermalblock(2) for basis generation')
    thermalblock(2);
  end;

  disp('demo_rb_gui:')
  descr            = thermalblock_model;
  [dmodel, rmodel] = gen_models(descr);
  load(fn);
  %dmodel = thermalblock_model;
  %rmodel = dmodel;
  rmodel.N = size(detailed_data.RB,2);
  %model_data       = gen_model_data(dmodel);
  %detailed_data    = gen_detailed_data(dmodel,model_data);
  plot_params.axis_tight = 1;
  plot_params.yscale_uicontrols = 0.5;
  demo_rb_gui(rmodel,detailed_data,[],plot_params);

 case 3.5
   fn=fullfile(rbmatlabresult, 'thermal_block_offline_data.mat');
   tmp = load(fn);
   detailed_data = tmp.detailed_data;
   descr = thermalblock_model(params);
   [dmodel, rmodel] = gen_models(descr);
   reduced_data = gen_reduced_data(rmodel, detailed_data);
   pr = ParameterSampling.Random(100);
   pr.init_sample(rmodel);
   M = pr.sample;

   nmus = size(M, 1);
   rtimes = zeros(nmus, 1);
   dtimes = zeros(nmus, 1);
   errs   = zeros(nmus, 1);
   for i = 1:size(M,1)
     rmodel = set_mu(rmodel, M(i,:));
     rntelapsed = 0;
     dntelapsed = 0;
     for j = 1:10
       tic;
       rbs = rb_simulation(rmodel, reduced_data);
       ttemp = toc;
       rrbs = rb_reconstruction(rmodel, detailed_data, rbs);
       rntelapsed = rntelapsed + ttemp;
       tic;
       dbs = detailed_simulation(dmodel, detailed_data.model_data);
       ttemp = toc;
       dntelapsed = dntelapsed + ttemp;
     end
     rtimes(i) = rntelapsed/10;
     dtimes(i) = dntelapsed/10;
     errs(i) = fem_l2_norm(rrbs.uh - dbs.uh);
   end

   save(fullfile(rbmatlabresult, 'thermalblock_timings'), 'errs', 'rtimes', 'dtimes');


 case 4 % compare with old version

  disp(['This step will only work after providing nice descr for' ...
       'thermalblock!']);
  return

  omodel         = thermalblock_model_old;
  omodel_data    = gen_model_data(omodel);
  odetailed_data = gen_detailed_data(omodel, omodel_data);
  oreduced_data  = gen_reduced_data(omodel, odetailed_data);
  pr = ParameterSampling.Random(100);
  pr.init_sample(omodel);
  M = pr.sample;

  otelapsed = 0;
  for i = 1:size(M,2)
    omodel = omodel.set_mu(omodel, M(i,:));
    for j = 1:10
      tic;
      rb_simulation(omodel, oreduced_data);
      ttemp = toc;
      otelapsed = otelapsed + ttemp;
    end
  end
  disp(num2str(otelapsed));

  descr           = thermalblock_model;
  [dmodel,rmodel] = gen_models(descr);
  model_data      = gen_model_data(dmodel);
  detailed_data   = gen_detailed_data(rmodel, model_data);
  reduced_data    = gen_reduced_data(rmodel, detailed_data);

  ntelapsed = 0;
  for i = 1:size(M,2)
    descr = omodel.set_mu(descr, M(i,:));
    rmodel = set_mu(rmodel, M(i,:));
    for j = 1:10
      tic;
      rb_simulation(rmodel, reduced_data);
      ttemp = toc;
      ntelapsed = ntelapsed + ttemp;
     end
  end
  disp(num2str(ntelapsed));

  ntelapsed2 = 0;
  for i = 1:size(M,2)
    descr = omodel.set_mu(descr, M(i,:));
    rmodel = set_mu(rmodel, M(i,:));
    descr = rmodel.descr;
    for j = 1:10
      tic;
      omodel.rb_simulation(descr, reduced_data);
      ttemp = toc;
      ntelapsed2 = ntelapsed2 + ttemp;
     end
  end
  disp(num2str(ntelapsed2));
  ntelapsed3 = 0;
  rd = struct(reduced_data);
  for i = 1:size(M,2)
    descr = omodel.set_mu(descr, M(i,:));
    rmodel = set_mu(rmodel, M(i,:));
    for j = 1:10
      tic;
      rb_simulation(rmodel, rd);
      ttemp = toc;
      ntelapsed3 = ntelapsed3 + ttemp;
     end
  end
  disp(num2str(ntelapsed3));

case 5

  % compute a bunch of error estimates and true errors...
  fn = fullfile(rbmatlabresult,'thermal_block_offline_data.mat');
  if ~exist(fn, 'file') % call step 2 first
    disp('calling thermalblock(2) for basis generation')
    thermalblock(2);
  end;

  descr            = thermalblock_model;
  [dmodel, rmodel] = gen_models(descr);
  model_data = gen_model_data(dmodel);
  load(fn);

  reduced_data = gen_reduced_data(rmodel, detailed_data);

  nr = 1000;

  ps=ParameterSampling.Random(nr);

  ps.init_sample(rmodel);

  deltas = zeros(1, nr);
  residuals = zeros(1, nr);
  errs = zeros(1, nr);
  h=waitbar(0,'Please wait..');
  for i=1:nr
    set_mu(rmodel, ps.sample(i, :));
    set_mu(dmodel, ps.sample(i, :));
    sim_data = detailed_simulation(dmodel, model_data);
    rb_sim_data = rb_simulation(rmodel, reduced_data);
    rb_sim_data = rb_reconstruction(rmodel, detailed_data, rb_sim_data);

    waitbar(i/nr);

    deltas(i) = rb_sim_data.Delta;
    residuals(i) = rb_sim_data.res_norm;
    df = sim_data.uh - rb_sim_data.uh;
    errs(i) = fem_l2_norm(df);
  end
  close(h);

  save(fullfile(rbmatlabresult, 'deltas_and_errs'), 'deltas', 'errs', 'residuals');

case 6
  dfn = fullfile(rbmatlabresult,'thermal_block_offline_data.mat');
  fn = fullfile(rbmatlabresult, 'deltas_and_errs.mat');
  if ~exist(fn, 'file') % call step 5 first
    disp('calling thermalblock(5) for estimator generation')
    thermalblock(5);
  end

  descr            = thermalblock_model;
  [dmodel, rmodel] = gen_models(descr);
  model_data = gen_model_data(dmodel);
  tmp = load(dfn);
  detailed_data = tmp.detailed_data;

  tmp = load(fn);
  deltas = tmp.deltas;
  errs = tmp.errs;
  residuals = tmp.residuals;

  figure
  % plot estimates and residuals against "true" errors
  subplot(2,1,1);
  plot(deltas, errs, 'g*', residuals, errs, 'b+'); set(gca, 'XScale', 'log');
  legend('estimates', 'residuals');
  ylabel('error');
  subplot(2,1,2);
  plot3(deltas, residuals, errs, 'g*');
  xlabel('estimates'); ylabel('residuals'); zlabel('errors');

  % register UQ classes
  cd ~/workspace/uq_rb_test/
  startup

  gFunc = @log10;
  gInvFunc = @(X) 10.^X;
  hFunc = @(X) log10(X);
  hInvFunc = @(X) 10.^(X);
  gPlotFunc = gInvFunc;
  hPlotFunc = hInvFunc;

  % models for Gaussian processes
  deltaModel = uqrb.base.ErrorModel(gFunc(deltas'), hFunc(errs'));
  resModel = uqrb.base.ErrorModel(gFunc(residuals'), hFunc(errs'));
  combiModel = uqrb.base.ErrorModel([gFunc(residuals'),gFunc(deltas')], hFunc(errs'));
  N(1) = 100; N(2) = 400;
  M = 2;
  deltaXns = cell(2,1);
  resXns = cell(2,1);
%  dbasis = cell(2,1);
%  rbasis = cell(2,1);
  dpb = cell(2,1);
  rpb = cell(2,1);
  cpb = cell(2,1);
  for i=1:2
    deltaXns{i} = deltaModel.samplePoints(N(i));
    resXns{i} = resModel.samplePoints(N(i));
    %dbasis{i} = uqrb.rvm.basisFunction.Radial(deltaXns{i}, [], deltaModel.domain);
    %rbasis{i} = uqrb.rvm.basisFunction.Radial(resXns{i}, [], resModel.domain);
    dpb{i} = uqrb.gp.kernel.Gaussian();
%    dpb{i} = uqrb.rvm.basisFunction.Polynomial(M, deltaModel.domain);
    rpb{i} = uqrb.rvm.basisFunction.Polynomial(M, resModel.domain);
    cpb{i} = uqrb.rvm.basisFunction.Mixed(4, combiModel.domain);
  end
%  keyboard
  GPs = {};
%  GPs = [ GPs, {uqrb.rvm.Inference(cpb{1}, combiModel)} ];
  GPs = [ GPs, {uqrb.gp.Inference(dpb{1}, deltaModel)} ];
  GPs = [ GPs, {uqrb.rvm.Inference(rpb{1}, deltaModel)} ];
  %GPs{1,1} = uqrb.rvm.Inference(dpb{1}, deltaModel);
  %GPs = [ GPs, {uqrb.rvm.Inference(rpb{1}, deltaModel)} ];
  %GPs = [ GPs, {uqrb.rvm.Inference(rpb{2}, deltaModel)} ];
  %GPs{1,4} = uqrb.rvm.Inference(dpb{2}, deltaModel);
  GPs = [ GPs, {uqrb.gp.Inference(dpb{2}, deltaModel)} ];
  GPs = [ GPs, {uqrb.rvm.Inference(rpb{2}, resModel)} ];
%  GPs = [ GPs, {uqrb.rvm.Inference(cpb{2}, combiModel)} ];
  %GPs = [ GPs, {uqrb.rvm.Inference(rpb{1}, resModel)} ];
  %GPs = [ GPs, {uqrb.rvm.Inference(rpb{2}, resModel)} ];
  %GPs = [ GPs, {uqrb.rvm.Inference(rpb{2}, resModel)} ];

  % declare variables for tikz plots
  overviewTitle = {'filename', 'abscissa', 'Mstart', 'M', 'sigma', 'basistype'};
  overviewData  = cell(6, numel(GPs));
  overviewFn    = fullfile(rbmatlabresult, 'overview.txt');
  validationTitle = {'vxns', 'muCurve', 'Sigma1', 'Sigma2', 'Sigma3', 'Sigma4', 'vns'};

  counter = 1;
%  alpGp   = cell(2,4);
%  s2s     = cell(2,4);

  for i=1:length(GPs)
    gp = GPs{i};
    dim = gp.model.dimension;
    if gp.model == deltaModel
      abscissa = 'estimates';
    elseif gp.model == resModel
      abscissa = 'residual';
    elseif gp.model == combiModel
      abscissa = 'combined';
    end

    stemFn = fullfile(rbmatlabresult, ['plots_', num2str(i)]);
    trainingFn   = [stemFn, '_training.dat'];

    gp.hyperMu.sigma2 = 1e-1;
    %      disp(num2str(gp.basis.M));
    n = N(1 + (i>2));
    [xns, tns] = gp.model.evaluateAtSamplePoints(1000);
    gp.model.domain = [ min(xns)-0.5; max(xns)+0.5 ]';
%    dwidth = gp.model.domain(2) - gp.model.domain(1);
    %      vxns = xns(n+1:1000);
    %      [vxns, si] = sort(vxns);
    NSP = [200, 20];
    vxns = gp.model.uniformSamplePoints(NSP(dim));
%    vxns = ((gp.model.domain(1) + 0.5):dwidth/200:(gp.model.domain(2) - 0.5))';
    xns = xns(1:n,:);
    %      vns = tns(n+1:1000);
    vns = zeros(size(vxns,1),1);
    %      vns = vns(si);
    tns = tns(1:n);

    trainingTitle = [arrayfun(@(X) sprintf('xns_%i', X), 1:dim, 'UniformOutput', false), {'tns'}];
    trainingData = zeros(dim+1, size(xns,1));
    trainingData(1:dim,:) = xns';
    trainingData(2,:) = tns;
%    keyboard
    print_datatable(trainingFn, trainingTitle, trainingData);

    validationData = zeros(5, size(vxns,1));

    validationFn = [stemFn, '_rvm_validation.dat'];

    %      keyboard

    Mstart = 0;
    if isa(gp, 'uqrb.rvm.Inference')
      Mstart = gp.basis.M;
      if isa(gp.basis, 'uqrb.rvm.basisFunction.Polynomial')
        disp(num2str(gp.basis.degrees));
        gp.alphaDiscardThreshold = 1e4;
        gp.stop_inferenceConvergenceRate = 1e-3;
      end
    end
    %      [alpGp{i,j}, s2s{i,j}] = gp.inferenceProcedure(tns);
    gp.inferHyperparameters(xns, tns);
    if isa(gp, 'uqrb.rvm.Inference') && isa(gp.basis, 'uqrb.rvm.baseFunction.Polynomial')
      disp(num2str(gp.basis.degrees));
    end
    if isa(gp, 'uqrb.gp.Inference')
      disp(num2str(gp.hyperMu.sigma2));
      disp(num2str(gp.hyperMu.radius));
      % gp.hyperMu.radius = 0.1;
    end
    %gp.alpha = alpGp{i,j};

    % gp.sigma2estimate = 20;
    %      gp.sigma2estimate = s2s{i,j};
    %      if isa(gp.basis, 'uqrb.rvm.basisFunction.Radial')
    %        gp.basis.radius = 0.9;
    %        gp.basis.xms = xns;
    %      end
    p1 = gp.inferPosterior(xns, tns);
    g1 = p1.predict(vxns);

    %      muCurve   = gp.predictForW(vxns,mu);
    %      bfs       = gp.regressionBasis.baseFunctionEval(vxns);
    figure
    subplot(2,1,1);
    [confX, confY, stddev] = g1.confidencePatch(0.95);
    g2=g1.copy();
    g2.Sigma = g1.Sigma + gp.hyperMu.sigma2*speye(size(g2.Sigma));
    [confX2, confY2, stddev2] = g2.confidencePatch(0.95);

    methodname = class(gp);
    if isa(gp, 'uqrb.rvm.Inference')
      methodname = class(gp.basis);
    end
    if dim == 1
      plot(xns, tns, ['*', 'g']);
      hold on
      plot(vxns, g1.mu, ['-', 'k']);
      hold on
      h=patch(vxns(confX2), confY2, ...
        [0.9,0.9,0.9], 'EdgeColor', [0.7,0.7,0.7]);
      alpha(h,0.25);

      h=patch(vxns(confX), confY, ...
        [0.7,0.7,0.7], 'EdgeColor', [0.6,0.6,0.6]);
      alpha(h,0.25);
%      set(gca, 'XScale', 'log');
      title(['k = ', num2str(n), ' abscissa: ', abscissa, ' bf: ', methodname]);
      subplot(2,1,2);
      %      [stdCurve1, stdCurve2] = gp.plot(xns, tns, vxns, [], p1.mu, p1.Sigma, gPlotFunc, hPlotFunc);
      plot(gPlotFunc(xns), hPlotFunc(tns), ['*', 'g']);
      hold on
      plot(gPlotFunc(vxns), hPlotFunc(g1.mu), ['-', 'k']);
      hold on
      h=patch(gPlotFunc(vxns(confX2)), hPlotFunc(confY2), ...
        [0.9,0.9,0.9], 'EdgeColor', [0.7,0.7,0.7]);
      alpha(h,0.25);
      h=patch(gPlotFunc(vxns(confX)), hPlotFunc(confY), ...
        [0.7,0.7,0.7], 'EdgeColor', [0.6,0.6,0.6]);
      alpha(h,0.25);
%      set(gca, 'YScale', 'log');
%      set(gca, 'XScale', 'log');
      title(['k = ', num2str(n), ' abscissa: ', abscissa, ' bf: ', methodname]);
    end

    %      stddev    = sqrt(diag(bfs*(Sigma)*bfs'));
    stdCurve1 = g1.mu + stddev;
    stdCurve2 = g1.mu - stddev;
    stdCurve3 = g1.mu + stddev2;
    stdCurve4 = g1.mu - stddev2;

    validationData(1,:)     = vxns;
    validationData(2,:)     = g1.mu;
    validationData(3,:)     = stdCurve1;
    validationData(4,:)     = stdCurve2;
    validationData(5,:)     = stdCurve3;
    validationData(6,:)     = stdCurve4;
    
    validationData(7,:) = vns;
    print_datatable(validationFn, validationTitle, validationData);
    overviewData{1,counter} = validationFn;
    overviewData{2,counter} = abscissa;
    overviewData{3,counter} = Mstart;
    overviewData{4,counter} = p1.getDimension();
    overviewData{5,counter} = gp.hyperMu.sigma2;
    overviewData{6,counter} = methodname;
    counter = counter + 1;
    %      keyboard;

  end
  print_datatable(overviewFn, overviewTitle, overviewData);

end

end