"""
Yakaboylu Weight Operator - Dirichlet Boundary Condition Projection
===================================================================

The weight function omega'(z,y) was shown to be indefinite on
the full L^2 space (6-23 negative eigenvalues depending on grid).

Escape route (a) from the original paper: Yakaboylu's construction
restricts to eigenfunctions satisfying f(0) = 0 (Dirichlet BC at
the origin). If the negative eigenvectors of the Gram matrix violate
this BC, the indefiniteness is irrelevant.

Test: project the Gram matrix onto functions vanishing at z = 0
and check if positive-definiteness is restored.

Three projection methods:
1. Grid exclusion: remove the z = 0 point from the grid
2. Sine basis: use sin(n*pi*z/z_max) which all vanish at z = 0
3. Chebyshev with BC: Chebyshev polynomials modified to vanish at z = 0

Requirements: numpy, scipy
"""

import numpy as np
from scipy.linalg import eigvalsh
import json


def omega_prime(z, y, eps=1e-15):
    """Weight function from equation (18)."""
    z = np.maximum(z, eps)
    y = np.maximum(y, eps)
    
    ratio1 = np.minimum(y / z, 500)
    term1 = (1.0 / z) / (1.0 + np.exp(ratio1))
    
    ratio2 = np.minimum(z / y, 500)
    term2 = (1.0 / y) / (1.0 + np.exp(ratio2))
    
    return term1 + term2


def test_method_1_grid_exclusion():
    """
    Method 1: Uniform grid on [epsilon, z_max].
    By starting at epsilon > 0, all basis functions automatically
    have f(0) undefined/excluded. This is the simplest test but
    doesn't strictly enforce f(0) = 0.
    """
    print("\n--- Method 1: Grid Exclusion (start at epsilon) ---")
    
    configs = [
        (50, 5.0, 0.01),
        (100, 5.0, 0.01),
        (100, 10.0, 0.01),
        (200, 10.0, 0.01),
        (50, 5.0, 0.001),
        (100, 5.0, 0.001),
        (200, 10.0, 0.001),
    ]
    
    for n, z_max, eps_start in configs:
        dz = (z_max - eps_start) / n
        z = np.linspace(eps_start + dz/2, z_max - dz/2, n)
        
        Z, Y = np.meshgrid(z, z)
        G = omega_prime(Z, Y) * dz * dz
        G = (G + G.T) / 2
        
        eigs = eigvalsh(G)
        n_neg = np.sum(eigs < -1e-14)
        min_e = eigs[0]
        
        status = "POSITIVE-DEFINITE" if n_neg == 0 else f"{n_neg} NEGATIVE"
        print(f"  N={n:>4d}, z_max={z_max:>4.0f}, eps={eps_start:.3f}: min_eig={min_e:>12.4e}  {status}")


def test_method_2_sine_basis():
    """
    Method 2: Project onto sine basis sin(n*pi*z/z_max).
    All basis functions vanish at z = 0 AND at z = z_max.
    This exactly enforces the Dirichlet BC.
    
    G_mn = integral_0^{z_max} integral_0^{z_max} 
           omega'(z,y) sin(m*pi*z/L) sin(n*pi*y/L) dz dy
    
    Computed by numerical quadrature on a fine grid.
    """
    print("\n--- Method 2: Sine Basis (exact Dirichlet BC) ---")
    
    configs = [
        (10, 5.0, 200),
        (20, 5.0, 200),
        (30, 5.0, 300),
        (10, 10.0, 200),
        (20, 10.0, 300),
        (30, 10.0, 400),
        (40, 10.0, 500),
        (50, 10.0, 600),
    ]
    
    for n_basis, z_max, n_quad in configs:
        # Quadrature grid (avoid z = 0)
        dz = z_max / n_quad
        z_quad = np.linspace(dz/2, z_max - dz/2, n_quad)
        
        # Build basis matrix: phi[q, m] = sin((m+1)*pi*z_q/z_max)
        phi = np.zeros((n_quad, n_basis))
        for m in range(n_basis):
            phi[:, m] = np.sin((m + 1) * np.pi * z_quad / z_max)
        
        # Build weight matrix on quadrature grid
        Z, Y = np.meshgrid(z_quad, z_quad)
        W = omega_prime(Z, Y) * dz * dz
        
        # Project: G_mn = phi^T W phi
        G = phi.T @ W @ phi
        G = (G + G.T) / 2
        
        eigs = eigvalsh(G)
        n_neg = np.sum(eigs < -1e-14)
        min_e = eigs[0]
        max_e = eigs[-1]
        
        status = "POSITIVE-DEFINITE" if n_neg == 0 else f"{n_neg} NEGATIVE"
        print(f"  basis={n_basis:>3d}, z_max={z_max:>4.0f}, quad={n_quad:>4d}: "
              f"min={min_e:>12.4e}, max={max_e:>10.4e}  {status}")


