Approx{1,2}_{v2,uniform}, af_gemm

arrayfire/blas.py

@@ -202,3 +202,109 @@ def dot(lhs, rhs, lhs_opts=MATPROP.NONE, rhs_opts=MATPROP.NONE, return_scalar =
         safe_call(backend.get().af_dot(c_pointer(out.arr), lhs.arr, rhs.arr,
                                        lhs_opts.value, rhs_opts.value))
         return out
+
+def gemm(lhs, rhs, alpha=1.0, beta=0.0, lhs_opts=MATPROP.NONE, rhs_opts=MATPROP.NONE, C=None):
+    """
+    BLAS general matrix multiply (GEMM) of two af_array objects.
+
+    This provides a general interface to the BLAS level 3 general matrix multiply (GEMM), which is generally defined as:
+
+    C = α ∗ opA(A) opB(B)+ β∗C
+
+    where α (alpha) and β (beta) are both scalars; A and B are the matrix multiply operands;
+    and opA and opB are noop (if AF_MAT_NONE) or transpose (if AF_MAT_TRANS) operations 
+    on A or B before the actual GEMM operation.
+    Batched GEMM is supported if at least either A or B have more than two dimensions
+    (see af::matmul for more details on broadcasting).
+    However, only one alpha and one beta can be used for all of the batched matrix operands.
+
+    Parameters
+    ----------
+
+    lhs : af.Array
+          A 2 dimensional, real or complex arrayfire array.
+
+    rhs : af.Array
+          A 2 dimensional, real or complex arrayfire array.
+
+    alpha : scalar
+
+    beta : scalar
+
+    lhs_opts: optional: af.MATPROP. default: af.MATPROP.NONE.
+              Can be one of
+               - af.MATPROP.NONE   - If no op should be done on `lhs`.
+               - af.MATPROP.TRANS  - If `lhs` has to be transposed before multiplying.
+               - af.MATPROP.CTRANS - If `lhs` has to be hermitian transposed before multiplying.
+
+    rhs_opts: optional: af.MATPROP. default: af.MATPROP.NONE.
+              Can be one of
+               - af.MATPROP.NONE   - If no op should be done on `rhs`.
+               - af.MATPROP.TRANS  - If `rhs` has to be transposed before multiplying.
+               - af.MATPROP.CTRANS - If `rhs` has to be hermitian transposed before multiplying.
+
+    Returns
+    -------
+
+    out : af.Array
+          Output of the matrix multiplication on `lhs` and `rhs`.
+
+    Note
+    -----
+
+    - The data types of `lhs` and `rhs` should be the same.
+    - Batches are not supported.
+
+    """
+    if C is None:
+        out = Array()
+    else:
+        out = C
+
+    ltype = lhs.dtype()
+
+    if ltype == Dtype.f32:
+        aptr = c_cast(c_pointer(c_float_t(alpha)),c_void_ptr_t)
+        bptr = c_cast(c_pointer(c_float_t(beta)), c_void_ptr_t)
+    elif ltype == Dtype.c32:
+        if isinstance(alpha, af_cfloat_t):
+            aptr = c_cast(c_pointer(alpha), c_void_ptr_t)
+        elif isinstance(alpha, tuple):
+            aptr = c_cast(c_pointer(af_cfloat_t(alpha[0], alpha[1])), c_void_ptr_t)
+        else:
+            aptr = c_cast(c_pointer(af_cfloat_t(alpha)), c_void_ptr_t)
+
+        if isinstance(beta, af_cfloat_t):
+            bptr = c_cast(c_pointer(beta), c_void_ptr_t)
+        elif isinstance(beta, tuple):
+            bptr = c_cast(c_pointer(af_cfloat_t(beta[0], beta[1])), c_void_ptr_t)
+        else:
+            bptr = c_cast(c_pointer(af_cfloat_t(beta)), c_void_ptr_t)
+
+    elif ltype == Dtype.f64:
+        aptr = c_cast(c_pointer(c_double_t(alpha)),c_void_ptr_t)
+        bptr = c_cast(c_pointer(c_double_t(beta)), c_void_ptr_t)
+    elif ltype == Dtype.c64:
+        if isinstance(alpha, af_cdouble_t):
+            aptr = c_cast(c_pointer(alpha), c_void_ptr_t)
+        elif isinstance(alpha, tuple):
+            aptr = c_cast(c_pointer(af_cdouble_t(alpha[0], alpha[1])), c_void_ptr_t)
+        else:
+            aptr = c_cast(c_pointer(af_cdouble_t(alpha)), c_void_ptr_t)
+
+        if isinstance(beta, af_cdouble_t):
+            bptr = c_cast(c_pointer(beta), c_void_ptr_t)
+        elif isinstance(beta, tuple):
+            bptr = c_cast(c_pointer(af_cdouble_t(beta[0], beta[1])), c_void_ptr_t)
+        else:
+            bptr = c_cast(c_pointer(af_cdouble_t(beta)), c_void_ptr_t)
+    elif ltype == Dtype.f16:
+        raise TypeError("fp16 currently unsupported gemm() input type")
+    else:
+        raise TypeError("unsupported input type")
+
+
+    safe_call(backend.get().af_gemm(c_pointer(out.arr),
+                                    lhs_opts.value, rhs_opts.value,
+                                    aptr, lhs.arr, rhs.arr, bptr))
+    return out

arrayfire/library.py

@@ -31,6 +31,13 @@
 c_void_ptr_t  = ct.c_void_p
 c_char_ptr_t  = ct.c_char_p
 c_size_t      = ct.c_size_t
+c_cast        = ct.cast
+
+class af_cfloat_t(ct.Structure):
+    _fields_ = [("real", ct.c_float), ("imag", ct.c_float)]
+
+class af_cdouble_t(ct.Structure):
+    _fields_ = [("real", ct.c_double), ("imag", ct.c_double)]
 
 
 AF_VER_MAJOR = '3'

arrayfire/signal.py

@@ -27,7 +27,7 @@ def _scale_pos_axis1(y_curr, y_orig):
     dy = y_orig[0, 1, 0, 0] - y0
     return((y_curr - y0) / dy)
 
-def approx1(signal, x, method=INTERP.LINEAR, off_grid=0.0, xp = None):
+def approx1(signal, x, method=INTERP.LINEAR, off_grid=0.0, xp = None, output = None):
     """
     Interpolate along a single dimension.Interpolation is performed along axis 0
     of the input array.
@@ -51,6 +51,10 @@ def approx1(signal, x, method=INTERP.LINEAR, off_grid=0.0, xp = None):
     xp : af.Array
          The x-coordinates of the input data points
 
+    output: None or af.Array
+        Optional preallocated output array. If it is a sub-array of an existing af_array,
+        only the corresponding portion of the af_array will be overwritten
+
     Returns
     -------
 
@@ -65,20 +69,86 @@ def approx1(signal, x, method=INTERP.LINEAR, off_grid=0.0, xp = None):
     where N is the length of the first dimension of `signal`.
     """
 
-    output = Array()
+    if output is None:
+        output = Array()
+
+        if(xp is not None):
+            pos0 = _scale_pos_axis0(x, xp)
+        else:
+            pos0 = x
+
+        safe_call(backend.get().af_approx1(c_pointer(output.arr), signal.arr, pos0.arr,
+                                           method.value, c_float_t(off_grid)))
 
-    if(xp is not None):
-        pos0 = _scale_pos_axis0(x, xp)
     else:
-        pos0 = x
+        if(xp is not None):
+            pos0 = _scale_pos_axis0(x, xp)
+        else:
+            pos0 = x
+        safe_call(backend.get().af_approx1_v2(c_pointer(output.arr), signal.arr, pos0.arr,
+                                              method.value, c_float_t(off_grid)))
+    return output
+
+
+def approx1_uniform(signal, x, interp_dim, idx_start, idx_step, method=INTERP.LINEAR, off_grid=0.0, output = None):
+    """
+    Interpolation on one dimensional signals along specified dimension.
+
+    af_approx1_uniform() accepts the dimension to perform the interpolation along the input.
+    It also accepts start and step values which define the uniform range of corresponding indices.
+
+    Parameters
+    ----------
+
+    signal: af.Array
+            Input signal array (signal = f(x))
+
+    x: af.Array
+       The x-coordinates of the interpolation points. The interpolation 
+       function is queried at these set of points.
+
+    interp_dim: scalar
+        is the dimension to perform interpolation across.
+
+    idx_start: scalar
+        is the first index value along interp_dim.
+
+    idx_step: scalar
+        is the uniform spacing value between subsequent indices along interp_dim.
 
-    safe_call(backend.get().af_approx1(c_pointer(output.arr), signal.arr, pos0.arr,
-                                       method.value, c_float_t(off_grid)))
+    method: optional: af.INTERP. default: af.INTERP.LINEAR.
+            Interpolation method.
+
+    off_grid: optional: scalar. default: 0.0.
+            The value used for positions outside the range.
+
+    output: None or af.Array
+        Optional preallocated output array. If it is a sub-array of an existing af_array,
+        only the corresponding portion of the af_array will be overwritten
+
+    Returns
+    -------
+
+    output: af.Array
+            Values calculated at interpolation points.
+
+    """
+
+    if output is None:
+        output = Array()
+
+        safe_call(backend.get().af_approx1_uniform(c_pointer(output.arr), signal.arr, x.arr,
+                                           c_dim_t(interp_dim), c_double_t(idx_start), c_double_t(idx_step),
+                                           method.value, c_float_t(off_grid)))
+    else:
+        safe_call(backend.get().af_approx1_uniform_v2(c_pointer(output.arr), signal.arr, x.arr,
+                                           c_dim_t(interp_dim), c_double_t(idx_start), c_double_t(idx_step),
+                                           method.value, c_float_t(off_grid)))
     return output
 
+
 def approx2(signal, x, y,
-            method=INTERP.LINEAR, off_grid=0.0, xp = None, yp = None 
-           ):
+            method=INTERP.LINEAR, off_grid=0.0, xp = None, yp = None, output = None):
     """
     Interpolate along a two dimension.Interpolation is performed along axes 0 and 1
     of the input array.
@@ -112,6 +182,10 @@ def approx2(signal, x, y,
          The y-coordinates of the input data points. The convention followed is that
          the y-coordinates vary along axis 1
 
+    output: None or af.Array
+        Optional preallocated output array. If it is a sub-array of an existing af_array,
+        only the corresponding portion of the af_array will be overwritten
+
     Returns
     -------
 
@@ -127,22 +201,114 @@ def approx2(signal, x, y,
     and N is the length of the second dimension of `signal`.
     """
 
