%
% file: a1068_cg_mcmc.sl
% 
% Commands to do CG fitting and MCMC exploration
% with the Abel 1068 ACIS data set.
% 
% 10/2/08 dd
% 


% - - - Read in the data:
.source a1068_data

% Includes both X-Y and RADIUS-ENERGY data sets.
% We'll fit just the X,Y Image, so
% ignore the RADIUS-ENERGY data:
e2d_ignore(2);

% Look at the X,Y data:
e2d_view_data(1); 
%     bins w/events have 50%, 90%, and max numbers of:  4,  22,  626


% - - - Setup the model:
.source a1068_model_simple
  
% Set some rough parameter values:
% first look at the available parameters and their values:
print(s3d_ps);

% Set these to values we might have guessed:
s3d_ps.rbeta = 15.0;
s3d_ps.ratio = 1.2;
s3d_ps.ang = -30.0;
s3d_ps.beta = 0.50;

% and look at a projection view of the components:
s3d_update;
s3d_proj_comps;



% - - - Setup the binning and science sigma

% Evaluate the model
e2d_pars2fom;
%      o Data-Model 1,  chi^2 = 8548.87  from 4717 cells;   chi^2/cell = 1.81235;  sum(D)/sum(M) = 1.03212  +/- 0.00446559
%      o Total FOM = 8548.87

% look at the residuals in 2D:
e2d_view_resid(1);

% change to a coarser binning for fitting:
e2d_view_data (1, [ -90 , 90 , 7 ] , [ -90 , 90 , 7 ] );
%     bins w/events have 50%, 90%, and max numbers of:  29,  159,  3126
e2d_list_fom;
%      o Data-Model 1,  chi^2 = 6912.19  from 676 cells;   chi^2/cell = 10.2251;  sum(D)/sum(M) = 1.03212  +/- 0.00446559
%      o Total FOM = 6912.19
% Look at the Data-Model residuals in 2D
e2d_view_resid(1);

% Check the "chi-percentage" plot
e2d_chipc_plot;
% Change to a 10% science-sigma limit:
e2d_meta[1].scisigma = 10.0e-2;
e2d_list_fom;
%      o Data-Model 1,  chi^2 = 2977.89  from 676 cells;   chi^2/cell = 4.40517;  sum(D)/sum(M) = 1.03212  +/- 0.00446559
%      o Total FOM = 2977.89
e2d_chipc_plot;

% Check that the norm is roughly correct and adjust it:
norm_adj = e2d_norm_optimum(NULL);
% Scale the norm (by ~ 0.80) for better fit:
s3d_comp[1].norm = s3d_comp[1].norm * norm_adj;

% update the model and look at our current residuals:
e2d_pars2fom;
%      o Data-Model 1,  chi^2 = 1837.93  from 676 cells;   chi^2/cell = 2.71884;  sum(D)/sum(M) = 1.29531  +/- 0.00560428
%      o Total FOM = 1837.93
e2d_view_resid(1);



% - - - Scan the parameters
% To get some idea of the scale over which each parameter changes,
% do 1-D scans of each parameter with the others fixed.

% Set the parameters close to where the scan has its
% minimum or low point and make a note of the 
% "should make a noticable difference" scale for each parameter.

% Each scan requires the "e2d_update_pars" routine be
% redefined for the appropriate scan, and then
% the e2d_scan_par routine can be run.

% reset the FOM-log values:
% (The scanned parameters and chi-squareds are saved for
%  viewing - this resets the save list.)
e2d_fomlog_i = 0;


