Separate tolerances for cont frac, cf_tol, & omega search, tol.

Starts to address #22

qnm/cached.py

@@ -349,13 +349,17 @@ def build_package_default_cache(ksc):
     QNMDict(init=True)
 
     a_max   = .9995
-    tol     = np.sqrt(np.finfo(float).eps)
+    tol     = 1e-10
+    cf_tol  = np.sqrt(np.finfo(float).eps)
     delta_a = 2.5e-3
     Nr_max  = 6000
 
-    ksc.compute_pars.update({'a_max': a_max, 'tol': tol,
+    ksc.compute_pars.update({'a_max': a_max, 'tol': tol, 'cf_tol': cf_tol,
                              'delta_a': delta_a, 'Nr_max': Nr_max})
 
+    # Known modes that become algebraically special
+    modes_to_avoid = [(-1, 1, 0, 6)]
+
     reruns = []
 
     ns=np.arange(0,7)
@@ -366,24 +370,49 @@ def build_package_default_cache(ksc):
             ms=np.arange(-l,l+1)
             for m in ms:
                 for n in ns:
+                    if ((s, l, m, n) in modes_to_avoid):
+                        print("Skipping known mode {}".format((s, l, m, n)))
+                        continue
                     try:
                         ksc(s, l, m, n)
                     except:
                         reruns.append((s,l,m,n))
 
+    print('{} modes in cache, going to rerun {} modes'.format(len(ksc.seq_dict),
+                                                              len(reruns)))
+
     # Tweak the params a little bit for the reruns
     re2runs = []
-    ksc.compute_pars.update({'delta_a': 1.9e-3})
+    ksc.compute_pars.update({'delta_a': 1.9e-3, 'cf_tol': cf_tol * 0.25, 'tol': tol * 3.})
     for s, l, m, n in reruns:
         try:
             ksc(s, l, m, n)
         except:
             re2runs.append((s, l, m, n))
 
-    # Experimentally determined we only need one more case of reruns
-    ksc.compute_pars.update({'a_max': .9992, 'delta_a': 1.8e-3})
+    print('{} modes in cache, going to rerun {} modes'.format(len(ksc.seq_dict),
+                                                              len(re2runs)))
+
+    ksc.compute_pars.update({'a_max': .9992, 'delta_a': 1.7e-3,
+                             'cf_tol': cf_tol * 0.05, 'tol': tol * 10.,
+                             'Nr_max': Nr_max * 2})
 
+    re3runs = []
     for s, l, m, n in re2runs:
+        try:
+            ksc(s, l, m, n)
+        except:
+            re3runs.append((s, l, m, n))
+
+    print('{} modes in cache, going to rerun {} modes'.format(len(ksc.seq_dict),
+                                                              len(re3runs)))
+
+    ksc.compute_pars.update({'a_max': .9993, 'delta_a': 0.24e-3,
+                             'cf_tol': cf_tol * 0.01, 'tol': tol * 70.,
+                             'Nr_max': Nr_max * 3})
+
+    # This is the last round of reruns we need to do
+    for s, l, m, n in re3runs:
         ksc(s, l, m, n)
 
     return ksc
@@ -464,7 +493,7 @@ def _clear_disk_cache(base_dir=None, delete_tarball=False):
     """Delete disk cache of precomputed spin sequences."""
 
     if base_dir is None:
-        base_dir = get_cachedir().parent.resolve()
+        base_dir = get_cachedir()
 
     data_dir = base_dir / 'data'
     pickle_files = data_dir.glob('*.pickle')

qnm/nearby.py

@@ -46,8 +46,11 @@ class NearbyRootFinder(object):
     omega_guess: complex [default: .5-.5j]
       Initial guess of omega for root-finding
 
-    tol: float [default: 1e-10]
-      Tolerance for root-finding
+    tol: float [default: sqrt(double epsilon)]
+      Tolerance for root-finding omega
+
+    cf_tol: float [defailt: 1e-10]
+      Tolerance for continued fraction calculation
 
     n_inv: int [default: 0]
       Inversion number of radial infinite continued fraction,
@@ -78,7 +81,8 @@ def __init__(self, *args, **kwargs):
         self.A0          = 4.+0.j
         self.l_max       = 20
         self.omega_guess = .5-.5j
-        self.tol         = 1e-10
+        self.tol         = np.sqrt(np.finfo(float).eps)
+        self.cf_tol      = 1e-10
         self.n_inv       = 0
         self.Nr          = 300
         self.Nr_min      = 300
@@ -101,6 +105,7 @@ def set_params(self, *args, **kwargs):
         self.l_max       = kwargs.get('l_max',        self.l_max)
         self.omega_guess = kwargs.get('omega_guess',  self.omega_guess)
         self.tol         = kwargs.get('tol',          self.tol)
+        self.cf_tol      = kwargs.get('cf_tol',       self.cf_tol)
         self.n_inv       = kwargs.get('n_inv',        self.n_inv)
         self.Nr          = kwargs.get('Nr',           self.Nr)
         self.Nr_min      = kwargs.get('Nr_min',       self.Nr_min)
@@ -151,15 +156,11 @@ def __call__(self, x):
         #                                            self.Nr, self.r_N)
 
         # TODO!
-        # The tolerance that was sent to optimize.root is the desired
-        # tolerance for omega.
-        # However, the tolerance that's sent to the Lentz method is
-        # the desired tolerance for the value of the error function
-        # E(\omega) whose root we desire.  These should be related by
-        # the Jacobian, Etol = |dE/d\omega| omegatol.
+        # Determine the value to use for cf_tol based on
+        # the Jacobian, cf_tol = |d cf(\omega)/d\omega| tol.
         inv_err, self.cf_err, self.n_frac = radial.leaver_cf_inv_lentz(omega, self.a,
                                                           self.s, self.m, A,
-                                                          self.n_inv, self.tol,
+                                                          self.n_inv, self.cf_tol,
                                                           self.Nr_min, self.Nr_max)
         # logging.info("Lentz terminated with cf_err={}, n_frac={}".format(self.cf_err, self.n_frac))
 

qnm/spinsequence.py

@@ -60,8 +60,11 @@ class KerrSpinSeq(object):
     omega_guess: complex [default: from schwarzschild.QNMDict]
       Initial guess of omega for root-finding
 
-    tol: float [default: 1e-10]
-      Tolerance for root-finding
+    tol: float [default: sqrt(double epsilon)]
+      Tolerance for root-finding omega
+
+    cf_tol: float [defailt: 1e-10]
+      Tolerance for continued fraction calculation
 
     n: int [default: 0]
       Overtone number of interest (sets the inversion number for
@@ -95,7 +98,8 @@ def __init__(self, *args, **kwargs):
         self.m           = kwargs.get('m',           2)
         self.l           = kwargs.get('l',           2)
         self.l_max       = kwargs.get('l_max',       20)
-        self.tol         = kwargs.get('tol',         1e-10)
+        self.tol         = kwargs.get('tol',         np.sqrt(np.finfo(float).eps))
+        self.cf_tol      = kwargs.get('cf_tol',      1e-10)
         self.n           = kwargs.get('n',           0)
 
         if ('omega_guess' in kwargs.keys()):
@@ -151,6 +155,9 @@ def do_find_sequence(self):
         i  = 0  # TODO Allow to start at other values
         _a = 0. # TODO Allow to start at other values
 
+        # Flag about whether we've warned re: imag axis
+        warned_imag_axis = False
+
         # Initializing the sequence, start with guesses
         A0 = swsphericalh_A(self.s, self.l, self.m)
         omega_guess = self.omega_guess
@@ -178,6 +185,17 @@ def do_find_sequence(self):
             # even if cf_err is larger than the desired error tol?
             cf_err, n_frac = self.solver.get_cf_err()
 
+            # Warn if we're very close to the imaginary axis
+            if ((np.abs(np.real(result)) < self.tol)
+                and not warned_imag_axis):
+                logging.warn("Danger! At a={}, found Re[omega]={}, "
+                             "twithin tol={} of the imaginary axis. "
+                             "this mode may become algebraically "
+                             "special, where Leaver's method would "
+                             "fail, or the mode may even disappear."
+                             .format(_a, result, self.tol))
+                warned_imag_axis = True
+
             # Done with this value of a
             self.a.append(_a)
             self.omega.append(result)