Python numba.float64() Examples

def csc_diagonal_from_array(m, array):
    """

    :param m:
    :param array:
    :return:
    """
    indptr = np.empty(m + 1, dtype=np.int32)
    indices = np.empty(m, dtype=np.int32)
    data = np.empty(m, dtype=np.float64)
    for i in range(m):
        indptr[i] = i
        indices[i] = i
        data[i] = array[i]
    indptr[m] = m

    return indices, indptr, data 

def csc_sprealloc_f(An, Aindptr, Aindices, Adata, nzmax):
    """
    Change the max # of entries a sparse matrix can hold.
    :param An: number of columns
    :param Aindptr: csc column pointers
    :param Aindices: csc row indices
    :param Adata: csc data
    :param nzmax:new maximum number of entries
    :return: indices, data, nzmax
    """

    if nzmax <= 0:
        nzmax = Aindptr[An]

    length = min(nzmax, len(Aindices))
    Ainew = np.empty(nzmax, dtype=nb.int32)
    for i in range(length):
        Ainew[i] = Aindices[i]

    length = min(nzmax, len(Adata))
    Axnew = np.empty(nzmax, dtype=nb.float64)
    for i in range(length):
        Axnew[i] = Adata[i]

    return Ainew, Axnew, nzmax 

def csc_mat_vec_ff(m, n, Ap, Ai, Ax, x):
    """
    Sparse matrix times dense column vector, y = A * x.
    :param m: number of rows
    :param n: number of columns
    :param Ap: pointers
    :param Ai: indices
    :param Ax: data
    :param x: vector x (n)
    :return: vector y (m)
    """

    assert n == x.shape[0]

    y = np.zeros(m, dtype=nb.float64)
    for j in range(n):
        for p in range(Ap[j], Ap[j + 1]):
            y[Ai[p]] += Ax[p] * x[j]
    return y 

def init_w(t, k, x, lx, w):
    tb = t[k]
    n = len(t)
    m = len(x)
    h = np.zeros(6, dtype=np.float64)#([0]*6 )
    hh = np.zeros(5, dtype=np.float64)##np.array([0]*5)
    te = t[n - k - 1]
    l1 = k + 1
    l2 = l1 + 1
    for i in range(m):
        arg = x[i]
        if arg < tb:
            arg = tb
        if arg > te:
            arg = te
        while not (arg < t[l1] or l1 == (n - k - 1)):
            l1 = l2
            l2 = l1 + 1
        h, hh = fpbspl(t, n, k, arg, l1, h, hh)

        lx[i] = l1 - k - 1
        for j in range(k + 1):
            w[i][j] = h[j]
    return w 

def csc_to_dense(m, n, indptr, indices, data):
    """
    Convert csc matrix to dense
    :param m:
    :param n:
    :param indptr:
    :param indices:
    :param data:
    :return: 2d numpy array
    """
    val = np.zeros((m, n), dtype=np.float64)

    for j in range(n):
        for p in range(indptr[j], indptr[j + 1]):
            val[indices[p], j] = data[p]
    return val 

def ij_bbox(self,
                xy_bbox: Tuple[float, float, float, float],
                xy_border: float = 0.0,
                ij_border: int = 0,
                gu: bool = False) -> Tuple[int, int, int, int]:
        """
        Compute bounding box in i,j pixel coordinates given a bounding box *xy_bbox* in x,y coordinates.

        :param xy_bbox: Bounding box (x_min, y_min, x_max, y_max) given in the same CS as x and y.
        :param xy_border: If non-zero, grows the bounding box *xy_bbox* before using it for comparisons. Defaults to 0.
        :param ij_border: If non-zero, grows the returned i,j bounding box and clips it to size. Defaults to 0.
        :param gu: Use generic ufunc for the computation (may be faster). Defaults to False.
        :return: Bounding box in (i_min, j_min, i_max, j_max) in pixel coordinates.
            Returns ``(-1, -1, -1, -1)`` if *xy_bbox* isn't intersecting any of the x,y coordinates.
        """
        xy_bboxes = np.array([xy_bbox], dtype=np.float64)
        ij_bboxes = np.full_like(xy_bboxes, -1, dtype=np.int64)
        self.ij_bboxes(xy_bboxes, xy_border=xy_border, ij_border=ij_border, ij_bboxes=ij_bboxes, gu=gu)
        # noinspection PyTypeChecker
        return tuple(map(int, ij_bboxes[0])) 

