#!/usr/bin/perl -w
#
# JaeSub Hong, 2003-2005, version 1.5
# Please report any problem or suggestion at jaesub@head.cfa.harvard.edu
# 
# quantile grid generator
# generate grid patterns for spectra modeled by any two parameters
# this is very similar to grid.pl, but instead of evaluating
# values along the grid lines, it maps out 2-d uniformly
# 
# Refer to 
#	J. Hong, E.M. Schlegel & J.E. Grindlay, 2004 ApJ 614 p508 
#		and references therein
# usage
#	map.pl 
#		-model  modelname 
#		-range 	x.x:y.y
# 		-par1 	name:value1,value2,value3,..:points
# 		-par2 	name:value1,value2,value3,..:points
#		-arf	arf_file
#		-rmf	rmf_file
#		[-etcp	nameA:valueA,nameB:valueB,...]
#		[-frac	x.x,y.y,...]
#		[-sherpa dir]
# 
# requires
#	sherpa (http://cxc.harvard.edu/sherpa/)
#	grid.sl (provided along with this file)
#
# examples
#
#	map.pl \
#		-model 	 "xszwabs[A]*xspowerlaw[B]" \
#		-range	 0.3:8.0 \
#		-par1	 A.1:0.1,1,10:100 \
#		-par2	 B.1:1,2,3:100 \
#		-etcp	 A.2:0.0,B.2:1	 	
#		-rmf	 test.rmf \
#		-arf	 test.arf 
#
#	map.pl \
#		-model 	 "xszwabs[A]*xspowerlaw[B]" \
#		-range	 0.3:8.0 \
#		-par1	 A.nH:0.01,0.1,1,10:200 \
#		-par2	 B.PhoIndx:0,1,2,3,4:200 \
#		-etcp	 A.Redshift:0.0,B.norm:1.0 \
#		-rmf	 test.rmf \
#		-arf	 test.arf \
#		> test_powerlaw.rdb							
#
#	map.pl \
#		-model 	 "xszwabs[A]*xsbremss[B]" \
#		-range	 0.3:8.0 \
#		-par1	 A.nH:0.01,0.1,1,10:100 \
#		-par2	 B.kT:0.5,1.0,2.0,4.0,10.0:400 \
#		-etcp	 A.Redshift:0.0,B.norm:1.0 \
#		-rmf	 test.rmf \
#		-arf	 test.arf \
#		> test_bremss.rdb				
#
#
# input
#	Run the program without option for basic options or see below.
#	For the model names and their parameters, refer Sherpa or
#	xspec manual. As for the parameter name, one can simply use numbers
#	(e.g. A.1 or B.1 instead of A.nh or B.PhoIndx), and the parameter 
#	names can be found in sherpa (load a model and type "show").
#
# output
#	The result will be dumped into the screen in RDB format.
#	One can save it by piping.
#		e.g.) grid.pl .... > result.rdb
#	For RDB format, refer 
#		http://hea-www.harvard.edu/MST/simul/software/docs/rdb.html
#	Basically, it is an ASCII text table file where 
#	each field is separated by <tab>.
#	The file will contains 2+n columns where n is the number of "frac"
#	you request.
#	First two columns are the two parameter values, and the rest will
#	be the quantiles (Ex% but not Qx).
#
# limitation
#	potentially many, but the accuracy of grid is limited by
#	the energy scale of rmf and arf files. Need to be re-visited at 
#	some point.
#
#----------------------------------------------------------------------
# get the parameter files from command line
use Getopt::Long;
GetOptions(
	"model=s"	=> \$model,
	"range=s"	=> \$range,
	"frac=s"	=> \$frac,
	"par1=s"	=> \$par1,
	"par2=s"	=> \$par2,
	"etcp=s"	=> \$etcp,
	"arf=s"		=> \$arf,
	"rmf=s"		=> \$rmf,
	"sherpa=s"	=> \$sherpa,
	"gridsl=s"	=> \$gridsl,
	"batch=s"	=> \$batch,
	"keep=s"	=> \$keep,
	"help" 		=> \$help,
	"debug:i" 	=> \$debug);

$debug = 0 unless defined $debug;

if (defined $help
	|| !defined $model
	|| !defined $par1
	|| !defined $par2
	|| !defined $arf
	|| !defined $rmf
	) {
	($usage = <<"	EOFHELP" ) =~ s/^\t\t//gm;
		Find quantiles for given fractions
		usage: $0  [-help] [-debug [X]] 
			-model	modelname 
			-range	X.X:Y.Y
			-par1	name:value1,value2,value3,..:points
			-par2	name:value1,value2,value3,..:points
			-arf	arf_file
			-rmf	rmf_file
			[-frac	X.X,Y.Y,...]
			[-etcp	nameA:valueA,nameB:valueB,...]
			[-sherpa dir]
			[-gridsl dir]
			[-batch  batch_file]
			[-keep]
		options
			-help     print this message
			-debug    set the level of debug
			-model    modelname 
			-range	  energy range
			-arf	  arf file 
			-rmf	  rmf file
			-frac	  fraction for quantile values
			-par1 	  parameter properties for grid lines
			-par2 	  parameter properties for grid lines
			-etcp	  the rest of the parameters
			-sherpa   sherpa directory
			-gridsl   location (directory) of grid.sl 
			-batch    temp batch file name for sherpa input
			-keep     if you want to keep the batch file

	EOFHELP
	print $usage;
	exit;
}

