From 2ef85d7169e40a5475b6cebc2182a7b6deb5bcb9 Mon Sep 17 00:00:00 2001
From: Gregory Ashton <gregory.ashton@ligo.org>
Date: Sat, 26 Aug 2017 12:42:54 +0200
Subject: [PATCH] Fix issues in the Nstar calculation

- Converts coordinates to more usual Alpha. Delta F0, F1 ordering
- Adds debug info on the Nstar for each dimension
- Returns maximum Nstar iterating over each set of dimensions (not
  completely, just [0], [0, 1], [0, 1, 2], etc.
---
 pyfstat/optimal_setup_functions.py | 20 ++++++++++++++------
 1 file changed, 14 insertions(+), 6 deletions(-)

diff --git a/pyfstat/optimal_setup_functions.py b/pyfstat/optimal_setup_functions.py
index dda36c2..dc0c5fb 100644
--- a/pyfstat/optimal_setup_functions.py
+++ b/pyfstat/optimal_setup_functions.py
@@ -226,13 +226,21 @@ def get_Nstar_estimate(
     lalpulsar.ConvertPhysicalToSuperskyPoints(
         out_rssky, in_phys, SSkyMetric.semi_rssky_transf)
 
-    parallelepiped = (out_rssky.data[i:, 1:].T - out_rssky.data[i:, 0]).T
+    d = out_rssky.data
 
-    sqrtdetG = np.sqrt(np.linalg.det(
-        SSkyMetric.semi_rssky_metric.data[i:, i:]))
+    g = SSkyMetric.semi_rssky_metric.data
 
-    dV = np.abs(np.linalg.det(parallelepiped))
+    d[2:] = d[2:][::-1]  # Convert to Alpha, Delta, F0, F1.. ordering
+    g[2:] = g[2:][::-1]  # Convert to Alpha, Delta, F0, F1.. ordering
 
-    Nstar = sqrtdetG * dV
+    g = g[i:, i:]  # Remove sky if required
+    parallelepiped = (out_rssky.data[i:, 1:].T - out_rssky.data[i:, 0]).T
 
-    return Nstar
+    Nstars = []
+    for j in range(1, len(g)+1):
+        dV = np.abs(np.linalg.det(parallelepiped[:j, :j]))
+        sqrtdetG = np.sqrt(np.abs(np.linalg.det(g[:j, :j])))
+        Nstars.append(sqrtdetG * dV)
+    logging.debug('Nstar for each dimension = {}'.format(
+        ', '.join(["{:1.1e}".format(n) for n in Nstars])))
+    return np.max(Nstars)
-- 
GitLab