Rewrite 'calc_entropy.py' from 'analyze_entropy_profile.py'

* XXX: 'calc_entropy.py' needs test
* XXX: check the uncertainty/error estimation/calculation
* Remove obsolete 'fit_{d,}beta_entropy.sh' and 'calc_all_entropy.sh'

5 files changed, 104 insertions, 392 deletions

diff --git a/mass_profile/analyze_entropy_profile.py b/mass_profile/analyze_entropy_profile.py
deleted file mode 100755
index c123681..0000000
--- a/mass_profile/analyze_entropy_profile.py
+++ /dev/null
@@ -1,54 +0,0 @@
-#!/usr/bin/env python
-
-import sys
-import numpy
-import scipy.interpolate
-
-center_entropy_file=open('entropy_center.qdp')
-entropy_file=open('summary_entropy.qdp')
-confidence_level=.68
-rout=float(sys.argv[1])
-
-center_s=0
-for i in center_entropy_file:
-    r,s=i.split()
-    r=float(r)
-    s=float(s)
-    if r>rout:
-        center_s=s
-        break
-
-new_data=True
-
-
-s_list=[]
-for i in entropy_file:
-    if i[0]=='n':
-        new_data=True
-        continue
-    if new_data==False:
-        continue
-    r,s=i.split()
-    r=float(r)
-    s=float(s)
-    if r>rout:
-        new_data=False
-        s_list.append(s)
-
-s_idx=-1
-
-s_list.sort()
-for i in range(len(s_list)-1):
-    if (center_s-s_list[i])*(center_s-s_list[i+1])<=0:
-        m_idx=i
-        break
-
-
-slidx=int(s_idx*(1-confidence_level))
-suidx=s_idx-1+int((len(s_list)-s_idx)*confidence_level)
-
-
-serr1=s_list[slidx]-center_s
-serr2=s_list[suidx]-center_s
-
-print("S=\t%e\t %e/+%e keV cm^2 (1 sigma)"%(center_s,serr1,serr2))
diff --git a/mass_profile/calc_all_entropy.sh b/mass_profile/calc_all_entropy.sh
deleted file mode 100755
index 5506ab8..0000000
--- a/mass_profile/calc_all_entropy.sh
+++ /dev/null
@@ -1,39 +0,0 @@
-#!/usr/bin/env bash
-
-if [ $# -eq 1 ]
-then
-    :
-else
-    echo "Usage: <global.cfg list file>"
-    exit
-fi
-
-bdir=`pwd`
-
-for i in `cat $1`
-do
-    dname=`dirname $i`
-    #dname=`dirname $dname`
-    echo $dname
-    cd $dname ||continue
-#    mkdir -p profile_entropy
-#    cp -rf profile/* profile_entropy/
-#    cd profile/
-    
-    [ -e result ] && r200file=result
-    [ -e result ] || r200file=result_checked
-    
-    r200=`grep r200 $r200file|awk -F = '{print $2}' |awk -F '-' '{print $1}'`
-    rout=`python -c "print($r200*.1)"`
-    echo $rout
-    if  grep n01 source.cfg  >/dev/null
-    then
-	echo "dbeta"
-	fit_dbeta_entropy.sh global.cfg $rout
-    else
-	echo "beta"
-	fit_beta_entropy.sh global.cfg $rout
-    fi
-    
-    cd $bdir
-done
diff --git a/mass_profile/calc_entropy.py b/mass_profile/calc_entropy.py
new file mode 100755
index 0000000..736fe7e
--- /dev/null
+++ b/mass_profile/calc_entropy.py
@@ -0,0 +1,104 @@
+#!/usr/bin/env python3
+#
+# Calculate the entropy within the specified radius.
+#
+# Junhua GU
+# Weitian LI
+# 2016-06-07
+#
+
+import argparse
+import re
+from itertools import groupby
+import numpy as np
+
+
+def isplit(iterable, splitters):
+    """
+    Credit: https://stackoverflow.com/a/4322780/4856091
+    """
+    return [list(g) for k, g in groupby(iterable,
+                                        lambda x:x in splitters) if not k]
+
+
+def get_entropy(data, r):
+    """
+    Get the entropy *at* the specified radius.
+
+    XXX: whether to interpolate first?
+    """
+    radius = data[:, 0]
+    entropy = data[:, 1]
+    s = np.min(entropy[radius > r])
+    return s
+
+
+def read_merged_qdp(infile):
+    """
+    Read merged QDP with multiple group of data separated by "no no no".
+    """
+    lines = map(lambda line: re.sub(r"^\s*no\s+no\s+no.*$", "X",
+                                    line.strip(), flags=re.I),
+                open(infile).readlines())
+    lines = isplit(lines, ("X",))
+    data_groups = []
+    for block in lines:
+        data = [list(map(float, l.split())) for l in block]
+        data.append(np.row_stack(data))
+    return data_groups
+
+
+def calc_error(center_value, mc_values, ci=0.683):
+    """
+    Calculate the uncertainties/errors.
+    """
+    data = np.concatenate([[center_value], mc_values])
+    median, q_lower, q_upper = np.percentile(data, q=(50, 50-50*ci, 50+50*ci))
+    mean = np.mean(data)
+    std = np.std(data)
+    return {
+        "mean":    mean,
+        "std":     std,
+        "median":  median,
+        "q_lower": q_lower,
+        "q_upper": q_upper,
+    }
+
+
+def main():
+    parser = argparse.ArgumentParser(
+            description="Calculate the entropy within the given radius")
+    parser.add_argument("-C", "--confidence-level", dest="ci",
+                        type=float, default=0.683,
+                        help="confidence level to estimate the errors")
+    parser.add_argument("center_data",
+                        help="calculate central entropy profile " +
+                             "(e.g., entropy_center.qdp)")
+    parser.add_argument("mc_data",
+                        help="Merged QDP file of all the Monte Carlo " +
+                             "simulated entropy profiles " +
+                             "(e.g., summary_entropy.qdp)")
+    parser.add_argument("rout", type=float, help="outer radius (kpc)")
+    args = parser.parse_args()
+
+    center_data = np.loadtxt(args.center_data)
+    center_s = get_entropy(center_data, r=args.rout)
+
+    data_groups = read_merged_qdp(args.mc_data)
+    entropy_list = []
+    for dg in data_groups:
+        s = get_entropy(dg, r=args.rout)
+        entropy_list.append(s)
+    results = calc_error(center_s, entropy_list, ci=args.ci)
+    s_err_lower = results["q_lower"] - center_s
+    s_err_upper = results["q_upper"] - center_s
+
+    print("entropy= %e %+e/%+e keV cm^2 (ci=%.1f%%)" %
+          (center_s, s_err_lower, s_err_upper, args.ci * 100))
+    print("entropy(mean)= %e" % results["mean"])
+    print("entropy(median)= %e" % results["median"])
+    print("entropy(std)= %e" % results["std"])
+
+
+if __name__ == "__main__":
+    main()
diff --git a/mass_profile/fit_beta_entropy.sh b/mass_profile/fit_beta_entropy.sh
deleted file mode 100755
index f52e4a1..0000000
--- a/mass_profile/fit_beta_entropy.sh
+++ /dev/null
@@ -1,150 +0,0 @@
-#!/bin/bash
-
-echo $#
-if [ $# -eq 2 ]
-then
-    :
-else
-    echo "Usage:$0 <cfg file> <radius in kpc>"
-    exit
-fi
-export PGPLOT_FONT=`locate grfont.dat|head -1`
-
-mkdir -p entropy_files
-
-cfg_file=$1
-base_path=`dirname $0`
-echo $base_path
-#initialize profile type name
-t_profile_type=`grep t_profile $cfg_file|awk '{print $2}'`
-#initialize data file name
-t_data_file=`grep t_data_file $cfg_file|awk '{print $2}'`
-#initialize sbp config file
-sbp_cfg=`grep sbp_cfg $cfg_file|awk '{print $2}'`
-#initialize the temperature profile file
-T_file=`grep '^T_file' $sbp_cfg|awk '{print $2}'`
-#echo $t_profile_type
-cm_per_pixel=`grep '^cm_per_pixel' $sbp_cfg|awk '{print $2}'`
-
-#determine which temperature profile to be used, and fit the T profile
-if [ $t_profile_type == zyy ]
-then
-    t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
-    $base_path/fit_zyy_model $t_data_file $t_param_file $cm_per_pixel
-    mv -f zyy_dump.qdp ${T_file}
-elif [ $t_profile_type == m0603246 ]
-then
-    $base_path/fit_m0603246 $t_data_file $cm_per_pixel
-    mv -f m0603246_dump.qdp ${T_file}
-elif [ $t_profile_type == wang2012 ]
-then
-    t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
-    $base_path/fit_wang2012_model $t_data_file $t_param_file $cm_per_pixel
-    mv -f wang2012_dump.qdp ${T_file}
-elif [ $t_profile_type == allen ]
-then
-    $base_path/fit_allen_model $t_data_file $cm_per_pixel
-    mv -f allen_dump.qdp ${T_file}
-elif [ $t_profile_type == zzl ]
-then
-    t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
-    $base_path/fit_zzl_model $t_data_file $t_param_file
-    mv -f zzl_dump.qdp ${T_file}
-else
-    echo temperature profile name invalid!
-    exit
-fi
-
-cfunc_file=`grep '^cfunc_file' ${sbp_cfg} |awk '{print $2}'`
-z=`grep  '^z' ${sbp_cfg}|awk '{print $2}'`
-abund=`grep '^abund' ${cfg_file} |awk '{print $2}'`
-nh=`grep '^nh' ${cfg_file} |awk '{print $2}'`
-$base_path/coolfunc_calc.sh ${T_file} $abund $nh $z $cfunc_file
-mv flux_cnt_ratio.txt flux_cnt_ratio_center_entropy.txt
-#fit sbp
-$base_path/fit_beta_sbp $sbp_cfg
-echo $cfunc_file
-#exit
-
-#store central valu
-mv sbp_fit.qdp sbp_fit_center_entropy.qdp
-mv mass_int.qdp mass_int_center_entropy.qdp
-mv overdensity.qdp overdensity_center_entropy.qdp
-mv gas_mass_int.qdp gas_mass_int_center_entropy.qdp
-mv entropy.qdp entropy_center.qdp
-sbp_data_file=`grep sbp_file $sbp_cfg|awk '{print $2}'`
-radius_sbp_file=`grep radius_sbp_file ${cfg_file}|awk '{print $2}'`
-
-if [ x"$radius_sbp_file" == x ]
-then
-    echo "Error, must have radius_sbp_file assigned, this file should be a 4-column file, which contains the radius, radius err, sbp, and sbp err"
-    exit
-fi
-
-cat ${radius_sbp_file} | sed 's/#.*$//' | grep -Ev '^\s*$' > .tmp.txt
-mv .tmp.txt ${radius_sbp_file}
-
-rm -f summary_entropy.qdp
-#100 times of Monte-carlo simulation to determine error
-#just repeat above steps
-for i in `seq 1 100`
-do
-    echo $t_data_file
-    $base_path/shuffle_T.py $t_data_file temp_shuffled_t.dat
-    $base_path/shuffle_sbp.py $sbp_data_file temp_shuffled_sbp.dat
-    
-#exit
-    if [ $t_profile_type == zyy ]
-    then
-	t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
-	$base_path/fit_zyy_model temp_shuffled_t.dat $t_param_file $cm_per_pixel
-	mv -f zyy_dump.qdp ${T_file}
-    elif [ $t_profile_type == m0603246 ]
-    then
-	$base_path/fit_m0603246 temp_shuffled_t.dat $cm_per_pixel
-	mv -f m0603246_dump.qdp ${T_file}
-    elif [ $t_profile_type == wang2012 ]
-    then
-	t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
-	$base_path/fit_wang2012_model temp_shuffled_t.dat $t_param_file $cm_per_pixel
-    mv -f wang2012_dump.qdp ${T_file}
-    elif [ $t_profile_type == allen ]
-    then
-	$base_path/fit_allen_model temp_shuffled_t.dat $cm_per_pixel
-	mv -f allen_dump.qdp ${T_file}
-    elif [ $t_profile_type == zzl ]
-    then
-	t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
-    $base_path/fit_zzl_model temp_shuffled_t.dat $t_param_file
-    mv -f zzl_dump.qdp ${T_file}
-    else
-	echo temperature profile name invalid!
-	exit
-    fi
-    
-#exit
-    echo >temp_sbp.cfg
-    
-    cat $sbp_cfg|while read l
-do
-    if echo $l|grep sbp_file >/dev/null
-    then
-	echo sbp_file temp_shuffled_sbp.dat >>temp_sbp.cfg
-    elif echo $l|grep T_file >/dev/null
-    then
-	echo T_file ${T_file} >>temp_sbp.cfg
-    else
-	echo $l >>temp_sbp.cfg
-    fi
-    
-done
-
-#$base_path/coolfunc_calc.sh ${T_file} $abund $nh $z $cfunc_file
-
-$base_path/fit_beta_sbp temp_sbp.cfg
-
-cat entropy.qdp >>summary_entropy.qdp
-echo no no no >>summary_entropy.qdp
-done
-
-$base_path/analyze_entropy_profile.py $2 |tee entropy_result.txt
diff --git a/mass_profile/fit_dbeta_entropy.sh b/mass_profile/fit_dbeta_entropy.sh
deleted file mode 100755
index 3f1bbda..0000000
--- a/mass_profile/fit_dbeta_entropy.sh
+++ /dev/null
@@ -1,149 +0,0 @@
-#!/bin/bash
-
-echo $#
-if [ $# -eq 2 ]
-then
-    :
-else
-    echo "Usage:$0 <cfg file> <radius in kpc>"
-    exit
-fi
-export PGPLOT_FONT=`locate grfont.dat|head -1`
-
-cfg_file=$1
-base_path=`dirname $0`
-echo $base_path
-#initialize profile type name
-t_profile_type=`grep t_profile $cfg_file|awk '{print $2}'`
-#initialize data file name
-t_data_file=`grep t_data_file $cfg_file|awk '{print $2}'`
-#initialize sbp config file
-sbp_cfg=`grep sbp_cfg $cfg_file|awk '{print $2}'`
-#initialize the temperature profile file
-T_file=`grep '^T_file' $sbp_cfg|awk '{print $2}'`
-#echo $t_profile_type
-cm_per_pixel=`grep '^cm_per_pixel' $sbp_cfg|awk '{print $2}'`
-#determine which temperature profile to be used, and fit the T profile
-if [ $t_profile_type == zyy ]
-then
-    t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
-    $base_path/fit_zyy_model $t_data_file $t_param_file $cm_per_pixel
-    mv -f zyy_dump.qdp ${T_file}
-elif [ $t_profile_type == m0603246 ]
-then
-    $base_path/fit_m0603246 $t_data_file $cm_per_pixel
-    mv -f m0603246_dump.qdp ${T_file}
-elif [ $t_profile_type == wang2012 ]
-then
-    t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
-    $base_path/fit_wang2012_model $t_data_file $t_param_file $cm_per_pixel
-    mv -f wang2012_dump.qdp ${T_file}
-elif [ $t_profile_type == allen ]
-then
-    $base_path/fit_allen_model $t_data_file $cm_per_pixel
-    mv -f allen_dump.qdp ${T_file}
-elif [ $t_profile_type == zzl ]
-then
-    t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
-    $base_path/fit_zzl_model $t_data_file $t_param_file
-    mv -f zzl_dump.qdp ${T_file}
-else
-    echo temperature profile name invalid!
-    exit
-fi
-
-cfunc_file=`grep '^cfunc_file' ${sbp_cfg} |awk '{print $2}'`
-z=`grep  '^z' ${sbp_cfg}|awk '{print $2}'`
-abund=`grep '^abund' ${cfg_file} |awk '{print $2}'`
-nh=`grep '^nh' ${cfg_file} |awk '{print $2}'`
-$base_path/coolfunc_calc.sh ${T_file} $abund $nh $z $cfunc_file
-mv flux_cnt_ratio.txt flux_cnt_ratio_center_entropy.txt
-#fit sbp
-$base_path/fit_dbeta_sbp $sbp_cfg
-echo $cfunc_file
-#exit
-
-#store central valu
-mv sbp_fit.qdp sbp_fit_center_entropy.qdp
-mv mass_int.qdp mass_int_center_entropy.qdp
-mv overdensity.qdp overdensity_center_entropy.qdp
-mv gas_mass_int.qdp gas_mass_int_center_entropy.qdp
-mv entropy.qdp entropy_center.qdp
-sbp_data_file=`grep sbp_file $sbp_cfg|awk '{print $2}'`
-radius_sbp_file=`grep radius_sbp_file ${cfg_file}|awk '{print $2}'`
-
-if [ x"$radius_sbp_file" == x ]
-then
-    echo "Error, must have radius_sbp_file assigned, this file should be a 4-column file, which contains the radius, radius err, sbp, and sbp err"
-    exit
-fi
-
-cat ${radius_sbp_file} | sed 's/#.*$//' | grep -Ev '^\s*$' > .tmp.txt
-mv .tmp.txt ${radius_sbp_file}
-
-rm -f summary_entropy.qdp
-
-#100 times of Monte-carlo simulation to determine error
-#just repeat above steps
-for i in `seq 1 100`
-do
-    echo $t_data_file
-    $base_path/shuffle_T.py $t_data_file temp_shuffled_t.dat
-    $base_path/shuffle_sbp.py $sbp_data_file temp_shuffled_sbp.dat
-
-#exit
-    if [ $t_profile_type == zyy ]
-    then
-	t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
-	$base_path/fit_zyy_model temp_shuffled_t.dat $t_param_file $cm_per_pixel
-	mv -f zyy_dump.qdp ${T_file}
-    elif [ $t_profile_type == m0603246 ]
-    then
-	$base_path/fit_m0603246 temp_shuffled_t.dat $cm_per_pixel
-	mv -f m0603246_dump.qdp ${T_file}
-    elif [ $t_profile_type == wang2012 ]
-    then
-	t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
-	$base_path/fit_wang2012_model temp_shuffled_t.dat $t_param_file $cm_per_pixel
-	mv -f wang2012_dump.qdp ${T_file}
-    elif [ $t_profile_type == allen ]
-    then
-	$base_path/fit_allen_model temp_shuffled_t.dat $cm_per_pixel
-	mv -f allen_dump.qdp ${T_file}
-    elif [ $t_profile_type == zzl ]
-    then
-	t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
-	$base_path/fit_zzl_model temp_shuffled_t.dat $t_param_file
-	mv -f zzl_dump.qdp ${T_file}
-    else
-	echo temperature profile name invalid!
-	exit
-    fi
-    
-#exit
-    echo >temp_sbp.cfg
-    
-    cat $sbp_cfg|while read l
-do
-    if echo $l|grep sbp_file >/dev/null
-    then
-	echo sbp_file temp_shuffled_sbp.dat >>temp_sbp.cfg
-    elif echo $l|grep T_file >/dev/null
-    then
-	echo T_file ${T_file} >>temp_sbp.cfg
-    else
-	echo $l >>temp_sbp.cfg
-    fi
-    
-done
-
-#$base_path/coolfunc_calc.sh ${T_file} $abund $nh $z $cfunc_file
-
-$base_path/fit_dbeta_sbp temp_sbp.cfg
-
-cat entropy.qdp >>summary_entropy.qdp
-echo no no no >>summary_entropy.qdp
-done
-
-$base_path/analyze_entropy_profile.py $2 |tee entropy_result.txt
-