Source code for pyretis.analysis.interface_placement

# -*- coding: utf-8 -*-
# Copyright (c) 2023, PyRETIS Development Team.
# Distributed under the LGPLv2.1+ License. See LICENSE for more info.
"""Automatic interface placement from a crossing-probability curve.

Ported from the upstream ``inftools`` ``infinit`` tool
(``misc/infinit_helper.py``): given a (binless) crossing-probability curve
-- order parameter ``x`` versus total crossing probability ``p`` -- place
TIS/RETIS interfaces so that every ensemble carries approximately the same
target local crossing probability ``pl_target`` (geometric spacing in ``p``).

This is the self-contained scientific core of ``infinit``; the iterative
"run a short simulation, re-estimate the crossing probability with WHAM,
re-place the interfaces, repeat" driver is a separate orchestration layer
(it runs ``pyretisrun`` and rewrites the input) and is not included here.
The crossing-probability curve consumed here is exactly the one
:func:`pyretis.analysis.path_weights.get_path_weights` can emit (its
binless ``Pcross``), or any monotonically-decreasing ``p(x)``.
"""
import numpy as np


__all__ = ['estimate_interface_positions']




def estimate_interface_positions(x, p, pl_target=0.3, num_ens=None):
    """Estimate interfaces equally spaced with respect to ``pl_target``.

    The interfaces are placed so the cumulative crossing probability drops
    by a constant factor (``pl_target``) from one interface to the next --
    i.e. each ensemble has a local crossing probability of about
    ``pl_target``. The final (state-B / cap) interface is not added, so the
    probability to reach it stays ``pl_target``.

    Parameters
    ----------
    x : array_like
        The order-parameter values of the crossing-probability curve
        (assumed to contain no duplicates, as produced by the binless
        ``Pcross``).
    p : array_like
        The total crossing probability at each ``x`` (monotonically
        decreasing, ``p[0] == 1``).
    pl_target : float, optional
        The target local crossing probability per ensemble (default 0.3).
    num_ens : int, optional
        The number of ensembles to place. If ``None`` (default) it is
        derived from ``pl_target`` and the final crossing probability
        ``p[-1]``.

    Returns
    -------
    interfaces : list of float
        The estimated interface positions (the first is ``x[0]``; the
        state-B interface is not included).
    pl_used : float
        The actual per-interface local crossing probability used (so the
        ``num_ens`` interfaces tile ``p[-1]`` exactly).

    """
    x = np.asarray(x, dtype=float)
    p = np.asarray(p, dtype=float)

    # Estimate how many ensembles we need (unless given). Floor to at
    # least 1: when ``p[-1]`` is not small enough relative to
    # ``pl_target`` the truncated estimate is 0, which would make the
    # ``p[-1] ** (1 / num_ens)`` below divide by zero.
    if not num_ens:
        num_ens = max(1, int(np.log(p[-1]) / np.log(pl_target)))

    # The actual per-interface local crossing probability so that the
    # num_ens steps tile p[-1] exactly.
    pl_used = p[-1] ** (1.0 / num_ens)

    # p decreases, so -log(p) increases: interpolate -log(p) vs x to find
    # the order parameter where the cumulative probability hits pl_used^(i+1).
    targets = [-np.log(pl_used ** (i + 1)) for i in range(num_ens - 1)]
    interior = np.interp(targets, -np.log(p), x)

    interfaces = [float(x[0])] + [float(v) for v in interior]
    return interfaces, float(pl_used)