points to points

binorm.py

import numpy as np
from numba import njit, prange
import numba

numba.config.NUMBA_NUM_THREADS = 20
numba.set_num_threads(20)

@njit(parallel=True, fastmath=True)
def binorm(xj, yj, x, y, xsig, ysig, nsig2_min=9.):
    r = np.zeros_like(x)

    for i in prange(x.shape[0]):
        for j in range(xj.shape[0]):
            dy2 = ((y[i] - yj[j]) / ysig[i])**2
            if dy2 > nsig2_min:
                continue
            dx2 = ((x[i] - xj[j]) / xsig[i])**2
            if dx2 > nsig2_min:
                continue
            r[i] += np.exp(-0.5 * (dx2 + dy2))
        r[i] *= 0.5 / (np.pi * xsig[i] * ysig[i])
    return r

binorm_weighted.py

import numpy as np
from numba import njit, prange
import numba

numba.config.NUMBA_NUM_THREADS = 20
# numba.set_num_threads(20)


@njit(parallel=True, fastmath=True)
def binorm(xj, yj, wj, x, y, xsig, ysig, nsig2_min=9.):
    r = np.zeros_like(x)

    for i in prange(x.shape[0]):
        for j in range(xj.shape[0]):
            dy2 = ((y[i] - yj[j]) / ysig[i])**2
            if dy2 > nsig2_min:
                continue
            dx2 = ((x[i] - xj[j]) / xsig[i])**2
            if dx2 > nsig2_min:
                continue
            r[i] += np.exp(-0.5 * (dx2 + dy2)) * wj[j]
        r[i] *= 0.5 / (np.pi * xsig[i] * ysig[i])
    return r

binorm_weighted_clip.py

import numpy as np
from numba import njit, prange
import numba

numba.config.NUMBA_NUM_THREADS = 20
# numba.set_num_threads(20)


@njit(parallel=True, fastmath=False)
def binorm_weighted_clip(xj, yj, wj, x, y, xsig, ysig, rcut=10.):
    """
    flat distribution beyond r_cut
    """
    r2_cut = rcut**2
    wj_sum = 0.0
    for j in prange(xj.shape[0]):
        wj_sum += wj[j]

    c = np.exp(-0.5 * r2_cut)
    r = np.full_like(x, c * wj_sum)

    for i in prange(x.shape[0]):
        for j in range(xj.shape[0]):
            r2 = ((y[i] - yj[j]) / ysig[i])**2
            if r2 >= r2_cut:
                continue
            r2 = ((x[i] - xj[j]) / xsig[i])**2 + r2
            if r2 >= r2_cut:
                continue
            r[i] += (np.exp(-0.5 * r2) - c) * wj[j]
        r[i] /= 2 * np.pi * xsig[i] * ysig[i]
    return r

binorm_weighted_clip_renorm.py

import numpy as np
from numba import njit, prange
import numba
import math

numba.config.NUMBA_NUM_THREADS = 20
# numba.set_num_threads(20)


@njit(parallel=True, fastmath=False)
def binorm_weighted_clip_renorm(xj, yj, wj, x, y, xsig, ysig, ymin, ymax, rcut=8.0):
    """
    flat distribution beyond r_cut, renormalized in [ymin, ymax]

    Parameters
    ----------
    xj, yj, wj : shape (nj,)
        Model points
    x, y, xsig, ysig : shape (n,)
        Observation points and errors
    ymin, ymax : float
        Observation window
    rcut : float
        rcut >= 7 is recommended

    Notes
    -----
    cdf(x1, x2) = 0.5 * (erf(x2 / sqrt(2)) - erf(x1 / sqrt(2)))
    """

    r2_cut = rcut**2
    wj_sum = 0.0
    for j in prange(xj.shape[0]):
        wj_sum += wj[j]

    c = math.exp(-0.5 * r2_cut)  # lower bound of normal
    p = np.full_like(x, c * wj_sum)
    s = math.sqrt(2)

    for i in prange(x.shape[0]):
        norm = 0.0
        for j in range(xj.shape[0]):
            # norm += (math.erf((ymax - yj[j]) / (ysig[i] * s)) -
            #          math.erf((ymin - yj[j]) / (ysig[i] * s))) * wj[j]  # *0.5

            d = (ymax - yj[j]) / ysig[i]
            if d < 6.0:
                pobs = math.erf(d / s)
            else:
                pobs = 1.0  # precision 1e-9
            d = (ymin - yj[j]) / ysig[i]
            if d > -6.0:
                pobs = pobs - math.erf(d / s)
            else:
                pobs = pobs + 1.0
            norm = norm + pobs * wj[j]

            r2 = ((y[i] - yj[j]) / ysig[i])**2
            if r2 >= r2_cut:
                continue
            r2 = ((x[i] - xj[j]) / xsig[i])**2 + r2
            if r2 >= r2_cut:
                continue

            p[i] += (math.exp(-0.5 * r2) - c) * wj[j]
        p[i] /= math.pi * xsig[i] * ysig[i] * norm  # * 2
    return p