#!/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
# 
# Refer to 
#	J. Hong, E.M. Schlegel & J.E. Grindlay, 2004 ApJ 614 p508 
#		and references therein
# usage
#	grid.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
#
#	grid.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 
#
#	grid.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							
#
#	grid.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,
	"relp=s"	=> \$relp,
	"arf=s"		=> \$arf,
	"rmf=s"		=> \$rmf,
	"full"		=> \$full,
	"sherpa=s"	=> \$sherpa,
	"gridsl=s"	=> \$gridsl,
	"batch=s"	=> \$batch,
	"keep=s"	=> \$keep,
	"raw"		=> \$raw,
	"time=s"	=> \$time,
	"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,...]
			[-relp	line_model_par:cont_model_name:low:high[:reln]]
			[-sherpa dir]
			[-gridsl dir]
			[-batch  batch_file]
			[-keep]
		options
			-help     print this message
			-debug    set the level of debug
			-model    modelname 
			-abs 	  abs_model_name
			-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
			-relp 	  set the parameter relative to other part of the model
			-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;


@model = split/;/, $model;
$model[$#model] = "source 1 = ". $model[$#model];
$rmodel = $model[$#model];
$model_ = join("\n",@model);

$rerun = "";
if (defined $relp) {
	($par,$cmod,$low,$high,$reln)  = split/:/, $relp;
	$reln = 1. unless defined $reln;
	($rerun .= << "	EOF") =~ s/^\t\t//gm;
		source 1 = $cmod
		fakeit 1
		K=get_pflux(1,[$low,$high])
		()=set_par("$par","value",K.value*$reln)
		$rmodel
	EOF
} 


if (defined $raw) { 
	$x_axe = "channels"; 
	$time = 1000000 unless defined $time;
	$timeset = "fakeit time = $time";
	$suffix = "_raw";
	$rerun .= "$timeset\nfakeit 1\n";
}
else { 
	$x_axe = "energy"; 
	$timeset = "";
	$suffix = "";
	$rerun .= "fakeit 1\n";
}

($input = <<"EOFINPUT" ) =~ s/^\t//gm;
	paramprompt off 
	INSTRUMENT 1 = RSP[instrumentA] 
	instrumentA.rmf = $rmf
	instrumentA.arf = $arf
	ignore $x_axe :${range}: 
	frac=[$frac]
	$model_
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

if (defined $full) {
	$gid = 0;
	foreach $j (1 .. $#p1v) {
		($p1_min, $p1_max) = ($p1v[$j-1], $p1v[$j]);
		# 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";

			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$rerun"
						."quant = mod_quantile${suffix}(1,frac)\n"
						."quant_ = pri_quantile(quant)\n"
						."print(\"==gx$gid\t\"+string($p1)+\"\t\"+string($p2)+\"\t\"+quant_)\n";
				}
			}
			$gid++;
		}
	}
} else {
	$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$rerun"
					."quant = mod_quantile${suffix}(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$rerun"
					."quant = mod_quantile${suffix}(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.

exit if $debug == 5;

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;