"""Synthetic instrument qualification: thermal corners, outliers and aliasing."""
import math, json
import numpy as np
from measurement import fit_calibration,convert,moving_average,uncertainty_budget,ZeroTracker
from completion_support import *
from support import plot_style,save_figure

def synthetic_calibration(seed=722,batch_bias_g=0.):
 p=json.loads((ROOT/'inputs.json').read_text());rng=np.random.default_rng(seed)
 masses=np.repeat(np.array(p['calibration_masses_g'],dtype=float),p['calibration_repeats'])
 v=p['synthetic_offset_uv']+p['synthetic_slope_uv_per_g']*(masses+(masses==8000)*batch_bias_g)+rng.normal(0,2,len(masses))
 fit=fit_calibration(masses,v)
 group_means=[float(np.mean(fit['residual_uv'][masses==m]))/fit['slope_uv_per_g'] for m in p['calibration_masses_g']]
 return fit,group_means

def alias_trace(analog_filter,frequency_hz=93.,sample_hz=100.,amplitude_g=20.):
 t=np.arange(10000)/sample_hz
 h=(1+1j*frequency_hz/10.)**-2 if analog_filter else 1+0j
 samples=amplitude_g*abs(h)*np.sin(2*np.pi*frequency_hz*t+np.angle(h))
 return samples,moving_average(samples,25),h

def selected_mass_error_bound(fit,p):
 # Exact extrema of the declared sinusoidal systematic-error model over 0..10 kg.
 a=p['synthetic_slope_uv_per_g']/fit['slope_uv_per_g'];b=(p['synthetic_offset_uv']-fit['offset_uv'])/fit['slope_uv_per_g'];candidates=[]
 for temperature in [-2.,2.]:
  derivative_constant=(a-1)+a*temperature/10*.0001
  cosine=-derivative_constant/(a*2*math.pi/10000)
  masses=[0.,10000.]+([math.acos(cosine)*10000/math.pi] if -1<=cosine<=1 else [])
  for direction in [-1,1]:
   for m in masses:
    mean=b+(a-1)*m+a*(temperature/10*(1.5+.0001*m)+2*math.sin(math.pi*m/10000)+direction+.6)
    candidates.append(dict(mass_g=m,temperature_change_c=temperature,direction=direction,mean_error_g=mean))
 worst=max(candidates,key=lambda r:abs(r['mean_error_g']));noise=2/fit['slope_uv_per_g']/math.sqrt(25)
 return dict(worst_systematic=worst,filtered_noise_std_g=noise,conservative_95pct_absolute_error_g=abs(worst['mean_error_g'])+1.959963984540054*noise)

