From 2ae719750d7e74666468a520478294dc9c57bb74 Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Mon, 8 Mar 2021 06:38:43 +0900
Subject: [PATCH] conditional pdf codes using fastKDE

---
 p_ddxc_given_c.py | 56 +++++++++++++++++++++++++++++++++++++++++++++++
 p_u_given_c.py    | 56 +++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 112 insertions(+)
 create mode 100644 p_ddxc_given_c.py
 create mode 100644 p_u_given_c.py

diff --git a/p_ddxc_given_c.py b/p_ddxc_given_c.py
new file mode 100644
index 0000000..a6b04a0
--- /dev/null
+++ b/p_ddxc_given_c.py
@@ -0,0 +1,56 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+import sys
+import argparse
+import numpy as np
+from fastkde import fastKDE
+
+
+# Commandline argument parser
+parser = argparse.ArgumentParser()
+parser.add_argument("-x", "--xindex", help="sampling x index", type=int, required=True)
+parser.add_argument("-p", "--point-sigma", help="number of pdf grid points per sigma", default=20, )
+args = parser.parse_args()
+params = vars(args)
+
+
+xidx = int(params["xindex"])
+numPointsPerSigma = int(params["point_sigma"])
+
+sigwidth = 1.
+
+print ("Computing pdf at {}".format(xidx))
+
+
+def sigmoid (y, a=1.0):
+    eps = np.finfo(np.float).eps
+    return 1. / (1. + np.exp(-y/a))
+
+def isigmoid(x, a=1.):
+    eps = np.finfo(np.float).eps
+    return - a * np.log((1.+2*eps)/(x+eps) - 1.)
+
+def disigmoid(x, a=1.):
+    eps = np.finfo(np.float).eps
+    return - a * (1./((1.+2*eps)/(x+eps) - 1.)) * (- 1. / ((x+eps)**2.))
+
+
+sc = np.memmap("c.dat", mode="r", dtype=np.double).reshape((512,-1,256,256))
+sabsk = np.memmap("ddxc.dat", mode="r", dtype=np.double).reshape((512,-1,256,256))
+
+c_absk_pdf, (ax1, ax2) = fastKDE.conditional(
+        sabsk[xidx].ravel(), 
+        isigmoid(sc[xidx].ravel(), a=sigwidth), 
+        peakFrac=0.001,
+        numPointsPerSigma=numPointsPerSigma
+        )
+
+finite_ax1 = sigmoid(ax1, a=sigwidth)
+
+arr_dict = {}
+arr_dict["cpdf"] = np.asarray(c_absk_pdf)
+arr_dict["ax1"] = np.asarray(ax1)
+arr_dict["ax2"] = np.asarray(ax2)
+
+np.savez("p_ddxc_given_c_{:03d}".format(xidx), **arr_dict)
diff --git a/p_u_given_c.py b/p_u_given_c.py
new file mode 100644
index 0000000..4a163b5
--- /dev/null
+++ b/p_u_given_c.py
@@ -0,0 +1,56 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+import sys
+import argparse
+import numpy as np
+from fastkde import fastKDE
+
+
+# Commandline argument parser
+parser = argparse.ArgumentParser()
+parser.add_argument("-x", "--xindex", help="sampling x index", type=int, required=True)
+parser.add_argument("-p", "--point-sigma", help="number of pdf grid points per sigma", default=20, )
+args = parser.parse_args()
+params = vars(args)
+
+
+xidx = int(params["xindex"])
+numPointsPerSigma = int(params["point_sigma"])
+
+sigwidth = 1.
+
+print ("Computing pdf at {}".format(xidx))
+
+
+def sigmoid (y, a=1.0):
+    eps = np.finfo(np.float).eps
+    return 1. / (1. + np.exp(-y/a))
+
+def isigmoid(x, a=1.):
+    eps = np.finfo(np.float).eps
+    return - a * np.log((1.+2*eps)/(x+eps) - 1.)
+
+def disigmoid(x, a=1.):
+    eps = np.finfo(np.float).eps
+    return - a * (1./((1.+2*eps)/(x+eps) - 1.)) * (- 1. / ((x+eps)**2.))
+
+
+sc = np.memmap("c.dat", mode="r", dtype=np.double).reshape((512,-1,256,256))
+sabsk = np.memmap("u.dat", mode="r", dtype=np.double).reshape((512,-1,256,256))
+
+c_absk_pdf, (ax1, ax2) = fastKDE.conditional(
+        sabsk[xidx].ravel(), 
+        isigmoid(sc[xidx].ravel(), a=sigwidth), 
+        peakFrac=0.001,
+        numPointsPerSigma=numPointsPerSigma
+        )
+
+finite_ax1 = sigmoid(ax1, a=sigwidth)
+
+arr_dict = {}
+arr_dict["cpdf"] = np.asarray(c_absk_pdf)
+arr_dict["ax1"] = np.asarray(ax1)
+arr_dict["ax2"] = np.asarray(ax2)
+
+np.savez("p_u_given_c_{:03d}".format(xidx), **arr_dict)