Skip to content

Commit

Permalink
multind: nested functions
Browse files Browse the repository at this point in the history
  • Loading branch information
@uecker
uecker committed Jan 29, 2018
1 parent a0f01db commit 3c1c2eb
Show file tree
Hide file tree
Showing 8 changed files with 92 additions and 144 deletions.
20 changes: 8 additions & 12 deletions src/homodyne.c
View file
Expand Up @@ -34,13 +34,10 @@ static const char help_str[] = "Perform homodyne reconstruction along dimension
struct wdata {

float frac;
float alpha;
int pfdim;
long wdims[DIMS];
long wstrs[DIMS];
complex float* weights;

bool clear;
};


Expand Down Expand Up @@ -74,12 +71,7 @@ static float homodyne_filter(long N, float frac, float alpha, bool clear, long p
return ret;
}

static void comp_weights(void* _data, const long pos[])
{
struct wdata* data = _data;
data->weights[md_calc_offset(DIMS, data->wstrs, pos) / CFL_SIZE]
= homodyne_filter(data->wdims[data->pfdim], data->frac, data->alpha, data->clear, pos[data->pfdim]);
}


static complex float* estimate_phase(struct wdata wdata, unsigned int flags,
unsigned int N, const long dims[N], const complex float* idata, bool center_fft)
Expand Down Expand Up @@ -167,10 +159,14 @@ int main_homodyne(int argc, char* argv[])
md_select_dims(N, MD_BIT(pfdim), wdata.wdims, dims);
md_calc_strides(N, wdata.wstrs, wdata.wdims, CFL_SIZE);
wdata.weights = md_alloc(N, wdata.wdims, CFL_SIZE);
wdata.alpha = alpha;
wdata.clear = clear;

md_loop(N, wdata.wdims, &wdata, comp_weights);
void comp_weights(const long pos[])
{
wdata.weights[md_calc_offset(DIMS, wdata.wstrs, pos) / CFL_SIZE]
= homodyne_filter(wdata.wdims[pfdim], frac, alpha, clear, pos[pfdim]);
}

md_loop(N, wdata.wdims, comp_weights);

long pstrs[N];
long pdims[N];
Expand Down
31 changes: 13 additions & 18 deletions src/num/filter.c
View file
Expand Up @@ -60,39 +60,34 @@ complex float median_complex_float(int N, const complex float ar[N])
}


struct median_s {

long length;
long stride;
};

static void nary_medianz(void* _data, void* ptr[])
{
struct median_s* data = (struct median_s*)_data;

complex float tmp[data->length];

for (long i = 0; i < data->length; i++)
tmp[i] = *((complex float*)(ptr[1] + i * data->stride));

*(complex float*)ptr[0] = median_complex_float(data->length, tmp);
}

void md_medianz2(int D, int M, const long dim[D], const long ostr[D], complex float* optr, const long istr[D], const complex float* iptr)
{
assert(M < D);
const long* nstr[2] = { ostr, istr };
void* nptr[2] = { optr, (void*)iptr };

struct median_s data = { dim[M], istr[M] };
long length = dim[M];
long stride = istr[M];

long dim2[D];
for (int i = 0; i < D; i++)
dim2[i] = dim[i];

dim2[M] = 1;

md_nary(2, D, dim2, nstr, nptr, (void*)&data, &nary_medianz);
void nary_medianz(void* ptr[])
{
complex float tmp[length];

for (long i = 0; i < length; i++)
tmp[i] = *((complex float*)(ptr[1] + i * stride));

*(complex float*)ptr[0] = median_complex_float(length, tmp);
}

md_nary(2, D, dim2, nstr, nptr, nary_medianz);
}

void md_medianz(int D, int M, const long dim[D], complex float* optr, const complex float* iptr)
Expand Down
44 changes: 11 additions & 33 deletions src/num/loop.c
View file
Expand Up @@ -19,53 +19,31 @@

// typedef complex float (*sample_fun_t)(void* _data, const long pos[]);

