7176301368

Committed 12 Dec 2023 02:54AM UTC coverage: 55.462% (+0.5%) from 54.987%

Build # 7176301368

Build Type

Pull #93

github

Committed by

web-flow

Commit Message

Merge e008ece66 into 04c61dd4e

Pull Request Pull Request #93: Galaxy WL reformulation

Run Details

7389 of 12715 branches covered (0.0%)

Branch coverage included in aggregate %.

1155 of 1436 new or added lines in 21 files covered. (80.43%)

2 existing lines in 2 files now uncovered.

48616 of 88264 relevant lines covered (55.08%)

1696139.76 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

77.58

/numcosmo/galaxy/nc_galaxy_wl_likelihood.c

/* -*- Mode: C; indent-tabs-mode: t; c-basic-offset: 2; tab-width: 2 -*-  */

/***************************************************************************
 *            nc_galaxy_wl_likelihood.c
 *
 *  Mon May 08 16:12:03 2023
 *  Copyright  2023  Caio Lima de Oliveira
 *  <caiolimadeoliveira@pm.me>
 ****************************************************************************/
/*
 * nc_galaxy_wl_likelihood.c
 * Copyright (C) 2023 Caio Lima de Oliveira <caiolimadeoliveira@pm.me>
 *
 * numcosmo is free software: you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * numcosmo is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 * See the GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License along
 * with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/**
 * SECTION: nc_galaxy_wl_likelihood
 * @title: NcGalaxyWLLikelihood
 * @short_description: Class describing galaxy weak lensing distributions.
 * @stability: Unstable
 *
 *
 * This class describes a galaxy weak lensing distribution.
 * It is composed by three distributions: a shape distribution $P(s)$, a proxy redshift distribution $P(z_p)$, and a position distribution $P(z)P(r)$.
 * The shape distribution is defined by the abstract class #NcGalaxySDShape.
 * The proxy redshift distribution is defined by the abstract class #NcGalaxySDZProxy.
 * The position distribution is defined by the abstract class #NcGalaxySDPosition.
 *
 */

#ifdef HAVE_CONFIG_H
#  include "config.h"
#endif /* HAVE_CONFIG_H */
#include "build_cfg.h"

#include "galaxy/nc_galaxy_wl_likelihood.h"
#include "galaxy/nc_galaxy_sd_shape.h"
#include "galaxy/nc_galaxy_sd_z_proxy.h"
#include "galaxy/nc_galaxy_sd_position.h"
#include "math/ncm_stats_dist.h"
#include "math/ncm_stats_dist_kde.h"
#include "math/ncm_stats_dist_vkde.h"
#include "math/ncm_stats_dist_kernel_gauss.h"
#include "math/ncm_stats_dist_kernel_st.h"
#include "math/ncm_integral_nd.h"
#include <math.h>
#include <gsl/gsl_math.h>

struct _NcGalaxyWLLikelihoodPrivate
{
  NcmMatrix *obs;
  NcGalaxySDShape *s_dist;
  NcGalaxySDZProxy *zp_dist;
  NcGalaxySDPosition *rz_dist;
  NcmStatsDist *kde;
  gdouble cut_fraction;
  gdouble zp_min;
  gdouble zp_max;
  gdouble r_min;
  gdouble r_max;
  gdouble scale_cut;
  gint ndata;
  gdouble prec;
  gboolean constructed;
  guint len;
};

enum
{
  PROP_0,
  PROP_OBS,
  PROP_S_DIST,
  PROP_ZP_DIST,
  PROP_RZ_DIST,
  PROP_KDE,
  PROP_R_MIN,
  PROP_R_MAX,
  PROP_ZP_MIN,
  PROP_ZP_MAX,
  PROP_SCALE_CUT,
  PROP_NDATA,
  PROP_PREC,
};

G_DEFINE_TYPE_WITH_PRIVATE (NcGalaxyWLLikelihood, nc_galaxy_wl_likelihood, G_TYPE_OBJECT);

static void
nc_galaxy_wl_likelihood_init (NcGalaxyWLLikelihood *gwl)
{
  NcGalaxyWLLikelihoodPrivate * const self = gwl->priv = nc_galaxy_wl_likelihood_get_instance_private (gwl);

  NcmStatsDistKernelGauss *kernel = ncm_stats_dist_kernel_gauss_new (4);

  self->obs         = NULL;
  self->s_dist      = NULL;
  self->zp_dist     = NULL;
  self->rz_dist     = NULL;
  self->kde         = NCM_STATS_DIST (ncm_stats_dist_vkde_new (NCM_STATS_DIST_KERNEL (kernel), NCM_STATS_DIST_CV_NONE));
  self->len         = 0;
  self->r_max       = 0.0;
  self->r_min       = 0.0;
  self->ndata       = 0.0;
  self->prec        = 1.0e-11;
  self->constructed = FALSE;
}

