diff --git a/src/toolbox_scs/routines/boz.py b/src/toolbox_scs/routines/boz.py index 1f231334d476cb429bea8cd7bd4d7b3e37d6bd1d..f118fb80e1b5cd49545d3f4b8429becc6aceea4f 100644 --- a/src/toolbox_scs/routines/boz.py +++ b/src/toolbox_scs/routines/boz.py @@ -1071,16 +1071,17 @@ def nl_crit(p, domain, alpha, arr_dark, arr, tid, rois, mask, flat_field, # drop saturated shots d = data.where(data['sat_sat'] == False, drop=True) - # calculated error from transmission of 1.0 - v = d['n'].values.flatten()/d['0'].values.flatten() - err1 = 1e8*np.nanmean((v - 1.0)**2) + v_1 = snr(d['n'].values.flatten(), d['0'].values.flatten(), + methods=['weighted']) + err_1 = 1e8*v_1['weighted']['s']**2 - err2 = np.sum((Fmodel-np.arange(2**9))**2) + v_2 = snr(d['p'].values.flatten(), d['0'].values.flatten(), + methods=['weighted']) + err_2 = 1e8*v_2['weighted']['s']**2 - # print(f'{err}: {p}') - # logging.info(f'{err}: {p}') + err_a = np.sum((Fmodel-np.arange(2**9))**2) - return (1.0 - alpha)*err1 + alpha*err2 + return (1.0 - alpha)*0.5*(err_1 + err_2) + alpha*err_a def nl_fit(params, domain): @@ -1173,10 +1174,27 @@ def inspect_nl_fit(res_fit): return f -def snr(sig, ref, verbose=False): - """ Compute mean, std and SNR with and without weight from transmitted signal sig - and I0 signal ref +def snr(sig, ref, methods=None, verbose=False): + """ Compute mean, std and SNR from transmitted signal sig and I0 signal ref. + + Inputs + ------ + sig: 1D signal samples + ref: 1D reference samples + methods: None by default or list of strings to select which methods to use. + Possible values are 'direct', 'weighted', 'diff'. In case of None, all + methods will be calculated. + verbose: booleand, if True prints calculated values + + Returns + ------- + dictionnary of [methods][value] where value is 'mu' for mean and 's' for + standard deviation. + """ + if methods is None: + methods = ['direct', 'weighted', 'diff'] + w = ref x = sig/ref @@ -1187,30 +1205,37 @@ def snr(sig, ref, verbose=False): ref = ref[mask] x = x[mask] - # direct mean and std - mu = np.mean(x) - s = np.std(x) - if verbose: - print(f'mu: {mu}, s: {s}, snr: {mu/s}') res = {} - res['direct'] = {'mu': mu, 's':s} + + # direct mean and std + if 'direct' in methods: + mu = np.mean(x) + s = np.std(x) + if verbose: + print(f'mu: {mu}, s: {s}, snr: {mu/s}') + + res['direct'] = {'mu': mu, 's':s} # weighted mean and std - wmu = np.sum(sig)/np.sum(ref) - v1 = np.sum(w) - v2 = np.sum(w**2) - ws = np.sqrt(np.sum(w*(x - wmu)**2)/(v1 - v2/v1)) + if 'weighted' in methods: + wmu = np.sum(sig)/np.sum(ref) + v1 = np.sum(w) + v2 = np.sum(w**2) + ws = np.sqrt(np.sum(w*(x - wmu)**2)/(v1 - v2/v1)) + + if verbose: + print(f'weighted mu: {wmu}, s: {ws}, snr: {wmu/ws}') - if verbose: - print(f'weighted mu: {wmu}, s: {ws}, snr: {wmu/ws}') - res['weighted'] = {'mu': wmu, 's':ws} + res['weighted'] = {'mu': wmu, 's':ws} # noise from diff - dmu = np.mean(x) - ds = np.std(np.diff(x))/np.sqrt(2) - if verbose: - print(f'diff mu: {dmu}, s: {ds}, snr: {dmu/ds}') - res['diff'] = {'mu': dmu, 's':ds} + if 'diff' in methods: + dmu = np.mean(x) + ds = np.std(np.diff(x))/np.sqrt(2) + if verbose: + print(f'diff mu: {dmu}, s: {ds}, snr: {dmu/ds}') + + res['diff'] = {'mu': dmu, 's':ds} return res