"""
Yakaboylu Weight Operator - Positive Definiteness Test

From equation (18) of the paper:
  omega'(z,y) = (1/z)/(1 + exp(y/z)) + (1/y)/(1 + exp(z/y))

The self-adjointness of the transformed Hamiltonian H' requires
this weight function to define a positive-definite inner product:

  <u, u> = integral integral omega'(z,y) u'*(y) u'(z) dz dy > 0

for all nonzero u' in the domain.

Test: discretize omega'(z,y) on a grid, form the Gram matrix,
check if all eigenvalues are positive.

If NOT positive-definite: the self-adjointness argument fails.
If positive-definite: numerical evidence it holds (on the grid).
"""

import numpy as np
from scipy.linalg import eigvalsh

def omega_prime(z, y, eps=1e-15):
    """Weight function from equation (18)."""
    # Avoid division by zero and overflow
    z = np.maximum(z, eps)
    y = np.maximum(y, eps)
    
    # Compute (1/z)/(1 + exp(y/z))
    ratio1 = y / z
    # Clip to avoid overflow in exp
    ratio1 = np.minimum(ratio1, 500)
    term1 = (1.0 / z) / (1.0 + np.exp(ratio1))
    
    # Compute (1/y)/(1 + exp(z/y))
    ratio2 = z / y
    ratio2 = np.minimum(ratio2, 500)
    term2 = (1.0 / y) / (1.0 + np.exp(ratio2))
    
    return term1 + term2


def build_gram_matrix(n_points, z_max, method='midpoint'):
    """
    Discretize the weight function on a grid and build the Gram matrix.
    
    G[i,j] = omega'(z_i, z_j) * dz * dy
    
    The inner product is:
    <u, v> = sum_{i,j} G[i,j] * u_i^* * v_j
    
    This is positive-definite iff all eigenvalues of G are positive.
    """
    if method == 'midpoint':
        dz = z_max / n_points
        z = np.linspace(dz/2, z_max - dz/2, n_points)
    elif method == 'log':
        # Log-spaced grid to better capture behavior near z=0
        z = np.logspace(-3, np.log10(z_max), n_points)
        dz = np.diff(z, prepend=0)
    
    # Build Gram matrix
    Z, Y = np.meshgrid(z, z)
    G = omega_prime(Z, Y)
    
    if method == 'midpoint':
        G *= dz * dz  # uniform grid spacing
    elif method == 'log':
        # Scale by grid spacing for non-uniform grid
        DZ, DY = np.meshgrid(dz, dz)
        G *= DZ * DY
    
    return G, z


def test_symmetry(G):
    """Check if G is symmetric (required by condition 14)."""
    asymm = np.max(np.abs(G - G.T))
    return asymm


def main():
    print("=" * 70)
    print("YAKABOYLU WEIGHT OPERATOR - POSITIVE DEFINITENESS TEST")
    print("Testing omega'(z,y) = (1/z)/(1+exp(y/z)) + (1/y)/(1+exp(z/y))")
    print("=" * 70)
    
    # First verify the symmetry property (condition 14):
    # omega'(z,y) = omega'*(y,z)
    # Since the function is real, this means omega'(z,y) = omega'(y,z)
    print("\n--- Symmetry Check ---")
    z_test = np.array([0.1, 0.5, 1.0, 2.0, 5.0])
    y_test = np.array([0.3, 0.7, 1.5, 3.0, 7.0])
    for zi in z_test:
        for yi in y_test:
            diff = abs(omega_prime(zi, yi) - omega_prime(yi, zi))
            if diff > 1e-10:
                print(f"  ASYMMETRIC at z={zi}, y={yi}: diff={diff:.2e}")
    print("  omega'(z,y) = omega'(y,z) verified (function is symmetric)")
    
    # Test on various grid sizes and domains
    configs = [
        (30, 5.0, 'midpoint'),
        (50, 5.0, 'midpoint'),
        (100, 5.0, 'midpoint'),
        (50, 10.0, 'midpoint'),
        (50, 20.0, 'midpoint'),
        (50, 5.0, 'log'),
        (100, 10.0, 'log'),
    ]
    
    for n, z_max, method in configs:
        print(f"\n--- Grid: {n} points, z_max={z_max}, method={method} ---")
        
        G, z = build_gram_matrix(n, z_max, method)
        
        # Symmetry of discretized matrix
        asymm = test_symmetry(G)
        print(f"  Matrix asymmetry: {asymm:.2e}")
        
        # Symmetrize (numerical noise)
        G = (G + G.T) / 2
        
        # Compute eigenvalues
        eigs = eigvalsh(G)
        eigs_sorted = np.sort(eigs)
        
        n_negative = np.sum(eigs < -1e-14)
        n_zero = np.sum(np.abs(eigs) < 1e-14)
        n_positive = np.sum(eigs > 1e-14)
        min_eig = eigs_sorted[0]
        max_eig = eigs_sorted[-1]
        
        print(f"  Eigenvalues: {n_negative} negative, {n_zero} zero, {n_positive} positive")
        print(f"  Min eigenvalue: {min_eig:.6e}")
        print(f"  Max eigenvalue: {max_eig:.6e}")
        print(f"  Condition number: {max_eig / max(abs(min_eig), 1e-30):.2e}")
        
        if n_negative > 0:
            print(f"  *** NOT POSITIVE-DEFINITE ***")
            print(f"  Smallest 5 eigenvalues: {eigs_sorted[:5]}")
        else:
            print(f"  POSITIVE-DEFINITE on this grid")
    
    # Detailed analysis: how eigenvalues behave as grid refines
    print(f"\n{'=' * 70}")
    print("CONVERGENCE: min eigenvalue as grid refines (z_max=10, midpoint)")
    print(f"{'=' * 70}")
    for n in [20, 40, 60, 80, 100, 150, 200, 300]:
        G, z = build_gram_matrix(n, 10.0, 'midpoint')
        G = (G + G.T) / 2
        eigs = eigvalsh(G)
        min_e = np.min(eigs)
        n_neg = np.sum(eigs < -1e-14)
        print(f"  N={n:>4}: min_eig = {min_e:>14.6e}  negative_count = {n_neg}")
    
    print(f"\n{'=' * 70}")
    print("INTERPRETATION:")
    print("If positive-definite at all grid sizes and the minimum eigenvalue")
    print("stays positive as the grid refines: strong evidence that the")
    print("weight operator defines a valid inner product, supporting")
    print("Yakaboylu's self-adjointness claim.")
    print("")
    print("If negative eigenvalues appear or min eigenvalue trends toward")
    print("zero/negative as grid refines: the weight operator fails to")
    print("define a positive-definite inner product, and the self-adjointness")
    print("argument collapses. This would identify the exact failure mode.")
    print(f"{'=' * 70}")


if __name__ == "__main__":
    main()