static void
_nc_galaxy_wl_likelihood_set_property (GObject *object, guint prop_id, const GValue *value, GParamSpec *pspec)
{
  NcGalaxyWLLikelihood *gwl                = NC_GALAXY_WL_LIKELIHOOD (object);
  NcGalaxyWLLikelihoodPrivate * const self = gwl->priv;

  g_return_if_fail (NC_IS_GALAXY_WL_LIKELIHOOD (object));

  switch (prop_id)
  {
    case PROP_OBS:
      nc_galaxy_wl_likelihood_set_obs (gwl, g_value_get_object (value));
      break;
    case PROP_S_DIST:
      self->s_dist = g_value_dup_object (value);
      break;
    case PROP_ZP_DIST:
      self->zp_dist = g_value_dup_object (value);
      break;
    case PROP_RZ_DIST:
      self->rz_dist = g_value_dup_object (value);
      break;
    case PROP_R_MIN:
      self->r_min = g_value_get_double (value);

      if (self->constructed)
        g_assert_cmpfloat (self->r_min, <, self->r_max);

      break;
    case PROP_R_MAX:
      self->r_max = g_value_get_double (value);

      if (self->constructed)
        g_assert_cmpfloat (self->r_min, <, self->r_max);

      break;
    case PROP_NDATA:
      nc_galaxy_wl_likelihood_set_ndata (gwl, g_value_get_int (value));
      break;
    case PROP_PREC:
      nc_galaxy_wl_likelihood_set_prec (gwl, g_value_get_double (value));
      break;
    default:
      G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
      break;
  }
}

static void
_nc_galaxy_wl_likelihood_get_property (GObject *object, guint prop_id, GValue *value, GParamSpec *pspec)
{
  NcGalaxyWLLikelihood *gwl                = NC_GALAXY_WL_LIKELIHOOD (object);
  NcGalaxyWLLikelihoodPrivate * const self = gwl->priv;

  g_return_if_fail (NC_IS_GALAXY_WL_LIKELIHOOD (object));

  switch (prop_id)
  {
    case PROP_OBS:
      g_value_set_object (value, nc_galaxy_wl_likelihood_peek_obs (gwl));
      break;
    case PROP_S_DIST:
      g_value_set_object (value, self->s_dist);
      break;
    case PROP_ZP_DIST:
      g_value_set_object (value, self->zp_dist);
      break;
    case PROP_RZ_DIST:
      g_value_set_object (value, self->rz_dist);
      break;
    case PROP_KDE:
      g_value_set_object (value, nc_galaxy_wl_likelihood_peek_kde (gwl));
      break;
    case PROP_R_MIN:
      g_value_set_double (value, self->r_min);
      break;
    case PROP_R_MAX:
      g_value_set_double (value, self->r_max);
      break;
    case PROP_NDATA:
      g_value_set_int (value, self->ndata);
      break;
    case PROP_PREC:
      g_value_set_double (value, self->prec);
      break;
    default:
      G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
      break;
  }
}

static void
_nc_galaxy_wl_likelihood_dispose (GObject *object)
{
  NcGalaxyWLLikelihood *gwl                = NC_GALAXY_WL_LIKELIHOOD (object);
  NcGalaxyWLLikelihoodPrivate * const self = gwl->priv;

  ncm_matrix_clear (&self->obs);
  nc_galaxy_sd_shape_clear (&self->s_dist);
  nc_galaxy_sd_z_proxy_clear (&self->zp_dist);
  nc_galaxy_sd_position_clear (&self->rz_dist);

  ncm_stats_dist_clear (&self->kde);

  G_OBJECT_CLASS (nc_galaxy_wl_likelihood_parent_class)->dispose (object);
}

static void
_nc_galaxy_wl_likelihood_finalize (GObject *object)
{
  G_OBJECT_CLASS (nc_galaxy_wl_likelihood_parent_class)->finalize (object);
}

static void
_nc_galaxy_wl_likelihood_constructed (GObject *object)
{
  /* Chain up : start */
  G_OBJECT_CLASS (nc_galaxy_wl_likelihood_parent_class)->constructed (object);
  {
    NcGalaxyWLLikelihood *gwl                = NC_GALAXY_WL_LIKELIHOOD (object);
    NcGalaxyWLLikelihoodPrivate * const self = gwl->priv;

    g_assert_cmpfloat (self->r_min, <, self->r_max);

    self->constructed = TRUE;
  }
}

