"""Small independent preflight reference, not a browser Engine or production fit."""
import numpy as np
from numpy.polynomial import Polynomial


def lens_positions(separation, mass_ratio):
    masses = np.array([1.0, mass_ratio]) / (1.0 + mass_ratio)
    return masses, np.array([-masses[1] * separation, masses[0] * separation])


def source_position(image, masses, positions):
    image = np.asarray(image, dtype=complex)
    return image - sum(mass / (image.conjugate() - position)
                       for mass, position in zip(masses, positions))


def image_jacobian(image, masses, positions):
    derivative = sum(mass / (image - position) ** 2
                     for mass, position in zip(masses, positions))
    return 1.0 - abs(derivative) ** 2


def binary_images(source, separation, mass_ratio):
    """Quintic in the planet-centered frame; polish/check the original equation."""
    masses, positions = lens_positions(separation, mass_ratio)
    origin = positions[1]
    local_positions = positions - origin
    local_source = source - origin
    coordinate = Polynomial([0.0j, 1.0])
    denominators = [coordinate - position for position in local_positions]
    product = denominators[0] * denominators[1]
    numerator = local_source.conjugate() * product + masses[0] * denominators[1] + masses[1] * denominators[0]
    conjugate_denominators = [numerator - position * product for position in local_positions]
    polynomial = (local_source - coordinate) * conjugate_denominators[0] * conjugate_denominators[1]
    polynomial += product * (masses[0] * conjugate_denominators[1] + masses[1] * conjugate_denominators[0])
    accepted = []
    for candidate in polynomial.roots():
        for _ in range(12):
            residual = source_position(candidate, masses, local_positions) - local_source
            derivative = sum(mass / (candidate.conjugate() - position) ** 2
                             for mass, position in zip(masses, local_positions))
            determinant = 1.0 - abs(derivative) ** 2
            if abs(determinant) < 1e-13:
                break
            step = (-residual + derivative * residual.conjugate()) / determinant
            if abs(step) > 0.1:
                break
            candidate += step
            if abs(residual) < 1e-13:
                break
        residual = abs(source_position(candidate, masses, local_positions) - local_source)
        if residual < 1e-9 and not any(abs(candidate - old) < 1e-8 for old in accepted):
            accepted.append(candidate)
    images = np.array(accepted) + origin
    jacobians = image_jacobian(images, masses, positions)
    if len(images) not in (3, 5) or np.sum(jacobians < 0) - np.sum(jacobians > 0) != 1:
        raise ArithmeticError('Incomplete/critical binary roots; never substitute a clipped magnification')
    return images, jacobians


