% file: mcmc_demo.sl    10/30/12 - dd
%

% - - - - -
% Setup Data, Model and Viewing parameters

% Load the data and the model
.source e0102_acis_data
;

.source e0102_model
;

% Use coarser cells and blur the data (and model) to be more comparable to the Einstein HRI:
e2d_set_view(1, [ -30 , 30 , 2.5 ] , [ -30 , 30 , 2.5 ] );
e2d_view[1].nsmoo=3;
% set a higher residual color-coding threshold:
e2d_view[1].chigreen=2.5;
% and set a minimum fractional sigma (10%):
e2d_meta[1].scisigma=0.1;

% The simple scans changed the model parameters to:
s3d_ps.norm1=0.62;
s3d_ps.norm2=0.42;
s3d_ps.norm3=0.028;
%
s3d_ps.az0=-125.0;

% and the Fitting changed these values:
s3d_ps.xc = -0.30;
s3d_ps.yc = 3.30;
s3d_ps.Ri = 13.2; 
s3d_ps.Ro = 16.60;
  s3d_ps.Si = 16.60;
s3d_ps.Ra = 53.0;
s3d_ps.norm1 = 0.62; 

% Compare model-data
e2d_pars2fom;
s3d_list_model;
e2d_norm_optimum(0.001);
%   - Minimum FOM of 814.606 at norm factor = 1.03587
e2d_view_resid;


#iffalse
% Commented out: manually cut/paste the following to command line in
% steps...

% - - - - -
% Demo MCMC
%

% From fitting the xc and yc are very well fixed,
% in contrast there seems less unique values for the other parameters.
% So, do MCMC exploration of those parameters: Ri, Ro, Ra, and norm1.
%
% The 'update_pars is setup in a similar ways as for fitting:
public define e2d_update_pars () 
{ 
   % this line is the same for all cg fits: 
   variable pvals = cg_p0 + cg_smand*cg_dir*e2d_scan_vals[e2d_iscan]; 
   % map the fitting parameter values to the desired model params: 
   s3d_ps.Ri = pvals[0]; 
   s3d_ps.Ro = pvals[1];  s3d_ps.Si = pvals[1];
   s3d_ps.Ra = pvals[2];
   s3d_ps.norm1 = pvals[3];
   % don't allow Ri to be too close to or greater than Ro:
   if ((s3d_ps.Ri+0.6) > s3d_ps.Ro) s3d_ps.Ri = s3d_ps.Ro - 0.6;
};

% The MCMC steps are monitored and saved in the "fomlog" arrays,
% can reset the fomlog values (=0) or leave as to append
% the MCMC points to existing values (from CG fitting or other MCMCs).
% Start anew:
e2d_fomlog_i=0;

% In the MCMC exploration that follows there are two constants
% that the user can adjust to control the MCMC behavior,
% typically their values are of order unity.
% In the code below they are set to 1.01 and 1.02 to distinquish them:
%  "1.01" - this value scales the MCMC STEP SIZE and determines
%           how big the jumps in parameter space are.
%  "1.02" - this value is a factor that scales the FOM value
%           and in this way change the amplitude of the FOM surface.

% define a variable to save the starting point,
% use the CG best fit:
bestfit = [13.2, 16.60, 53.0, 0.62];

% Set the MCMC starting point and step size:
cg_p0 = 1.0*bestfit;
cg_smand = 1.01*[1.0, 1.0, 5.0, 0.05];

% Do MCMC exploration of 200 accepted points
e2d_mcmc_explore(200, 1.02);
% and start there again two more times:
cg_p0 = 1.0*bestfit;
e2d_mcmc_explore(200, 1.02);
%
cg_p0 = 1.0*bestfit;
e2d_mcmc_explore(200, 1.02);

% Reset the parameter update:
e2d_update_pars_NULL;

% Make the plots to show what happened (same as in fitting):
e2d_fomlog_plot("Ri","Ro", 10.0);
e2d_fomlog_plot("Ra","norm1", 10.0);
e2d_fomlog_plot("Ra","Ro", 10.0);
e2d_fomlog_history(["Ri","Ro","Ra","norm1"]);

% Looking at these, ~ best FOM region is around:
bestfit = [14.8, 15.6, 55.0, 0.51];

% So, update the parameters to:
s3d_ps.xc = -0.30;
s3d_ps.yc = 3.30;
s3d_ps.Ri = 14.8; 
s3d_ps.Ro = 15.6;
  s3d_ps.Si = 15.6;
s3d_ps.Ra = 55.0;
s3d_ps.norm1 = 0.51; 

% Compare model-data
e2d_pars2fom;
s3d_list_model;
e2d_norm_optimum(0.001);
%   - Minimum FOM of 827.886 at norm factor = 1.15332
e2d_view_resid;

% Scale the norms by the factor:
s3d_ps.norm1*=1.15;
s3d_ps.norm2*=1.15;
s3d_ps.norm3*=1.15;

% and now the norm factor is ~ 1.00:
e2d_pars2fom;
s3d_list_model;
e2d_norm_optimum(0.001);
%   - Minimum FOM of 811.564 at norm factor = 1.00196
e2d_view_resid;

#endif