static void
nc_galaxy_wl_likelihood_class_init (NcGalaxyWLLikelihoodClass *klass)
{
  GObjectClass *object_class = G_OBJECT_CLASS (klass);

  object_class->set_property = &_nc_galaxy_wl_likelihood_set_property;
  object_class->get_property = &_nc_galaxy_wl_likelihood_get_property;
  object_class->dispose      = &_nc_galaxy_wl_likelihood_dispose;
  object_class->finalize     = &_nc_galaxy_wl_likelihood_finalize;
  object_class->constructed  = &_nc_galaxy_wl_likelihood_constructed;

  /**
   * NcGalaxyWLLikelihood:obs:
   *
   * Galaxy weak lensing observables matrix.
   *
   */
  g_object_class_install_property (object_class,
                                   PROP_OBS,
                                   g_param_spec_object ("obs",
                                                        NULL,
                                                        "Galaxy weak lensing observables",
                                                        NCM_TYPE_MATRIX,
                                                        G_PARAM_READWRITE | G_PARAM_STATIC_NAME | G_PARAM_STATIC_BLURB));

  /**
   * NcGalaxyWLLikelihood:s-dist:
   *
   * A #NcGalaxySDShape object.
   *
   */

  g_object_class_install_property (object_class,
                                   PROP_S_DIST,
                                   g_param_spec_object ("s-dist",
                                                        NULL,
                                                        "Galaxy sample shape distribution",
                                                        NC_TYPE_GALAXY_SD_SHAPE,
                                                        G_PARAM_READWRITE | G_PARAM_CONSTRUCT_ONLY | G_PARAM_STATIC_NAME | G_PARAM_STATIC_BLURB));

  /**
   * NcGalaxyWLLikelihood:zp-dist:
   *
   * A #NcGalaxySDZProxy object.
   *
   */

  g_object_class_install_property (object_class,
                                   PROP_ZP_DIST,
                                   g_param_spec_object ("zp-dist",
                                                        NULL,
                                                        "Galaxy sample proxy redshift distribution",
                                                        NC_TYPE_GALAXY_SD_Z_PROXY,
                                                        G_PARAM_READWRITE | G_PARAM_CONSTRUCT_ONLY | G_PARAM_STATIC_NAME | G_PARAM_STATIC_BLURB));

  /**
   * NcGalaxyWLLikelihood:rz-dist:
   *
   * A #NcGalaxySDZPosition object.
   *
   */

  g_object_class_install_property (object_class,
                                   PROP_RZ_DIST,
                                   g_param_spec_object ("rz-dist",
                                                        NULL,
                                                        "Galaxy sample position distribution",
                                                        NC_TYPE_GALAXY_SD_POSITION,
                                                        G_PARAM_READWRITE | G_PARAM_CONSTRUCT_ONLY | G_PARAM_STATIC_NAME | G_PARAM_STATIC_BLURB));

  /**
   * NcGalaxyWLLikelihood:r-min:
   *
   * Minimum radius of the weak lensing observables.
   *
   */

  g_object_class_install_property (object_class,
                                   PROP_R_MIN,
                                   g_param_spec_double ("r-min",
                                                        NULL,
                                                        "Minimum radius of the weak lensing observables",
                                                        0.0, G_MAXDOUBLE, 0.0,
                                                        G_PARAM_READWRITE | G_PARAM_CONSTRUCT | G_PARAM_STATIC_NAME | G_PARAM_STATIC_BLURB));

  /**
   * NcGalaxyWLLikelihood:r-max:
   *
   * Maximum radius of the weak lensing observables.
   *
   */

  g_object_class_install_property (object_class,
                                   PROP_R_MAX,
                                   g_param_spec_double ("r-max",
                                                        NULL,
                                                        "Maximum radius of the weak lensing observables",
                                                        0.0, G_MAXDOUBLE, 10.0,
                                                        G_PARAM_READWRITE | G_PARAM_CONSTRUCT | G_PARAM_STATIC_NAME | G_PARAM_STATIC_BLURB));

  /**
   * NcGalaxyWLLikelihood:ndata:
   *
   * Number of data points to sample for KDE.
   *
   */

  g_object_class_install_property (object_class,
                                   PROP_NDATA,
                                   g_param_spec_int ("ndata",
                                                     NULL,
                                                     "Number of data points to sample for KDE",
                                                     0, G_MAXINT, 10000,
                                                     G_PARAM_READWRITE | G_PARAM_CONSTRUCT | G_PARAM_STATIC_NAME | G_PARAM_STATIC_BLURB));

  /**
   * NcGalaxyWLLikelihood:prec:
   *
   * Precision for integral.
   *
   */

  g_object_class_install_property (object_class,
                                   PROP_PREC,
                                   g_param_spec_double ("prec",
                                                        NULL,
                                                        "Precision for integral",
                                                        0.0, G_MAXDOUBLE, 1.e-11,
                                                        G_PARAM_READWRITE | G_PARAM_CONSTRUCT | G_PARAM_STATIC_NAME | G_PARAM_STATIC_BLURB));
}

