"""Original duplicate-unit readings and the same-unit averaging variance floor."""
from pathlib import Path
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np


def main():
    readings = np.array([[99.4,99.6],[100.0,100.2],[99.8,100.0],[100.6,100.8]])
    m,q = readings.shape
    means,grand = readings.mean(axis=1),readings.mean()
    msw = np.sum((readings-means[:,None])**2)/(m*(q-1))
    msb = q*np.sum((means-grand)**2)/(m-1)
    between = max(0.,(msb-msw)/q)
    assert np.allclose([grand,msw,msb,between],[100.05,.02,.5,.24])
    # Independently check the leave-one-out covariance derivation in the article.
    n,sigma = 4,2.
    transform = (n*np.eye(n)-np.ones((n,n)))/(n-1)
    covariance = sigma**2*transform@transform.T
    assert np.allclose(np.diag(covariance), n*sigma**2/(n-1))
    assert np.allclose(covariance[~np.eye(n,dtype=bool)],-n*sigma**2/(n-1)**2)
    assert np.allclose(covariance@np.ones(n),0.)
    repeats = np.arange(1,21)
    scale = .2
    total = np.sqrt(between+msw/repeats)/scale
    analytical_only = np.sqrt(msw/repeats)/scale
    floor = np.sqrt(between)/scale
    fig,(left,right) = plt.subplots(1,2,figsize=(10,4.1),dpi=100,facecolor='white')
    fig.subplots_adjust(left=.065,right=.98,bottom=.20,top=.76,wspace=.34)
    for j,pair in enumerate(readings):
        left.plot(pair,[j+1,j+1],color='#7e8d9b',lw=1.6)
        left.scatter(pair,[j+1,j+1],s=38,color='#196b93',zorder=3)
    left.scatter(means,np.arange(1,m+1),marker='|',s=180,color='black',label='Unit mean',zorder=4)
    left.axvline(grand,color='#a44e29',ls='--',lw=1.4,label='Grand mean: 100.05')
    left.set(xlim=(99.25,100.9),ylim=(.5,4.5),yticks=[1,2,3,4],xlabel='Concentration (mg/L)',ylabel='Supplied unit')
    left.set_title('Two test portions from each unit',fontsize=11)
    left.legend(fontsize=8,loc='lower right',framealpha=1)
    right.plot(repeats,total,color='#196b93',lw=2.3,label='Item + analytical variation')
    right.axhline(floor,color='#a44e29',ls='--',lw=1.5,label=f'Item floor: {floor:.2f}')
    right.plot(repeats,analytical_only,color='#43894a',lw=1.8,label='Analytical variation alone')
    right.set(xlim=(1,20),ylim=(0,2.85),xticks=[1,5,10,15,20],xlabel='Analyses averaged from the same unit',ylabel='Score standard deviation')
    right.set_title('Assessment scale: 0.20 mg/L',fontsize=11)
    right.legend(fontsize=8,loc='center right',framealpha=1)
    for ax in [left,right]:ax.grid(alpha=.20)
    fig.suptitle('Item heterogeneity persists under repeated analysis',fontsize=14,y=.97)
    fig.text(.5,.84,'Original illustrative data; between-unit variance 0.24 and within-unit variance 0.02 (mg/L)²',ha='center',fontsize=10)
    fig.text(.5,.035,'Conditional model: independent analytical repeats share one unchanged unit effect.',ha='center',fontsize=10)
    fig.savefig(Path.cwd()/'proficiency-item-heterogeneity-and-averaging.png',dpi=100,facecolor='white')
    plt.close(fig)


if __name__=='__main__':
    main()