#-----------------------------------------------------------------------
# here we go
# there's got to be a way to do all these in sherpa-slang, but
# my kung-fu is not there yet.

unless (defined $sherpa)
	{ $sherpa = "sherpa"; } 
	else { $sherpa .= "/sherpa"; }

unless (defined $gridsl)
	{ $gridsl = "grid.sl"; } 
	else { $gridsl .= "/grid.sl"; }

$batch = "sherpa.in" unless defined $batch;

#-----------------------------------------------------------------------
# write the sherpa batch input file

$range =~ s/:/,/;
$frac = "0.25,0.5,0.75" unless defined $frac;
$rmf = "ideal.rmf" unless defined $rmf;
$arf = "ideal.arf" unless defined $arf;

($input = <<"EOFINPUT" ) =~ s/^\t//gm;
	paramprompt off 
	INSTRUMENT 1 = RSP[instrumentA] 
	instrumentA.rmf = $rmf
	instrumentA.arf = $arf
	ignore energy :${range}: 
	source 1 = $model
	fakeit 1
	frac=[$frac]
EOFINPUT

if (defined $etcp) {
	@etcp = split/,/, $etcp;
	foreach (@etcp) {
		($name, $value) = split/:/, $_;
		$input .= "$name = $value\n";
	}
}

($p1n,$p1v,$p1p) = split/:/, $par1;
($p2n,$p2v,$p2p) = split/:/, $par2;

@p1v = split/,/, $p1v;
@p2v = split/,/, $p2v;

#($p1_min, $p1_max) = ($p1v[0], $p1v[$#p1v]);
#($p2_min, $p2_max) = ($p2v[0], $p2v[$#p2v]);

# yes, switching (1<->2) is right
#$p1s = ($p1_max-$p1_min)/$p2p;
#$p2s = ($p2_max-$p2_min)/$p1p;

# go for the par1 first

$gid = 0;
foreach (@p1v) {
	$p1 = $_;
	$input .="$p1n = $p1\n";

	foreach $i (1 .. $#p2v) {
		($p2_min, $p2_max) = ($p2v[$i-1], $p2v[$i]);
		# yes, switching (1<->2) is right
		$p2s = ($p2_max-$p2_min)/$p1p * ($#p2v);
		foreach ($p2=$p2_min;$p2<=$p2_max;$p2 +=$p2s) {
			$input .="$p2n = $p2\n"
				."quant = mod_quantile(1,frac)\n"
				."quant_ = pri_quantile(quant)\n"
				."print(\"==gx$gid\t\"+string($p1)+\"\t\"+string($p2)+\"\t\"+quant_)\n";
		}
	}
	$gid++;
}

$gid = 0;
foreach (@p2v) {
	$p2 = $_;
	$input .="$p2n = $p2\n";

	foreach $i (1 .. $#p1v) {
		($p1_min, $p1_max) = ($p1v[$i-1], $p1v[$i]);
		# yes, switching (1<->2) is right
		$p1s = ($p1_max-$p1_min)/$p2p * ($#p1v);
		foreach ($p1=$p1_min;$p1<=$p1_max;$p1 +=$p1s) {
			$input .="$p1n = $p1\n"
				."quant = mod_quantile(1,frac)\n"
				."quant_ = pri_quantile(quant)\n"
				."print(\"==gy$gid\t\"+string($p1)+\"\t\"+string($p2)+\"\t\"+quant_)\n";
		}
	}
	$gid++;
}


open(BAT,"> $batch") or die "can't find $batch\n";
print BAT $input;
close(BAT);


#-----------------------------------------------------------------------
# now go, really.

die "can't find $gridsl\n" unless -e $gridsl;
@output = `$sherpa --slscript $gridsl --batch $batch`;
`rm -rf $batch` unless defined $keep;

#-----------------------------------------------------------------------
# let's try to understand the results

@frac = split/,/,$frac;

$tit ="gridID\t$p1n\t$p2n";
$unit ="S\tN\tN";
for ($i=0;$i<=$#frac;$i++){ 
	$ifrac = sprintf("%d",$frac[$i]*100.);
	$tit .= "\tE$ifrac";
	$unit .= "\tN";
}

$tit =~ s/\./_/g;

print "#\n";
print "# $0 \\ \n";
print "# 	-model \"$model\" \\ \n";
print "# 	-range $range \\ \n";
print "# 	-frac $frac \\ \n";
print "# 	-par1 $par1  \\ \n";
print "# 	-par2 $par2  \\ \n";
print "# 	-etcp $etcp  \\ \n" if defined $etcp;
print "# 	-rmf $rmf  \\ \n";
print "# 	-arf $arf  \\ \n";
print "#\n";
print "$tit\n";
print "$unit\n";
foreach (@output){
	next unless /^==(.*)/;
	print "$1\n";
}

exit;