1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
|
/*
* This file is part of the micropython-ulab project,
*
* https://github.com/v923z/micropython-ulab
*
* The MIT License (MIT)
*
* Copyright (c) 2019-2021 Zoltán Vörös
*
*/
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "py/obj.h"
#include "py/runtime.h"
#include "py/misc.h"
#include "../ulab.h"
#include "../ulab_tools.h"
#include "transform.h"
#if ULAB_MAX_DIMS > 1
#if ULAB_NUMPY_HAS_DOT
//| def dot(m1: ulab.numpy.ndarray, m2: ulab.numpy.ndarray) -> Union[ulab.numpy.ndarray, _float]:
//| """
//| :param ~ulab.numpy.ndarray m1: a matrix, or a vector
//| :param ~ulab.numpy.ndarray m2: a matrix, or a vector
//|
//| Computes the product of two matrices, or two vectors. In the letter case, the inner product is returned."""
//| ...
//|
mp_obj_t transform_dot(mp_obj_t _m1, mp_obj_t _m2) {
// TODO: should the results be upcast?
// This implements 2D operations only!
if(!mp_obj_is_type(_m1, &ulab_ndarray_type) || !mp_obj_is_type(_m2, &ulab_ndarray_type)) {
mp_raise_TypeError(translate("arguments must be ndarrays"));
}
ndarray_obj_t *m1 = MP_OBJ_TO_PTR(_m1);
ndarray_obj_t *m2 = MP_OBJ_TO_PTR(_m2);
uint8_t *array1 = (uint8_t *)m1->array;
uint8_t *array2 = (uint8_t *)m2->array;
mp_float_t (*func1)(void *) = ndarray_get_float_function(m1->dtype);
mp_float_t (*func2)(void *) = ndarray_get_float_function(m2->dtype);
if(m1->shape[ULAB_MAX_DIMS - 1] != m2->shape[ULAB_MAX_DIMS - m2->ndim]) {
mp_raise_ValueError(translate("dimensions do not match"));
}
uint8_t ndim = MIN(m1->ndim, m2->ndim);
size_t shape1 = m1->ndim == 2 ? m1->shape[ULAB_MAX_DIMS - m1->ndim] : 1;
size_t shape2 = m2->ndim == 2 ? m2->shape[ULAB_MAX_DIMS - 1] : 1;
size_t *shape = NULL;
if(ndim == 2) { // matrix times matrix -> matrix
shape = ndarray_shape_vector(0, 0, shape1, shape2);
} else { // matrix times vector -> vector, vector times vector -> vector (size 1)
shape = ndarray_shape_vector(0, 0, 0, shape1);
}
ndarray_obj_t *results = ndarray_new_dense_ndarray(ndim, shape, NDARRAY_FLOAT);
mp_float_t *rarray = (mp_float_t *)results->array;
for(size_t i=0; i < shape1; i++) { // rows of m1
for(size_t j=0; j < shape2; j++) { // columns of m2
mp_float_t dot = 0.0;
for(size_t k=0; k < m1->shape[ULAB_MAX_DIMS - 1]; k++) {
// (i, k) * (k, j)
dot += func1(array1) * func2(array2);
array1 += m1->strides[ULAB_MAX_DIMS - 1];
array2 += m2->strides[ULAB_MAX_DIMS - m2->ndim];
}
*rarray++ = dot;
array1 -= m1->strides[ULAB_MAX_DIMS - 1] * m1->shape[ULAB_MAX_DIMS - 1];
array2 -= m2->strides[ULAB_MAX_DIMS - m2->ndim] * m2->shape[ULAB_MAX_DIMS - m2->ndim];
array2 += m2->strides[ULAB_MAX_DIMS - 1];
}
array1 += m1->strides[ULAB_MAX_DIMS - m1->ndim];
array2 = m2->array;
}
if((m1->ndim * m2->ndim) == 1) { // return a scalar, if product of two vectors
return mp_obj_new_float(*(--rarray));
} else {
return MP_OBJ_FROM_PTR(results);
}
}
MP_DEFINE_CONST_FUN_OBJ_2(transform_dot_obj, transform_dot);
#endif
#endif
|