def finite_source_image_area(source, source_radius, separation, mass_ratio,
                             relative_width=0.01, maximum_depth=20, maximum_cells=3000000,
                             limb_gamma=0.0, use_curvature_bounds=False):
    """Image-plane inverse integration with conservative geometric cell bounds.

    The finite box is guaranteed to contain every preimage of the source disk.
    A Lipschitz bound certifies cells wholly inside/outside the source. Unresolved
    boundary cells retain their entire area in the upper bound. This is slow but
    independent of polynomial roots and contour reconnection.
    """
    masses, positions = lens_positions(separation, mass_ratio)
    farthest_lens = float(max(abs(positions)))
    source_extent = abs(source) + source_radius
    bound = (source_extent + farthest_lens + np.sqrt((source_extent - farthest_lens) ** 2 + 4.0)) / 2.0
    bound = float(np.nextafter(bound, np.inf))
    centers = np.array([0.0j])
    half_width = bound
    if not 0 <= limb_gamma <= 1:
        raise ValueError('The mean-normalized source intensity must remain nonnegative')
    lower_retired = upper_retired = midpoint_retired = 0.0
    examined = 0
    last = None
    for depth in range(maximum_depth + 1):
        examined += len(centers)
        radius = np.sqrt(2.0) * half_width
        with np.errstate(divide='ignore', invalid='ignore'):
            mapped = source_position(centers, masses, positions)
        minimum_distances = [np.hypot(np.maximum(abs(centers.real - position) - half_width, 0.0),
                                      np.maximum(abs(centers.imag) - half_width, 0.0)) for position in positions]
        with np.errstate(divide='ignore'):
            variation = radius * (1.0 + sum(mass / distance ** 2 for mass, distance in zip(masses, minimum_distances)))
        source_distances = abs(mapped - source)
        with np.errstate(invalid='ignore'):
            minimum_source_radius = np.maximum(source_distances - variation, 0.0)
        maximum_source_radius = source_distances + variation
        singular = ~np.isfinite(variation) | ~np.isfinite(source_distances)
        minimum_source_radius[singular] = 0.0
        maximum_source_radius[singular] = np.inf
        def intensity(radius_values):
            inside_disk = radius_values <= source_radius
            intensity_values = (1.0 - limb_gamma + 1.5 * limb_gamma
                                * np.sqrt(np.maximum(1.0 - (radius_values / source_radius) ** 2, 0.0)))
            return np.where(inside_disk, intensity_values, 0.0)
        minimum_intensity = intensity(maximum_source_radius)
        maximum_intensity = intensity(minimum_source_radius)
        midpoint_intensity = intensity(source_distances)
        if use_curvature_bounds and limb_gamma:
            # On a complete interior cell, the symmetric square's average
            # differs from its central value by at most H*h^2/3. H bounds the
            # spectral norm of the composed intensity Hessian throughout it.
            interior = maximum_source_radius < source_radius
            safe_radius = np.minimum(maximum_source_radius[interior], source_radius)
            mu = np.sqrt(1.0 - (safe_radius / source_radius) ** 2)
            radial_hessian = 1.5 * limb_gamma / (source_radius ** 2 * mu ** 3)
            radial_gradient = 1.5 * limb_gamma * safe_radius / (source_radius ** 2 * mu)
            mapping_gradient = 1.0 + sum(mass / distance[interior] ** 2
                                          for mass, distance in zip(masses, minimum_distances))
            mapping_hessian = 2.0 * sum(mass / distance[interior] ** 3
                                        for mass, distance in zip(masses, minimum_distances))
            hessian = radial_hessian * mapping_gradient ** 2 + radial_gradient * mapping_hessian
            average_error = hessian * half_width ** 2 / 3.0
            minimum_intensity[interior] = np.maximum(minimum_intensity[interior],
                                                     midpoint_intensity[interior] - average_error)
            maximum_intensity[interior] = np.minimum(maximum_intensity[interior],
                                                     midpoint_intensity[interior] + average_error)
        estimate_intensity = midpoint_intensity
        if use_curvature_bounds:
            # Positive tensor two-point Gauss rule improves the estimate only;
            # it does not replace or shrink the independent enclosure above.
            estimate_intensity = np.zeros_like(midpoint_intensity)
            for offset in (-1-1j, -1+1j, 1-1j, 1+1j):
                with np.errstate(divide='ignore', invalid='ignore'):
                    quadrature_source = source_position(centers + half_width / np.sqrt(3.0) * offset,
                                                        masses, positions)
                estimate_intensity += .25 * intensity(abs(quadrature_source - source))
        outside = maximum_intensity == 0
        # Interior limb-darkened cells retire only with a small intensity
        # interval; keep its nonzero lower/upper uncertainty in the totals.
        inside = ((maximum_source_radius < source_radius)
                  & (maximum_intensity - minimum_intensity <= .0025 * minimum_intensity))
        # Cells containing a lens can have a nonfinite central map. They remain
        # unresolved and are refined; a lens point itself has zero area.
        unresolved = ~(inside | outside)
        cell_area = 4.0 * half_width ** 2
        lower_retired += float(np.sum(minimum_intensity[inside])) * cell_area
        upper_retired += float(np.sum(maximum_intensity[inside])) * cell_area
        midpoint_retired += float(np.sum(estimate_intensity[inside])) * cell_area
        lower_area = lower_retired + float(np.sum(minimum_intensity[unresolved])) * cell_area
        upper_area = upper_retired + float(np.sum(maximum_intensity[unresolved])) * cell_area
        midpoint_area = midpoint_retired + float(np.sum(estimate_intensity[unresolved])) * cell_area
        source_area = np.pi * source_radius ** 2
        last = {'lower': lower_area / source_area, 'upper': upper_area / source_area,
                ('estimate' if use_curvature_bounds else 'midpoint'): midpoint_area / source_area, 'depth': depth,
                'examined_cells': examined, 'remaining_boundary_cells': int(np.count_nonzero(unresolved)),
                'box_half_width': bound, 'status': 'budget-limited',
                'estimate_rule': 'positive tensor Gauss2' if use_curvature_bounds else 'central midpoint',
                'bounds_rule': 'geometric range plus interior Hessian remainder' if use_curvature_bounds else 'geometric intensity range'}
        if lower_area > 0 and (upper_area - lower_area) / lower_area <= relative_width:
            last['status'] = 'bounded'
            break
        if examined + 4 * np.count_nonzero(unresolved) > maximum_cells:
            break
        half_width /= 2.0
        centers = (centers[unresolved, None] + half_width * np.array([-1-1j, -1+1j, 1-1j, 1+1j])).ravel()
    return last