"""Original synthetic illustration, not a plagiarism classifier. Run from project root."""
from pathlib import Path
import json,csv
import numpy as np
from scipy.signal import stft
from scipy.io.wavfile import write
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
OUT=Path('dist/assets/music');OUT.mkdir(parents=True,exist_ok=True)
SR=8000
notes=np.array([60,64,67,69,67,64,62,60, 65,66,70,63,68,61,66,70, 72,71,67,64,62,64,67,60])
def synth(seq,dur=.5,shift=0,variant=False):
 parts=[]
 for m in seq:
  t=np.arange(round(SR*dur))/SR;f=440*2**((m+shift-69)/12)
  harmonics=[(1,1),(2,.35),(3,.14)] if not variant else [(1,1),(3,.42),(5,.17)]
  y=sum(a*np.sin(2*np.pi*f*k*t+.04*k) for k,a in harmonics)
  env=np.minimum(t/.025,1)*np.minimum((dur-t)/.055,1)
  parts.append(.48*y*env/max(1,sum(a for _,a in harmonics)))
 return np.concatenate(parts)
a=synth(notes)
b=np.concatenate([synth([53,56,58,54],variant=True),synth(notes[:8],.625,3,True),synth([58,61,55,56],variant=True),synth(notes[16:],.5/1.2,3,True),synth([54,57,59,53],variant=True)])
for name,y in [('a',a),('b',b)]:write(OUT/f'synthetic-{name}.wav',SR,np.int16(y*32767))
def features(y):
 f,t,z=stft(y,SR,nperseg=4096,noverlap=3896,boundary=None);mag=abs(z);valid=(f>=100)&(f<=1800)
 midi=np.rint(69+12*np.log2(f[valid]/440)).astype(int);c=np.zeros((12,len(t)))
 for k,pc in enumerate(midi%12):c[pc]+=mag[valid][k]
 c/=np.maximum(np.linalg.norm(c,axis=0,keepdims=True),1e-12)
 return f,t,mag,c
fa,ta,ma,ca=features(a);fb,tb,mb,cb=features(b);cb=np.roll(cb,-3,axis=0)
regions=[('R1',0,4,2,7,'#d47827'),('R2',8,12,9,9+4/1.2,'#198b98')]
rows=[];paths=[]
def dtw(x,y):
 cost=np.linalg.norm(x.T[:,None,:]-y.T[None,:,:],axis=2);n,m=cost.shape;d=np.full((n+1,m+1),np.inf);d[0,0]=0
 for i in range(1,n+1):
  for j in range(1,m+1):d[i,j]=cost[i-1,j-1]+min(d[i-1,j-1],d[i-1,j],d[i,j-1])
 i,j=n,m;p=[]
 while i>0 and j>0:
  p.append((i-1,j-1));s=int(np.argmin([d[i-1,j-1],d[i-1,j],d[i,j-1]]));i-=s in [0,1];j-=s in [0,2]
 p=np.array(p[::-1]);return cost,p
for rid,s,e,u,v,color in regions:
 ia=np.where((ta>s+.2)&(ta<e-.2))[0];ib=np.where((tb>u+.2)&(tb<v-.2))[0];cost,p=dtw(ca[:,ia],cb[:,ib]);dist=cost[p[:,0],p[:,1]];mean=float(dist.mean());score=100*(1-mean/np.sqrt(2));paths.append((ta[ia[p[:,0]]],tb[ib[p[:,1]]],color))
 rows.append(dict(region=rid,a_start=s,a_end=e,b_start=u,b_end=v,pitch_correction_semitones=-3,b_duration_over_a=(v-u)/(e-s),mean_euclidean=mean,feature_proximity_percent=score,path_pairs=len(p)))
 with (OUT/f'{rid.lower()}-alignment.csv').open('w') as f:
  writer=csv.writer(f);writer.writerow(['a_seconds','b_seconds','euclidean']);writer.writerows(zip(ta[ia[p[:,0]]],tb[ib[p[:,1]]],dist))
report={'kind':'synthetic demonstration; known inserted windows; DTW inside those windows','sample_rate':SR,'stft_samples':4096,'hop_samples':200,'pitch_correction_B':-3,'score':'100*(1-mean_path_distance/sqrt(2)); not calibrated probability','coverage_a':8/12,'coverage_b':sum(r['b_end']-r['b_start'] for r in rows)/(len(b)/SR),'regions':rows}
(OUT/'measurements.json').write_text(json.dumps(report,indent=2))
with (OUT/'regions.csv').open('w') as f:
 w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False,'figure.facecolor':'#f7f6ef','axes.facecolor':'#f7f6ef'})
