"""Independent derivative and Gaussian-posterior checks; no Engine dependency."""
import json
from pathlib import Path
import numpy as np
from scipy.integrate import quad
from rv_kernel import latent_covariance, activity_covariance, linear_gaussian_posterior


def main():
    times = np.array([0., .31, 1.27, 3.83, 8.1])
    evolution, period, width = 17., 23.81, .47
    covariance, derivative, mixed = latent_covariance(times, times, evolution, period, width)
    delta = 1e-3
    kernel = lambda first, second: latent_covariance(first, second, evolution, period, width)[0]
    derivative_numeric = (kernel(times + delta, times) - kernel(times - delta, times)) / (2 * delta)
    mixed_numeric = (kernel(times + delta, times + delta) - kernel(times + delta, times - delta)
                     - kernel(times - delta, times + delta) + kernel(times - delta, times - delta)) / (4 * delta ** 2)
    derivative_error = float(np.max(np.abs(derivative - derivative_numeric)))
    mixed_error = float(np.max(np.abs(mixed - mixed_numeric)))
    assert derivative_error < 3e-8 and mixed_error < 3e-8
    assert np.max(np.abs(derivative + derivative.T)) < 1e-14
    assert np.min(np.linalg.eigvalsh(mixed)) > -1e-12
    parameters = dict(evolution_days=evolution, rotation_days=period, periodic_scale=width,
                      rv_coupling_m_s=8., rv_derivative_m_s_day=-45., calcium_amplitude_scaled=1.3)
    block = activity_covariance(times, times, parameters)
    assert np.max(np.abs(block - block.T)) < 1e-14
    assert np.min(np.linalg.eigvalsh(block)) > -1e-10
    cross = activity_covariance(times + .17, times, parameters)
    reverse = activity_covariance(times, times + .17, parameters)
    assert np.max(np.abs(cross - reverse.T)) < 1e-13
    values = np.array([1.1, .8, -.2]); design = np.array([[1.], [.4], [-.7]])
    noise = np.array([[.3,.04,0],[.04,.5,.03],[0,.03,.2]])
    log_density, mean, posterior_covariance = linear_gaussian_posterior(values, noise, design, [2.])
    inverse_noise = np.linalg.inv(noise)
    def density(coefficient):
        residual = values - design[:,0] * coefficient
        return np.exp(-.5 * (residual @ inverse_noise @ residual + coefficient ** 2 / 4))
    normalization = quad(density,-20,20,epsabs=1e-12)[0]
    integral_mean = quad(lambda coefficient: coefficient * density(coefficient),-20,20,epsabs=1e-12)[0] / normalization
    integral_variance = quad(lambda coefficient: (coefficient-integral_mean)**2 * density(coefficient),-20,20,epsabs=1e-12)[0] / normalization
    assert abs(mean[0]-integral_mean)<1e-12
    assert abs(posterior_covariance[0,0]-integral_variance)<1e-12
    result = dict(kind="independent finite-difference derivatives and numerical scalar posterior integral",
                  derivative_max_error=derivative_error, mixed_derivative_max_error=mixed_error,
                  block_min_eigenvalue=float(np.linalg.eigvalsh(block).min()),
                  posterior_mean_error=float(abs(mean[0]-integral_mean)),
                  posterior_variance_error=float(abs(posterior_covariance[0,0]-integral_variance)),
                  scalar_log_marginal_density=log_density)
    print(json.dumps(result,indent=2))
    return result


if __name__ == "__main__":
    main()