Added test for normal distributions.

include/ccs/distribution.h

@@ -35,16 +35,16 @@ ccs_create_distribution(ccs_distribution_type_t distribution_type,
 
 extern ccs_error_t
 _ccs_create_normal_distribution(ccs_data_type_t     data_type,
-                                ccs_value_t         mu,
+                                ccs_float_t         mu,
                                 ccs_float_t         sigma,
                                 ccs_scale_type_t    scale,
                                 ccs_value_t         quantization,
                                 ccs_distribution_t *distribution_ret);
 #define ccs_create_normal_distribution(t, m, s, sc, q, d) \
-	_ccs_create_normal_distribution(t, (ccs_value_t)(m), s, sc, (ccs_value_t)(q), d)
+	_ccs_create_normal_distribution(t, m, s, sc, (ccs_value_t)(q), d)
 
 extern ccs_error_t
-ccs_create_normal_int_distribution(ccs_int_t           mu,
+ccs_create_normal_int_distribution(ccs_float_t         mu,
                                    ccs_float_t         sigma,
                                    ccs_scale_type_t    scale,
                                    ccs_int_t           quantization,

src/distribution_normal.c

@@ -6,7 +6,7 @@
 
 struct _ccs_distribution_normal_data_s {
 	_ccs_distribution_common_data_t common_data;
-	ccs_value_t                     mu;
+	ccs_float_t                     mu;
 	ccs_float_t                     sigma;
 	int                             quantize;
 };
@@ -61,9 +61,9 @@ _ccs_distribution_normal_get_parameters(_ccs_distribution_data_t *data,
 		*num_parameters_ret = 2;
 
 	if (parameters) {
-		parameters[0].type = d->common_data.data_type;
-		parameters[0].value = d->mu;
-		parameters[1].type = d->common_data.data_type;
+		parameters[0].type = CCS_FLOAT;
+		parameters[0].value.f = d->mu;
+		parameters[1].type = CCS_FLOAT;
 		parameters[1].value.f = d->sigma;
 	}
 	return CCS_SUCCESS;
@@ -79,35 +79,66 @@ _ccs_distribution_normal_samples(_ccs_distribution_data_t *data,
 	const ccs_data_type_t  data_type      = d->common_data.data_type;
 	const ccs_scale_type_t scale_type     = d->common_data.scale_type;
 	const ccs_value_t      quantization   = d->common_data.quantization;
-	const ccs_value_t      mu             = d->mu;
+	const ccs_float_t      mu             = d->mu;
 	const ccs_float_t      sigma          = d->sigma;
 	const int              quantize       = d->quantize;
 	gsl_rng *grng;
 	ccs_error_t err = ccs_rng_get_gsl_rng(rng, &grng);
 	if (err)
 		return err;
-	if (data_type == CCS_FLOAT) {
-		for (i = 0; i < num_values; i++) {
-			values[i].f = gsl_ran_gaussian(grng, sigma) + mu.f;
+	if (data_type == CCS_FLOAT)
+		if (scale_type == CCS_LOGARITHMIC && quantize) {
+			double lq = log(quantization.f*0.5);
+			if (mu - lq >= 0.0)
+				//at least 50% chance to get a valid value
+				for (i = 0; i < num_values; i++)
+					do {
+						values[i].f = gsl_ran_gaussian(grng, sigma) + mu;
+					} while (values[i].f < lq);
+			else
+				//use tail distribution
+				for (i = 0; i < num_values; i++)
+					values[i].f = gsl_ran_gaussian_tail(grng, lq - mu, sigma) + mu;
 		}
-	} else {
-		for (i = 0; i < num_values; i++) {
-			values[i].f = gsl_ran_gaussian(grng, sigma) + mu.i;
+		else
+			for (i = 0; i < num_values; i++)
+				values[i].f = gsl_ran_gaussian(grng, sigma) + mu;
+	else {
+		if (scale_type == CCS_LOGARITHMIC) {
+			double lq;
+			if (quantize)
+				lq = log(quantization.i*0.5);
+			else
+				lq = log(0.5);
+			if (mu - lq >= 0.0)
+				for (i = 0; i < num_values; i++)
+					do {
+						values[i].f = gsl_ran_gaussian(grng, sigma) + mu;
+					} while (values[i].f < lq);
+			else
+				for (i = 0; i < num_values; i++)
+					values[i].f = gsl_ran_gaussian_tail(grng, lq - mu, sigma) + mu;
 		}
+		else
+			for (i = 0; i < num_values; i++)
+				values[i].f = gsl_ran_gaussian(grng, sigma) + mu;
 	}
 
 	if (scale_type == CCS_LOGARITHMIC)
 		for (i = 0; i < num_values; i++)
 			values[i].f = exp(values[i].f);
 	if (data_type == CCS_FLOAT) {
-		if (quantize)
+		if (quantize) {
+			ccs_float_t rquantization = 1.0 / quantization.f;
 			for (i = 0; i < num_values; i++)
-				values[i].f = round(values[i].f/quantization.f) * quantization.f;
+				values[i].f = round(values[i].f * rquantization) * quantization.f;
+		}
 	} else {
-		if (quantize)
+		if (quantize) {
+			ccs_float_t rquantization = 1.0 / quantization.i;
 			for (i = 0; i < num_values; i++)
-				values[i].i = (ccs_int_t)round(values[i].f/quantization.i) * quantization.i;
-		else
+				values[i].i = (ccs_int_t)round(values[i].f * rquantization) * quantization.i;
+		} else
 			for (i = 0; i < num_values; i++)
 				values[i].i = round(values[i].f);
 	}
@@ -116,17 +147,21 @@ _ccs_distribution_normal_samples(_ccs_distribution_data_t *data,
 
 extern ccs_error_t
 _ccs_create_normal_distribution(ccs_data_type_t     data_type,
-                               ccs_value_t         mu,
-                               ccs_float_t         sigma,
-                               ccs_scale_type_t    scale_type,
-                               ccs_value_t         quantization,
-                               ccs_distribution_t *distribution_ret) {
+                                ccs_float_t         mu,
+                                ccs_float_t         sigma,
+                                ccs_scale_type_t    scale_type,
+                                ccs_value_t         quantization,
+                                ccs_distribution_t *distribution_ret) {
 	if (!distribution_ret)
 		return -CCS_INVALID_VALUE;
 	if (data_type != CCS_FLOAT && data_type != CCS_INTEGER)
 		return -CCS_INVALID_TYPE;
 	if (scale_type != CCS_LINEAR && scale_type != CCS_LOGARITHMIC)
 		return -CCS_INVALID_SCALE;
+	if (data_type == CCS_INTEGER && quantization.i < 0 )
+		return -CCS_INVALID_VALUE;
+	if (data_type == CCS_FLOAT && quantization.f < 0.0 )
+		return -CCS_INVALID_VALUE;
 	uintptr_t mem = (uintptr_t)calloc(1, sizeof(struct _ccs_distribution_s) + sizeof(_ccs_distribution_normal_data_t));
 
 	if (!mem)
@@ -152,4 +187,27 @@ _ccs_create_normal_distribution(ccs_data_type_t     data_type,
 	return CCS_SUCCESS;
 }
 
+extern ccs_error_t
+ccs_normal_distribution_get_parameters(ccs_distribution_t  distribution,
+                                       ccs_datum_t        *mu,
+                                       ccs_datum_t        *sigma) {
+	if (!distribution || distribution->obj.type != CCS_DISTRIBUTION)
+		return -CCS_INVALID_OBJECT;
+	if (!distribution->data || ((_ccs_distribution_common_data_t*)distribution->data)->type != CCS_NORMAL)
+		return -CCS_INVALID_OBJECT;
+	if (!mu && !sigma)
+		return -CCS_INVALID_VALUE;
+	_ccs_distribution_normal_data_t * data = (_ccs_distribution_normal_data_t *)distribution->data;
+
+	if (mu) {
+		mu->value.f = data->mu;
+		mu->type = CCS_FLOAT;
+	}
+	if (sigma) {
+		sigma->value.f = data->sigma;
+		sigma->type = CCS_FLOAT;
+	}
+	return CCS_SUCCESS;
+}
+
 

src/distribution_uniform.c

@@ -204,7 +204,9 @@ extern ccs_error_t
 ccs_uniform_distribution_get_parameters(ccs_distribution_t  distribution,
                                         ccs_datum_t        *lower,
                                         ccs_datum_t        *upper) {
-	if (!distribution || !distribution->data)
+	if (!distribution || distribution->obj.type != CCS_DISTRIBUTION)
+		return -CCS_INVALID_OBJECT;
+	if (!distribution->data || ((_ccs_distribution_common_data_t*)distribution->data)->type != CCS_UNIFORM)
 		return -CCS_INVALID_OBJECT;
 	if (!lower && !upper)
 		return -CCS_INVALID_VALUE;

tests/Makefile.am

@@ -5,7 +5,8 @@ AM_LDFLAGS = ../src/libcconfigspace.la
 
 RNG_TESTS = \
 		test_rng \
-		test_distribution
+		test_uniform_distribution \
+		test_normal_distribution
 
 # unit tests
 UNIT_TESTS = \

tests/test_normal_distribution.c

+#include <math.h>
+#include <stdlib.h>
+#include <assert.h>
+#include <cconfigspace.h>
+#include <gsl/gsl_statistics.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_cdf.h>
+
+static void test_create_normal_distribution() {
+	ccs_distribution_t      distrib = NULL;
+	ccs_error_t             err = CCS_SUCCESS;
+	int32_t                 refcount;
+	size_t                  num_parameters;
+	ccs_object_type_t       otype;
+	ccs_distribution_type_t dtype;
+	ccs_scale_type_t        stype;
+	ccs_data_type_t         data_type;
+	ccs_datum_t             parameters[2], quantization, mu, sigma;
+
+	err = ccs_create_normal_distribution(
+		CCS_FLOAT,
+		1.0,
+		2.0,
+		CCS_LINEAR,
+		0.0,
+		&distrib);
+	assert( err == CCS_SUCCESS );
+
+	err = ccs_object_get_type(distrib, &otype);
+	assert( err == CCS_SUCCESS );
+	assert( otype == CCS_DISTRIBUTION );
+
+	err = ccs_distribution_get_type(distrib, &dtype);
+	assert( err == CCS_SUCCESS );
+	assert( dtype == CCS_NORMAL );
+
+	err = ccs_distribution_get_data_type(distrib, &data_type);
+	assert( err == CCS_SUCCESS );
+	assert( data_type == CCS_FLOAT );
+
+	err = ccs_distribution_get_scale_type(distrib, &stype);
+	assert( err == CCS_SUCCESS );
+	assert( stype == CCS_LINEAR );
+
+	err = ccs_distribution_get_quantization(distrib, &quantization);
+	assert( err == CCS_SUCCESS );
+	assert( quantization.type == CCS_FLOAT );
+	assert( quantization.value.f == 0.0 );
+
+	err = ccs_distribution_get_num_parameters(distrib, &num_parameters);
+	assert( err == CCS_SUCCESS );
+	assert( num_parameters == 2 );
+
+	err = ccs_distribution_get_parameters(distrib, 2, parameters, &num_parameters);
+	assert( err == CCS_SUCCESS );
+	assert( num_parameters == 2 );
+	assert( parameters[0].type == CCS_FLOAT );
+	assert( parameters[0].value.f == 1.0 );
+	assert( parameters[1].type == CCS_FLOAT );
+	assert( parameters[1].value.f == 2.0 );
+
+	err = ccs_normal_distribution_get_parameters(distrib, &mu, &sigma);
+	assert( err == CCS_SUCCESS );
+	assert( mu.type == CCS_FLOAT );
+	assert( mu.value.f == 1.0 );
+	assert( sigma.type == CCS_FLOAT );
+	assert( sigma.value.f == 2.0 );
+
+	err = ccs_object_get_refcount(distrib, &refcount);
+	assert( err == CCS_SUCCESS );
+	assert( refcount == 1 );
+
+	err = ccs_release_object(distrib);
+	assert( err == CCS_SUCCESS );
+}
+
+static void test_create_normal_distribution_errors() {
+	ccs_distribution_t      distrib = NULL;
+	ccs_error_t             err = CCS_SUCCESS;
+
+	// check wrong data_type
+	err = ccs_create_normal_distribution(
+		CCS_OBJECT,
+		1.0,
+		2.0,
+		CCS_LINEAR,
+		0.0,
+		&distrib);
+	assert( err == -CCS_INVALID_TYPE );
+
+	// check wrong data_type
+	err = ccs_create_normal_distribution(
+		CCS_FLOAT,
+		1.0,
+		2.0,
+		0xdeadbeef,
+		0.0,
+		&distrib);
+	assert( err == -CCS_INVALID_SCALE );
+
+	// check wrong quantization
+	err = ccs_create_normal_distribution(
+		CCS_FLOAT,
+		1.0,
+		2.0,
+		CCS_LINEAR,
+		-1.0,
+		&distrib);
+	assert( err == -CCS_INVALID_VALUE );
+
+	// check wrong pointer
+	err = ccs_create_normal_distribution(
+		CCS_FLOAT,
+		1.0,
+		2.0,
+		CCS_LINEAR,
+		0.0,
+		NULL);
+	assert( err == -CCS_INVALID_VALUE );
+
+
+}
+
+static void to_float(ccs_value_t * values, size_t length) {
+	for (size_t i = 0; i < length; i++) {
+		values[i].f = values[i].i;
+	}
+}
+
+static void to_log(ccs_value_t * values, size_t length) {
+	for (size_t i = 0; i < length; i++) {
+		values[i].f = log(values[i].f);
+	}
+}
+
+static void test_normal_distribution_int() {
+	ccs_distribution_t distrib = NULL;
+	ccs_rng_t          rng = NULL;
+	ccs_error_t        err = CCS_SUCCESS;
+	const size_t       num_samples = 10000;
+	const ccs_float_t  mu = 1;
+	const ccs_float_t  sigma = 2;
+	ccs_value_t        samples[num_samples];
+	double mean, sig;
+
+	err = ccs_rng_create(&rng);
+	assert( err == CCS_SUCCESS );
+	err = ccs_create_normal_distribution(
+		CCS_INTEGER,
+		mu,
+		sigma,
+		CCS_LINEAR,
+		0L,
+		&distrib);
+	assert( err == CCS_SUCCESS );
+
+	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
+	assert( err == CCS_SUCCESS );
+
+	to_float(samples, num_samples);
+	mean = gsl_stats_mean((double*)samples, 1, num_samples);
+	//fprintf(stderr, "mu: %lf \n", mean);
+	assert( mean < mu + 0.1 );
+	assert( mean > mu - 0.1 );
+	sig    = gsl_stats_sd_m((double*)samples, 1, num_samples, mu);
+	//fprintf(stderr, "sig: %lf \n", sig);
+	assert( sig < sigma + 0.1 );
+	assert( sig > sigma - 0.1 );
+
+	err = ccs_release_object(distrib);
+	assert( err == CCS_SUCCESS );
+}
+
+static void test_normal_distribution_float() {
+	ccs_distribution_t distrib = NULL;
+	ccs_rng_t          rng = NULL;
+	ccs_error_t        err = CCS_SUCCESS;
+	const size_t       num_samples = 10000;
+	const ccs_float_t  mu = 1;
+	const ccs_float_t  sigma = 2;
+	ccs_value_t        samples[num_samples];
+	double mean, sig;
+
+	err = ccs_rng_create(&rng);
+	assert( err == CCS_SUCCESS );
+	err = ccs_create_normal_distribution(
+		CCS_FLOAT,
+		mu,
+		sigma,
+		CCS_LINEAR,
+		0.0,
+		&distrib);
+	assert( err == CCS_SUCCESS );
+
+	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
+	assert( err == CCS_SUCCESS );
+
+	mean = gsl_stats_mean((double*)samples, 1, num_samples);
+	//fprintf(stderr, "mu: %lf \n", mean);
+	assert( mean < mu + 0.1 );
+	assert( mean > mu - 0.1 );
+	sig    = gsl_stats_sd_m((double*)samples, 1, num_samples, mu);
+	//fprintf(stderr, "sig: %lf \n", sig);
+	assert( sig < sigma + 0.1 );
+	assert( sig > sigma - 0.1 );
+
+	err = ccs_release_object(distrib);
+	assert( err == CCS_SUCCESS );
+}
+
+static void test_normal_distribution_int_log() {
+	ccs_distribution_t distrib = NULL;
+	ccs_rng_t          rng = NULL;
+	ccs_error_t        err = CCS_SUCCESS;
+	const size_t       num_samples = 10000;
+	const ccs_float_t  mu = 1;
+	const ccs_float_t  sigma = 2;
+	ccs_value_t        samples[num_samples];
+	double mean, sig;
+	double tmean, tsigma, alpha, zee, pdfa;
+
+	err = ccs_rng_create(&rng);
+	assert( err == CCS_SUCCESS );
+	err = ccs_create_normal_distribution(
+		CCS_INTEGER,
+		mu,
+		sigma,
+		CCS_LOGARITHMIC,
+		0L,
+		&distrib);
+	assert( err == CCS_SUCCESS );
+
+	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
+	assert( err == CCS_SUCCESS );
+
+	to_float(samples, num_samples);
+	to_log(samples, num_samples);
+	mean = gsl_stats_mean((double*)samples, 1, num_samples);
+	//fprintf(stderr, "mean: %lf \n", mean);
+	// cutoff at 0.0 to have exp(v) >= 1
+	// see https://en.wikipedia.org/wiki/Truncated_normal_distribution
+	alpha = (log(0.5) - mu)/sigma;
+	zee = (1.0 - gsl_cdf_ugaussian_P(alpha));
+	pdfa = gsl_ran_ugaussian_pdf(alpha);
+	tmean = mu + pdfa * sigma / zee;
+	//fprintf(stderr, "tmean: %lf \n", tmean);
+	assert( mean < tmean + 0.1 );
+	assert( mean > tmean - 0.1 );
+	sig    = gsl_stats_sd_m((double*)samples, 1, num_samples, tmean);
+	//fprintf(stderr, "sig: %lf \n", sig);
+	tsigma = sqrt( sigma * sigma * ( 1.0 + alpha * pdfa / zee - ( pdfa * pdfa )/( zee * zee ) ) );
+	//fprintf(stderr, "tsig: %lf \n", tsigma);
+	assert( sig < tsigma + 0.1 );
+	assert( sig > tsigma - 1.1 );
+
+	err = ccs_release_object(distrib);
+	assert( err == CCS_SUCCESS );
+}
+
+static void test_normal_distribution_float_log() {
+	ccs_distribution_t distrib = NULL;
+	ccs_rng_t          rng = NULL;
+	ccs_error_t        err = CCS_SUCCESS;
+	const size_t       num_samples = 10000;
+	const ccs_float_t  mu = 1;
+	const ccs_float_t  sigma = 2;
+	ccs_value_t        samples[num_samples];
+	double mean, sig;
+
+	err = ccs_rng_create(&rng);
+	assert( err == CCS_SUCCESS );
+	err = ccs_create_normal_distribution(
+		CCS_FLOAT,
+		mu,
+		sigma,
+		CCS_LOGARITHMIC,
+		0.0,
+		&distrib);
+	assert( err == CCS_SUCCESS );
+
+	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
+	assert( err == CCS_SUCCESS );
+
+	to_log(samples, num_samples);
+	mean = gsl_stats_mean((double*)samples, 1, num_samples);
+	//fprintf(stderr, "mu: %lf \n", mean);
+	assert( mean < mu + 0.1 );
+	assert( mean > mu - 0.1 );
+	sig    = gsl_stats_sd_m((double*)samples, 1, num_samples, mu);
+	//fprintf(stderr, "sig: %lf \n", sig);
+	assert( sig < sigma + 0.1 );
+	assert( sig > sigma - 0.1 );
+
+	err = ccs_release_object(distrib);
+	assert( err == CCS_SUCCESS );
+}
+
+static void test_normal_distribution_int_quantize() {
+	ccs_distribution_t distrib = NULL;
+	ccs_rng_t          rng = NULL;
+	ccs_error_t        err = CCS_SUCCESS;
+	const size_t       num_samples = 10000;
+	const ccs_float_t  mu = 1;
+	const ccs_float_t  sigma = 2;
+	ccs_value_t        samples[num_samples];
+	double mean, sig;
+
+	err = ccs_rng_create(&rng);
+	assert( err == CCS_SUCCESS );
+	err = ccs_create_normal_distribution(
+		CCS_INTEGER,
+		mu,
+		sigma,
+		CCS_LINEAR,
+		2L,
+		&distrib);
+	assert( err == CCS_SUCCESS );
+
+	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
+	assert( err == CCS_SUCCESS );
+
+	to_float(samples, num_samples);
+	mean = gsl_stats_mean((double*)samples, 1, num_samples);
+	//fprintf(stderr, "mu: %lf \n", mean);
+	assert( mean < mu + 0.1 );
+	assert( mean > mu - 0.1 );
+	sig    = gsl_stats_sd_m((double*)samples, 1, num_samples, mu);
+	//fprintf(stderr, "sig: %lf \n", sig);
+	assert( sig < sigma + 0.1 );
+	assert( sig > sigma - 0.1 );
+
+	err = ccs_release_object(distrib);
+	assert( err == CCS_SUCCESS );
+}
+
+static void test_normal_distribution_float_quantize() {
+	ccs_distribution_t distrib = NULL;
+	ccs_rng_t          rng = NULL;
+	ccs_error_t        err = CCS_SUCCESS;
+	const size_t       num_samples = 10000;
+	const ccs_float_t  mu = 1;
+	const ccs_float_t  sigma = 2;
+	ccs_value_t        samples[num_samples];
+	double mean, sig;
+
+	err = ccs_rng_create(&rng);
+	assert( err == CCS_SUCCESS );
+	err = ccs_create_normal_distribution(
+		CCS_FLOAT,
+		mu,
+		sigma,
+		CCS_LINEAR,
+		0.2,
+		&distrib);
+	assert( err == CCS_SUCCESS );
+
+	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
+	assert( err == CCS_SUCCESS );
+
+	mean = gsl_stats_mean((double*)samples, 1, num_samples);
+	//fprintf(stderr, "mu: %lf \n", mean);
+	assert( mean < mu + 0.1 );
+	assert( mean > mu - 0.1 );
+	sig    = gsl_stats_sd_m((double*)samples, 1, num_samples, mu);
+	//fprintf(stderr, "sig: %lf \n", sig);
+	assert( sig < sigma + 0.1 );
+	assert( sig > sigma - 0.1 );
+
+	err = ccs_release_object(distrib);
+	assert( err == CCS_SUCCESS );
+}
+
+static void test_normal_distribution_int_log_quantize() {
+	ccs_distribution_t distrib = NULL;
+	ccs_rng_t          rng = NULL;
+	ccs_error_t        err = CCS_SUCCESS;
+	const size_t       num_samples = 10000;
+	const ccs_float_t  mu = 3;
+	const ccs_float_t  sigma = 2;
+	const ccs_int_t    quantize = 2;
+	ccs_value_t        samples[num_samples];
+	double mean, sig;
+	double tmean, tsigma, alpha, zee, pdfa;
+
+	err = ccs_rng_create(&rng);
+	assert( err == CCS_SUCCESS );
+	err = ccs_create_normal_distribution(
+		CCS_INTEGER,
+		mu,
+		sigma,
+		CCS_LOGARITHMIC,
+		quantize,
+		&distrib);
+	assert( err == CCS_SUCCESS );
+
+	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
+	assert( err == CCS_SUCCESS );
+
+	to_float(samples, num_samples);
+	to_log(samples, num_samples);
+	mean = gsl_stats_mean((double*)samples, 1, num_samples);
+	//fprintf(stderr, "mean: %lf \n", mean);
+	// cutoff at 0.0 to have exp(v) >= 1
+	// see https://en.wikipedia.org/wiki/Truncated_normal_distribution
+	alpha = (log(0.5*quantize) - mu)/sigma;