-    output = Array()
-
-    if(xp is not None):
-        pos0 = _scale_pos_axis0(x, xp)
+    if output is None:
+        output = Array()
+
+        if(xp is not None):
+            pos0 = _scale_pos_axis0(x, xp)
+        else:
+            pos0 = x
+
+        if(yp is not None):
+            pos1 = _scale_pos_axis1(y, yp)
+        else:
+            pos1 = y
+
+        safe_call(backend.get().af_approx2(c_pointer(output.arr), signal.arr,
+                                           pos0.arr, pos1.arr, method.value, c_float_t(off_grid)))
     else:
-        pos0 = x
+        if(xp is not None):
+            pos0 = _scale_pos_axis0(x, xp)
+        else:
+            pos0 = x
+
+        if(yp is not None):
+            pos1 = _scale_pos_axis1(y, yp)
+        else:
+            pos1 = y
+
+        safe_call(backend.get().af_approx2_v2(c_pointer(output.arr), signal.arr,
+                                              pos0.arr, pos1.arr, method.value, c_float_t(off_grid)))
+
+    return output
+
+def approx2_uniform(signal, pos0, interp_dim0, idx_start0, idx_step0, pos1, interp_dim1, idx_start1, idx_step1,
+            method=INTERP.LINEAR, off_grid=0.0, output = None):
+    """
+    Interpolate along two uniformly spaced dimensions of the input array.
+    af_approx2_uniform() accepts two dimensions to perform the interpolation along the input.
+    It also accepts start and step values which define the uniform range of corresponding indices.
+
+    Parameters
+    ----------
+
+    signal: af.Array
+            Input signal array (signal = f(x, y))
+
+    pos0 : af.Array
+        positions of the interpolation points along interp_dim0.
+
+    interp_dim0: scalar
+        is the first dimension to perform interpolation across.
+
+    idx_start0: scalar
+        is the first index value along interp_dim0.
+
+    idx_step0: scalar
+        is the uniform spacing value between subsequent indices along interp_dim0.
+
+    pos1 : af.Array
+        positions of the interpolation points along interp_dim1.
+
+    interp_dim1: scalar
+        is the second dimension to perform interpolation across.
+
+    idx_start1: scalar
+        is the first index value along interp_dim1.
+
+    idx_step1: scalar
+        is the uniform spacing value between subsequent indices along interp_dim1.
+
+    method: optional: af.INTERP. default: af.INTERP.LINEAR.
+            Interpolation method.
+
+    off_grid: optional: scalar. default: 0.0.
+            The value used for positions outside the range.
+
+    output: None or af.Array
+        Optional preallocated output array. If it is a sub-array of an existing af_array,
+        only the corresponding portion of the af_array will be overwritten
+
+    Returns
+    -------
+
+    output: af.Array
+            Values calculated at interpolation points.
+
+    Note
+    -----
+    This holds applicable when x_input/y_input isn't provided:
 
-    if(yp is not None):
-        pos1 = _scale_pos_axis1(y, yp)
+    The initial measurements are assumed to have taken place at equal steps between [(0,0) - [M - 1, N - 1]]
+    where M is the length of the first dimension of `signal`,
+    and N is the length of the second dimension of `signal`.
+    """
+
+    if output is None:
+        output = Array()
+        safe_call(backend.get().af_approx2_uniform(c_pointer(output.arr), signal.arr,
+                                           pos0.arr, c_dim_t(interp_dim0), c_double_t(idx_start0), c_double_t(idx_step0),
+                                           pos1.arr, c_dim_t(interp_dim1), c_double_t(idx_start1), c_double_t(idx_step1),
+                                           method.value, c_float_t(off_grid)))
     else:
-        pos1 = y
+        safe_call(backend.get().af_approx2_uniform_v2(c_pointer(output.arr), signal.arr,
+                                           pos0.arr, c_dim_t(interp_dim0), c_double_t(idx_start0), c_double_t(idx_step0),
+                                           pos1.arr, c_dim_t(interp_dim1), c_double_t(idx_start1), c_double_t(idx_step1),
+                                           method.value, c_float_t(off_grid)))
 
-    safe_call(backend.get().af_approx2(c_pointer(output.arr), signal.arr,
-                                       pos0.arr, pos1.arr, method.value, c_float_t(off_grid)))
     return output
 
+
 def fft(signal, dim0 = None , scale = None):
     """
     Fast Fourier Transform: 1D