/**
 * nc_galaxy_wl_likelihood_new:
 * @s_dist: a #NcGalaxySDShape
 * @zp_dist: a #NcGalaxySDZProxy
 * @rz_dist: a #NcGalaxySDPosition
 *
 * Creates a new galaxy weak lensing object.
 * Requires an instance of #NcGalaxySDShape, #NcGalaxySDZProxy, and #NcGalaxySDPosition.
 *
 * Returns: (transfer full): a new NcGalaxyWLLikelihood.
 */
NcGalaxyWLLikelihood *
nc_galaxy_wl_likelihood_new (NcGalaxySDShape *s_dist, NcGalaxySDZProxy *zp_dist, NcGalaxySDPosition *rz_dist)
{
  NcGalaxyWLLikelihood *gwl = g_object_new (NC_TYPE_GALAXY_WL_LIKELIHOOD,
                                            "s-dist", s_dist,
                                            "zp-dist", zp_dist,
                                            "rz-dist", rz_dist,
                                            NULL);

  return gwl;
}

/**
 * nc_galaxy_wl_likelihood_ref:
 * @gwl: a #NcGalaxyWLLikelihood
 *
 * Increase the reference of @gwl by one.
 *
 * Returns: (transfer full): @gwl.
 */
NcGalaxyWLLikelihood *
nc_galaxy_wl_likelihood_ref (NcGalaxyWLLikelihood *gwl)
{
  return g_object_ref (gwl);
}

/**
 * nc_galaxy_wl_likelihood_free:
 * @gwl: a #NcGalaxyWLLikelihood
 *
 * Decrease the reference count of @gwl by one.
 *
 */
void
nc_galaxy_wl_likelihood_free (NcGalaxyWLLikelihood *gwl)
{
  g_object_unref (gwl);
}

/**
 * nc_galaxy_wl_likelihood_clear:
 * @gwl: a #NcGalaxyWLLikelihood
 *
 * Decrease the reference count of @gwl by one, and sets the pointer *@gwl to
 * NULL.
 *
 */
void
nc_galaxy_wl_likelihood_clear (NcGalaxyWLLikelihood **gwl)
{
  g_clear_object (gwl);
}

/**
 * nc_galaxy_wl_likelihood_set_obs:
 * @gwl: a #NcGalaxyWLLikelihood
 * @obs: a #NcmMatrix
 *
 * Sets the observables matrix @obs.
 */
void
nc_galaxy_wl_likelihood_set_obs (NcGalaxyWLLikelihood *gwl, NcmMatrix *obs)
{
  NcGalaxyWLLikelihoodPrivate * const self = gwl->priv;

  g_assert_cmpuint (ncm_matrix_ncols (obs), ==, 4);
  g_assert_cmpuint (ncm_matrix_nrows (obs), >, 0);

  ncm_matrix_clear (&self->obs);

  self->len = ncm_matrix_nrows (obs);
  self->obs = ncm_matrix_ref (obs);
}

/**
 * nc_galaxy_wl_likelihood_gen_obs:
 * @gwl: a #NcGalaxyWLLikelihood
 * @nobs: number of observables to generate
 * @rng: a #NcmRNG
 *
 * Generates @nobs observables.
 */
