bisection method + associated fancy charts

Author: Exr0n

calculate.py

@@ -7,8 +7,11 @@
 from operator import itemgetter # https://stackoverflow.com/a/52083390/10372825
 from ErrorProp import ErroredValue as EV
 
+from config import PRINT_PRECISION
+
 CM_PER_INCH = 2.54
 MAX_COUNTRATE = 3500
+GOOD_ENOUGH = 1e-9
 
 def process(dataitem):
     meta, data, background = itemgetter('meta', 'data', 'background')(dataitem)
@@ -46,31 +49,25 @@ def SSE(indicies, prediction, logits, logits_error, function):
     ri = function(prediction, indicies)
     return ((logits-ri)**2/logits_error**2).sum()
 
-def sMinFit(datatable, function, param=1, lr=1e-4, epsilon=1e-8, epochs=10000):
-    logits = datatable.normalized_count_rate.apply(lambda x:x.value)
-    logits_err = datatable.normalized_count_rate.apply(lambda x:x.delta)
-    inches  = datatable.inches
+def unwrap(datatable):
+    return datatable.inches, datatable.normalized_count_rate.apply(lambda x:x.value), datatable.normalized_count_rate.apply(lambda x:x.delta)
+
+def sMinFit(datatable, function, param=1, lr=1e-4, epsilon=1e-8, epochs=100000):
+    inches, logits, logits_err = unwrap(datatable)
     dydx = epsilon
 
     dydx = epsilon
     bar = tqdm.tqdm(range(epochs))
 
-    minParam = param
-    minS = int(1e10)
-
     for _ in bar:
         param_next = param-(dydx*lr+epsilon)
         loss = SSE(inches, param+epsilon, logits, logits_err, function)
         dydx = (loss-SSE(inches, param_next, logits, logits_err, function))/(param-param_next)
         param = param_next
 
-        bar.set_description(f'Current fit: {param:2f}, Best fit: {minParam:2f}, Best loss: {minS:2f}')
+        bar.set_description(f'Current fit: {param:2f}, Update: {dydx:2f}, Loss: {loss:2f}')
 
-        if loss < minS and param > 0:
-            minParam = param
-            minS = loss
-
-    return minParam, minS
+    return param, SSE(inches, param+epsilon, logits, logits_err, function)
 
     # breakpoint()
 
@@ -79,3 +76,43 @@ def sMinFit(datatable, function, param=1, lr=1e-4, epsilon=1e-8, epochs=10000):
     #     breakpoint()
     #     print(dydx)
     #     # weight update
+
+def calculateSfitUncert(bestx, besty, targety, function, ax=None, low=0, high=100):
+    def bisect(low, high, target, function):
+        triedx, triedy = [], []
+        if (low > high): low, high = high, low
+        while (high-low > GOOD_ENOUGH):
+            mid = (low+high)/2
+            triedx.append(mid)
+            triedy.append(function(mid))
+            if (function(mid) < target) == (function(low) < target):
+                low = mid
+            else:
+                high = mid
+        return low, (triedx, triedy)
+
+    # binary search on min and max
+    param_min, tries_min = bisect(low, bestx, targety, function)
+    param_max, tries_max = bisect(bestx, high, targety, function)
+
+    # plot everything
+    if ax is not None:
+        # plot neighborhood
+        abs_err = max(bestx - param_min, param_max-bestx)
+        neighborhood = np.arange(bestx-abs_err*2, bestx+abs_err*2, abs_err*4/200)   # 200 evenly spaced points
+        ax.scatter(neighborhood, list(map(function, neighborhood)), color='black', label='S(T)')
+
+        # # plot binary search tries
+        # ax.scatter(*tries_min, color='grey')
+        # ax.scatter(*tries_max, color='grey')
+
+        # plot given information
+        ax.axhline(y=besty, label=f"S min = {besty:.6f}", color='green')
+        ax.axhline(y=targety, label=f"S min + chi^2 = {targety:.6f}", color='lightgreen')
+        ax.axvline(x=bestx, label=f"T_best = {bestx:.6f}", color='green')
+
+        # plot found information
+        ax.axvline(x=param_min, label=f"T_min = {param_min:.6f}", color='red')
+        ax.axvline(x=param_max, label=f"T_min = {param_max:.6f}", color='blue')
+
+    return param_min, param_max

config.py

@@ -0,0 +1,2 @@
+PRINT_PRECISION = 6
+

drive.py

@@ -1,6 +1,6 @@
 from datacleaning import globit
 from multiprocessing import Pool
-from calculate import process, sMinFit, RelativeIntersity
+from calculate import SSE, process, unwrap, sMinFit, RelativeIntersity, calculateSfitUncert
 
 import matplotlib.pyplot as plt
 
@@ -62,3 +62,16 @@ def process(indx):
     # breakpoint()
 
 # >>>>>>> 68cb3cc3a83d6f09391e99a4e1cc04d712bebe16
+
+if __name__ == '__main__':
+    # plot(7)
+    fig, ax = plt.subplots()
+    inches, logits, logits_err = unwrap(results[7])
+    param = 0.946
+    smin = 82
+    calculateSfitUncert(param, smin, smin+1, lambda T: SSE(inches, T, logits, logits_err, RelativeIntersity), ax=ax, low=0.1, high=30)
+
+    ax.set_xlabel(f"T ({results[7].attrs['material']})")
+    ax.set_ylabel("S(T)")
+    ax.legend()
+    plt.show()

fits.org

@@ -0,0 +1,11 @@
+| job number | fit |
++------------+-----+
+|  2 | 0.426406 |
+|  3 | 0.016489 |
+|  4 | 0.347708 |
+|  5 | 0.019284 |
+|  6 |      nan |
+|  7 | 7.097304 |
+|  8 |      nan |
+|  9 | 0.945550 |
+| 10 | 0.019839 |