Re: is this possible?

On Thu, Sep 04, 2003 at 13:23 -0400, John E. Davis wrote:
> Hi,
>
>    I created a custom kernel for isis that has parameters.  I would
> like to find the best set of parameters that minimize the difference
> between two flux-corrected spectra.  To make this more clear, suppose
> the data set ids are 1 and 2.  The function that I would like to
> minimize is
>
>     define my_chisqr ()
>     {
>         flux_corr (1);
>         flux_corr (2);
>         f1 = get_data_flux(1);
>         f2 = get_data_flux(2);
>
> 	return sum ((f1.value-f2.value)^2);
>     }
>
> Note that since the kernel has parameters, the value returned by the
> get_data_flux function will depend upon those parameters.
>
> What is the best way to accomplish this in isis?

Hi John,

My first attempt to do this triggered an internal isis bug.  
I had not anticipated that flux_corr() might be called from
inside a fit-function -- doing so corrupted an internal data
structure.  I've fixed this in my version so it will be
corrected in the next release.

For the fitting problem, I think you were on the right track.
I've appended a script which should do what you want.  The
script is relatively untested however, because I don't have a
good test case.

Please try it out and let me know if you run into problems.

Needless to say that at the moment it works only with my
version of isis but, since you have access to that, it
shouldn't be a problem.  If no further problems crop up, I'll
release an updated version of isis.

>
> Also, if that wasn't enough, I have a feature request: allow an array
> as the first parameter of set_par, e.g.,
>
>     set_par ([1:3], 0.5);
>
> would set parameters 1,2, and 3 to 0.5.
>

I've added this in my version.  For example, given
   fit_fun("poly(1)");
the following now work:

  set_par ([1:3], 0.5);
  set_par ("poly(1)", [-1,1,1], 0, [-2, 0.5, 0], 10);

The 2nd one frees all parameters and sets 10 as an upper limit
for all of them.  The parameter values and lower limits are set
individually.

Thanks,

John

--
John C. Houck                       MIT Center for Space Research
NE80-6005: 617-253-3849                   77 Massachusetts Avenue
42:21:55.105N, 71:05:28.122W             Cambridge, MA 02139-4307



% load data and assign responses
variable dir, phas, arfs, rmfs;
dir  = "/nfs/cxc/h1/houck/src/isis-examples/data";
phas = dir + "/acisf01318N003_pha2.fits";
arfs = dir + ["/acisf01318_000N001MEG_-1_garf.fits",
	      "/acisf01318_000N001MEG_1_garf.fits"];
rmfs = dir + ["/acismeg1D1999-07-22rmfN0002.fits.gz",
	      "/acismeg1D1999-07-22rmfN0002.fits.gz"];

phas = load_data (phas, [9,10]);
arfs = array_map (Integer_Type, &load_arf, arfs);
rmfs = array_map (Integer_Type, &load_rmf, rmfs);

assign_rmf (rmfs, phas);
assign_arf (arfs, phas);

% matching wavelength grids simplifies the problem:
match_dataset_grids (phas);

% This fit-statistic minimizes the difference
% between the flux-corrected versions of the
% two datasets

static define stat_function (y, fx, w)
{
   variable f1, f2, v;
   flux_corr (phas);
   f1 = get_data_flux(1);
   f2 = get_data_flux(2);
   v = f1.value - f2.value;
   return (y-fx, sum (v^2));
}
static define stat_report (stat, npts, nvpars)
{
   variable s = sprintf ("  stat = %0.4g\n", stat);
   return s;
}
add_slang_statistic ("stat", &stat_function, &stat_report);
set_fit_statistic ("stat");

% The fit-function is irrelevant here because
% the fit-statistic is determined by combining
% the assigned kernel and the user-defined fit-statistic.

fit_fun ("bin_width(1)");

% Assign the fit-kernels
%set_kernel (phas, "kernel");

() = fit_counts ();

% Note that, because of the unusual fit-function definition,
% rplot_counts() and rplot_flux() won't produce useful plots.
% For what its worth, here's a cute trick to plot the difference
% in the flux-corrected datasets.

variable weights, gid, c;

weights = [1.0, -1.0];
gid = combine_datasets (phas, weights);
c = get_combined (gid, &get_data_flux);
uncombine_datasets (gid);

hplot(c);