void
nc_galaxy_wl_likelihood_gen_obs (NcGalaxyWLLikelihood *gwl, NcHICosmo *cosmo, NcHaloDensityProfile *dp, NcWLSurfaceMassDensity *smd, const gdouble z_cluster, guint nobs, NcmRNG *rng)
{
  NcGalaxyWLLikelihoodPrivate * const self = gwl->priv;
  NcmMatrix *obs                           = ncm_matrix_new (nobs, 4);
  NcmVector *sample                        = ncm_vector_new (4);
  guint i;

  for (i = 0; i < nobs; i++)
  {
    gdouble r  = 0.0;
    gdouble z  = 0.0;
    gdouble zp = 0.0;

    do {
      r = nc_galaxy_sd_position_gen_r (self->rz_dist, rng);
    } while ((self->r_min > r) || (r > self->r_max));

    while (TRUE)
    {
      z = nc_galaxy_sd_position_gen_z (self->rz_dist, rng);

      if (nc_galaxy_sd_z_proxy_gen (self->zp_dist, rng, z, &zp))
        break;
    }

    ncm_vector_set (sample, 0, r);
    ncm_vector_set (sample, 1, zp);
    nc_galaxy_sd_shape_gen (self->s_dist, cosmo, dp, smd, z_cluster, rng, r, z, ncm_vector_ptr (sample, 2), ncm_vector_ptr (sample, 3));

    ncm_matrix_set_row (obs, i, sample);
  }

  nc_galaxy_wl_likelihood_set_obs (gwl, obs);

  ncm_vector_free (sample);
  ncm_matrix_free (obs);
}

/**
 * nc_galaxy_wl_likelihood_peek_obs:
 * @gwl: a #NcGalaxyWLLikelihood
 *
 * Gets the observables matrix.
 *
 * Returns: (transfer none): the observables matrix.
 */
NcmMatrix *
nc_galaxy_wl_likelihood_peek_obs (NcGalaxyWLLikelihood *gwl)
{
  NcGalaxyWLLikelihoodPrivate * const self = gwl->priv;

  return self->obs;
}

/**
 * nc_galaxy_wl_likelihood_peek_kde:
 * @gwl: a #NcGalaxyWLLikelihood
 *
 * Gets the observables matrix.
 *
 * Returns: (transfer none): the observables matrix.
 */
NcmStatsDist *
nc_galaxy_wl_likelihood_peek_kde (NcGalaxyWLLikelihood *gwl)
{
  NcGalaxyWLLikelihoodPrivate * const self = gwl->priv;

  return self->kde;
}

struct int_data
{
  NcGalaxySDPosition *rz_dist;
  NcGalaxySDZProxy *zp_dist;
  NcGalaxySDShape *s_dist;
  NcHICosmo *cosmo;
  NcHaloDensityProfile *dp;
  NcWLSurfaceMassDensity *smd;
  gdouble z_cluster;
  gdouble zp;
  gdouble et;
  gdouble ex;
  gdouble r;
};

static void nc_galaxy_wl_likelihood_integ (NcmIntegralND *intnd, NcmVector *x, guint dim, guint npoints, guint fdim, NcmVector *fval);
static void nc_galaxy_wl_likelihood_int_dim (NcmIntegralND *intnd, guint *dim, guint *fdim);

NCM_INTEGRAL_ND_DEFINE_TYPE (NC, GALAXY_WL_LIKELIHOOD_INT, NcGalaxyWLLikelihoodInt, nc_galaxy_wl_likelihood_int, nc_galaxy_wl_likelihood_int_dim, nc_galaxy_wl_likelihood_integ, struct int_data);

static void
nc_galaxy_wl_likelihood_integ (NcmIntegralND *intnd, NcmVector *x, guint dim, guint npoints, guint fdim, NcmVector *fval)
{
  NcGalaxyWLLikelihoodInt *lh_int = NC_GALAXY_WL_LIKELIHOOD_INT (intnd);
  guint i;

  for (i = 0; i < npoints; i++)
  {
    const gdouble z = ncm_vector_get (x, i);
    gdouble res     = nc_galaxy_sd_position_integ (lh_int->data.rz_dist, lh_int->data.r, z) *
                      nc_galaxy_sd_z_proxy_integ (lh_int->data.zp_dist, z, lh_int->data.zp) *
                      nc_galaxy_sd_shape_integ (lh_int->data.s_dist,
                                                lh_int->data.cosmo,
                                                lh_int->data.dp,
                                                lh_int->data.smd,
                                                lh_int->data.z_cluster,
                                                lh_int->data.r,
                                                z,
                                                lh_int->data.et,
                                                lh_int->data.ex);

    ncm_vector_set (fval, i, res);
  }
}

static void
nc_galaxy_wl_likelihood_int_dim (NcmIntegralND *intnd, guint *dim, guint *fdim)
{
  *dim  = 1;
  *fdim = 1;
}