struct sample_data {

unsigned int N;
const long* strs;
complex float* out;
void* data;
sample_fun_t fun;
};

static void sample_kernel(void* _data, const long pos[])
{
struct sample_data* data = _data;
data->out[md_calc_offset(data->N, data->strs, pos)] = data->fun(data->data, pos);
}

void md_zsample(unsigned int N, const long dims[N], complex float* out, void* data, sample_fun_t fun)
{
struct sample_data sdata;

sdata.N = N;

long strs[N];
md_calc_strides(N, strs, dims, 1); // we use size = 1 here
sdata.strs = strs;

sdata.out = out;
sdata.data = data;
sdata.fun = fun;
void sample_kernel(const long pos[])
{
out[md_calc_offset(N, strs, pos)] = fun(data, pos);
}

md_loop(N, dims, &sdata, sample_kernel);
md_loop(N, dims, sample_kernel);
}


void md_parallel_zsample(unsigned int N, const long dims[N], complex float* out, void* data, sample_fun_t fun)
{
struct sample_data sdata;

sdata.N = N;

long strs[N];
md_calc_strides(N, strs, dims, 1); // we use size = 1 here
sdata.strs = strs;

sdata.out = out;
sdata.data = data;
sdata.fun = fun;
void sample_kernel(const long pos[])
{
out[md_calc_offset(N, strs, pos)] = fun(data, pos);
}

md_parallel_loop(N, dims, ~0u, &sdata, sample_kernel);
md_parallel_loop(N, dims, ~0u, sample_kernel);
}


