"""Alternating-tangent routines for binary binodals, spinodals,
and LCST estimates."""

import numpy as np
from scipy.optimize import brentq, minimize_scalar, minimize
from scipy.special import expit, logit  # pylint: disable=no-name-in-module

from .interface_feos import get_feos_npt_state


_EPS = 1e-290


def _sigmoid(zeta):
    return expit(zeta)


def _composition_from_w1(x1=None, w1=None):
    if x1 is not None:
        return [x1, 1.0 - x1], None
    if w1 is None:
        raise ValueError(
            "Provide either x1 (mole fraction) or w1 (mass fraction)."
        )
    return None, [w1, 1.0 - w1]


def spinodal_condition(temperature, pressure, params, x1=None, w1=None):
    """Return the spinodal condition value for a binary mixture."""
    xi_lst, wi_lst = _composition_from_w1(x1=x1, w1=w1)
    mix_state = get_feos_npt_state(
        temperature, pressure, params, xi_lst, wi_lst
    )
    return (mix_state.total_moles * mix_state.dln_phi_dnj())[1][0] - 1


def find_spinodal_roots(
    temperature, pressure, params, eps=_EPS, verbose=False
):
    """Find the two spinodal points (roots) with respect to w1 in (0, 1)."""

    def f(w1):
        return spinodal_condition(temperature, pressure, params, w1=w1)

    res = minimize_scalar(
        lambda w: -f(w), bounds=(0.0 + eps, 1.0 - eps), method="bounded"
    )

    if -res.fun > 0.0:
        return [
            brentq(f, 0.0, res.x, xtol=1e-12, rtol=1e-10),
            brentq(f, res.x, 1.0, xtol=1e-12, rtol=1e-10),
        ]

    if verbose:
        print("No spinodal found")
    return None


def find_lcst_temperature(
    temperature_low, temperature_high, pressure, params, xtol=1e-2
):
    """Find the LCST temperature where spinodal roots appear."""

    def has_spinodal(temperature):
        return (
            1.0 if find_spinodal_roots(temperature, pressure, params) else -1.0
        )

    return brentq(
        has_spinodal, temperature_low, temperature_high, xtol=xtol, rtol=1e-10
    )


# pylint: disable=too-many-arguments,too-many-positional-arguments


def find_binodal_by_alternating_tangents(
    temperature, pressure, params, max_iter=20, tol=1e-10, verbose=False
):
    """Find binodal points using the alternating tangents method."""
    w_spinodal = find_spinodal_roots(temperature, pressure, params)
    if not w_spinodal:
        if verbose:
            print("No spinodal found")
        return None

    x_spinodal = [
        get_feos_npt_state(
            temperature, pressure, params, wi_lst=[w, 1.0 - w]
        ).molefracs[0]
        for w in w_spinodal
    ]

    x_candidates = list(x_spinodal)
    change = tol + 1.0

    for _ in range(max_iter):
        x_new = solve_tangent_newton(
            temperature,
            pressure,
            params,
            x_fixed=x_candidates[-1],
            x_init=x_candidates[-2],
        )

        for inner in range(max_iter):
            if x_candidates[-1] < x_candidates[-2]:
                if x_new > max(x_spinodal):
                    break
                x_init = x_candidates[-2] + (inner + 1) * (
                    1.0 - x_candidates[-2]
                ) / (max_iter + 1)
            else:
                if x_new < min(x_spinodal):
                    break
                x_init = x_candidates[-2] - (inner + 1) * (
                    x_candidates[-2]
                ) / (max_iter + 1)

            x_new = solve_tangent_newton(
                temperature,
                pressure,
                params,
                x_fixed=x_candidates[-1],
                x_init=x_init,
            )

        if inner == max_iter - 1:
            print("NR failed")

        x_candidates.append(x_new)

        if len(x_candidates) > 3:
            change = abs(x_candidates[-1] - x_candidates[-3]) + abs(
                x_candidates[-2] - x_candidates[-4]
            )
        if change <= tol:
            break

    return [
        get_feos_npt_state(
            temperature, pressure, params, xi_lst=[x, 1.0 - x]
        ).massfracs()[0]
        for x in x_candidates[-2:]
    ]


# pylint: disable=too-many-arguments,too-many-positional-arguments


def _tangent_objective(
    x1, temperature, pressure, params, g_red_fixed, x1_fixed
):
    g_red, d_g_red_dx1, _ = reduced_gibbs_and_derivatives(
        x1, temperature, pressure, params
    )
    return (g_red - g_red_fixed) / (x1 - x1_fixed) - d_g_red_dx1


# pylint: disable=too-many-arguments,too-many-positional-arguments


def _tangent_objective_dlnx1(
    x1, temperature, pressure, params, g_red_fixed, x1_fixed
):
    g_red, d_g_red_dx1, x1_dd_g_red_ddx1 = reduced_gibbs_and_derivatives(
        x1, temperature, pressure, params
    )

    return (
        (x1 / (x1 - x1_fixed)) * d_g_red_dx1
        - (x1 / (x1 - x1_fixed) ** 2) * (g_red - g_red_fixed)
        - x1_dd_g_red_ddx1
    )


# pylint: disable=too-many-arguments,too-many-positional-arguments


def solve_tangent_newton(
    temperature, pressure, params, x_fixed, x_init, eps=_EPS
):
    """Solve the one-tangent condition using Newton-CG on sigmoid(zeta)."""
    g_red_fixed, _, _ = reduced_gibbs_and_derivatives(
        x_fixed, temperature, pressure, params
    )

    zeta_0 = logit(float(np.clip(x_init, eps, 1.0 - eps)))

    def f_zeta(zeta):
        x1 = _sigmoid(zeta)
        return _tangent_objective(
            x1, temperature, pressure, params, g_red_fixed, x_fixed
        )

    def sq_f_zeta(zeta):
        f_val = f_zeta(zeta[0])
        return f_val * f_val

    def df_dzeta(zeta):
        x1 = _sigmoid(zeta)
        df_dlnx1_val = _tangent_objective_dlnx1(
            x1, temperature, pressure, params, g_red_fixed, x_fixed
        )
        return (1 - x1) * df_dlnx1_val

    def sq_f_dzeta(zeta):
        f_val = f_zeta(zeta[0])
        df_dzeta_val = df_dzeta(zeta[0])
        return np.array([2 * f_val * df_dzeta_val])

    res = minimize(
        sq_f_zeta,
        x0=[zeta_0],
        jac=sq_f_dzeta,
        method="Newton-CG",
        options={"maxiter": 200, "xtol": 1e-10},
    )

    return _sigmoid(res.x[0])


# pylint: disable=too-many-arguments,too-many-positional-arguments


def reduced_gibbs_and_derivatives(x1, temperature, pressure, params, eps=_EPS):
    """Return reduced Gibbs energy and derivatives for a binary mixture."""
    x1 = float(np.clip(x1, eps, 1.0 - eps))
    x2 = float(np.clip(1.0 - x1, eps, 1.0 - eps))
    xi_lst = [x1, x2]

    mix_state = get_feos_npt_state(
        temperature, pressure, params, xi_lst=xi_lst
    )
    ln_phi = mix_state.ln_phi()

    g_red = x1 * (np.log(x1) + ln_phi[0]) + x2 * (np.log(x2) + ln_phi[1])
    d_g_red_dx1 = np.log(x1) - np.log(x2) + ln_phi[0] - ln_phi[1]
    x1_dd_g_red_ddx1 = (
        1 - (mix_state.total_moles * mix_state.dln_phi_dnj())[1][0]
    ) / x2

    return g_red, d_g_red_dx1, x1_dd_g_red_ddx1
