Skip to content

Commit

Permalink
Merge pull request #5 from v923z/testing
Browse files Browse the repository at this point in the history 
Testing
  • Loading branch information
@v923z
v923z committed Nov 1, 2019
2 parents 74a3b17 + 2af876e commit 43cda00
Show file tree
Hide file tree
Showing 16 changed files with 946 additions and 395 deletions.
6 changes: 3 additions & 3 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,9 +8,9 @@ Documentation can be found under https://micropython-ulab.readthedocs.io/en/late
The source for the manual is in https://github.com/v923z/micropython-ulab/blob/master/docs/ulab-manual.ipynb,
while developer help is in https://github.com/v923z/micropython-ulab/blob/master/docs/ulab.ipynb.

Firmware for pyboard.v.1.1 is updated once in a while, and can be downloaded
from https://github.com/v923z/micropython-ulab/releases, otherwise, it can
be compiled from the source by following the steps
Firmware for pyboard.v.1.1, as well as for PYBD_SF6 is updated once
in a while, and can be downloaded from https://github.com/v923z/micropython-ulab/releases,
otherwise, it can be compiled from the source by following the steps
https://micropython-usermod.readthedocs.io/en/latest/usermods_05.html#compiling-our-module.


Expand Down
20 changes: 10 additions & 10 deletions code/fft.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,13 +26,13 @@ enum FFT_TYPE {
FFT_SPECTRUM,
};

