"""Owned finite linear-kernel GP joint covariance; exact Fraction, no optimizer."""
from fractions import Fraction as F
import json
def dot(a,b):return sum(x*y for x,y in zip(a,b))
def kernel(x,z):return 1+x*z
def solve(A,b):
    n=len(A)
    if n not in (1,2) or any(len(r)!=n for r in A) or len(b)!=n:raise ValueError('one/two training points subset')
    M=[r[:]+[v] for r,v in zip(A,b)]
    for j in range(n):
        pivot=next((i for i in range(j,n) if M[i][j]!=0),None)
        if pivot is None:raise ValueError('singular ordinary linear solve; not a general GP invalidity verdict')
        M[j],M[pivot]=M[pivot],M[j];scale=M[j][j];M[j]=[x/scale for x in M[j]]
        for i in range(n):
            if i!=j:
                scale=M[i][j];M[i]=[x-scale*y for x,y in zip(M[i],M[j])]
    x=[r[-1] for r in M];assert [dot(row,x) for row in A]==b;return x
def posterior(X,y,noise_variance,queries,future_noise_variance):
    if noise_variance<0 or future_noise_variance<0 or len(queries)!=2:raise ValueError('variance and two-query contract')
    K=[[kernel(x,z)+(noise_variance if i==j else 0) for j,z in enumerate(X)] for i,x in enumerate(X)]
    cross=[[kernel(q,x) for x in X] for q in queries];prior=[[kernel(a,b) for b in queries] for a in queries]
    alpha=solve(K,y);solved=[solve(K,row) for row in cross];mean=[dot(row,alpha) for row in cross]
    C=[[prior[i][j]-dot(cross[i],solved[j]) for j in range(2)] for i in range(2)]
    Cy=[[C[i][j]+(future_noise_variance if i==j else 0) for j in range(2)] for i in range(2)]
    assert C[0][1]==C[1][0] and C[0][0]>=0 and C[1][1]>=0 and C[0][0]*C[1][1]-C[0][1]**2>=0
    return {'train_X':X,'train_y':y,'training_noise_variance':noise_variance,'query_X':queries,'K_train_with_noise':K,'cross_query_train':cross,'K_query_prior':prior,'alpha_solved_from_K_y':alpha,'solves_K_cross_columns':solved,'posterior_mean':mean,'latent_joint_covariance':C,'independent_future_measurement_noise_variance':future_noise_variance,'future_observation_joint_covariance':Cy,'latent_sum_variance_full':sum(map(sum,C)),'latent_sum_variance_if_independence_wrongly_assumed':C[0][0]+C[1][1]}
def demo():
    q=[F(0),F(1)];one=posterior([F(0)],[F(1)],F(1),q,F(1));changed=posterior([F(0)],[F(3)],F(1),q,F(1))
    assert one['posterior_mean']==[F(1,2),F(1,2)] and one['latent_joint_covariance']==[[F(1,2),F(1,2)],[F(1,2),F(3,2)]]
    assert one['future_observation_joint_covariance']==[[F(3,2),F(1,2)],[F(1,2),F(5,2)]]
    assert one['latent_sum_variance_full']==3 and one['latent_sum_variance_if_independence_wrongly_assumed']==2
    # Independent coefficient posterior: prior(a,b)=N(0,I), observe a+noise with variance1.
    coefficient_cov=[[F(1,2),F(0)],[F(0),F(1)]];Phi=[[F(1),x] for x in q]
    independent_cov=[[sum(Phi[i][a]*coefficient_cov[a][b]*Phi[j][b] for a in range(2) for b in range(2)) for j in range(2)] for i in range(2)]
    assert independent_cov==one['latent_joint_covariance'] and changed['latent_joint_covariance']==one['latent_joint_covariance'] and changed['posterior_mean']==[F(3,2),F(3,2)]
    repeated=posterior([F(0),F(0)],[F(1),F(3)],F(1),q,F(1));assert repeated['posterior_mean']==[F(4,3),F(4,3)]
    try:posterior([F(0),F(0)],[F(1),F(3)],F(0),q,F(1))
    except ValueError as e:singular=str(e)
    else:raise AssertionError('singular ordinary solve accepted')
    return {'kernel':'1+x*z, prior f(x)=a+b*x with independent standard-normal a,b','one_training_point':one,'independent_coefficient_posterior_covariance_check':independent_cov,'changed_y_fixed_kernel_and_noise':changed,'repeated_X_positive_noise':repeated,'repeated_X_zero_noise_ordinary_solve_rejected':singular,'hyperparameters_or_noise_estimated':False,'acquisition_optimization_or_sklearn_executed':False}
if __name__=='__main__':print(json.dumps(demo(),default=str,indent=2))