fig,axes=plt.subplots(4,1,figsize=(12,10),layout='constrained')
for idx,(y,t,mag,label) in enumerate([(a,ta,ma,'A · original synthetic score'),(b,tb,mb,'B · different harmonics; inserted phrases +3 semitones')]):
 ax=axes[idx*2];ax.plot(np.arange(len(y))[::20]/SR,y[::20],color='#334b5a',lw=.4);ax.set(ylabel='Amplitude',title=label,xlim=(0,len(y)/SR),ylim=(-.7,.7))
 sx=axes[idx*2+1];sx.pcolormesh(t,fa,20*np.log10(np.maximum(mag,1e-6)),cmap='magma',vmin=-75,vmax=-15,shading='auto',rasterized=True);sx.set(ylim=(100,1800),ylabel='Frequency (Hz)',xlabel='Time (seconds)',xlim=(0,len(y)/SR))
 for row,reg in zip(rows,regions):
  start,end=(reg[1],reg[2]) if idx==0 else (reg[3],reg[4]);color=reg[5]
  ax.add_patch(Rectangle((start,-.67),end-start,1.34,fill=False,lw=2.5,edgecolor=color));ax.text((start+end)/2,.47,f"{reg[0]} · {row['feature_proximity_percent']:.1f}%",ha='center',color=color,weight='bold',bbox={'facecolor':'#f7f6ef','edgecolor':'none','pad':2})
  sx.add_patch(Rectangle((start,110),end-start,1670,fill=False,lw=2.5,edgecolor=color))
fig.suptitle('Known inserted regions · percentages measure chroma proximity after correction, not plagiarism',fontsize=13);fig.savefig(OUT/'paired-regions.png',dpi=170);plt.close(fig)
fig,axes=plt.subplots(1,2,figsize=(12,4.8),layout='constrained')
# Chroma histogram is an average profile, not amplitude bins or an identification proof.
h1=ca.mean(axis=1);h2=cb.mean(axis=1);h1/=h1.sum();h2/=h2.sum();x=np.arange(12);axes[0].bar(x-.18,h1,width=.36,label='A',color='#d47827');axes[0].bar(x+.18,h2,width=.36,label='B corrected -3 st',color='#198b98');axes[0].set(xticks=x,xticklabels=['C','C#','D','D#','E','F','F#','G','G#','A','A#','B'],ylabel='Share of mean chroma energy',title='Pitch-class histogram · temporal order discarded');axes[0].legend(fontsize=8)
for row,reg in zip(rows,regions):
 ia=(ta>reg[1]+.2)&(ta<reg[2]-.2);ib=(tb>reg[3]+.2)&(tb<reg[4]-.2);c,p=dtw(ca[:,ia],cb[:,ib]);axes[1].hist(c[p[:,0],p[:,1]],bins=np.linspace(0,np.sqrt(2),30),alpha=.5,label=f"{reg[0]} · mean {row['mean_euclidean']:.3f}",color=reg[5])
axes[1].set(xlabel='Euclidean distance · nonnegative unit chroma vectors',ylabel='Aligned frame pairs',title='Distance histograms inside the colored regions');axes[1].legend();fig.savefig(OUT/'histograms.png',dpi=170);plt.close(fig)
fig,ax=plt.subplots(figsize=(10,6),layout='constrained');cost=np.linalg.norm(ca.T[:,None,:]-cb.T[None,:,:],axis=2);im=ax.imshow(cost,origin='lower',aspect='auto',extent=[tb[0],tb[-1],ta[0],ta[-1]],cmap='viridis_r',vmin=0,vmax=np.sqrt(2))
for xa,xb,color in paths:ax.plot(xb,xa,color=color,lw=2)
ax.set(xlabel='B time (seconds)',ylabel='A time (seconds)',title='Cross-distance matrix · local DTW paths in known inserted windows');fig.colorbar(im,ax=ax,label='Euclidean chroma distance');fig.savefig(OUT/'distance-matrix.png',dpi=170);plt.close(fig)
print(json.dumps(report,indent=2))