"""Analytic correlated estimates, clock drift and delayed acceleration subtraction."""
import math,csv,json
from pathlib import Path
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
P=Path(__file__).resolve().parent
fusion=[dict(correlation=-1+i/100,shared_bias_units=b,single_noise_variance_units2=1,average_variance_units2=(1+(-1+i/100))/2,average_MSE_units2=(1+(-1+i/100))/2+b*b) for b in [0,.5,1] for i in range(201)]
timing=[dict(elapsed_s=i*60,clock_rate_error_ppm=10,clock_offset_s=i*60*10e-6,constant_velocity_m_s=2,position_mismatch_m=2*i*60*10e-6) for i in range(181)]
motion=[]
for delay in [0,.01,.1]:
 for i in range(401):
  t=i/200;a=math.sin(2*math.pi*t);old=math.sin(2*math.pi*(t-delay));res=a-old
  motion.append(dict(time_s=t,delay_s=delay,acceleration_m_s2=a,specific_force_up_m_s2=9.81+a,delayed_acceleration_m_s2=old,gravity_estimate_m_s2=9.81+res,residual_m_s2=res))
for name,rows in [('correlated-average.csv',fusion),('clock-drift.csv',timing),('motion-subtraction.csv',motion)]:
 with (P/name).open('w',newline='') as f:w=csv.DictWriter(f,fieldnames=rows[0]);w.writeheader();w.writerows(rows)
err=0
for r in fusion:
 rho=r['correlation'];independent=(1+1+2*rho)/4;assert abs(independent-r['average_variance_units2'])<1e-14
 assert r['average_MSE_units2']>=r['average_variance_units2']>=0
 assert abs(r['average_MSE_units2']-r['average_variance_units2']-r['shared_bias_units']**2)<1e-14
for r in motion:
 analytic=2*math.cos(2*math.pi*r['time_s']-math.pi*r['delay_s'])*math.sin(math.pi*r['delay_s']);err=max(err,abs(analytic-r['residual_m_s2']))
 assert abs(r['gravity_estimate_m_s2']-(r['specific_force_up_m_s2']-r['delayed_acceleration_m_s2']))<1e-14
assert err<1e-14
assert all(r['gravity_estimate_m_s2']==9.81 for r in motion if r['delay_s']==0)
assert math.isclose(timing[-1]['clock_offset_s'],.108) and math.isclose(timing[-1]['position_mismatch_m'],.216)
assert timing[0]['clock_offset_s']==0
assert 0+9.81==9.81 and -9.81+9.81==0 # supported rest and freefall specific force
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False,'svg.fonttype':'none','pdf.fonttype':42})
fig,ax=plt.subplots(figsize=(8,4.5))
for b,c in [(0,'#087e91'),(.5,'#a14786'),(1,'#7451a0')]:
 rr=[r for r in fusion if r['shared_bias_units']==b];ax.plot([r['correlation'] for r in rr],[r['average_MSE_units2'] for r in rr],color=c,label=f'Shared bias {b:g} unit')
ax.set(xlabel='Noise correlation',ylabel='Average mean squared error (units squared)',title='Averaging does not remove a shared bias');ax.legend(fontsize=9);ax.grid(alpha=.2);fig.text(.5,.035,'Two estimates of the same scalar; each noise standard deviation is 1 unit; equal weights. Deterministic assigned bias, no measurements.',ha='center',fontsize=8);fig.text(.5,.01,'Perfect anticorrelation is an ideal algebraic cancellation of noise, not proof of real certainty.',ha='center',fontsize=8);fig.tight_layout(rect=[0,.075,1,1])
for ext in ['png','svg','pdf']:fig.savefig(P/f'correlation-and-bias.{ext}',dpi=190)
fig,ax=plt.subplots(1,2,figsize=(10.5,4.5));ax[0].plot([r['elapsed_s']/3600 for r in timing],[r['clock_offset_s']*1000 for r in timing],color='#087e91');ax[0].set(xlabel='Elapsed time (hours)',ylabel='Clock offset (ms)',title='A  |  Constant 10 ppm rate error accumulates')
for delay,c in [(0,'#777'),(.01,'#087e91'),(.1,'#a14786')]:
 rr=[r for r in motion if r['delay_s']==delay];ax[1].plot([r['time_s'] for r in rr],[r['residual_m_s2'] for r in rr],color=c,label=f'Delay {delay:g} s')
ax[1].set(xlabel='Time (s)',ylabel='Residual acceleration (m/s²)',title='B  |  Delayed subtraction leaves motion behind');ax[1].legend(fontsize=8)
for a in ax:a.grid(alpha=.2)
fig.suptitle('Measurements must refer to compatible times',fontsize=14);fig.text(.5,.04,'A: no initial offset or resynchronisation; at 3 hours the offset is 108 ms, giving 0.216 m mismatch at assigned constant speed 2 m/s.',ha='center',fontsize=8);fig.text(.5,.012,'B: assigned upward acceleration sin(2πt) m/s²; ideal g = 9.81 m/s². Delay example is separate from the clock-drift scenario.',ha='center',fontsize=8);fig.tight_layout(rect=[0,.085,1,.92])
for ext in ['png','svg','pdf']:fig.savefig(P/f'timing-and-motion.{ext}',dpi=190)
meta={'title':'Correlated averages and synchronised motion correction','status':'Deterministic analytical teaching study; no sensor measurements or accuracy claims','rows':{'fusion':len(fusion),'timing':len(timing),'motion':len(motion)},'parameters':{'noise_sd_units':1,'correlation_range':[-1,1],'correlation_step':.01,'shared_bias_units':[0,.5,1],'clock_rate_error_ppm':10,'clock_duration_hours':3,'constant_velocity_m_s':2,'gravity_m_s2':9.81,'acceleration_amplitude_m_s2':1,'acceleration_frequency_Hz':1,'delays_s':[0,.01,.1]},'equations':{'average_variance':'sigma²*(1+rho)/2','MSE':'average_variance+bias²','clock_offset':'fractional_rate_error*elapsed_time','position_mismatch':'constant_velocity*clock_offset','specific_force_up':'a(t)+g','gravity_estimate':'g+a(t)-a(t-delay)'},'checks':{'all_pass':True,'covariance_expansion':True,'shared_bias_MSE_identity':True,'delay_trig_identity_max_error_m_s2':err,'zero_delay':True,'clock_108ms_at3h':True,'zero_clock_offset':True,'stationary_and_freefall_convention':True},'limits':['Only equal-weight estimates of the same scalar in compatible units are averaged. Noise covariance must be physically meaningful; no independence is inferred from sensor count.','MSE includes an assigned deterministic common bias. A variance reduction is not a guarantee of improved real accuracy.','Perfect correlation endpoints are singular ideal covariance cases; no covariance inversion is performed.','Clock example omits initial offset, jitter, network delay, synchronisation correction and packet loss. Constant-velocity mismatch is a separate kinematic illustration.','Vertical specific force is expressed in a known upward inertial axis, with downward gravity. No body-axis orientation error, Earth rotation, vibration, lever arm, drift or gravity field model.','The sinusoid is defined before t=0 so delayed samples need no invented zero padding.','No classifier, unknown-class detection, altitude estimate, physiological validation or physical sensor performance is established.']}
(P/'study.json').write_text(json.dumps(meta,indent=2)+'\n');(P/'validation.json').write_text(json.dumps(meta['checks'],indent=2)+'\n');print(meta['rows'])