Merge branch 'posmoves-recent-rehome' into prepare_posparams

bin/fit_posparams

@@ -41,9 +41,9 @@ parser.add_argument('-i', '--infiles', type=str, required=True, nargs='*',
                          'Multiple file args also ok (like M0001.csv M0002.csv M01*.csv)')
 parser.add_argument('-o', '--outdir', type=str, required=True,
                     help='path to directory where to save output files')
-parser.add_argument('-dw', '--data_window', type=int, default=100,
+parser.add_argument('-dw', '--data_window', type=int, default=1000000,
                     help='minimum width (num pts) for window swept through historical data')
-parser.add_argument('-pd', '--period_days', type=int, default=10000,
+parser.add_argument('-pd', '--period_days', type=int, default=1000000,
                     help='spacing of datum dates at which to run the fit')
 parser.add_argument('-np', '--n_processes_max', type=int, default=None,
                     help='max number of processors to use. Note: can argue 1 to avoid multiprocessing pool (useful for debugging)')

py/desimeter/circles.py

@@ -94,7 +94,11 @@ def _fast_fit_circle(x,y,use_median=False) :
     # 2 equations:
     # mx1 + nx1*s1 = mx2 + nx2*s2
     # my1 + ny1*s1 = my2 + ny2*s2
-    s1 = (ny[i2]*mx[i2]-nx[i2]*my[i2]-ny[i2]*mx[i1]+nx[i2]*my[i1])/(ny[i2]*nx[i1]-nx[i2]*ny[i1])
+    num   = (ny[i2]*mx[i2]-nx[i2]*my[i2]-ny[i2]*mx[i1]+nx[i2]*my[i1])
+    denom = (ny[i2]*nx[i1]-nx[i2]*ny[i1])
+    ok = (denom!=0)
+    s1 = np.zeros(nn)
+    s1[ok] = num[ok]/denom[ok]
 
     # coordinates of intersections are estimates of center of circle
     xc=mx[i1]+nx[i1]*s1[i1]

py/desimeter/posparams/fithandler.py

@@ -77,18 +77,57 @@ def run_best_fits(posid, path, period_days, data_window, savedir,
     datum_dates = datum_dates.tolist()
     cases = _define_cases(table, datum_dates, data_window, printf=printf)
 
-    # FIRST-PASS:  STATIC PARAMETERS
+
+    if not np.all(table['RECENT_REHOME']) :
+        print("{}: has exposure(s) with RECENT_REHOME=False. Will first initiate the fit without those".format(posid))
+        recent_rehome = table['RECENT_REHOME'].astype(int)
+        selection = np.where(recent_rehome>0)[0]
+        print("{}: number of points with RECENT_REHOME=True = {}".format(posid,selection.size))
+        tmp_table = table[selection]
+        tmp_cases = _define_cases(tmp_table, datum_dates, data_window, printf=printf)
+        static_out = _process_cases(tmp_table, tmp_cases, printf=printf, mode='static',
+                                    param_nominals=fitter.default_values.copy())
+        flags=static_out["FLAGS"]
+        if np.any((flags&movemask["FAILED_FIT"])>0) :
+            flag=0
+            for f in flags : flag |=f
+            write_failed_fit(posid,savedir,flag)
+            return "{}: (in run_best_fits) has FAILED_FIT {}".format(posid,flags.tolist())
+
+        # DECIDE ON BEST STATIC PARAMS
+        best_static = fitter.default_values.copy()
+        errors = static_out['FIT_ERROR_STATIC']
+        printf('{}: rms of residuals (init static) = {}'.format(posid,list(errors)))
+        quantile = np.percentile(errors, static_quantile * 100)
+        selection = static_out[errors <= quantile]
+        these_best = {key:np.median(selection[key + _mode_suffix('static')]) for key in best_static}
+        best_static.update(these_best)
+        printf(f'{posid}: Selected best static params = {_dict_str(best_static)}')
+
+    else :
+        best_static = fitter.default_values.copy()
+
+
+    # FIRST-PASS (or second) :  STATIC PARAMETERS
     static_out = _process_cases(table, cases, printf=printf, mode='static',
-                                param_nominals=fitter.default_values.copy(),
+                                param_nominals=best_static,
                                 )
 
+    flags=static_out["FLAGS"]
+    if np.any((flags&movemask["FAILED_FIT"])>0) :
+        flag=0
+        for f in flags : flag |=f
+        write_failed_fit(posid,savedir,flag)
+        return "{}: (in run_best_fits) has FAILED_FIT {}".format(posid,flags.tolist())
+
     # DECIDE ON BEST STATIC PARAMS
     best_static = fitter.default_values.copy()
     # add EPSILONs parameters if exist
     for k in static_out.keys() :
         if k.find("EPSILON")>=0 :
             best_static[k.replace("_STATIC","")]=0.
     errors = static_out['FIT_ERROR_STATIC']
+    printf('{}: rms of residuals (static) = {}'.format(posid,list(errors)))
     quantile = np.percentile(errors, static_quantile * 100)
     selection = static_out[errors <= quantile]
     these_best = {key:np.median(selection[key + _mode_suffix('static')]) for key in best_static}
@@ -99,6 +138,8 @@ def run_best_fits(posid, path, period_days, data_window, savedir,
     dynamic_out = _process_cases(table, cases, printf=printf, mode='dynamic',
                                  param_nominals=best_static)
 
+    printf('{}: rms of residuals (dynamic) = {}'.format(posid,list(dynamic_out['FIT_ERROR_DYNAMIC'])))
+
     # MERGED STATIC + DYNAMIC
     same_keys = [key for key, is_variable in output_keys.items() if not is_variable]
     merged = join(static_out, dynamic_out, keys=same_keys)
@@ -317,6 +358,7 @@ def _process_cases(table, cases, mode, param_nominals, printf=print):
                                                                       nominals=param_nominals,
                                                                       bounds=fitter.default_bounds,
                                                                       keep_fixed=[])
+
         output['ANALYSIS_DATE'].append(Time.now().iso)
         output['POS_ID'].append(posid)
         for suffix in {'', '_SEC'}:
@@ -329,7 +371,7 @@ def _process_cases(table, cases, mode, param_nominals, printf=print):
 
         for key in params.keys() :
             if key not in fitted_keys and key.find("EPSILON")>=0 :
-                print("warning adding key={} which was not in fitted_keys={}".format(key,fitted_keys))
+                #print("warning adding key={} which was not in fitted_keys={}".format(key,fitted_keys))
                 fitted_keys.append(key)
                 if key not in output.keys() :
                     output[key] = list()
@@ -362,8 +404,11 @@ def _process_cases(table, cases, mode, param_nominals, printf=print):
                     output[okey].append(0)
 
         printf(f'{posid}: best params = {_dict_str(params)}')
-        printf(f'  fit error = {rms_of_residuals:.3f}')
+        printf(f'{posid}: fit error = {rms_of_residuals:.3f}')
 
+    keys=list(output.keys())
+    for k in keys :
+        if len(output[k])==0 : output.pop(k)
     table = Table(output)
 
     if True : # drop columns with empty covariance to make the output file a bit smaller

py/desimeter/posparams/fitter.py

@@ -194,6 +194,16 @@ def fit_params(posintT, posintP, ptlX, ptlY, gearT, gearP, recent_rehome, sequen
     xerr_flat = np.abs(d_x_flat_d_x_ptl) * xerr_ptl
     yerr_flat = np.abs(d_y_flat_d_y_ptl) * yerr_ptl
 
+    # we need to set one new pair of params EPSILON_P,T for each sequence_id except
+    # for the one(s) with recent_rehome = true
+    # values are fit independently for the DYNAMIC case
+    sequences_without_recent_rehome = np.unique(sequence_id[recent_rehome==False])
+    if mode == 'static' :
+        for s in sequences_without_recent_rehome :
+            for k in ["EPSILON_T_{}".format(s) , "EPSILON_P_{}".format(s)] :
+                nominals[k]=0.
+                bounds[k]=(-200.,200)
+
     if mode == 'dynamic':
         n_pts -= 1 # since in dynamic mode we are using delta values
 
@@ -208,6 +218,13 @@ def fit_params(posintT, posintP, ptlX, ptlY, gearT, gearP, recent_rehome, sequen
             nominals['LENGTH_R1'], nominals['LENGTH_R2'],
             nominals['OFFSET_T'], nominals['OFFSET_P']
             )
+
+        # need to subtract epsilons if exist to measured_t_int_0, measured_p_int_0
+        for s in sequences_without_recent_rehome :
+            ii=(sequence_id==s)
+            measured_t_int_0[ii] -= nominals["EPSILON_T_{}".format(s)]
+            measured_p_int_0[ii] -= nominals["EPSILON_P_{}".format(s)]
+
         del unreachable
         dt_int = pos2ptl.delta_angle(u_start=t_int[:-1], u_final=t_int[1:], direction=0)
         dp_int = pos2ptl.delta_angle(u_start=p_int[:-1], u_final=p_int[1:], direction=0)
@@ -230,17 +247,18 @@ def fit_params(posintT, posintP, ptlX, ptlY, gearT, gearP, recent_rehome, sequen
         measured_p_int_0 = measured_p_int_0[:-1]
         sequence_id   = sequence_id[1:]
         recent_rehome = recent_rehome[1:]
+        # also remove first point for t_int,p_int
+        t_int = t_int[1:]
+        p_int = p_int[1:]
 
 
-    # we need to set one new pair of params EPSILON_P,T for each sequence_id except
-    # for the one(s) with recent_rehome = true
-    # values are fit independently for the DYNAMIC case
-    sequences_without_recent_rehome = np.unique(sequence_id[recent_rehome==False])
-    if mode == 'static' :
-        for s in sequences_without_recent_rehome :
-            for k in ["EPSILON_T_{}".format(s) , "EPSILON_P_{}".format(s)] :
-                nominals[k]=0.
-                bounds[k]=(-200.,200)
+
+    if mode == 'dynamic' : # keep epsilons fixed in dynamic mode
+        keep_fixed = list(keep_fixed)
+        for k in list(nominals.keys()) :
+            if k.find("EPSILON")>=0 and k not in keep_fixed :
+                keep_fixed.append(k)
+        keep_fixed = set(keep_fixed)
 
     params_to_fit = set(nominals).difference(keep_fixed)
     params_to_fit = sorted(params_to_fit) # because optimizer expects a vector with consistent ordering
@@ -260,9 +278,11 @@ def compute_offset_and_delta_tp(sequence_id,p0) :
         return offset_and_delta_t,offset_and_delta_p
 
     def expected_xy(params):
+
+        for key, idx in param_idx.items():
+            p0[key] = params[idx]
+
         if mode =='static':
-            for key, idx in param_idx.items():
-                p0[key] = params[idx]
             t_int_expected = t_int
             p_int_expected = p_int
         else:
@@ -325,6 +345,7 @@ def compute_chi2(params):
         optimizer_result = scipy.optimize.minimize(fun=compute_chi2,
                                                    x0=initial_params,
                                                    bounds=bounds_vector)
+
         x_exp, y_exp = expected_xy(optimizer_result.x)
         dist = np.sqrt( (x_exp - x_flat)**2 + (y_exp - y_flat)**2 )
         meddist = np.median(dist)
@@ -357,6 +378,7 @@ def compute_chi2(params):
             dp_int = np.delete(dp_int,worst)
             measured_t_int_0 = np.delete(measured_t_int_0,worst)
             measured_p_int_0 = np.delete(measured_p_int_0,worst)
+
         npts=x_flat.size
         npar=len(initial_params)
         if npts<npar :