void fft_kernel(float *real, float *imag, int n, int isign) {
void fft_kernel(mp_float_t *real, mp_float_t *imag, int n, int isign) {
// This is basically a modification of four1 from Numerical Recipes
// The main difference is that this function takes two arrays, one
// for the real, and one for the imaginary parts.
int j, m, mmax, istep;
float tempr, tempi;
float wtemp, wr, wpr, wpi, wi, theta;
mp_float_t tempr, tempi;
mp_float_t wtemp, wr, wpr, wpi, wi, theta;

j = 0;
for(int i = 0; i < n; i++) {
Expand All @@ -52,9 +52,9 @@ void fft_kernel(float *real, float *imag, int n, int isign) {
while (n > mmax) {
istep = mmax << 1;
theta = -1.0*isign*6.28318530717959/istep;
wtemp = sinf(0.5 * theta);
wtemp = MICROPY_FLOAT_C_FUN(sin)(0.5 * theta);
wpr = -2.0 * wtemp * wtemp;
wpi = sinf(theta);
wpi = MICROPY_FLOAT_C_FUN(sin)(theta);
wr = 1.0;
wi = 0.0;
for(m = 0; m < mmax; m++) {
Expand Down Expand Up @@ -92,27 +92,27 @@ mp_obj_t fft_fft_ifft_spectrum(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_im,
}

ndarray_obj_t *out_re = create_new_ndarray(1, len, NDARRAY_FLOAT);
float *data_re = (float *)out_re->array->items;
mp_float_t *data_re = (mp_float_t *)out_re->array->items;

if(re->array->typecode == NDARRAY_FLOAT) {
// By treating this case separately, we can save a bit of time.
// I don't know if it is worthwhile, though...
memcpy((float *)out_re->array->items, (float *)re->array->items, re->bytes);
memcpy((mp_float_t *)out_re->array->items, (mp_float_t *)re->array->items, re->bytes);
} else {
for(size_t i=0; i < len; i++) {
data_re[i] = ndarray_get_float_value(re->array->items, re->array->typecode, i);
}
}
ndarray_obj_t *out_im = create_new_ndarray(1, len, NDARRAY_FLOAT);
float *data_im = (float *)out_im->array->items;
mp_float_t *data_im = (mp_float_t *)out_im->array->items;

if(n_args == 2) {
ndarray_obj_t *im = MP_OBJ_TO_PTR(arg_im);
if (re->array->len != im->array->len) {
mp_raise_ValueError("real and imaginary parts must be of equal length");
}
if(im->array->typecode == NDARRAY_FLOAT) {
memcpy((float *)out_im->array->items, (float *)im->array->items, im->bytes);
memcpy((mp_float_t *)out_im->array->items, (mp_float_t *)im->array->items, im->bytes);
} else {
for(size_t i=0; i < len; i++) {
data_im[i] = ndarray_get_float_value(im->array->items, im->array->typecode, i);
Expand All @@ -123,7 +123,7 @@ mp_obj_t fft_fft_ifft_spectrum(size_t n_args, mp_obj_t arg_re, mp_obj_t arg_im,
fft_kernel(data_re, data_im, len, 1);
if(type == FFT_SPECTRUM) {
for(size_t i=0; i < len; i++) {
data_re[i] = sqrtf(data_re[i]*data_re[i] + data_im[i]*data_im[i]);
data_re[i] = MICROPY_FLOAT_C_FUN(sqrt)(data_re[i]*data_re[i] + data_im[i]*data_im[i]);
}
}
} else { // inverse transform
Expand Down
57 changes: 29 additions & 28 deletions code/linalg.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,24 +108,24 @@ mp_obj_t linalg_size(size_t n_args, const mp_obj_t *pos_args, mp_map_t *kw_args)
}
}

bool linalg_invert_matrix(float *data, size_t N) {
bool linalg_invert_matrix(mp_float_t *data, size_t N) {
// returns true, of the inversion was successful,
// false, if the matrix is singular

// initially, this is the unit matrix: the contents of this matrix is what
// will be returned after all the transformations
float *unit = m_new(float, N*N);
mp_float_t *unit = m_new(mp_float_t, N*N);

float elem = 1.0;
mp_float_t elem = 1.0;
// initialise the unit matrix
memset(unit, 0, sizeof(float)*N*N);
memset(unit, 0, sizeof(mp_float_t)*N*N);
for(size_t m=0; m < N; m++) {
memcpy(&unit[m*(N+1)], &elem, sizeof(float));
memcpy(&unit[m*(N+1)], &elem, sizeof(mp_float_t));
}
for(size_t m=0; m < N; m++){
// this could be faster with ((c < epsilon) && (c > -epsilon))
if(abs(data[m*(N+1)]) < epsilon) {
m_del(float, unit, N*N);
m_del(mp_float_t, unit, N*N);
return false;
}
for(size_t n=0; n < N; n++){
Expand All @@ -145,8 +145,8 @@ bool linalg_invert_matrix(float *data, size_t N) {
unit[N*m+n] /= elem;
}
}
memcpy(data, unit, sizeof(float)*N*N);
m_del(float, unit, N*N);
memcpy(data, unit, sizeof(mp_float_t)*N*N);
m_del(mp_float_t, unit, N*N);
return true;
}

Expand All @@ -163,21 +163,21 @@ mp_obj_t linalg_inv(mp_obj_t o_in) {
mp_raise_ValueError("only square matrices can be inverted");
}
ndarray_obj_t *inverted = create_new_ndarray(o->m, o->n, NDARRAY_FLOAT);
float *data = (float *)inverted->array->items;
mp_float_t *data = (mp_float_t *)inverted->array->items;
mp_obj_t elem;
for(size_t m=0; m < o->m; m++) { // rows first
for(size_t n=0; n < o->n; n++) { // columns next
// this could, perhaps, be done in single line...
// On the other hand, we probably spend little time here
elem = mp_binary_get_val_array(o->array->typecode, o->array->items, m*o->n+n);
data[m*o->n+n] = (float)mp_obj_get_float(elem);
data[m*o->n+n] = (mp_float_t)mp_obj_get_float(elem);
}
}

if(!linalg_invert_matrix(data, o->m)) {
// TODO: I am not sure this is needed here. Otherwise,
// how should we free up the unused RAM of inverted?
m_del(float, inverted->array->items, o->n*o->n);
m_del(mp_float_t, inverted->array->items, o->n*o->n);
mp_raise_ValueError("input matrix is singular");
}
return MP_OBJ_FROM_PTR(inverted);
Expand All @@ -192,8 +192,8 @@ mp_obj_t linalg_dot(mp_obj_t _m1, mp_obj_t _m2) {
}
// TODO: numpy uses upcasting here
ndarray_obj_t *out = create_new_ndarray(m1->m, m2->n, NDARRAY_FLOAT);
float *outdata = (float *)out->array->items;
float sum, v1, v2;
mp_float_t *outdata = (mp_float_t *)out->array->items;
mp_float_t sum, v1, v2;
for(size_t i=0; i < m1->n; i++) {
for(size_t j=0; j < m2->m; j++) {
sum = 0.0;
Expand Down Expand Up @@ -301,14 +301,14 @@ mp_obj_t linalg_det(mp_obj_t oin) {
mp_raise_ValueError("input must be square matrix");
}

float *tmp = m_new(float, in->n*in->n);
mp_float_t *tmp = m_new(mp_float_t, in->n*in->n);
for(size_t i=0; i < in->array->len; i++){
tmp[i] = ndarray_get_float_value(in->array->items, in->array->typecode, i);
}
float c;
mp_float_t c;
for(size_t m=0; m < in->m-1; m++){
if(abs(tmp[m*(in->n+1)]) < epsilon) {
m_del(float, tmp, in->n*in->n);
m_del(mp_float_t, tmp, in->n*in->n);
return mp_obj_new_float(0.0);
}
for(size_t n=0; n < in->n; n++){
Expand All @@ -320,12 +320,12 @@ mp_obj_t linalg_det(mp_obj_t oin) {
}
}
}
float det = 1.0;
mp_float_t det = 1.0;

for(size_t m=0; m < in->m; m++){
det *= tmp[m*(in->n+1)];
}
m_del(float, tmp, in->n*in->n);
m_del(mp_float_t, tmp, in->n*in->n);
return mp_obj_new_float(det);
}

Expand All @@ -337,14 +337,15 @@ mp_obj_t linalg_eig(mp_obj_t oin) {
if(in->m != in->n) {
mp_raise_ValueError("input must be square matrix");
}
float *array = m_new(float, in->array->len);
mp_float_t *array = m_new(mp_float_t, in->array->len);
for(size_t i=0; i < in->array->len; i++) {
array[i] = ndarray_get_float_value(in->array->items, in->array->typecode, i);
}
// make sure the matrix is symmetric
for(size_t m=0; m < in->m; m++) {
for(size_t n=m+1; n < in->n; n++) {
// compare entry (m, n) to (n, m)
// TODO: this must probably be scaled!
if(epsilon < abs(array[m*in->n + n] - array[n*in->n + m])) {
mp_raise_ValueError("input matrix is asymmetric");
}
Expand All @@ -354,12 +355,12 @@ mp_obj_t linalg_eig(mp_obj_t oin) {
// if we got this far, then the matrix will be symmetric

ndarray_obj_t *eigenvectors = create_new_ndarray(in->m, in->n, NDARRAY_FLOAT);
float *eigvectors = (float *)eigenvectors->array->items;
mp_float_t *eigvectors = (mp_float_t *)eigenvectors->array->items;
// start out with the unit matrix
for(size_t m=0; m < in->m; m++) {
eigvectors[m*(in->n+1)] = 1.0;
}
float largest, w, t, c, s, tau, aMk, aNk, vm, vn;
mp_float_t largest, w, t, c, s, tau, aMk, aNk, vm, vn;
size_t M, N;
size_t iterations = JACOBI_MAX*in->n*in->n;
do {
Expand Down Expand Up @@ -387,12 +388,12 @@ mp_obj_t linalg_eig(mp_obj_t oin) {
// The following if/else chooses the smaller absolute value for the tangent
// of the rotation angle. Going with the smaller should be numerically stabler.
if(w > 0) {
t = sqrtf(w*w + 1.0) - w;
t = MICROPY_FLOAT_C_FUN(sqrt)(w*w + 1.0) - w;
} else {
t = -1.0*(sqrtf(w*w + 1.0) + w);
t = -1.0*(MICROPY_FLOAT_C_FUN(sqrt)(w*w + 1.0) + w);
}
s = t / sqrtf(t*t + 1.0); // the sine of the rotation angle
c = 1.0 / sqrtf(t*t + 1.0); // the cosine of the rotation angle
s = t / MICROPY_FLOAT_C_FUN(sqrt)(t*t + 1.0); // the sine of the rotation angle
c = 1.0 / MICROPY_FLOAT_C_FUN(sqrt)(t*t + 1.0); // the cosine of the rotation angle
tau = (1.0-c)/s; // this is equal to the tangent of the half of the rotation angle

// at this point, we have the rotation angles, so we can transform the matrix
Expand Down Expand Up @@ -438,15 +439,15 @@ mp_obj_t linalg_eig(mp_obj_t oin) {

if(iterations == 0) {
// the computation did not converge; numpy raises LinAlgError
m_del(float, array, in->array->len);
m_del(mp_float_t, array, in->array->len);
mp_raise_ValueError("iterations did not converge");
}
ndarray_obj_t *eigenvalues = create_new_ndarray(1, in->n, NDARRAY_FLOAT);
float *eigvalues = (float *)eigenvalues->array->items;
mp_float_t *eigvalues = (mp_float_t *)eigenvalues->array->items;
for(size_t i=0; i < in->n; i++) {
eigvalues[i] = array[i*(in->n+1)];
}
m_del(float, array, in->array->len);
m_del(mp_float_t, array, in->array->len);

mp_obj_tuple_t *tuple = MP_OBJ_TO_PTR(mp_obj_new_tuple(2, NULL));
tuple->items[0] = MP_OBJ_FROM_PTR(eigenvalues);
Expand Down
10 changes: 8 additions & 2 deletions code/linalg.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,13 +14,19 @@
#include "ndarray.h"

#define SWAP(t, a, b) { t tmp = a; a = b; b = tmp; }
#define epsilon 1e-6

#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT
#define epsilon 1.2e-7
#elif MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE
#define epsilon 2.3e-16
#endif

#define JACOBI_MAX 20

mp_obj_t linalg_transpose(mp_obj_t );
mp_obj_t linalg_reshape(mp_obj_t , mp_obj_t );
mp_obj_t linalg_size(size_t , const mp_obj_t *, mp_map_t *);
bool linalg_invert_matrix(float *, size_t );
bool linalg_invert_matrix(mp_float_t *, size_t );
mp_obj_t linalg_inv(mp_obj_t );
mp_obj_t linalg_dot(mp_obj_t , mp_obj_t );
mp_obj_t linalg_zeros(size_t , const mp_obj_t *, mp_map_t *);
Expand Down

0 comments on commit 43cda00

Please sign in to comment.