def csc_diagonal(m, value=1.0):
    """
    Build CSC diagonal matrix of the given value
    :param m: size
    :param value: value
    :return: CSC matrix
    """
    indptr = np.empty(m + 1, dtype=np.int32)
    indices = np.empty(m, dtype=np.int32)
    data = np.empty(m, dtype=np.float64)
    for i in range(m):
        indptr[i] = i
        indices[i] = i
        data[i] = value
    indptr[m] = m

    return indices, indptr, data 

def rotate_(a, theta):
    """
    rotate a theta from vector a
    """
    d0 = np.cos(theta[2]) * (np.cos(theta[1]) * a[0]
                             + np.sin(theta[1]) * (
                                     np.sin(theta[0]) * a[1] + np.cos(
                                 theta[0]) * a[2])) \
         - np.sin(theta[2]) * (
                 np.cos(theta[0]) * a[1] - np.sin(theta[0]) * a[2])
    d1 = np.sin(theta[2]) * (np.cos(theta[1]) * a[0]
                             + np.sin(theta[1]) * (
                                     np.sin(theta[0]) * a[1] + np.cos(
                                 theta[0]) * a[2])) \
         + np.cos(theta[2]) * (
                 np.cos(theta[0]) * a[1] - np.sin(theta[0]) * a[2])
    d2 = -np.sin(theta[1]) * a[0] + np.cos(theta[1]) * \
         (np.sin(theta[0]) * a[1] + np.cos(theta[0]) * a[2])

    vector = np.zeros(3, dtype=numba.float64)
    vector[0] = d0
    vector[1] = d1
    vector[2] = d2

    return vector 

def csc_sub_matrix_cols(Am, Anz, Ap, Ai, Ax, cols):
    """
    Get SCS arbitrary sub-matrix with all the rows
    :param Am: number of rows
    :param Anz: number of non-zero entries
    :param Ap: Column pointers
    :param Ai: Row indices
    :param Ax: Data
    :param cols: column indices to keep
    :return: CSC sub-matrix (n, new_col_ptr, new_row_ind, new_val)
    """

    n_cols = len(cols)
    n = 0
    p = 0
    Bx = np.empty(Anz, dtype=nb.float64)
    Bi = np.empty(Anz, dtype=nb.int32)
    Bp = np.empty(n_cols + 1, dtype=nb.int32)

    Bp[p] = 0

    for j in cols:  # for each column selected ...
        for k in range(Ap[j], Ap[j + 1]):  # for each row of the column j of A...
            # store the values if the row was found in rows
            Bx[n] = Ax[k]  # store the value
            Bi[n] = Ai[k]  # store the row index
            n += 1
        p += 1
        Bp[p] = n

    Bp[p] = n

    return n, Bp, Bi[:n], Bx[:n]


# @nb.njit("Tuple((i8, i4[:], i4[:], f8[:]))(i8, i8, i4[:], i4[:], f8[:], i4[:])") 

def csc_sub_matrix(Am, Anz, Ap, Ai, Ax, rows, cols):
    """
    Get SCS arbitrary sub-matrix
    :param Am: number of rows
    :param Anz: number of non-zero entries
    :param Ap: Column pointers
    :param Ai: Row indices
    :param Ax: Data
    :param rows: row indices to keep
    :param cols: column indices to keep
    :return: CSC sub-matrix (n, new_col_ptr, new_row_ind, new_val)
    """
    n_cols = len(cols)

    Bx = np.zeros(Anz, dtype=np.float64)
    Bi = np.empty(Anz, dtype=np.int32)
    Bp = np.empty(n_cols + 1, dtype=np.int32)

    n = 0
    p = 0
    Bp[p] = 0

    for j in cols:  # for each column selected ...
        i = 0
        for r in rows:
            for k in range(Ap[j], Ap[j + 1]):  # for each row of the column j of A...
                if Ai[k] == r:
                    Bx[n] = Ax[k]  # store the value
                    Bi[n] = i  # row index in the new matrix
                    i += 1
                    n += 1
            if i == 0:
                i += 1
        p += 1
        Bp[p] = n

    Bp[p] = n

    return n, Bp, Bi[:n], Bx[:n] 