def test_method_3_vanishing_polynomials():
    """
    Method 3: Use basis functions z * L_n(z) where L_n are Laguerre
    polynomials. The factor z ensures f(0) = 0. The Laguerre functions
    provide good spectral resolution on [0, infinity).
    
    This is the natural basis for the Yakaboylu Hamiltonian because
    the eigenfunctions are of the form z^s * (confluent hypergeometric)
    where s >= 0 ensures vanishing at z = 0.
    """
    print("\n--- Method 3: z * Laguerre Basis (natural BC) ---")
    
    from numpy.polynomial.laguerre import lagval
    
    configs = [
        (10, 10.0, 200),
        (20, 10.0, 300),
        (30, 10.0, 400),
        (20, 20.0, 400),
        (30, 20.0, 500),
        (40, 20.0, 600),
        (50, 20.0, 800),
    ]
    
    for n_basis, z_max, n_quad in configs:
        dz = z_max / n_quad
        z_quad = np.linspace(dz/2, z_max - dz/2, n_quad)
        
        # Basis: phi_m(z) = z * L_m(z) * exp(-z/2) (normalized Laguerre)
        phi = np.zeros((n_quad, n_basis))
        for m in range(n_basis):
            coeffs = np.zeros(m + 1)
            coeffs[m] = 1.0
            L_m = lagval(z_quad, coeffs)
            phi[:, m] = z_quad * L_m * np.exp(-z_quad / 2)
        
        # Normalize each basis function
        for m in range(n_basis):
            norm = np.sqrt(np.sum(phi[:, m]**2) * dz)
            if norm > 1e-15:
                phi[:, m] /= norm
        
        # Weight matrix
        Z, Y = np.meshgrid(z_quad, z_quad)
        W = omega_prime(Z, Y) * dz * dz
        
        # Project
        G = phi.T @ W @ phi
        G = (G + G.T) / 2
        
        eigs = eigvalsh(G)
        n_neg = np.sum(eigs < -1e-14)
        min_e = eigs[0]
        
        status = "POSITIVE-DEFINITE" if n_neg == 0 else f"{n_neg} NEGATIVE"
        print(f"  basis={n_basis:>3d}, z_max={z_max:>4.0f}, quad={n_quad:>4d}: "
              f"min={min_e:>12.4e}  {status}")


def test_negative_eigenvector_bc():
    """
    Direct test: compute the negative eigenvectors of the full
    Gram matrix and check whether they satisfy f(0) = 0.
    
    If the negative eigenvectors have f(0) != 0, the BC projection
    removes them and the construction might survive.
    If they have f(0) = 0, they're in the domain and the construction fails.
    """
    print("\n--- Direct Test: Do Negative Eigenvectors Satisfy f(0) = 0? ---")
    
    from scipy.linalg import eigh
    
    for n in [50, 100, 200]:
        z_max = 10.0
        dz = z_max / n
        z = np.linspace(dz/2, z_max - dz/2, n)
        
        Z, Y = np.meshgrid(z, z)
        G = omega_prime(Z, Y) * dz * dz
        G = (G + G.T) / 2
        
        eigs, vecs = eigh(G)
        
        n_neg = np.sum(eigs < -1e-14)
        
        print(f"\n  N={n}, z_max={z_max}")
        print(f"  Negative eigenvalues: {n_neg}")
        
        if n_neg > 0:
            print(f"  {'Eig idx':>8s}  {'Eigenvalue':>14s}  {'v[0] (z~0)':>12s}  {'v[1]':>12s}  {'v[N-1] (z~L)':>12s}  {'|v[0]|/max':>10s}")
            
            for k in range(min(n_neg, 10)):
                ev = eigs[k]
                vec = vecs[:, k]
                v0 = vec[0]        # value near z = 0
                v1 = vec[1]        # next point
                vend = vec[-1]     # value near z = z_max
                vmax = np.max(np.abs(vec))
                ratio = abs(v0) / vmax if vmax > 0 else 0
                
                vanishes = "YES (< 0.01)" if ratio < 0.01 else "NO"
                print(f"  {k:>8d}  {ev:>14.6e}  {v0:>12.6f}  {v1:>12.6f}  {vend:>12.6f}  {ratio:>10.4f}  {vanishes}")


def main():
    print("=" * 70)
    print("YAKABOYLU WEIGHT OPERATOR - DIRICHLET BC PROJECTION")
    print("Does f(0) = 0 eliminate the negative eigenvalues?")
    print("=" * 70)
    
    # First: check if negative eigenvectors vanish at z = 0
    test_negative_eigenvector_bc()
    
    # Method 1: Grid exclusion
    test_method_1_grid_exclusion()
    
    # Method 2: Sine basis (exact BC)
    test_method_2_sine_basis()
    
    # Method 3: Laguerre basis (natural BC)
    test_method_3_vanishing_polynomials()
    
    print(f"\n{'='*70}")
    print("INTERPRETATION")
    print(f"{'='*70}")
    print("If negative eigenvectors have f(0) != 0 AND projected")
    print("  Gram matrix is positive-definite:")
    print("  -> The BC projection saves the construction.")
    print("  -> Yakaboylu's self-adjointness may hold on the domain.")
    print()
    print("If negative eigenvectors already have f(0) ~ 0 OR projected")
    print("  Gram matrix is still indefinite:")
    print("  -> The construction fails on its own terms.")
    print("  -> The BC does not help. Escape route (a) is closed.")
    print(f"{'='*70}")


if __name__ == '__main__':
    main()