Expand Down
73 changes: 32 additions & 41 deletions src/num/multind.c
View file
Expand Up @@ -52,11 +52,11 @@
* Generic functions which loops over all dimensions of a set of
* multi-dimensional arrays and calls a given function for each position.
*/
void md_nary(unsigned int C, unsigned int D, const long dim[D], const long* str[C], void* ptr[C], void* data, md_nary_fun_t fun)
void md_nary(unsigned int C, unsigned int D, const long dim[D], const long* str[C], void* ptr[C], md_nary_fun_t fun)
{
if (0 == D) {

fun(data, ptr);
fun(ptr);
return;
}

Expand All @@ -67,7 +67,7 @@ void md_nary(unsigned int C, unsigned int D, const long dim[D], const long* str[
for (unsigned int j = 0; j < C; j++)
moving_ptr[j] = ptr[j] + i * str[j][D - 1];

md_nary(C, D - 1, dim, str, moving_ptr, data, fun);
md_nary(C, D - 1, dim, str, moving_ptr, fun);
}
}

Expand All @@ -79,11 +79,11 @@ void md_nary(unsigned int C, unsigned int D, const long dim[D], const long* str[
* This functions tries to parallelize over the dimensions indicated
* with flags.
*/
void md_parallel_nary(unsigned int C, unsigned int D, const long dim[D], unsigned long flags, const long* str[C], void* ptr[C], void* data, md_nary_fun_t fun)
void md_parallel_nary(unsigned int C, unsigned int D, const long dim[D], unsigned long flags, const long* str[C], void* ptr[C], md_nary_fun_t fun)
{
if (0 == flags) {

md_nary(C, D, dim, str, ptr, data, fun);
md_nary(C, D, dim, str, ptr, fun);
return;
}

Expand Down Expand Up @@ -138,17 +138,17 @@ void md_parallel_nary(unsigned int C, unsigned int D, const long dim[D], unsigne
moving_ptr[j] += iter_i[p] * str[j][parallel_b[p]];
}

md_nary(C, D, dimc, str, moving_ptr, data, fun);
md_nary(C, D, dimc, str, moving_ptr, fun);
}
}



static void md_parallel_loop_r(unsigned int D, unsigned int N, const long dim[static N], unsigned int flags, const long pos[static N], void* data, md_loop_fun_t fun)
static void md_parallel_loop_r(unsigned int D, unsigned int N, const long dim[static N], unsigned int flags, const long pos[static N], md_loop_fun_t fun)
{
if (0 == D) {

fun(data, pos);
fun(pos);
return;
}

Expand All @@ -165,7 +165,7 @@ static void md_parallel_loop_r(unsigned int D, unsigned int N, const long dim[st

pos_copy[D] = i;

md_parallel_loop_r(D, N, dim, flags, pos_copy, data, fun);
md_parallel_loop_r(D, N, dim, flags, pos_copy, fun);
}
}

Expand All @@ -176,39 +176,39 @@ static void md_parallel_loop_r(unsigned int D, unsigned int N, const long dim[st
* Runs fun(data, position) for all position in dim
*
*/
void md_parallel_loop(unsigned int D, const long dim[static D], unsigned long flags, void* data, md_loop_fun_t fun)
void md_parallel_loop(unsigned int D, const long dim[static D], unsigned long flags, md_loop_fun_t fun)
{
long pos[D];
md_parallel_loop_r(D, D, dim, flags, pos, data, fun);
md_parallel_loop_r(D, D, dim, flags, pos, fun);
}



static void md_loop_r(unsigned int D, const long dim[D], long pos[D], void* data, md_loop_fun_t fun)
static void md_loop_r(unsigned int D, const long dim[D], long pos[D], md_loop_fun_t fun)
{
if (0 == D) {

fun(data, pos);
fun(pos);
return;
}

D--;

for (pos[D] = 0; pos[D] < dim[D]; pos[D]++)
md_loop_r(D, dim, pos, data, fun);
md_loop_r(D, dim, pos, fun);
}

/**
* Generic function which loops over all dimensions and calls a given
* function passing the current indices as argument.
*
* Runs fun( data, position ) for all position in dim
* Runs fun( position ) for all position in dim
*
*/
void md_loop(unsigned int D, const long dim[D], void* data, md_loop_fun_t fun)
void md_loop(unsigned int D, const long dim[D], md_loop_fun_t fun)
{
long pos[D];
md_loop_r(D, dim, pos, data, fun);
md_loop_r(D, dim, pos, fun);
}


Expand Down Expand Up @@ -1328,51 +1328,42 @@ bool md_compare(unsigned int D, const long dims[D], const void* src1, const void



struct septrafo_s {

long N;
long str;
void* data;
md_trafo_fun_t fun;
};

static void nary_septrafo(void* _data, void* ptr[])
{
struct septrafo_s* data = (struct septrafo_s*)_data;

data->fun(data->data, data->N, data->str, ptr[0]);
}

static void md_septrafo_r(unsigned int D, unsigned int R, long dimensions[D], unsigned long flags, const long strides[D], void* ptr, md_trafo_fun_t fun, void* _data)
static void md_septrafo_r(unsigned int D, unsigned int R, long dimensions[D], unsigned long flags, const long strides[D], void* ptr, md_trafo_fun_t fun)
{
if (0 == R--)
return;

md_septrafo_r(D, R, dimensions, flags, strides, ptr, fun, _data);
md_septrafo_r(D, R, dimensions, flags, strides, ptr, fun);

if (MD_IS_SET(flags, R)) {

struct septrafo_s data = { dimensions[R], strides[R], _data, fun };
void* nptr[1] = { ptr };
const long* nstrides[1] = { strides };

dimensions[R] = 1; // we made a copy in md_septrafo2
long dimsR = dimensions[R];

void nary_septrafo(void* ptr[])
{
fun(dimsR, strides[R], ptr[0]);
}

//md_nary_parallel(1, D, dimensions, nstrides, nptr, &data, nary_septrafo);
md_nary(1, D, dimensions, nstrides, nptr, &data, nary_septrafo);
dimensions[R] = data.N;
md_nary(1, D, dimensions, nstrides, nptr, nary_septrafo);
dimensions[R] = dimsR;
}
}

/**
* Apply a separable transformation along selected dimensions.
*
*/
void md_septrafo2(unsigned int D, const long dimensions[D], unsigned long flags, const long strides[D], void* ptr, md_trafo_fun_t fun, void* _data)
void md_septrafo2(unsigned int D, const long dimensions[D], unsigned long flags, const long strides[D], void* ptr, md_trafo_fun_t fun)
{
long dimcopy[D];
md_copy_dims(D, dimcopy, dimensions);

md_septrafo_r(D, D, dimcopy, flags, strides, ptr, fun, _data);
md_septrafo_r(D, D, dimcopy, flags, strides, ptr, fun);
}


Expand All @@ -1381,9 +1372,9 @@ void md_septrafo2(unsigned int D, const long dimensions[D], unsigned long flags,
* Apply a separable transformation along selected dimensions.
*
*/
void md_septrafo(unsigned int D, const long dims[D], unsigned long flags, void* ptr, size_t size, md_trafo_fun_t fun, void* _data)
void md_septrafo(unsigned int D, const long dims[D], unsigned long flags, void* ptr, size_t size, md_trafo_fun_t fun)
{
md_septrafo2(D, dims, flags, MD_STRIDES(D, dims, size), ptr, fun, _data);
md_septrafo2(D, dims, flags, MD_STRIDES(D, dims, size), ptr, fun);
}


Expand Down
18 changes: 9 additions & 9 deletions src/num/multind.h
View file
Expand Up @@ -15,20 +15,20 @@

#include "misc/cppwrap.h"

typedef void (*md_nary_fun_t)(void* data, void* ptr[]);
typedef void (*md_trafo_fun_t)(void* data, long N, long str, void* ptr);
typedef void (*md_loop_fun_t)(void* data, const long* pos);
typedef void (*md_nary_fun_t)(void* ptr[]);
typedef void (*md_trafo_fun_t)(long N, long str, void* ptr);
typedef void (*md_loop_fun_t)(const long* pos);


extern void md_nary(unsigned int C, unsigned int D, const long dim[__VLA(D)], const long* str[__VLA(C)], void* ptr[__VLA(C)], void* data, md_nary_fun_t fun);
extern void md_nary(unsigned int C, unsigned int D, const long dim[__VLA(D)], const long* str[__VLA(C)], void* ptr[__VLA(C)], md_nary_fun_t fun);

extern void md_parallel_nary(unsigned int C, unsigned int D, const long dim[__VLA(D)], unsigned long flags, const long* str[__VLA(C)], void* ptr[__VLA(C)], void* data, md_nary_fun_t fun);
extern void md_parallel_loop(unsigned int D, const long dim[__VLA(D)], unsigned long flags, void* data, md_loop_fun_t fun);
extern void md_parallel_nary(unsigned int C, unsigned int D, const long dim[__VLA(D)], unsigned long flags, const long* str[__VLA(C)], void* ptr[__VLA(C)], md_nary_fun_t fun);
extern void md_parallel_loop(unsigned int D, const long dim[__VLA(D)], unsigned long flags, md_loop_fun_t fun);

extern void md_loop(unsigned int D, const long dim[__VLA(D)], void* data, md_loop_fun_t fun);
extern void md_loop(unsigned int D, const long dim[__VLA(D)], md_loop_fun_t fun);

extern void md_septrafo2(unsigned int D, const long dimensions[__VLA(D)], unsigned long flags, const long strides[__VLA(D)], void* ptr, md_trafo_fun_t fun, void* _data);
extern void md_septrafo(unsigned int D, const long dimensions[__VLA(D)], unsigned long flags, void* ptr, size_t size, md_trafo_fun_t fun, void* _data);
extern void md_septrafo2(unsigned int D, const long dimensions[__VLA(D)], unsigned long flags, const long strides[__VLA(D)], void* ptr, md_trafo_fun_t fun);
extern void md_septrafo(unsigned int D, const long dimensions[__VLA(D)], unsigned long flags, void* ptr, size_t size, md_trafo_fun_t fun);


extern void md_clear2(unsigned int D, const long dim[__VLA(D)], const long str[__VLA(D)], void* ptr, size_t size);
Expand Down

0 comments on commit 3c1c2eb

Please sign in to comment.