"""Case 009: normal-ML covariance CFA, NumPy + SciPy.
Controlled covariance matrix; nominal N=400, not sampled observations.
Compact positive-loading example, not a general SEM package.
"""
import numpy as np
import scipy
from scipy.optimize import minimize
S=np.full((6,6),.4)
S[:3,:3]=1
S[3:,3:]=1
np.fill_diagonal(S,2)
N=400
logdetS=np.linalg.slogdet(S)[1]
def implied(p,two):
    loadings=np.exp(p[:6])
    residual=np.exp(p[6:12])
    if two:
        L=np.zeros((6,2))
        L[:3,0]=loadings[:3]
        L[3:,1]=loadings[3:]
        rho=np.tanh(p[12])
        phi=np.array([[1,rho],[rho,1]])
    else:
        L=loadings[:,None]
        phi=np.ones((1,1))
    return L@phi@L.T+np.diag(residual)
def discrepancy(p,two):
    sigma=implied(p,two)
    return np.linalg.slogdet(sigma)[1]+np.trace(np.linalg.solve(sigma,S))-logdetS-6
baseline=N*(np.log(np.diag(S)).sum()-logdetS)
results=[]
for two in [False,True]:
    start=np.r_[np.log(np.repeat(.9,6)),np.zeros(6)]
    if two: start=np.r_[start,np.arctanh(.3)]
    starts=[start]
    if not two:
        starts += [np.r_[np.log([1,1,1,.6,.6,.6]),np.zeros(6)],
                   np.r_[np.log([.6,.6,.6,1,1,1]),np.zeros(6)]]
    candidates=[minimize(discrepancy,initial,args=(two,),method="L-BFGS-B",jac="3-point",
        options=dict(ftol=1e-14,gtol=1e-8,maxiter=2000,maxls=50)) for initial in starts]
    fit=min(candidates,key=lambda result:result.fun)
    if not two:
        print("One-factor multistart discrepancies:",np.round([result.fun for result in candidates],9))
    sigma=implied(fit.x,two)
    chi=max(N*discrepancy(fit.x,two),0)
    df=21-len(start)
    cfi=1-max(chi-df,0)/max(chi-df,baseline-15,0)
    rmsea=np.sqrt(max(chi/df-1,0)/N)
    residual=(S-sigma)/np.sqrt(np.outer(np.diag(S),np.diag(S)))
    srmr=np.sqrt(np.mean(residual[np.tril_indices(6)]**2))
    print("two_factor" if two else "one_factor", "chisq df CFI RMSEA SRMR:",
          np.round([chi,df,cfi,rmsea,srmr],6))
    print("Optimizer success:",fit.success,"max numerical gradient:",np.max(np.abs(fit.jac)))
    assert fit.success and np.max(np.abs(fit.jac))<1e-5
    assert np.linalg.eigvalsh(sigma).min()>0
    results.append((chi,sigma))
print("One-factor raw covariance residuals:\n",np.round(S-results[0][1],6))
assert results[0][0]>100 and results[1][0]<1e-6
assert np.max(np.abs(S-results[1][1]))<1e-6
print("All assertions passed; NumPy",np.__version__,"SciPy",scipy.__version__)
print("Does not calculate lavaan modification indices or robust fit statistics.")