void
nc_galaxy_wl_likelihood_prepare (NcGalaxyWLLikelihood *gwl, NcHICosmo *cosmo, NcHaloDensityProfile *dp, NcWLSurfaceMassDensity *smd, const gdouble z_cluster)
{
  NcGalaxyWLLikelihoodPrivate * const self = gwl->priv;
  NcmVector *sample                        = ncm_vector_new (4);
  NcmRNG *rng                              = ncm_rng_new (NULL);
  gdouble in_cut                           = 0;
  gdouble out_cut                          = 0;
  gint i                                   = 0;

  ncm_stats_dist_reset (self->kde);

  while (i < self->ndata)
  {
    const gdouble r = nc_galaxy_sd_position_gen_r (self->rz_dist, rng);
    gdouble z       = 0.0;
    gdouble zp      = 0.0;

    while (TRUE)
    {
      z = nc_galaxy_sd_position_gen_z (self->rz_dist, rng);

      if (nc_galaxy_sd_z_proxy_gen (self->zp_dist, rng, z, &zp))
        break;
    }

    if ((self->r_min <= r) && (r <= self->r_max))
      in_cut++;
    else
      out_cut++;

    i++;

    ncm_vector_set (sample, 0, r);
    ncm_vector_set (sample, 1, zp);
    nc_galaxy_sd_shape_gen (self->s_dist, cosmo, dp, smd, z_cluster, rng, r, z, ncm_vector_ptr (sample, 2), ncm_vector_ptr (sample, 3));

    ncm_stats_dist_add_obs (self->kde, sample);
  }

  printf ("# in_cut: %g\n", in_cut);

  self->cut_fraction = (gdouble) in_cut / (gdouble) (in_cut + out_cut);

  ncm_stats_dist_prepare (self->kde);

  ncm_rng_clear (&rng);
  ncm_vector_clear (&sample);
}

/**
 * nc_galaxy_wl_likelihood_eval_m2lnP:
 * @gwl: a #NcGalaxyWLLikelihood
 * @cosmo: a #NcHICosmo
 * @dp: a #NcHaloDensityProfile
 * @smd: a #NcWLSurfaceMassDensity
 * @z_cluster: cluster redshift $z_\mathrm{cl}$
 * @gal_obs: a #NcmMatrix
 * @m2lnP_gal: (out) (optional): a #NcmVector
 *
 * Computes the observables probability given the theoretical modeling using
 * integration method.
 *
 *
 * Returns: $-2\ln(P)$.
 */
gdouble
nc_galaxy_wl_likelihood_eval_m2lnP (NcGalaxyWLLikelihood *gwl, NcHICosmo *cosmo, NcHaloDensityProfile *dp, NcWLSurfaceMassDensity *smd, const gdouble z_cluster, NcmMatrix *gal_obs, NcmVector *m2lnP_gal)
{
  NcGalaxyWLLikelihoodPrivate * const self = gwl->priv;
  gdouble verr, vzpi, vzpf, vres, vtmp;
  NcmVector *err  = ncm_vector_new_data_static (&verr, 1, 1);
  NcmVector *zpi  = ncm_vector_new_data_static (&vzpi, 1, 1);
  NcmVector *zpf  = ncm_vector_new_data_static (&vzpf, 1, 1);
  NcmVector *res  = ncm_vector_new_data_static (&vres, 1, 1);
  NcmVector *tmp  = ncm_vector_new_data_static (&vtmp, 1, 1);
  const guint len = ncm_matrix_nrows (gal_obs);
  gdouble zp_i    = 1.0e-11;
  gdouble zp_f    = 100.0;
  gdouble result  = 0;
  guint gal_i;

  ncm_vector_set (zpi, 0, zp_i);
  ncm_vector_set (zpf, 0, zp_f);

  g_assert_cmpuint (ncm_matrix_ncols (gal_obs), ==, 4);

  if (m2lnP_gal != NULL)
    g_assert_cmpuint (ncm_vector_len (m2lnP_gal), ==, len);

  for (gal_i = 0; gal_i < len; gal_i++)
  {
    NcGalaxyWLLikelihoodInt *likelihood_integral = g_object_new (nc_galaxy_wl_likelihood_int_get_type (), NULL);
    NcmIntegralND *lh_int                        = NCM_INTEGRAL_ND (likelihood_integral);
    gdouble m2lnP_gal_i;
    gdouble z_min, z_max;

    likelihood_integral->data.rz_dist   = NC_GALAXY_SD_POSITION (self->rz_dist);
    likelihood_integral->data.zp_dist   = NC_GALAXY_SD_Z_PROXY (self->zp_dist);
    likelihood_integral->data.s_dist    = NC_GALAXY_SD_SHAPE (self->s_dist);
    likelihood_integral->data.cosmo     = cosmo;
    likelihood_integral->data.dp        = dp;
    likelihood_integral->data.smd       = smd;
    likelihood_integral->data.z_cluster = z_cluster;
    likelihood_integral->data.r         = ncm_matrix_get (gal_obs, gal_i, 0);
    likelihood_integral->data.zp        = ncm_matrix_get (gal_obs, gal_i, 1);
    likelihood_integral->data.et        = ncm_matrix_get (gal_obs, gal_i, 2);
    likelihood_integral->data.ex        = ncm_matrix_get (gal_obs, gal_i, 3);

    ncm_integral_nd_set_method (lh_int, NCM_INTEGRAL_ND_METHOD_CUBATURE_P);

    nc_galaxy_sd_z_proxy_get_true_z_lim (self->zp_dist, likelihood_integral->data.zp, &z_min, &z_max);
    ncm_integral_nd_set_reltol (lh_int, self->prec);
    ncm_integral_nd_set_abstol (lh_int, 0.0);

    if (z_min < z_cluster)
    {
      vzpi = z_min;
      vzpf = z_cluster;
      ncm_integral_nd_eval (lh_int, zpi, zpf, tmp, err);

      vzpi = z_cluster;
      vzpf = z_max;
      ncm_integral_nd_eval (lh_int, zpi, zpf, res, err);
      vres += vtmp;
    }
    else
    {
      vzpi = z_min;
      vzpf = z_max;
      ncm_integral_nd_eval (lh_int, zpi, zpf, res, err);
    }

    m2lnP_gal_i = -2.0 * log (vres);

    if (m2lnP_gal != NULL)
      ncm_vector_set (m2lnP_gal, gal_i, m2lnP_gal_i);

    result += m2lnP_gal_i;

    ncm_integral_nd_clear (&lh_int);
  }

  ncm_vector_free (err);
  ncm_vector_free (zpi);
  ncm_vector_free (zpf);
  ncm_vector_free (res);
  ncm_vector_free (tmp);

  return result;
}

