Merge pull request #90 from jmfee-usgs/change-sqdist-defaults

Set b0 to 0 and s0 to list of 0s by default

geomagio/algorithm/SqDistAlgorithm.py

@@ -330,8 +330,7 @@ class SqDistAlgorithm(Algorithm):
                 raise AlgorithmException("l0 must be a scalar")
 
         if b0 is None:
-            b = (np.nanmean(yobs[m:2 * m]) - np.nanmean(yobs[0:m])) / m
-            b = 0 if np.isnan(b) else b  # replace NaN with 0
+            b = 0
         else:
             b = b0
             if not np.isscalar(b0):
@@ -345,8 +344,7 @@ class SqDistAlgorithm(Algorithm):
                 raise AlgorithmException("yhat0 must have length %d" % hstep)
 
         if s0 is None:
-            s = [yobs[i] - l for i in range(m)]
-            s = [i if ~np.isnan(i) else 0 for i in s]  # replace NaNs with 0s
+            s = [0 for i in range(m)]
         else:
             s = list(s0)
             if len(s) != m:

test/algorithm_test/SQDistAlgorithm_test.py

@@ -146,6 +146,11 @@ def test_sqdistalgorithm_additive2():
     t000to050 = np.arange(5001)
     syn000to050 = 10.0 * np.sin(t000to050 * (2 * np.pi) / 100.0)
 
+    # these are the old "defaults" computed when l0, b0, and s0 were None
+    l0 = np.nanmean(syn000to050[0:1440])
+    b0 = (np.nanmean(syn000to050[m:2 * m]) - np.nanmean(syn000to050[0:m])) / m
+    s0 = [syn000to050[i] - l0 for i in range(m)]
+
     # run additive method on first 50 "days"
     (synHat000to050, sHat000to050, sigma000to050,
         syn050, s050, l050, b050, sigma050) = sq.additive(