% rbeta:
public define e2d_update_pars () 
{         
   s3d_ps.rbeta =  e2d_scan_vals[e2d_iscan];
}                 
e2d_scan_par([8.0:18.0:#9],3);
%  Set it to ~ min and note that smand ~ 3.0
s3d_ps.rbeta=10.5;


% ratio:
public define e2d_update_pars () 
{         
   s3d_ps.ratio = e2d_scan_vals[e2d_iscan]; 
}                 
e2d_scan_par([1.0:1.8:#9],3);
%  Set it to ~ min and note that smand ~ 0.2
s3d_ps.ratio = 1.35;


% angle:   * Note: nominal value is negative ! *
public define e2d_update_pars () 
{         
   s3d_ps.ang =  e2d_scan_vals[e2d_iscan];
}                 
e2d_scan_par([-60.0:0.0:#9],3);
%  Set it to ~ min and note that smand ~ 15.0
s3d_ps.ang = -40.0;

 
% beta:
public define e2d_update_pars () 
{         
   s3d_ps.beta =  e2d_scan_vals[e2d_iscan];
}                 
e2d_scan_par([0.40:0.60:#9],3);
%  Set it to ~ min and note that smand ~ 0.04
s3d_ps.beta = 0.50;


% Reset/cancel the definition of e2d_update_pars :
% (  It is good to be sure this is done habitually 
%    after scanning or fitting to cancel the
%    e2d_update_pars definition.   )
e2d_update_pars_NULL;

% Look at our new scan-fit parameter values
% and evaluate the model etc as before:
s3d_ps.rbeta=10.5;
s3d_ps.ratio = 1.35;
s3d_ps.ang = -40.0;
s3d_ps.beta = 0.50;

e2d_pars2fom;
norm_adj = e2d_norm_optimum(NULL);
% Scale the norm (by ~ 1.09) for better fit:
s3d_comp[1].norm = s3d_comp[1].norm * norm_adj;
e2d_pars2fom;
% That's an improvement:
%      o Data-Model 1,  chi^2 = 1301.51  from 676 cells;   chi^2/cell = 1.92532;  sum(D)/sum(M) = 1.1782  +/- 0.00509763
%      o Total FOM = 1301.51
e2d_view_resid(1);

% From these residuals it's clear that the center is
% unlikely to be well-fit with a simple beta model only, 
% e.g., we should add a central point source, etc.
% However, let's go with this model and fit it and
% MCMC it ...



% - - - Fit the paramters

% Each fit requires the "e2d_update_pars" routine be
% redefined for the appropriate fit, and then
% the e2d_cg_fit routine can be run.
% (this looks similar to the scan case because the
%  cg fit routine calls the scan routine.)

% Set up to fit 4 parameters:
% Assign them to the cg_p0 array in the order that
% we scanned them above for convenience:

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.rbeta = pvals[0]; 
   s3d_ps.ratio = pvals[1]; 
   s3d_ps.ang = pvals[2]; 
   s3d_ps.beta = pvals[3]; 
   % can also clamp parameters, etc. here:
   if (s3d_ps.rbeta < 3.0) s3d_ps.rbeta = 3.0;
} 

% Of course it would be better to start the fitting at the scan
% values we found above, but let's make the cg routine
% show us how it works.  (We could have skipped the scanning
% if we had some idea of the "smand" scales to use.)

% Start at the same very rough guess set of values as
% we used above before the scanning.
cg_p0 = [15.0, 1.20, -30.0, 0.50];

% and set the "should make a noticable difference" values
% for the parameters according to what the scans suggested:
cg_smand = [3.0, 0.2, 15.0, 0.04];

% === and do the fit:  (the "1" is a dummy argument)
e2d_cg_fit(1);
% ===


% Plot the FOM-log values
% Show the delta chi-squared that corresponds to
% Data-Model agreement =  sqrt(2*N_dof)
variable delta_chi = sqrt(2.0*e2d_data[1].nfom);
delta_chi;
e2d_fomlog_plot("rbeta","beta", delta_chi);
% output is in fomlog_plot.ps: rename the plot
()=system("cp fomlog_plot.ps fomlog_rbeta_beta_plot.ps");
e2d_fomlog_plot("ang","ratio", delta_chi);
()=system("cp fomlog_plot.ps fomlog_ang_ratio_plot.ps");

% Look at the data and model and residual
e2d_update_pars_NULL;
s3d_ps.rbeta = 7.832;
s3d_ps.ratio = 1.419;
s3d_ps.ang = -41.73;
s3d_ps.beta = 0.4897;
e2d_pars2fom;
norm_adj = e2d_norm_optimum(NULL);
s3d_comp[1].norm = s3d_comp[1].norm * norm_adj;
e2d_pars2fom;
% 
%      o Data-Model 1,  chi^2 = 1257.12  from 676 cells;   chi^2/cell = 1.85965;  sum(D)/sum(M) = 1.12413  +/- 0.0048637
%      o Total FOM = 1257.12
e2d_view_resid(1);



% - - - Do MCMC exploration

% Uses the same e2d_update_pars routine and the 
% cg_p0 and cg_smand parameters already defined above -
% this makes it easy to do MCMC and CG fitting together/alternately.

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.rbeta = pvals[0]; 
   s3d_ps.ratio = pvals[1]; 
   s3d_ps.ang = pvals[2]; 
   s3d_ps.beta = pvals[3]; 
   % can also clamp parameters, etc. here:
   if (s3d_ps.rbeta < 3.0) s3d_ps.rbeta = 3.0;
} 

% Start the MCMC at the same point as the cg_fit did:
cg_p0 = [15.0, 1.20, -30.0, 0.50];
cg_smand = [3.0, 0.2, 15.0, 0.04];

% Put that same delta-chi here to help scale delta-chisquared
% to a probability, like simulated annealing with "T" = sqrt(2*N_dof)
e2d_mcmc_explore(100, 1.0/delta_chi);

% start at the fit min:
cg_p0 = [7.832, 1.419, -41.73, 0.4897];
e2d_mcmc_explore(100, 1.0/delta_chi);

% start between them:  :-)
cg_p0 = 0.5*( [7.832, 1.419, -41.73, 0.4897] + [15.0, 1.20, -30.0, 0.50] );
% and "lower the temp" by half to emphasize the minimums:
e2d_mcmc_explore(100, 2.0/delta_chi);

% start at the fit min: and T/3
cg_p0 = [7.832, 1.419, -41.73, 0.4897];
e2d_mcmc_explore(100, 3.0/delta_chi);


% Make FOM-log plots with the cg fit and MCMC points together:
e2d_fomlog_plot("rbeta","beta", delta_chi);
()=system("cp fomlog_plot.ps fomlog_rbeta_beta_wMCMC_plot.ps");
e2d_fomlog_plot("ang","ratio",delta_chi);
()=system("cp fomlog_plot.ps fomlog_ang_ratio_wMCMC_plot.ps");


% Make a plot of the parameter values vs iteration:
PSwin=open_plot("fomlog_ps_vs_iter.ps/VCPS",1,4);
title(""); xlabel("Iteration");
variable str_tag = "rbeta";
ylabel(str_tag); yrange();
plot( [0:e2d_fomlog_i-1], array_struct_field(e2d_fomlog_ps[[0:e2d_fomlog_i-1]], str_tag));
variable str_tag = "ratio";
ylabel(str_tag); yrange();
plot( [0:e2d_fomlog_i-1], array_struct_field(e2d_fomlog_ps[[0:e2d_fomlog_i-1]], str_tag));
variable str_tag = "ang";
ylabel(str_tag); yrange();
plot( [0:e2d_fomlog_i-1], array_struct_field(e2d_fomlog_ps[[0:e2d_fomlog_i-1]], str_tag));
variable str_tag = "beta";
ylabel(str_tag); yrange();
plot( [0:e2d_fomlog_i-1], array_struct_field(e2d_fomlog_ps[[0:e2d_fomlog_i-1]], str_tag));
close_plot(PSwin);


% restart MCMC at 196:
e2d_fomlog_i=196;

% Do MCMC with rbeta fixed:
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.ratio = pvals[0]; 
   s3d_ps.ang = pvals[1]; 
   s3d_ps.beta = pvals[2]; 
} 

% fix it at the cg fit value:
s3d_ps.rbeta=7.832;
% the other three parameters:
cg_p0 = [1.20, -30.0, 0.50];
cg_smand = [0.2, 15.0, 0.04];
variable delta_chi = sqrt(2.0*e2d_data[1].nfom);
e2d_mcmc_explore(400, 1.0/delta_chi);

e2d_fomlog_plot("rbeta","beta", delta_chi);
()=system("cp fomlog_plot.ps fomlog_rbeta_beta_wMCMCrbfixed_plot.ps");
e2d_fomlog_plot("ang","ratio",delta_chi);
()=system("cp fomlog_plot.ps fomlog_ang_ratio_wMCMCrbfixed_plot.ps");
e2d_fomlog_plot("ratio","beta", delta_chi);
()=system("cp fomlog_plot.ps fomlog_ratio_beta_wMCMCrbfixed_plot.ps");
e2d_fomlog_plot("ang","beta", delta_chi);
()=system("cp fomlog_plot.ps fomlog_ang_beta_wMCMCrbfixed_plot.ps");

% Remake the ps vs iter plots for this case too.