/**
 * nc_galaxy_wl_likelihood_kde_eval_m2lnP:
 * @gwl: a #NcGalaxyWLLikelihood
 * @cosmo: a #NcHICosmo
 * @dp: a #NcHaloDensityProfile
 * @smd: a #NcWLSurfaceMassDensity
 * @z_cluster: cluster redshift $z_\mathrm{cl}$
 * @gal_obs: a #NcmMatrix
 * @m2lnP_gal: (out) (optional): a #NcmVector
 *
 * Computes the observables probability given the theoretical modeling using
 * kernel density estimation method.
 *
 *
 * Returns: $-2\ln(P)$.
 */
gdouble
nc_galaxy_wl_likelihood_kde_eval_m2lnP (NcGalaxyWLLikelihood *gwl, NcHICosmo *cosmo, NcHaloDensityProfile *dp, NcWLSurfaceMassDensity *smd, const gdouble z_cluster, NcmMatrix *gal_obs, NcmVector *m2lnP_gal)
{
  NcGalaxyWLLikelihoodPrivate * const self = gwl->priv;
  NcmVector *data_vec                      = ncm_vector_new (4);
  const guint len                          = ncm_matrix_nrows (gal_obs);
  gdouble res                              = 0.0;
  glong in_cut                             = 0;

  guint gal_i;

  g_assert_cmpuint (ncm_matrix_ncols (gal_obs), ==, 4);

  if (m2lnP_gal != NULL)
    g_assert_cmpuint (ncm_vector_len (m2lnP_gal), ==, len);

  nc_galaxy_wl_likelihood_prepare (gwl, cosmo, dp, smd, z_cluster);

  if (m2lnP_gal != NULL)
  {
    for (gal_i = 0; gal_i < len; gal_i++)
    {
      const gdouble r_i  = ncm_matrix_get (gal_obs, gal_i, 0);
      const gdouble z_i  = ncm_matrix_get (gal_obs, gal_i, 1);
      const gdouble et_i = ncm_matrix_get (gal_obs, gal_i, 2);
      const gdouble ex_i = ncm_matrix_get (gal_obs, gal_i, 3);

      ncm_vector_set (data_vec, 0, r_i);
      ncm_vector_set (data_vec, 1, z_i);
      ncm_vector_set (data_vec, 2, et_i);
      ncm_vector_set (data_vec, 3, ex_i);

      if ((self->r_min <= r_i) && (r_i <= self->r_max))
      {
        const gdouble m2lnP_gal_i = ncm_stats_dist_eval_m2lnp (self->kde, data_vec);

        ncm_vector_set (m2lnP_gal, gal_i, m2lnP_gal_i);

        in_cut++;
        res += m2lnP_gal_i;
      }
    }
  }
  else
  {
    for (gal_i = 0; gal_i < self->len; gal_i++)
    {
      const gdouble r_i = ncm_matrix_get (gal_obs, gal_i, 0);
      const gdouble z_i = ncm_matrix_get (gal_obs, gal_i, 1);
      const gdouble s_i = ncm_matrix_get (gal_obs, gal_i, 2);

      ncm_vector_set (data_vec, 0, r_i);
      ncm_vector_set (data_vec, 1, z_i);
      ncm_vector_set (data_vec, 2, s_i);

      if ((self->r_min <= r_i) && (r_i <= self->r_max))
      {
        const gdouble m2lnP_gal_i = ncm_stats_dist_eval_m2lnp (self->kde, data_vec);

        in_cut++;
        res += m2lnP_gal_i;
      }
    }
  }

  res += 2.0 * in_cut * log (self->cut_fraction);

  ncm_vector_clear (&data_vec);

  return res;
}

