Find the intersection points of a half-closed simplex
    (pt0, pt1] and IC constraints.

    Parameters
    ----------
    C : ndarray(float, ndim=2)
        The 4 by 2 array for storing the generated points of
        one action profile. One action profile can only
        generate at most 4 points.

    n : scalar(int)
        The number of intersection points that have been found.

    weights : ndarray(float, ndim=1)
        The size 2 array for storing the weights when calculate
        the intersection point of simplex and IC constraints.

    IC : ndarray(float, ndim=1)
        The minimum IC continuation values.

    pt0 : ndarray(float, ndim=1)
        Coordinates of the starting point of the simplex.

    pt1 : ndarray(float, ndim=1)
        Coordinates of the ending point of the simplex.

    tol : scalar(float)
        The tolerance for checking if two values are equal.

    Returns
    -------
    n : scalar(int)
        The updated number of found intersection points.

def _intersect(C, n, weights, IC, pt0, pt1, tol):
    """
    Find the intersection points of a half-closed simplex
    (pt0, pt1] and IC constraints.

    Parameters
    ----------
    C : ndarray(float, ndim=2)
        The 4 by 2 array for storing the generated points of
        one action profile. One action profile can only
        generate at most 4 points.

    n : scalar(int)
        The number of intersection points that have been found.

    weights : ndarray(float, ndim=1)
        The size 2 array for storing the weights when calculate
        the intersection point of simplex and IC constraints.

    IC : ndarray(float, ndim=1)
        The minimum IC continuation values.

    pt0 : ndarray(float, ndim=1)
        Coordinates of the starting point of the simplex.

    pt1 : ndarray(float, ndim=1)
        Coordinates of the ending point of the simplex.

    tol : scalar(float)
        The tolerance for checking if two values are equal.

    Returns
    -------
    n : scalar(int)
        The updated number of found intersection points.
    """
    for i in range(2):
        if (abs(pt0[i] - pt1[i]) < tol):
            if (abs(pt1[i] - IC[i]) < tol):
                x = pt1[1-i]
            else:
                continue
        else:
            weights[i] = (pt0[i] - IC[i]) / (pt0[i] - pt1[i])
            # pt0 is not included to avoid duplication
            # weights in (0, 1]
            if (0 < weights[i] <= 1):
                x = (1 - weights[i]) * pt0[1-i] + weights[i] * pt1[1-i]
            else:
                continue
        # x has to be strictly higher than IC[1-j]
        # if it is equal, then it means IC is one of the vertex
        # it will be added to C in below
        if x - IC[1-i] > tol:
            C[n, i] = IC[i]
            C[n, 1-i] = x
            n += 1
        elif x - IC[1-i] > -tol:
            # to avoid duplication when IC is a vertex
            break

    return n