def annhilate_k_device(K_val, C_val, output):
    k_4 = numba.cuda.local.array(32, dtype=numba.float64)
    project_val_cuda(K_val, k_4, 4)
    for i in range(32):
        k_4[i] = -k_4[i]
    k_4[0] += K_val[0]
    gp_device(k_4, C_val, output)
    normalise_mv_device(output) 

def xalloc(n):
    return np.zeros(n, dtype=nb.float64) 

def ray_trace(x, y, poly):
    """
    Determines for some set of x and y coordinates, which of those coordinates is within `poly`. Ray trace is \
    generally called as an internal function, see :func:`.poly_clip`

    :param x: A 1D numpy array of x coordinates.
    :param y: A 1D numpy array of y coordinates.
    :param poly: The coordinates of a polygon as a numpy array (i.e. from geo_json['coordinates']
    :return: A 1D boolean numpy array, true values are those points that are within `poly`.
    """

    @vectorize([bool_(float64, float64)])
    def ray(x, y):
        # where xy is a coordinate
        n = len(poly)
        inside = False
        p2x = 0.0
        p2y = 0.0
        xints = 0.0
        p1x, p1y = poly[0]
        for i in range(n + 1):
            p2x, p2y = poly[i % n]
            if y > min(p1y, p2y):
                if y <= max(p1y, p2y):
                    if x <= max(p1x, p2x):
                        if p1y != p2y:
                            xints = (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x
                        if p1x == p2x or x <= xints:
                            inside = not inside
            p1x, p1y = p2x, p2y
        return inside

    return ray(x, y) 

def cost_line_to_line_device(L1, L2):
    R_val = numba.cuda.local.array(32, dtype=numba.float64)
    rotor_between_lines_device(L1, L2, R_val)
    return rotor_cost_device(R_val) 

def rotor_cost_device(R_val):
    translation_val = numba.cuda.local.array(32, dtype=numba.float64)
    rotation_val = numba.cuda.local.array(32, dtype=numba.float64)
    ep_val = numba.cuda.local.array(32, dtype=numba.float64)
    for i in range(32):
        ep_val[i] = 0.0
    ep_val[4] = 1.0
    ip_device(R_val, ep_val, translation_val)
    for i in range(32):
        rotation_val[i] = R_val[i]
    rotation_val[0] -= 1
    a = abs(gp_mult_with_adjoint_to_scalar(rotation_val))
    b = abs(gp_mult_with_adjoint_to_scalar(translation_val))
    return a + b 

def gp_mult_with_adjoint(value, output):
    other_value = numba.cuda.local.array(32, dtype=numba.float64)
    adjoint_device(value, other_value)
    gp_device(value, other_value, output) 

def rotor_between_lines_device(L1, L2, rotor):
    L21_val = numba.cuda.local.array(32, dtype=numba.float64)
    L12_val = numba.cuda.local.array(32, dtype=numba.float64)

    gp_device(L2, L1, L21_val)
    gp_device(L1, L2, L12_val)

    beta_val = numba.cuda.local.array(32, dtype=numba.float64)
    K_val = numba.cuda.local.array(32, dtype=numba.float64)
    for i in range(32):
        K_val[i] = L21_val[i] + L12_val[i]
        beta_val[i] = 0.0
    K_val[0] += 2.0

    project_val_cuda(K_val, beta_val, 4)

    alpha = 2.0 * K_val[0]

    denominator = math.sqrt(alpha / 2)

    normalisation_val = numba.cuda.local.array(32, dtype=numba.float64)
    output_val = numba.cuda.local.array(32, dtype=numba.float64)
    for i in range(32):
        if i == 0:
            numerator_val = 1.0 - beta_val[i] / alpha
        else:
            numerator_val = -beta_val[i] / alpha
        normalisation_val[i] = numerator_val / denominator
        output_val[i] = L21_val[i]

    output_val[0] += 1
    gp_device(normalisation_val, output_val, rotor) 

def cost_between_objects_device(L1, L2):
    R_val = numba.cuda.local.array(32, dtype=numba.float64)
    rotor_between_objects_device(L1, L2, R_val)
    return rotor_cost_device(R_val) 

def rotor_between_objects_device(L1, L2, rotor):
    L1sqrd_val = numba.cuda.local.array(32, dtype=numba.float64)
    gp_device(L1, L1, L1sqrd_val)
    if L1sqrd_val[0] > 0:
        C_val = numba.cuda.local.array(32, dtype=numba.float64)
        sigma_val = numba.cuda.local.array(32, dtype=numba.float64)
        k_value = numba.cuda.local.array(32, dtype=numba.float64)
        gp_device(L2, L1, C_val)
        C_val[0] += 1.0
        gp_mult_with_adjoint(C_val, sigma_val)
        positive_root_device(sigma_val, k_value)
        annhilate_k_device(k_value, C_val, rotor)
    else:
        L21 = numba.cuda.local.array(32, dtype=numba.float64)
        L12 = numba.cuda.local.array(32, dtype=numba.float64)
        gp_device(L2, L1, L21)
        gp_device(L1, L2, L12)
        sumval = 0.0
        for i in range(32):
            if i == 0:
                sumval += abs(L12[i] + L21[i] - 2.0)
            else:
                sumval += abs(L12[i] + L21[i])
            rotor[i] = -L21[i]
        if sumval < 0.0000001:
            rotor[0] = rotor[0] - 1.0
        else:
            rotor[0] = rotor[0] + 1.0
        normalise_mv_device(rotor) 

def dorst_norm_val_device(sigma_val):
    """ Square Root of Rotors - Implements the norm of a rotor"""
    s_4 = numba.cuda.local.array(32, dtype=numba.float64)
    s_4_sqrd = numba.cuda.local.array(32, dtype=numba.float64)
    project_val_cuda(sigma_val, s_4, 4)
    gp_device(s_4, s_4, s_4_sqrd)
    sqrd_ans = sigma_val[0]*sigma_val[0] - s_4_sqrd[0]
    return math.sqrt(abs(sqrd_ans)) 

def square_root_of_rotor_device(rotor, rotor_root):
    k_value = numba.cuda.local.array(32, dtype=numba.float64)
    sigma_val = numba.cuda.local.array(32, dtype=numba.float64)
    C_val = numba.cuda.local.array(32, dtype=numba.float64)
    for i in range(32):
        C_val[i] = rotor[i]
    C_val[0] += 1.0
    gp_mult_with_adjoint(C_val, sigma_val)
    positive_root_device(sigma_val, k_value)
    annhilate_k_device(k_value, C_val, rotor_root) 

def apply_rotor_device(mv, rotor, output):
    rotor_adjoint = numba.cuda.local.array(32, dtype=numba.float64)
    temp = numba.cuda.local.array(32, dtype=numba.float64)
    adjoint_device(rotor, rotor_adjoint)
    gp_device(mv, rotor_adjoint, temp)
    gp_device(rotor, temp, output) 

def sequential_rotor_estimation_cuda(reference_model_array, query_model_array, n_samples=None, n_objects_per_sample=None, mutation_probability=None):

    if n_samples is None:
        n_samples = int(len(query_model_array)/2)
    if n_objects_per_sample is None:
        n_objects_per_sample = max(int(len(query_model_array)/10), 5)

    # Stack up a list of numbers
    total_matches = n_samples*n_objects_per_sample
    sample_indices = random.sample(range(total_matches), total_matches)

    n_mvs = reference_model_array.shape[0]
    sample_indices = [i % n_mvs for i in sample_indices]

    if mutation_probability is not None:
        reference_model_array_new = []
        mutation_flag = np.random.binomial(1, mutation_probability, total_matches)
        for mut, i in zip(mutation_flag, sample_indices):
            if mut:
                ref_ind = random.sample(range(len(reference_model_array)), 1)[0]
            else:
                ref_ind = i
            reference_model_array_new.append(reference_model_array[ref_ind, :])
        reference_model_array_new = np.array(reference_model_array_new)
    else:
        reference_model_array_new = np.array([reference_model_array[i, :] for i in sample_indices], dtype=np.float64)
    query_model_array_new = np.array([query_model_array[i, :] for i in sample_indices], dtype=np.float64)

    output = np.zeros((n_samples, 32), dtype=np.float64)
    cost_array = np.zeros(n_samples, dtype=np.float64)

    blockdim = 64
    griddim = int(math.ceil(reference_model_array_new.shape[0] / blockdim))

    sequential_rotor_estimation_kernel[griddim, blockdim](reference_model_array_new, query_model_array_new, output, cost_array)

    return output, cost_array 

def sequential_rotor_estimation_chunks(reference_model_array, query_model_array, n_samples, n_objects_per_sample, mutation_probability=None):

    # Stack up a list of numbers
    total_matches = n_samples*n_objects_per_sample
    sample_indices = random.sample(range(total_matches), total_matches)

    n_mvs = reference_model_array.shape[0]
    sample_indices = [i % n_mvs for i in sample_indices]

    if mutation_probability is not None:
        reference_model_array_new = []
        mutation_flag = np.random.binomial(1, mutation_probability, total_matches)
        for mut, i in zip(mutation_flag, sample_indices):
            if mut:
                ref_ind = random.sample(range(len(reference_model_array)), 1)[0]
            else:
                ref_ind = i
            reference_model_array_new.append(reference_model_array[ref_ind, :])
        reference_model_array_new = np.array(reference_model_array_new)
    else:
        reference_model_array_new = np.array([reference_model_array[i, :] for i in sample_indices], dtype=np.float64)
    query_model_array_new = np.array([query_model_array[i, :] for i in sample_indices], dtype=np.float64)

    output = np.zeros((n_samples, 32), dtype=np.float64)
    cost_array = np.zeros(n_samples, dtype=np.float64)

    sequential_rotor_estimation_chunks_jit(reference_model_array_new, query_model_array_new, output, cost_array)

    return output, cost_array 

def transform_lists_to_arrays(module, to_change, __funcs, vec=False, cache_blacklist=set([])):
    if vec:
        conv_fun = numba.vectorize
        extra_args = extra_args_vec
    else:
        conv_fun = numba.njit
        extra_args = extra_args_std
    for s in to_change:
        func = s.split('.')[-1]
        mod = '.'.join(s.split('.')[:-1])
        fake_mod = __funcs[mod]

        try:
            real_mod = getattr(module, mod)
        except:
            real_mod = module
            for s in mod.split('.'):
                real_mod = getattr(real_mod, s)

        orig_func = getattr(real_mod, func)
        source = inspect.getsource(orig_func)
        source = remove_for_numba(source) # do before anything else
        source = return_value_numpy(source)
        source = re.sub(list_mult_expr, numpy_not_list_expr, source)
#        if 'longitude_obliquity_nutation' in s:
#        print(source)
        numba_exec_cacheable(source, fake_mod.__dict__, fake_mod.__dict__)
        new_func = fake_mod.__dict__[func]
        do_cache = caching and func not in cache_blacklist
        obj = conv_fun(cache=do_cache, **extra_args)(new_func)
        __funcs[func] = obj
        fake_mod.__dict__[func] = obj
        obj.__doc__ = ''


#set_signatures = {'Clamond': [numba.float64(numba.float64, numba.float64, numba.boolean),
#                              numba.float64(numba.float64, numba.float64, numba.optional(numba.boolean))
#                              ]
#                    } 

def points_to_angles(points):
    angles = np.empty(len(points) - 1, dtype=np.float64)
    for i, (p1, p2) in enumerate(zip(points[0:-1], points[1:])):
        angles[i] = np.arctan2(p2[1] - p1[1], p2[0] - p1[0])
    return angles 

def __init__(self, j2_max, j3_max):
        self._size = j2_max + j3_max + 1
        self.workspace = np.empty(4 * self._size, dtype=np.float64) 

def _rotation_matrix(axis, angle):
    """
    This method produces a rotation matrix given an axis and an angle.
    """
    axis_norm = _norm(axis)
    for k in range(3):
        axis[k] = axis[k] / axis_norm
    axis_squared = axis**2
    cos_angle = np.cos(angle)
    sin_angle = np.sin(angle)

    rotation_matrix = np.zeros((3,3), dtype=float64)

    rotation_matrix[0, 0] = cos_angle + axis_squared[0]*(1.0-cos_angle)
    rotation_matrix[0, 1] = axis[0]*axis[1]*(1.0-cos_angle) - axis[2]*sin_angle
    rotation_matrix[0, 2] = axis[0]*axis[2]*(1.0-cos_angle) + axis[1]*sin_angle

    rotation_matrix[1, 0] = axis[1]*axis[0]*(1.0-cos_angle) + axis[2]*sin_angle
    rotation_matrix[1, 1] = cos_angle + axis_squared[1]*(1.0-cos_angle)
    rotation_matrix[1, 2] = axis[1]*axis[2]*(1.0-cos_angle) - axis[0]*sin_angle

    rotation_matrix[2, 0] = axis[2]*axis[0]*(1.0-cos_angle) - axis[1]*sin_angle
    rotation_matrix[2, 1] = axis[2]*axis[1]*(1.0-cos_angle) + axis[0]*sin_angle
    rotation_matrix[2, 2] = cos_angle + axis_squared[2]*(1.0-cos_angle)

    return rotation_matrix 

def rot_axis(angle, axis):
    # RX = np.array([ [1,             0,              0],
    #                 [0, np.cos(gamma), -np.sin(gamma)],
    #                 [0, np.sin(gamma),  np.cos(gamma)]])
    #
    # RY = np.array([ [ np.cos(beta), 0, np.sin(beta)],
    #                 [            0, 1,            0],
    #                 [-np.sin(beta), 0, np.cos(beta)]])
    #
    # RZ = np.array([ [np.cos(alpha), -np.sin(alpha), 0],
    #                 [np.sin(alpha),  np.cos(alpha), 0],
    #                 [            0,              0, 1]])
    cg = np.cos(angle)
    sg = np.sin(angle)
    if axis == 0:  # X
        v = [0, 4, 5, 7, 8]
    elif axis == 1:  # Y
        v = [4, 0, 6, 2, 8]
    else:  # Z
        v = [8, 0, 1, 3, 4]
    RX = np.zeros(9, dtype=numba.float64)
    RX[v[0]] = 1.0
    RX[v[1]] = cg
    RX[v[2]] = -sg
    RX[v[3]] = sg
    RX[v[4]] = cg
    return RX.reshape(3, 3) 

def calc_nrst_dist_gpu(gids,xs,ys,densities):   
    from numba import cuda, float64,float32
    @cuda.jit
    def calc_nrst_dist_cuda(gids,xs,ys,densities,nrst_dists,parent_gids,n):
        '''
        Numba Cuda code for calculate nearest point with higher density.
        gids: identifier of geometry
        xs: lons' array
        ys: lats' array
        densities: densities' array
        nrst_dists: results of the nearest distance
        parent_gids: results of gid of the nearest point with higher density
        n: array size
        '''
        i=cuda.grid(1)
        if i<n:
            xi=xs[i]
            yi=ys[i]
            density=densities[i]
            nrst_dist=float32(1e100)
            parent_gid=np.int(-1)
            for j in range(n):
                xd=xs[j]-xi
                yd=ys[j]-yi
                gidd=gids[j]
                distpow2=xd**2+yd**2
                if densities[j]>density and distpow2<nrst_dist:
                    nrst_dist=distpow2
                    parent_gid=gidd
            nrst_dists[i]=math.sqrt(float64(nrst_dist))
            parent_gids[i]=parent_gid
    
    n=xs.shape[0]
    xs=np.ascontiguousarray((xs-xs.min()).astype(np.float32))
    ys=np.ascontiguousarray((ys-ys.min()).astype(np.float32))
    threadsperblock = 1024
    blockspergrid = np.int((n + (threadsperblock - 1)) / threadsperblock)
    dev_nrst_dists=cuda.device_array(n)
    dev_parent_gids=cuda.device_array_like(gids)
    calc_nrst_dist_cuda[blockspergrid,threadsperblock](cuda.to_device(np.ascontiguousarray(gids))\
                       ,cuda.to_device(xs),cuda.to_device(ys),cuda.to_device(np.ascontiguousarray(densities))\
                       ,dev_nrst_dists,dev_parent_gids,n)
    return (dev_nrst_dists.copy_to_host(),dev_parent_gids.copy_to_host())