import numpy as np


#def project(a, b):  # projects b onto the unit cube, using a as an offset origin
#    a = np.asanyarray(a)
#    b = np.asanyarray(b)
#    d = b - a
#    if all(np.abs(d) <= 1e-8):
#        return a
#    inner = -1 / np.where(np.abs(d) > 1e-8, d, np.where(d >= 0, 1e-8, -1e-8))
#    small = np.minimum(1 - a, a - 1)
#    large = np.maximum(0 - a, a - 0)
#    rescale = np.min(np.maximum(inner * small, inner * large))
#    assert np.isfinite(rescale), rescale
#    if rescale < 1:
#        return a + (rescale - 1e-6) * d
#    else:
#        return b


def project(p, a, eps=1e-8):
    # https://www.desmos.com/calculator/gdcu0ivk0i
    p = np.asanyarray(p)
    a = np.asanyarray(a)
    d = p - a
    if all(np.abs(d) <= eps):
        # we might still be inching out of bounds, so just to be sure:
        a[a <= 0] = 0
        a[a >= 1] = 1
        return a

    inner = 1 / np.where(np.abs(d) > eps, d, np.where(d >= 0, eps, -eps))
    small = -np.abs(p - 1)  # np.minimum(1 - p, p - 1)
    large = np.abs(p)  # np.maximum(0 - p, p - 0)
    rescale = np.min(np.maximum(inner * small, inner * large))

    if rescale <= 1:
        b = p - max(0, rescale - 1e-8) * d
        return b
    else:
        return a