/**
 * nc_galaxy_wl_likelihood_len:
 * @gwll: a #NcGalaxyWL
 *
 * Returns: the number of galaxies in @gwl.
 */
guint
nc_galaxy_wl_likelihood_len (NcGalaxyWLLikelihood *gwll)
{
  NcGalaxyWLLikelihoodPrivate * const self = gwll->priv;

  return self->len;
}

/**
 * nc_galaxy_wl_likelihood_set_cut:
 * @gwl: a #NcGalaxyWL
 * @r_min: minimum projected radius $r_\mathrm{min}$
 * @r_max: maximum projected radius $r_\mathrm{max}$
 *
 * Sets the cut in the observables.
 *
 */
void
nc_galaxy_wl_likelihood_set_cut (NcGalaxyWLLikelihood *gwl, const gdouble r_min, const gdouble r_max)
{
  NcGalaxyWLLikelihoodPrivate * const self = gwl->priv;

  g_assert_cmpfloat (r_min, <, r_max);

  self->r_min = r_min;
  self->r_max = r_max;
}

/**
 * nc_galaxy_wl_likelihood_set_ndata:
 * @gwl: a #NcGalaxyWL
 * @ndata: number of samples to take for KDE
 *
 * Sets the number of samples ndata.
 *
 */
void
nc_galaxy_wl_likelihood_set_ndata (NcGalaxyWLLikelihood *gwl, gdouble ndata)
{
  NcGalaxyWLLikelihoodPrivate * const self = gwl->priv;

  self->ndata = ndata;
}

/**
 * nc_galaxy_wl_likelihood_set_prec:
 * @gwl: a #NcGalaxyWL
 * @prec: precision for integral
 *
 * Sets the number of samples ndata.
 *
 */
void
nc_galaxy_wl_likelihood_set_prec (NcGalaxyWLLikelihood *gwl, gdouble prec)
{
  NcGalaxyWLLikelihoodPrivate * const self = gwl->priv;

  self->prec = prec;
}

/**
 * nc_galaxy_wl_likelihood_gen:
 * @gwl: a #NcGalaxyWL
 * @cosmo: a #NcHICosmo
 * @dp: a #NcHaloDensityProfile
 * @smd: a #NcWLSurfaceMassDensity
 * @z_cluster: cluster redshift $z_\mathrm{cl}$
 * @rng: a #NcmRNG
 *
 * Returns: (transfer full): a #NcmVector with a randomly generated observation.
 *
 */
NcmVector *
nc_galaxy_wl_likelihood_gen (NcGalaxyWLLikelihood *gwl, NcHICosmo *cosmo, NcHaloDensityProfile *dp, NcWLSurfaceMassDensity *smd, const gdouble z_cluster, NcmRNG *rng)
{
  NcGalaxyWLLikelihoodPrivate * const self = gwl->priv;
  NcmVector *sample                        = ncm_vector_new (3);
  gdouble r                                = 0.0;
  gdouble z                                = 0.0;
  gdouble zp                               = 0.0;
  gdouble s                                = 0.0;

  r = nc_galaxy_sd_position_gen_r (self->rz_dist, rng);

  while (TRUE)
  {
    z = nc_galaxy_sd_position_gen_z (self->rz_dist, rng);

    if (nc_galaxy_sd_z_proxy_gen (self->zp_dist, rng, z, &zp))
      break;
  }

  ncm_vector_set (sample, 0, r);
  ncm_vector_set (sample, 1, zp);
  nc_galaxy_sd_shape_gen (self->s_dist, cosmo, dp, smd, z_cluster, rng, r, z, ncm_vector_ptr (sample, 2), ncm_vector_ptr (sample, 3));

  return sample;
}


NumCosmo / NumCosmo / 7176301368

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous