Convert PDF spectrum to CDF by default

appletree/config.py

@@ -2,9 +2,10 @@
 
 from immutabledict import immutabledict
 from jax import numpy as jnp
+from warnings import warn
 
 from appletree.share import _cached_configs
-from appletree.utils import exporter, load_json, get_file_path
+from appletree.utils import exporter, load_json, get_file_path, integrate_midpoint, cum_integrate_midpoint
 
 export, __all__ = exporter()
 
@@ -125,6 +126,12 @@ def build(self):
     def build_point(self, data):
         """Cache the map to jnp.array if bins_type is point"""
 
+        if data['coordinate_name'] == 'pdf':
+            warn(f'Convert {self.name} from (x, pdf) to (cdf, x).')
+            x, cdf = self.pdf_to_cdf(data['coordinate_system'], data['map'])
+            data['coordinate_name'] = 'cdf'
+            data['coordinate_system'] = cdf
+            data['map'] = x
         self.coordinate_name = data['coordinate_name']
         self.coordinate_system = jnp.asarray(data['coordinate_system'], dtype=float)
         self.map = jnp.asarray(data['map'], dtype=float)
@@ -136,3 +143,10 @@ def build_regbin(self, data):
         self.coordinate_lowers = jnp.asarray(data['coordinate_lowers'], dtype=float)
         self.coordinate_uppers = jnp.asarray(data['coordinate_uppers'], dtype=float)
         self.map = jnp.asarray(data['map'], dtype=float)
+
+    def pdf_to_cdf(self, x, pdf):
+        """Convert pdf map to cdf map"""
+        norm = integrate_midpoint(x, pdf)
+        x, cdf = cum_integrate_midpoint(x, pdf)
+        cdf /= norm
+        return x, cdf

appletree/utils.py

@@ -498,3 +498,22 @@ def _add_extension(module, subclass, base):
             )
         else:
             setattr(module, subclass.__name__, subclass)
+
+
+def integrate_midpoint(x, y):
+    """Calculate the integral using midpoint method.
+    :param x: 1D array-like 
+    :param y: 1D array-like, with the same length as x.
+    """
+    _, res = cum_integrate(x, y)
+    return res[-1]
+
+
+def cum_integrate_midpoint(x, y):
+    """Calculate the cumulative integral using midpoint method.
+    :param x: 1D array-like 
+    :param y: 1D array-like, with the same length as x."""
+    dx = x[1:] - x[:-1]
+    x_mid = 0.5 * (x[1:] + x[:-1])
+    y_mid = 0.5 * (y[1:] + y[:-1])
+    return x_mid, np.cumsum(dx * y_mid)