diff --git a/bin/fit_posparams b/bin/fit_posparams index b119dad8..4390c2c5 100755 --- a/bin/fit_posparams +++ b/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)') diff --git a/py/desimeter/circles.py b/py/desimeter/circles.py index 03e8028d..a3ed8c8e 100644 --- a/py/desimeter/circles.py +++ b/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] diff --git a/py/desimeter/posparams/fithandler.py b/py/desimeter/posparams/fithandler.py index c0da6352..9b5e7357 100644 --- a/py/desimeter/posparams/fithandler.py +++ b/py/desimeter/posparams/fithandler.py @@ -77,11 +77,49 @@ 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 @@ -89,6 +127,7 @@ def run_best_fits(posid, path, period_days, data_window, savedir, 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 diff --git a/py/desimeter/posparams/fitter.py b/py/desimeter/posparams/fitter.py index 7f6f77ba..ae9bc350 100644 --- a/py/desimeter/posparams/fitter.py +++ b/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