smile_cloud.cpp

/*
 * smile_cloud.cpp
 */

#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <iostream>

#include "common/log.h"
#include "smile_cloud.hpp"

namespace ug {
namespace smile {

void Cloud::init(int _d, double *_x0, double *_dx, int *_n, int _m, int *_log)
{
	d = _d; m = _m;
	memcpy(x0,     _x0,   sizeof(double)*CLOUD_MAX_DIM);
    memcpy(dx,     _dx,   sizeof(double)*CLOUD_MAX_DIM);
    memcpy(n,       _n,   sizeof(int)*CLOUD_MAX_DIM);
    memcpy(log,   _log,   sizeof(int)*CLOUD_MAX_DIM);
    int s = _n[0];
    for (int i = 1; i < _d; i++)
        s *= _n[i];
    a = (float *)malloc(s * _m * sizeof(float));
}

Cloud::Cloud(int _d, double *_x0, double *_dx, int *_n, int _m, int *_log)
{
	init(_d, _x0, _dx, _n, _m, _log);
}

Cloud::~Cloud()
{
	free(a);
}

float Cloud::read(int *i, int k)
{
	int j, l;
	
	l = i[d-1];
	for (j = d-2; j >= 0; j--)
		l = l * n[j] + i[j]; // how to save the data? 
	l *= m;
	
	return a[l+k];
}

float Cloud::read_o(int *i, int *o, int k)
{
	int j[CLOUD_MAX_DIM], l;
	
	for (l = 0; l < d; l++)
		j[l] = i[l] + o[l];
	
	return read(j, k);
}

/*
 * See section 'Trilinear Interpolation' of
 * http://local.wasp.uwa.edu.au/~pbourke/miscellaneous/interpolation/index.html
 * to get an idea how this playing with bits works even in n dimensions.
 */

double Cloud::interpolate(double *x, int k)
{
	int t, b, i[CLOUD_MAX_DIM], o[CLOUD_MAX_DIM];
	double M, w[CLOUD_MAX_DIM], z, g;
	
	for (t = 0; t < d; t++) {
		/* log if necessary */
		double origin=x[t];
		if (log[t]) {
			if (x[t] > 0.0)
				x[t] = log10(x[t]);
			else
				x[t] = -1e10;
		}
		/* check coordinates and project if necessary */
		M = x0[t] + (n[t] - 1.0) * dx[t];
		if (x[t] < x0[t] || x[t] > M) {
			double xnew=std::min(std::max(x[t], x0[t]), M);
			std::cout << "(W) interpolate_cloud: out of bounds => projecting (" << t << "): x="<<x[t]<<", Min="<<x0[t]<<", Max="<<M<<")!\n";
			//x[t] = std::min(std::max(x[t], x0[t]), M);
			x[t]=xnew;
		}
	}
	
	/* indices & weights */
	for (t = 0; t < d; t++) {
		z = x[t] - x0[t];
		if (dx[t] > 0.0) z /= dx[t];
		i[t] = (int)z;
		w[t] = z - i[t];
	}
	
	/* n-linear interpolation */
	z = 0.0;
	for (b = 0; b < (1 << d); b++) {
		g = 1.0;
		for (t = 0; t < d; t++) {
			o[t] = (b >> (d-1-t)) & 1;
			if (o[t] == 1)
				g *= w[t];
			else
				g *= (1.0-w[t]);
		}
		if (g > 0.0)z+=read_o(i, o, k)*g;
	}

	return z;
}

} // end namespace smile
} // end namespace ug