27#include <gsl/gsl_randist.h>
28#include <gsl/gsl_vector.h>
29#include <gsl/gsl_matrix.h>
30#include <gsl/gsl_linalg.h>
31#include <gsl/gsl_blas.h>
33#include <lal/LALConstants.h>
34#include <lal/LALDatatypes.h>
35#include <lal/LALStdlib.h>
37#include <lal/LALInference.h>
38#include <lal/LALInferenceKDE.h>
72 for (
i = 0;
i < npts;
i++) {
83 for (
i = 0;
i < npts;
i++) {
84 if (mask == NULL || mask[
i]) {
85 for (
j = 0;
j < dim;
j++)
86 gsl_matrix_set(kde->
data, idx,
j, pts[
i*dim +
j]);
117 for (
i = 0;
i < npts;
i++) {
128 for (
i = 0;
i < npts;
i++) {
129 if (mask == NULL || mask[
i]) {
130 for (
j = 0;
j < dim;
j++)
131 gsl_matrix_set(kde->
data, idx,
j, gsl_matrix_get(
data,
i,
j));
160 kde->
mean = gsl_vector_calloc(dim);
163 kde->
cov = gsl_matrix_calloc(dim, dim);
170 for (
p = 0;
p < dim;
p++) {
176 kde->
data = gsl_matrix_alloc(npts, dim);
191 gsl_vector_free(kde->
mean);
194 gsl_matrix_free(kde->
cov);
196 if (kde->
npts > 0) gsl_matrix_free(kde->
data);
222 gsl_error_handler_t *default_gsl_error_handler =
223 gsl_set_error_handler_off();
225 status = gsl_linalg_cholesky_decomp(mat);
228 fprintf(stderr,
"ERROR: Unexpected problem \
229 Cholesky-decomposing matrix.\n");
235 gsl_set_error_handler(default_gsl_error_handler);
249 REAL8 det_cov, cov_factor;
254 if (kde->
npts == 0) {
264 gsl_matrix_scale(kde->
cov, cov_factor*cov_factor);
278 for (
i = 0 ;
i < kde->
dim;
i++) {
279 for (
j =
i+1 ;
j < kde->
dim;
j++)
311 gsl_vector_view
x = gsl_vector_view_array(point, dim);
314 for (
p = 0;
p < dim;
p++) {
315 val = gsl_vector_get(&
x.vector,
p);
327 for (
p = 0;
p < dim;
p++) {
337 gsl_matrix *points = gsl_matrix_alloc(n_evals, dim);
338 gsl_matrix_set_row(points, 0, &
x.vector);
341 for (
p = 0;
p < dim;
p++) {
345 val = gsl_vector_get(&
x.vector,
p);
348 gsl_matrix_set_row(points,
i, &
x.vector);
349 gsl_matrix_set(points,
i,
p, min - (val - min));
354 gsl_matrix_set_row(points,
i, &
x.vector);
355 gsl_matrix_set(points,
i,
p,
max + (
max - val));
361 gsl_matrix_set_row(points,
i, &
x.vector);
362 gsl_matrix_set(points,
i,
p, val +
width);
365 gsl_matrix_set_row(points,
i, &
x.vector);
366 gsl_matrix_set(points,
i,
p, val -
width);
376 for (
i = 0;
i < n_evals;
i++) {
377 gsl_vector_view pt = gsl_matrix_row(points,
i);
384 gsl_vector *diff = gsl_vector_alloc(dim);
385 gsl_vector *tdiff = gsl_vector_alloc(dim);
387 #pragma omp for schedule(static)
388 for (
j = 0;
j < npts;
j++) {
389 gsl_vector_view d = gsl_matrix_row(kde->
data,
j);
390 gsl_vector_memcpy(diff, &d.vector);
392 gsl_vector_sub(diff, &pt.vector);
394 gsl_vector_mul(diff, tdiff);
398 energy += gsl_vector_get(diff,
k);
399 results[
j] = -energy/2.;
402 gsl_vector_free(diff);
403 gsl_vector_free(tdiff);
413 gsl_matrix_free(points);
434 INT4 within_bounds = 0;
436 gsl_vector *unit_draw = gsl_vector_alloc(dim);
439 gsl_vector_view pt = gsl_vector_view_array(point, dim);
442 while (!within_bounds) {
446 INT4 ind = gsl_rng_uniform_int(rng, kde->
npts);
447 gsl_vector_view d = gsl_matrix_row(kde->
data, ind);
448 gsl_vector_memcpy(&pt.vector, &d.vector);
451 for (
j = 0;
j < dim;
j++) {
452 val = gsl_ran_ugaussian(rng);
453 gsl_vector_set(unit_draw,
j, val);
457 gsl_blas_dgemv(CblasNoTrans, 1.0,
462 for (
p = 0;
p < dim;
p++) {
467 if (point[
p] < min) {
468 offset = min - point[
p];
472 point[
p] = min + offset;
476 }
else if (point[
p] >
max) {
477 offset = point[
p] -
max;
481 point[
p] =
max - offset;
488 gsl_vector_free(unit_draw);
503 INT4 * signum = &sign;
504 INT4 dim = mat->size1;
506 gsl_permutation *
p = gsl_permutation_calloc(dim);
507 gsl_matrix * tmp_ptr = gsl_matrix_calloc(dim, dim);
508 gsl_matrix_memcpy(tmp_ptr, mat);
510 gsl_linalg_LU_decomp(tmp_ptr,
p, signum);
511 det = gsl_linalg_LU_det(tmp_ptr, *signum);
513 gsl_permutation_free(
p);
514 gsl_matrix_free(tmp_ptr);
534 gsl_vector_set(
mean,
i, 0.0);
537 for (
i = 0;
i < npts;
i++) {
538 x = gsl_matrix_row(
data,
i);
539 gsl_vector_add(
mean, &
x.vector);
561 gsl_vector *adiff = gsl_vector_alloc(npts);
562 gsl_vector *bdiff = gsl_vector_alloc(npts);
564 gsl_vector *
mean = gsl_vector_alloc(dim);
566 for (
i = 0;
i < dim;
i++) {
567 gsl_vector_view
a = gsl_matrix_column(
data,
i);
568 gsl_vector_memcpy(adiff, &
a.vector);
569 gsl_vector_add_constant(adiff, -gsl_vector_get(
mean,
i));
571 for (
j = 0;
j < dim;
j++) {
572 gsl_vector_view b = gsl_matrix_column(
data,
j);
573 gsl_vector_memcpy(bdiff, &b.vector);
574 gsl_vector_add_constant(bdiff, -gsl_vector_get(
mean,
j));
576 gsl_vector_mul(bdiff, adiff);
578 for (
k = 0;
k < npts;
k++)
579 var += gsl_vector_get(bdiff,
k);
580 var /= (
REAL8)(npts-1);
582 gsl_matrix_set(cov,
i,
j, var);
586 gsl_vector_free(adiff);
587 gsl_vector_free(bdiff);
588 gsl_vector_free(
mean);
606 if (max_comp <
vals[
i])
612 result += exp(
vals[
i]-max_comp);
613 result = max_comp + log(result);
INT4 LALInferenceCholeskyDecompose(gsl_matrix *mat)
Compute the Cholesky decomposition of a matrix.
void LALInferenceSetKDEBandwidth(LALInferenceKDE *kde)
Calculate the bandwidth and normalization factor for a KDE.
LALInferenceKDE * LALInferenceNewKDE(REAL8 *pts, INT4 npts, INT4 dim, INT4 *mask)
Allocate, fill, and tune a Gaussian kernel density estimate from an array of samples.
REAL8 * LALInferenceDrawKDESample(LALInferenceKDE *kde, gsl_rng *rng)
Draw a sample from a kernel density estimate.
void LALInferenceComputeMean(gsl_vector *mean, gsl_matrix *data)
Calculate the mean of a data set.
REAL8 LALInferenceKDEEvaluatePoint(LALInferenceKDE *kde, REAL8 *point)
Evaluate the (log) PDF from a KDE at a single point.
void LALInferenceDestroyKDE(LALInferenceKDE *kde)
Free an allocated KDE structure.
REAL8 log_add_exps(REAL8 *vals, INT4 size)
Determine the log of the sum of an array of exponentials.
void LALInferenceComputeCovariance(gsl_matrix *cov, gsl_matrix *data)
Calculate the covariance matrix of a data set.
LALInferenceKDE * LALInferenceInitKDE(INT4 npts, INT4 dim)
Construct an empty KDE structure.
REAL8 LALInferenceMatrixDet(gsl_matrix *mat)
Calculate the determinant of a matrix.
LALInferenceKDE * LALInferenceNewKDEfromMat(gsl_matrix *data, INT4 *mask)
Allocate, fill, and tune a Gaussian kernel density estimate from a matrix of points.
static REAL8 mean(REAL8 *array, int N)
LALInferenceParamVaryType
An enumerated type for denoting the topology of a parameter.
@ LALINFERENCE_PARAM_OUTPUT
A parameter that never changes, functions should respect this.
@ LALINFERENCE_PARAM_FIXED
A parameter that is cyclic, such as an angle between 0 and 2pi.
@ LALINFERENCE_PARAM_CIRCULAR
A parameter that simply has a maximum and a minimum.
@ LALINFERENCE_PARAM_LINEAR
void * XLALMalloc(size_t n)
void * XLALCalloc(size_t m, size_t n)
Structure containing the Guassian kernel density of a set of samples.
REAL8 log_norm_factor
Normalization factor of the KDE.
REAL8 * upper_bounds
Upper param bounds.
LALInferenceParamVaryType * lower_bound_types
Array of param boundary types.
INT4 npts
Number of points in data.
gsl_matrix * cov
The covariance matrix of data.
LALInferenceParamVaryType * upper_bound_types
Array of param boundary types.
gsl_matrix * cholesky_decomp_cov
The Cholesky decomposition of the covariance matrix, containing both terms as returned by gsl_linalg_...
gsl_matrix * data
Data to estimate the underlying distribution of.
INT4 dim
Dimension of points in data.
REAL8 * lower_bounds
Lower param bounds.
gsl_vector * mean
The mean of data.
gsl_matrix * cholesky_decomp_cov_lower
Just the lower portion of cholesky_decomp_cov.