def main():
 p=json.loads((ROOT/'inputs.json').read_text());fit,_=synthetic_calibration();rng=np.random.default_rng(919);rows=[]
 for mass in [0.,2000.,5000.,10000.]:
  for delta_t in [-10.,0.,10.]:
   for direction in [-1,1]:
    thermal=delta_t/10*(p['temperature_zero_bound_g']+mass*p['temperature_span_bound_fraction'])
    nonlinearity=p['nonlinearity_bound_g']*math.sin(math.pi*mass/p['capacity_g'])
    hysteresis=direction*p['hysteresis_bound_g'];drift=60*p['drift_g_per_s']
    true_signal=p['synthetic_offset_uv']+p['synthetic_slope_uv_per_g']*(mass+thermal+nonlinearity+hysteresis+drift)
    raw=convert(true_signal+rng.normal(0,p['noise_standard_uv'],20000),fit)
    error=moving_average(raw,25)[100:]-mass
    mean=float(np.mean(error));q95=float(np.quantile(np.abs(error),.95))
    rows.append(dict(mass_g=mass,temperature_change_c=delta_t,direction='loading' if direction==1 else 'unloading',drift_at_60_s_g=drift,systematic_mean_error_g=mean,absolute_error_95th_g=q95,meets_reference_5_g=q95<=5.))
 write_csv('results/virtual-temperature-corners.csv',rows)
 selected=[];rng_selected=np.random.default_rng(921)
 for mass in [0.,2000.,5000.,10000.]:
  for delta_t in [-2.,2.]:
   for direction in [-1,1]:
    thermal=delta_t/10*(1.5+mass*.0001);nonlinearity=2*math.sin(math.pi*mass/10000);systematic=thermal+nonlinearity+direction+.6
    signal=p['synthetic_offset_uv']+p['synthetic_slope_uv_per_g']*(mass+systematic)
    errors=moving_average(convert(signal+rng_selected.normal(0,2,20000),fit),25)[100:]-mass
    q=float(np.quantile(np.abs(errors),.95));selected.append(dict(mass_g=mass,temperature_change_c=delta_t,hysteresis_direction=direction,max_elapsed_s=30,absolute_error_95th_g=q,meets_5g=q<=5.))
 write_csv('results/virtual-selected-envelope.csv',selected)
 selected_bound=selected_mass_error_bound(fit,p)
 write_json('selected-operating-envelope.json',dict(calibration_offset_uv=fit['offset_uv'],calibration_slope_uv_per_g=fit['slope_uv_per_g'],continuous_mass_error_bound=selected_bound,evidence_type='synthetic/theoretical',mass_range_g=[0,10000],temperature_change_from_calibration_max_c=2.,time_since_zero_max_s=30.,sample_hz=100,filter='MA 25',settled_before_read_s=.25,drift_bound_g=.6,calibration_gate='Each calibration-mass group mean residual <= 2 g; reference accuracy remains an assumption',screen='Empirical 95th absolute error <= 5 g',executed_corners=len(selected),maximum_corner_error_95th_g=max(r['absolute_error_95th_g'] for r in selected),all_corners_meet_screen=all(r['meets_5g'] for r in selected),physical_validation=False))
 calibration=[]
 for bias in [0.,10.,50.]:
  candidate,means=synthetic_calibration(batch_bias_g=bias);res=max(abs(x) for x in means)
  calibration.append(dict(injected_8000g_batch_bias_g=bias,max_group_residual_g=res,accepted_by_2g_group_gate=res<=2.,slope_uv_per_g=candidate['slope_uv_per_g']))
 write_csv('results/virtual-calibration-rejection.csv',calibration)
 noise_sd=2/fit['slope_uv_per_g']/math.sqrt(25);adc_step=(2*p['adc_reference_v']/p['adc_pga']*1e6/2**p['adc_bits'])/fit['slope_uv_per_g']
 _,independent=uncertainty_budget(fit,noise_sd,5000,adc_step,1.2,p)
 n=50000;rng=np.random.default_rng(920);coeff=rng.multivariate_normal([0,0],fit['covariance'],n);fit_error=(-coeff[:,0]-5000*coeff[:,1])/fit['slope_uv_per_g']
 temperature=rng.uniform(-1,1,n)*(1.5+.5)
 error=(fit_error+rng.normal(0,noise_sd,n)+rng.normal(0,.05,n)+rng.uniform(-2,2,n)+rng.uniform(-1,1,n)+rng.uniform(-1.2,1.2,n)+rng.uniform(-adc_step/2,adc_step/2,n)+temperature)
 correlated=math.sqrt(independent**2+2*1.5*.5/3)
 mc_std=float(np.std(error,ddof=1));coverage=float(np.mean(np.abs(error)<=2*correlated))
 uncertainty=dict(samples=n,independent_budget_standard_g=independent,temperature_correlated_standard_g=correlated,temperature_correlated_expanded_k2_g=2*correlated,monte_carlo_standard_g=mc_std,monte_carlo_2_5_percentile_g=float(np.quantile(error,.025)),monte_carlo_97_5_percentile_g=float(np.quantile(error,.975)),fraction_within_k2=coverage,physical_observations=False)
 write_json('results/virtual-correlated-uncertainty.json',uncertainty)
 alias=[];alias_samples={}
 for active in [False,True]:
  raw,filtered,h=alias_trace(active);alias.append(dict(analog_two_pole_10hz=active,source_hz=93,alias_hz=7,analog_amplitude_gain=abs(h),sampled_rms_g=float(np.std(raw[100:],ddof=0)),after_ma25_rms_g=float(np.std(filtered[100:],ddof=0))))
  alias_samples[active]=(raw,filtered)
 write_csv('results/virtual-aliasing-summary.csv',alias)
 write_csv('results/virtual-aliasing-trace.csv',[dict(time_s=i/100,raw_without_analog_g=alias_samples[False][0][i],ma25_without_analog_g=alias_samples[False][1][i],raw_with_analog_g=alias_samples[True][0][i],ma25_with_analog_g=alias_samples[True][1][i]) for i in range(400)])
 correct=ZeroTracker();wrong=ZeroTracker()
 for _ in range(300):correct.update(np.full(25,2.),empty_confirmed=False);wrong.update(np.full(25,2.),empty_confirmed=True)
 checks=[
 check('thermal-count','All declared thermal/load/direction corners executed',24,len(rows),len(rows)==24,'4 masses x 3 temperatures x 2 loading directions'),
 check('thermal-failures','Aligned thermal/systematic error can exceed 5 g','At least one adverse corner',sum(not x['meets_reference_5_g'] for x in rows),any(not x['meets_reference_5_g'] for x in rows),'Declared +/-5 g screening criterion on 95th absolute error; failures retained'),
 check('selected-envelope','Selected temperature/time envelope meets the error screen','All 16 corners <= 5 g',max(r['absolute_error_95th_g'] for r in selected),all(r['meets_5g'] for r in selected),'Independent seed; +/-2 C since calibration and <=30 s drift; both hysteresis directions'),
 check('continuous-envelope','Selected mass range satisfies a continuous theoretical bound','Conservative 95% absolute-error bound <= 5 g',selected_bound,selected_bound['conservative_95pct_absolute_error_g']<=5.,'Solve exact extrema of the declared sinusoidal error curve; add 1.96 filtered-noise standard deviations'),
 check('calibration-valid','Uncontaminated batch passes residual gate',True,calibration[0]['accepted_by_2g_group_gate'],calibration[0]['accepted_by_2g_group_gate'],'2 g maximum mean residual at each calibration mass'),
 check('calibration-corrupt','Biased mass group is rejected','Both 10 g and 50 g injected biases rejected',[x['accepted_by_2g_group_gate'] for x in calibration[1:]],all(not x['accepted_by_2g_group_gate'] for x in calibration[1:]),'Inject group-wide offset into 30 readings at 8 kg; fit with original OLS implementation'),
 check('covariance','Shared temperature source increases variance','Correlated budget greater than independent',dict(independent_g=independent,correlated_g=correlated),correlated>independent,'Add 2*u_zero*u_span for the shared rectangular temperature variable'),
 check('monte-carlo','Synthetic correlated uncertainty agrees with analytic variance','Standard deviation difference < 2%',abs(mc_std/correlated-1),abs(mc_std/correlated-1)<.02,'50000 seeded draws; all bounds and correlation explicit'),
 check('alias-frequency','Sampled 93 Hz tone aliases to 7 Hz','Dominant discrete frequency 7 Hz',float(np.fft.rfftfreq(10000,1/100)[np.argmax(np.abs(np.fft.rfft(alias_samples[False][0])))]),abs(np.fft.rfftfreq(10000,1/100)[np.argmax(np.abs(np.fft.rfft(alias_samples[False][0])))]-7)<1e-9,'100 Hz acquisition, frequency folding f_alias = abs(93 - 100)'),
 check('analog-suppression','Declared analog filter reduces the aliased amplitude','RMS ratio agrees with transfer magnitude within 1e-8',alias[1]['sampled_rms_g']/alias[0]['sampled_rms_g'],abs(alias[1]['sampled_rms_g']/alias[0]['sampled_rms_g']-alias[1]['analog_amplitude_gain'])<1e-8,'H(f)=(1+j*f/10)^-2; two identical 10 Hz RC poles'),
 check('zero-confirmation','Correct empty-pan gate preserves a real 2 g load',0.,correct.offset_g,correct.offset_g==0.,'300 updates with empty confirmation false'),
 check('wrong-confirmation','False empty confirmation absorbs a real load','Offset exceeds 1.9 g',wrong.offset_g,wrong.offset_g>1.9,'Same updates with deliberately incorrect external empty indication; adverse result'),
 ]
 plt=plot_style();fig,axes=plt.subplots(1,2,figsize=(10,4))
 for temperature in [-10.,0.,10.]:
  chosen=[r for r in rows if r['temperature_change_c']==temperature and r['direction']=='loading'];axes[0].plot([r['mass_g']/1000 for r in chosen],[r['absolute_error_95th_g'] for r in chosen],marker='o',label=f'{temperature:+g} C')
 axes[0].axhline(5,color='#e7a387',ls='--',label='Declared 5 g screen');axes[0].set(xlabel='Applied mass (kg)',ylabel='95th absolute error (g)',title='Loading direction, 60 s drift');axes[0].legend(fontsize=8);axes[0].grid()
 t=np.arange(100)/100;axes[1].plot(t,alias_samples[False][1][:100],label='MA 25, no analog filter');axes[1].plot(t,alias_samples[True][1][:100],label='MA 25 after two analog poles');axes[1].set(xlabel='Time (s)',ylabel='Indicated interference (g)',title='93 Hz tone sampled at 100 Hz');axes[1].legend(fontsize=8);axes[1].grid()
 fig.tight_layout();save_figure(fig,'virtual-instrument');plt.close(fig)
 finish('Load-cell virtual qualification','Executed thermal/load corners, contaminated-calibration rejection, correlated uncertainty, aliasing and zero-tracking fault cases using synthetic data.',[
 'The original 10 kg bridge, 100 Hz sample rate and 25-sample moving average are retained. The corner study adds temperature changes of -10, 0 and +10 C, masses 0/2/5/10 kg, both hysteresis directions and 1.2 g of accumulated drift at 60 s.',
 'The corner criterion is an explicitly selected +/-5 g screen: the empirical 95th percentile of absolute error must be <=5 g. It is not a legal metrology class or a promise of device accuracy. Each corner contains 20,000 generated observations; the first 100 filtered samples are discarded.',
 'Selected synthetic operating envelope: temperature remains within +/-2 C of calibration, zero is established at most 30 s before the reading, and the 25-sample filter is allowed 0.25 s to settle. Sixteen independent-seed load/temperature/direction corners meet the same 5 g screen. Exact extrema of the assumed sinusoidal systematic-error curve, plus 1.96 filtered-noise standard deviations, also bound the entire 0-10 kg interval below 5 g for this generated calibration. The wider +/-10 C, 60 s cases are rejected where they exceed it.',
 'The calibration diagnostic checks mean residual at each mass, rejecting a group exceeding 2 g. This is a completed synthetic acceptance gate for the declared dataset, not a universal detector of bad calibration. Common-mode reference-mass error or a perfectly linear span error can pass residual tests.',
 f'The shared-temperature study gives {correlated:.3f} g standard uncertainty and {2*correlated:.3f} g expanded uncertainty at k=2 and 5 kg. It uses a separate generated calibration and analytical filter noise. The uncertainty experiment draws 50,000 independent synthetic realizations of all budget terms. Zero and span temperature terms share one rectangular random variable, so their covariance is explicit. k=2 coverage is reported empirically for this assumed error distribution only.',
 'The aliasing example is a 20 g-amplitude, 93 Hz interference tone sampled at 100 Hz. A hypothetical analog filter has two identical 10 Hz first-order poles. Its exact complex response is applied before sampling. Digital averaging acts after aliasing and cannot identify the original frequency.'
 ],rows,checks,[
 'This completes the theoretical and synthetic qualification run-through. No real bridge, ADC, temperature chamber, weights or calibration certificate was used.',
 'The analog filter is an ideal transfer function, not a released circuit. ADC common-mode limits, reference tracking, component tolerances, EMI and physical overload survival are not established by this model.',
 'Systematic effects can align rather than cancel. The failed thermal corners remain visible; a small white-noise standard deviation does not guarantee a small total error.',
 'Zero tracking still depends on a trustworthy external empty-pan indication. The simulation shows why a false confirmation can remove a genuine small load.'
 ],[('selected-operating-envelope.json','Selected measurement envelope'),('results/virtual-selected-envelope.csv','16 verified selected-envelope corners'),('figures/virtual-instrument.svg','Thermal screening and synthetic interference after sampling'),('results/virtual-temperature-corners.csv','All 24 temperature/load corners'),('results/virtual-calibration-rejection.csv','Calibration contamination cases'),('results/virtual-correlated-uncertainty.json','Correlated uncertainty calculation'),('results/virtual-aliasing-summary.csv','Analog filtering and aliasing'),('results/virtual-aliasing-trace.csv','Synthetic interference trace')],[dict(label='Thermal/load/direction corners',value=24),dict(label='Correlated standard uncertainty at 5 kg',value=round(correlated,3),unit='g'),dict(label='Shared-temperature study at 5 kg, k = 2',value=round(2*correlated,3),unit='g')])
if __name__=='__main__':main()