6817476894

Committed 09 Nov 2023 09:11PM UTC coverage: 53.294% (-0.2%) from 53.454%

Build # 6817476894

Build Type

Pull #117

github

Committed by

web-flow

Commit Message

Merge 23da33f81 into 7269172d1

Pull Request Pull Request #117: Adding interface to nc_cluster_mass_ascaso to compute mean and variances

Run Details

8770 of 20841 branches covered (0.0%)

Branch coverage included in aggregate %.

62 of 333 new or added lines in 6 files covered. (18.62%)

19 existing lines in 7 files now uncovered.

44664 of 79421 relevant lines covered (56.24%)

2680942.57 hits per line

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

76.92

/numcosmo/lss/nc_cluster_mass_ascaso.c

/***************************************************************************
 *            nc_cluster_mass_ascaso.c
 *
 *  Thu Jan 26 18:25:11 2017
 *  Copyright  2017  Mariana Penna Lima and Begoña Ascaso
 *  <pennalima@gmail.com>, <bego.ascaso.work@gmail.com>
 ****************************************************************************/
/*
 * numcosmo
 * Copyright (C) Mariana Penna Lima and Begoña Ascaso 2017 <pennalima@gmail.com>
 *
 * 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_cluster_mass_ascaso
 * @title: NcClusterMassAscaso
 * @short_description: FIXME
 *
 * FIXME
 */

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

#include "lss/nc_cluster_mass_ascaso.h"
#include "math/ncm_c.h"
#include "math/ncm_cfg.h"

#ifndef NUMCOSMO_GIR_SCAN
#include <gsl/gsl_randist.h>
#endif /* NUMCOSMO_GIR_SCAN */

struct _NcClusterMassAscasoPrivate
{
  gdouble M0;
  gdouble z0;
  gdouble lnM0;
  gdouble ln1pz0;
  gdouble lnR_max;
  gdouble lnR_min;
};

G_DEFINE_TYPE_WITH_PRIVATE (NcClusterMassAscaso, nc_cluster_mass_ascaso, NC_TYPE_CLUSTER_MASS);

#define VECTOR (NCM_MODEL (ascaso)->params)
#define MU_P0 (ncm_vector_get (VECTOR, NC_CLUSTER_MASS_ASCASO_MU_P0))
#define MU_P1 (ncm_vector_get (VECTOR, NC_CLUSTER_MASS_ASCASO_MU_P1))
#define MU_P2 (ncm_vector_get (VECTOR, NC_CLUSTER_MASS_ASCASO_MU_P2))
#define SIGMA_P0  (ncm_vector_get (VECTOR, NC_CLUSTER_MASS_ASCASO_SIGMA_P0))
#define SIGMA_P1  (ncm_vector_get (VECTOR, NC_CLUSTER_MASS_ASCASO_SIGMA_P1))
#define SIGMA_P2  (ncm_vector_get (VECTOR, NC_CLUSTER_MASS_ASCASO_SIGMA_P2))

enum
{
  PROP_0,
  PROP_M0,
  PROP_Z0,
  PROP_LNRICHNESS_MIN,
  PROP_LNRICHNESS_MAX,
  PROP_SIZE,
};

static void
nc_cluster_mass_ascaso_init (NcClusterMassAscaso *ascaso)
{
  NcClusterMassAscasoPrivate * const self = ascaso->priv = nc_cluster_mass_ascaso_get_instance_private (ascaso);

  self->M0      = 0.0;
  self->z0      = 0.0;
  self->lnM0    = 0.0;
  self->ln1pz0  = 0.0;
  self->lnR_min = GSL_NEGINF;
  self->lnR_max = GSL_POSINF;
}

static void
_nc_cluster_mass_ascaso_set_property (GObject *object, guint prop_id, const GValue *value, GParamSpec *pspec)
{
  NcClusterMassAscaso *ascaso             = NC_CLUSTER_MASS_ASCASO (object);
  NcClusterMassAscasoPrivate * const self = ascaso->priv;

  g_return_if_fail (NC_IS_CLUSTER_MASS_ASCASO (object));

  switch (prop_id)
  {
    case PROP_M0:
      self->M0   = g_value_get_double (value);
      self->lnM0 = log (self->M0);
      break;
    case PROP_Z0:
      self->z0     = g_value_get_double (value);
      self->ln1pz0 = log1p (self->z0);
      break;
    case PROP_LNRICHNESS_MIN:
      self->lnR_min = g_value_get_double (value);
      g_assert (self->lnR_min < self->lnR_max);
      break;
    case PROP_LNRICHNESS_MAX:
      self->lnR_max = g_value_get_double (value);
      g_assert (self->lnR_min < self->lnR_max);
      break;
    default:
      G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
      break;
  }
}

static void
_nc_cluster_mass_ascaso_get_property (GObject *object, guint prop_id, GValue *value, GParamSpec *pspec)
{
  NcClusterMassAscaso *ascaso             = NC_CLUSTER_MASS_ASCASO (object);
  NcClusterMassAscasoPrivate * const self = ascaso->priv;

  g_return_if_fail (NC_IS_CLUSTER_MASS_ASCASO (object));

  switch (prop_id)
  {
    case PROP_M0:
      g_value_set_double (value, self->M0);
      break;
    case PROP_Z0:
      g_value_set_double (value, self->z0);
      break;
    case PROP_LNRICHNESS_MIN:
      g_value_set_double (value, self->lnR_min);
      break;
    case PROP_LNRICHNESS_MAX:
      g_value_set_double (value, self->lnR_max);
      break;
    default:
      G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
      break;
  }
}

static void
_nc_cluster_mass_ascaso_finalize (GObject *object)
{
  /* Chain up : end */
  G_OBJECT_CLASS (nc_cluster_mass_ascaso_parent_class)->finalize (object);
}

static gdouble _nc_cluster_mass_ascaso_p (NcClusterMass *clusterm,  NcHICosmo *cosmo, gdouble lnM, gdouble z, const gdouble *lnM_obs, const gdouble *lnM_obs_params);
static gdouble _nc_cluster_mass_ascaso_intp (NcClusterMass *clusterm,  NcHICosmo *cosmo, gdouble lnM, gdouble z);
static gdouble _nc_cluster_mass_ascaso_intp_bin (NcClusterMass *clusterm, NcHICosmo *cosmo, gdouble lnM, gdouble z, const gdouble *lnM_obs_lower, const gdouble *lnM_obs_upper, const gdouble *lnM_obs_params);
static gboolean _nc_cluster_mass_ascaso_resample (NcClusterMass *clusterm,  NcHICosmo *cosmo, gdouble lnM, gdouble z, gdouble *lnM_obs, const gdouble *lnM_obs_params, NcmRNG *rng);
static void _nc_cluster_mass_ascaso_p_limits (NcClusterMass *clusterm,  NcHICosmo *cosmo, const gdouble *lnM_obs, const gdouble *lnM_obs_params, gdouble *lnM_lower, gdouble *lnM_upper);
static void _nc_cluster_mass_ascaso_p_bin_limits (NcClusterMass *clusterm, NcHICosmo *cosmo, const gdouble *lnM_obs_lower, const gdouble *lnM_obs_upper, const gdouble *lnM_obs_params, gdouble *lnM_lower, gdouble *lnM_upper);
static void _nc_cluster_mass_ascaso_n_limits (NcClusterMass *clusterm,  NcHICosmo *cosmo, gdouble *lnM_lower, gdouble *lnM_upper);
static gdouble _nc_cluster_mass_ascaso_volume (NcClusterMass *clusterm);

static void _nc_cluster_mass_ascaso_p_vec_z_lnMobs (NcClusterMass *clusterm, NcHICosmo *cosmo, const gdouble lnM, const NcmVector *z, const NcmMatrix *lnM_obs, const NcmMatrix *lnM_obs_params, GArray *res);

static void
nc_cluster_mass_ascaso_class_init (NcClusterMassAscasoClass *klass)
{
  GObjectClass *object_class       = G_OBJECT_CLASS (klass);
  NcClusterMassClass *parent_class = NC_CLUSTER_MASS_CLASS (klass);
  NcmModelClass *model_class       = NCM_MODEL_CLASS (klass);

  model_class->set_property = &_nc_cluster_mass_ascaso_set_property;
  model_class->get_property = &_nc_cluster_mass_ascaso_get_property;
  object_class->finalize    = &_nc_cluster_mass_ascaso_finalize;

  ncm_model_class_set_name_nick (model_class, "Ascaso Ln-normal richness distribution", "Ascaso");
  ncm_model_class_add_params (model_class, 6, 0, PROP_SIZE);


  /**
   * NcClusterMassAscaso:M0:
   *
   * Pivot mass FIXME Set correct values (limits)
   */
  g_object_class_install_property (object_class,
                                   PROP_M0,
                                   g_param_spec_double ("M0",
                                                        NULL,
                                                        "Pivot mass",
                                                        11.0 * M_LN10, G_MAXDOUBLE, 3.0e14 / 0.71,
                                                        G_PARAM_READWRITE | G_PARAM_CONSTRUCT_ONLY | G_PARAM_STATIC_NAME | G_PARAM_STATIC_BLURB));

/**
 * NcClusterMassAscaso:Z0:
 *
 * Pivot redshift FIXME Set correct values (limits)
 */
  g_object_class_install_property (object_class,
                                   PROP_Z0,
                                   g_param_spec_double ("z0",
                                                        NULL,
                                                        "Pivot redshift",
                                                        0.0, G_MAXDOUBLE, 0.6,
                                                        G_PARAM_READWRITE | G_PARAM_CONSTRUCT_ONLY | G_PARAM_STATIC_NAME | G_PARAM_STATIC_BLURB));



  /**
   * NcClusterMassAscaso:lnRichness_min:
   *
   * FIXME Set correct values (limits)
   */
  g_object_class_install_property (object_class,
                                   PROP_LNRICHNESS_MIN,
                                   g_param_spec_double ("lnRichness-min",
                                                        NULL,
                                                        "Minimum LnRichness",
                                                        0.0, G_MAXDOUBLE, M_LN10 * 1.0,
                                                        G_PARAM_READWRITE | G_PARAM_CONSTRUCT_ONLY | G_PARAM_STATIC_NAME | G_PARAM_STATIC_BLURB));

  /**
   * NcClusterMassAscaso:lnRichness_max:
   *
   * FIXME Set correct values (limits)
   */
  g_object_class_install_property (object_class,
                                   PROP_LNRICHNESS_MAX,
                                   g_param_spec_double ("lnRichness-max",
                                                        NULL,
                                                        "Maximum LnRichness",
                                                        0.0, G_MAXDOUBLE,  M_LN10 * 2.0,
                                                        G_PARAM_READWRITE | G_PARAM_CONSTRUCT_ONLY | G_PARAM_STATIC_NAME | G_PARAM_STATIC_BLURB));

  /**
   * NcClusterMassAscaso:MU_P0:
   *
   * Distribution's  bias in the mean.
   * FIXME Set correct values (limits)
   */
  ncm_model_class_set_sparam (model_class, NC_CLUSTER_MASS_ASCASO_MU_P0, "mu_p0", "mup0",
                              0.0,  6.0, 1.0e-1,
                              NC_CLUSTER_MASS_ASCASO_DEFAULT_PARAMS_ABSTOL, NC_CLUSTER_MASS_ASCASO_DEFAULT_MU_P0,
                              NCM_PARAM_TYPE_FIXED);

  /**
   * NcClusterMassAscaso:MU_P1:
   *
   * Distribution's slope with respect to the mass in the mean.
   * FIXME Set correct values (limits)
   */
  ncm_model_class_set_sparam (model_class, NC_CLUSTER_MASS_ASCASO_MU_P1, "mu_p1", "mup1",
                              -10.0,  10.0, 1.0e-2,
                              NC_CLUSTER_MASS_ASCASO_DEFAULT_PARAMS_ABSTOL, NC_CLUSTER_MASS_ASCASO_DEFAULT_MU_P1,
                              NCM_PARAM_TYPE_FIXED);

  /**
   * NcClusterMassAscaso:MU_P2:
   *
   * Distribution's slope with respect to the redshift in the mean.
   * FIXME Set correct values (limits)
   */
  ncm_model_class_set_sparam (model_class, NC_CLUSTER_MASS_ASCASO_MU_P2, "mu_p2", "mup2",
                              -10.0,  10.0, 1.0e-2,
                              NC_CLUSTER_MASS_ASCASO_DEFAULT_PARAMS_ABSTOL, NC_CLUSTER_MASS_ASCASO_DEFAULT_MU_P2,
                              NCM_PARAM_TYPE_FIXED);

  /**
   * NcClusterMassAscaso:sigma_P0:
   *
   * Distribution's bias in the standard deviation, $\sigma \in [10^{-4}, 10]$.
   *
   */
  ncm_model_class_set_sparam (model_class, NC_CLUSTER_MASS_ASCASO_SIGMA_P0, "\\sigma_p0", "sigmap0",
                              1.0e-4, 10.0, 1.0e-2,
                              NC_CLUSTER_MASS_ASCASO_DEFAULT_PARAMS_ABSTOL, NC_CLUSTER_MASS_ASCASO_DEFAULT_SIGMA_P0,
                              NCM_PARAM_TYPE_FIXED);

  /**
   * NcClusterMassAscaso:sigma_P1:
   *
   * Distribution's slope with respect to the mass in the standard deviation.
   *
   */
  ncm_model_class_set_sparam (model_class, NC_CLUSTER_MASS_ASCASO_SIGMA_P1, "\\sigma_p1", "sigmap1",
                              -10.0, 10.0, 1.0e-2,
                              NC_CLUSTER_MASS_ASCASO_DEFAULT_PARAMS_ABSTOL, NC_CLUSTER_MASS_ASCASO_DEFAULT_SIGMA_P1,
                              NCM_PARAM_TYPE_FIXED);


/**
 * NcClusterMassAscaso:sigma_P2:
 *
 * Distribution's slope with respect to the redshift in the standard deviation.
 *
 */
  ncm_model_class_set_sparam (model_class, NC_CLUSTER_MASS_ASCASO_SIGMA_P2, "\\sigma_p2", "sigmap2",
                              -10.0,  10.0, 1.0e-2,
                              NC_CLUSTER_MASS_ASCASO_DEFAULT_PARAMS_ABSTOL, NC_CLUSTER_MASS_ASCASO_DEFAULT_SIGMA_P2,
                              NCM_PARAM_TYPE_FIXED);

  /* Check for errors in parameters initialization */
  ncm_model_class_check_params_info (model_class);

  parent_class->P              = &_nc_cluster_mass_ascaso_p;
  parent_class->intP           = &_nc_cluster_mass_ascaso_intp;
  parent_class->intP_bin       = &_nc_cluster_mass_ascaso_intp_bin;
  parent_class->resample       = &_nc_cluster_mass_ascaso_resample;
  parent_class->P_limits       = &_nc_cluster_mass_ascaso_p_limits;
  parent_class->P_bin_limits   = &_nc_cluster_mass_ascaso_p_bin_limits;
  parent_class->N_limits       = &_nc_cluster_mass_ascaso_n_limits;
  parent_class->volume         = &_nc_cluster_mass_ascaso_volume;
  parent_class->P_vec_z_lnMobs = &_nc_cluster_mass_ascaso_p_vec_z_lnMobs;
  parent_class->obs_len        = 1;
  parent_class->obs_params_len = 0;

  ncm_model_class_add_impl_flag (model_class, NC_CLUSTER_MASS_IMPL_ALL);
}

static void
_nc_cluster_mass_ascaso_lnR_sigma (NcClusterMass *clusterm, const gdouble lnM, const gdouble z, gdouble *lnR, gdouble *sigma)
{
  NcClusterMassAscaso *ascaso             = NC_CLUSTER_MASS_ASCASO (clusterm);
  NcClusterMassAscasoPrivate * const self = ascaso->priv;
  const gdouble DlnM                      = lnM - self->lnM0;
  const gdouble Dln1pz                    = log1p (z) - self->ln1pz0;

  lnR[0]   = MU_P0    + MU_P1    * DlnM + MU_P2    * Dln1pz;
  sigma[0] = SIGMA_P0 + SIGMA_P1 * DlnM + SIGMA_P2 * Dln1pz;
}

static gdouble
_nc_cluster_mass_ascaso_p (NcClusterMass *clusterm,  NcHICosmo *cosmo, gdouble lnM, gdouble z, const gdouble *lnM_obs, const gdouble *lnM_obs_params)
{
  /*NcClusterMassAscaso *ascaso    = NC_CLUSTER_MASS_ASCASO (clusterm);*/
  gdouble lnR_true, sigma;

  _nc_cluster_mass_ascaso_lnR_sigma (clusterm, lnM, z, &lnR_true, &sigma);

  {
    const gdouble x = (lnM_obs[0] - lnR_true) / sigma;

    return 1.0 / (ncm_c_sqrt_2pi () * sigma) * exp (-0.5 * x * x);
  }
}

static gdouble
_nc_cluster_mass_ascaso_intp (NcClusterMass *clusterm,  NcHICosmo *cosmo, gdouble lnM, gdouble z)
{
  NcClusterMassAscaso *ascaso             = NC_CLUSTER_MASS_ASCASO (clusterm);
  NcClusterMassAscasoPrivate * const self = ascaso->priv;
  gdouble lnR_true, sigma;

  _nc_cluster_mass_ascaso_lnR_sigma (clusterm, lnM, z, &lnR_true, &sigma);

  {
    const gdouble x_min = (lnR_true - self->lnR_min) / (M_SQRT2 * sigma);
    const gdouble x_max = (lnR_true - self->lnR_max) / (M_SQRT2 * sigma);

    if (x_max > 4.0)
      return -(erfc (x_min) - erfc (x_max)) / 2.0;
    else
      return (erf (x_min) - erf (x_max)) / 2.0;
  }
}

static gdouble
_nc_cluster_mass_ascaso_intp_bin (NcClusterMass *clusterm, NcHICosmo *cosmo, gdouble lnM, gdouble z, const gdouble *lnM_obs_lower, const gdouble *lnM_obs_upper, const gdouble *lnM_obs_params)
{
  /*NcClusterMassAscaso *ascaso = NC_CLUSTER_MASS_ASCASO (clusterm);*/
  gdouble lnR_true, sigma;

  _nc_cluster_mass_ascaso_lnR_sigma (clusterm, lnM, z, &lnR_true, &sigma);

  {
    const gdouble x_min = (lnR_true - lnM_obs_lower[0]) / (M_SQRT2 * sigma);
    const gdouble x_max = (lnR_true - lnM_obs_upper[0]) / (M_SQRT2 * sigma);

    if (x_max > 4.0)
      return -(erfc (x_min) - erfc (x_max)) / 2.0;
    else
      return (erf (x_min) - erf (x_max)) / 2.0;
  }
}

static gboolean
_nc_cluster_mass_ascaso_resample (NcClusterMass *clusterm,  NcHICosmo *cosmo, gdouble lnM, gdouble z, gdouble *lnM_obs, const gdouble *lnM_obs_params, NcmRNG *rng)
{
  NcClusterMassAscaso *ascaso             = NC_CLUSTER_MASS_ASCASO (clusterm);
  NcClusterMassAscasoPrivate * const self = ascaso->priv;
  gdouble lnR_true, sigma;

  _nc_cluster_mass_ascaso_lnR_sigma (clusterm, lnM, z, &lnR_true, &sigma);

  ncm_rng_lock (rng);
  lnM_obs[0] = lnR_true + gsl_ran_gaussian (rng->r, sigma);
  ncm_rng_unlock (rng);

  return (lnM_obs[0] <= self->lnR_max) && (lnM_obs[0] >= self->lnR_min);
}

static void
_nc_cluster_mass_ascaso_p_limits (NcClusterMass *clusterm,  NcHICosmo *cosmo, const gdouble *lnM_obs, const gdouble *lnM_obs_params, gdouble *lnM_lower, gdouble *lnM_upper)
{
  NcClusterMassAscaso *ascaso = NC_CLUSTER_MASS_ASCASO (clusterm);
  const gdouble mean          = lnM_obs[0] - MU_P0; /* - P2 * log10(1.0 + z);  FIX This!!!! What is the mean richeness? */
  const gdouble logRichnessl  = M_LN10 * log10 (1e13);
  const gdouble logRichnessu  = M_LN10 * log10 (1e15);

  *lnM_lower = logRichnessl;
  *lnM_upper = logRichnessu;

  return;
}

static void
_nc_cluster_mass_ascaso_p_bin_limits (NcClusterMass *clusterm, NcHICosmo *cosmo, const gdouble *lnM_obs_lower, const gdouble *lnM_obs_upper, const gdouble *lnM_obs_params, gdouble *lnM_lower, gdouble *lnM_upper)
{
  NcClusterMassAscaso *ascaso = NC_CLUSTER_MASS_ASCASO (clusterm);
  const gdouble lnMl          = M_LN10 * log10 (1e13);
  const gdouble lnMu          = M_LN10 * log10 (1e15);

  *lnM_lower = lnMl;
  *lnM_upper = lnMu;
}

static void
_nc_cluster_mass_ascaso_n_limits (NcClusterMass *clusterm,  NcHICosmo *cosmo, gdouble *lnM_lower, gdouble *lnM_upper)
{
  NcClusterMassAscaso *ascaso = NC_CLUSTER_MASS_ASCASO (clusterm);
  const gdouble lnMl          =  M_LN10 * log10 (1e13);
  const gdouble lnMu          =  M_LN10 * log10 (1e15);

  *lnM_lower = lnMl;
  *lnM_upper = lnMu;

  return;
}

static gdouble
_nc_cluster_mass_ascaso_volume (NcClusterMass *clusterm)
{
  NcClusterMassAscaso *ascaso             = NC_CLUSTER_MASS_ASCASO (clusterm);
  NcClusterMassAscasoPrivate * const self = ascaso->priv;

  return (self->lnR_max - self->lnR_min);
}

static void
_nc_cluster_mass_ascaso_p_vec_z_lnMobs (NcClusterMass *clusterm, NcHICosmo *cosmo, const gdouble lnM, const NcmVector *z, const NcmMatrix *lnM_obs, const NcmMatrix *lnM_obs_params, GArray *res)
{
  NcClusterMassAscaso *ascaso             = NC_CLUSTER_MASS_ASCASO (clusterm);
  NcClusterMassAscasoPrivate * const self = ascaso->priv;

  const gdouble *lnM_obs_ptr = ncm_matrix_const_data (lnM_obs);
  const gdouble *z_ptr       = ncm_vector_const_data (z);
  const guint tda            = ncm_matrix_tda (lnM_obs);
  const guint sz             = ncm_vector_stride (z);
  const guint len            = ncm_vector_len (z);
  const gdouble DlnM         = lnM - self->lnM0;
  const gdouble sqrt_2pi     = ncm_c_sqrt_2pi ();
  const gdouble lnR_pre      = MU_P0    + MU_P1    * DlnM;
  const gdouble sigma_pre    = SIGMA_P0 + SIGMA_P1 * DlnM;
  const gdouble mu_p2        = MU_P2;
  const gdouble sigma_p2     = SIGMA_P2;
  gdouble *res_ptr           = &g_array_index (res, gdouble, 0);
  gint i;

  if ((tda == 1) && (sz == 1))
  {
    for (i = 0; i < len; i++)
    {
      const gdouble Dln1pz = log1p (z_ptr[i]) - self->ln1pz0;
      const gdouble lnR    = lnR_pre + mu_p2 * Dln1pz;
      const gdouble sigma  = sigma_pre + sigma_p2 * Dln1pz;

      const gdouble x = (lnM_obs_ptr[i] - lnR) / sigma;

      res_ptr[i] = exp (-0.5 * x * x) / (sqrt_2pi * sigma);
    }
  }
  else
  {
    for (i = 0; i < len; i++)
    {
      const gdouble Dln1pz = log1p (z_ptr[i * sz]) - self->ln1pz0;
      const gdouble lnR    = lnR_pre + mu_p2 * Dln1pz;
      const gdouble sigma  = sigma_pre + sigma_p2 * Dln1pz;

      const gdouble x = (lnM_obs_ptr[i * tda] - lnR) / sigma;

      res_ptr[i] = exp (-0.5 * x * x) / (sqrt_2pi * sigma);
    }
  }
}

/**
 * nc_cluster_mass_ascaso_get_mean_richness:
 * @ascaso: a #NcClusterMassAscaso
 * @lnM: ln of the mass
 * @z: redshift
 *
 * Computes the mean of the richness distribution.
 *
 */
gdouble
nc_cluster_mass_ascaso_get_mean_richness (NcClusterMassAscaso *ascaso, gdouble lnM, gdouble z)
{
  NcClusterMassAscasoPrivate * const self = ascaso->priv;
  const gdouble DlnM                      = lnM - self->lnM0;
  const gdouble Dln1pz                    = log1p (z) - self->ln1pz0;

  return MU_P0    + MU_P1    * DlnM + MU_P2    * Dln1pz;
}

/**
 * nc_cluster_mass_ascaso_get_std_richness:
 * @ascaso: a #NcClusterMassAscaso
 * @lnM: ln of the mass
 * @z: redshift
 *
 * Computes the standard deviation of the richness distribution.
 *
 */
gdouble
nc_cluster_mass_ascaso_get_std_richness (NcClusterMassAscaso *ascaso, gdouble lnM, gdouble z)
{
  NcClusterMassAscasoPrivate * const self = ascaso->priv;
  const gdouble DlnM                      = lnM - self->lnM0;
  const gdouble Dln1pz                    = log1p (z) - self->ln1pz0;

  return SIGMA_P0 + SIGMA_P1 * DlnM + SIGMA_P2 * Dln1pz;
}


1	/***************************************************************************
2	* nc_cluster_mass_ascaso.c
3	*
4	* Thu Jan 26 18:25:11 2017
5	* Copyright 2017 Mariana Penna Lima and Begoña Ascaso
6	* <pennalima@gmail.com>, <bego.ascaso.work@gmail.com>
7	****************************************************************************/
8	/*
9	* numcosmo
10	* Copyright (C) Mariana Penna Lima and Begoña Ascaso 2017 <pennalima@gmail.com>
11	*
12	* numcosmo is free software: you can redistribute it and/or modify it
13	* under the terms of the GNU General Public License as published by the
14	* Free Software Foundation, either version 3 of the License, or
15	* (at your option) any later version.
16	*
17	* numcosmo is distributed in the hope that it will be useful, but
18	* WITHOUT ANY WARRANTY; without even the implied warranty of
19	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
20	* See the GNU General Public License for more details.
21	*
22	* You should have received a copy of the GNU General Public License along
23	* with this program. If not, see <http://www.gnu.org/licenses/>.
24	*/
25
26	/**
27	* SECTION:nc_cluster_mass_ascaso
28	* @title: NcClusterMassAscaso
29	* @short_description: FIXME
30	*
31	* FIXME
32	*/
33
34	#ifdef HAVE_CONFIG_H
35	#include "config.h"
36	#endif /* HAVE_CONFIG_H */
37	#include "build_cfg.h"
38
39	#include "lss/nc_cluster_mass_ascaso.h"
40	#include "math/ncm_c.h"
41	#include "math/ncm_cfg.h"
42
43	#ifndef NUMCOSMO_GIR_SCAN
44	#include <gsl/gsl_randist.h>
45	#endif /* NUMCOSMO_GIR_SCAN */
46
47	struct _NcClusterMassAscasoPrivate
48	{
49	gdouble M0;
50	gdouble z0;
51	gdouble lnM0;
52	gdouble ln1pz0;
53	gdouble lnR_max;
54	gdouble lnR_min;
55	};
56
57	G_DEFINE_TYPE_WITH_PRIVATE (NcClusterMassAscaso, nc_cluster_mass_ascaso, NC_TYPE_CLUSTER_MASS);	110!
58
59	#define VECTOR (NCM_MODEL (ascaso)->params)
60	#define MU_P0 (ncm_vector_get (VECTOR, NC_CLUSTER_MASS_ASCASO_MU_P0))
61	#define MU_P1 (ncm_vector_get (VECTOR, NC_CLUSTER_MASS_ASCASO_MU_P1))
62	#define MU_P2 (ncm_vector_get (VECTOR, NC_CLUSTER_MASS_ASCASO_MU_P2))
63	#define SIGMA_P0 (ncm_vector_get (VECTOR, NC_CLUSTER_MASS_ASCASO_SIGMA_P0))
64	#define SIGMA_P1 (ncm_vector_get (VECTOR, NC_CLUSTER_MASS_ASCASO_SIGMA_P1))
65	#define SIGMA_P2 (ncm_vector_get (VECTOR, NC_CLUSTER_MASS_ASCASO_SIGMA_P2))
66
67	enum
68	{
69	PROP_0,
70	PROP_M0,
71	PROP_Z0,
72	PROP_LNRICHNESS_MIN,
73	PROP_LNRICHNESS_MAX,
74	PROP_SIZE,
75	};
76
77	static void
78	nc_cluster_mass_ascaso_init (NcClusterMassAscaso *ascaso)	23✔
79	{
80	NcClusterMassAscasoPrivate * const self = ascaso->priv = nc_cluster_mass_ascaso_get_instance_private (ascaso);	23✔
81
82	self->M0 = 0.0;	23✔
83	self->z0 = 0.0;	23✔
84	self->lnM0 = 0.0;	23✔
85	self->ln1pz0 = 0.0;	23✔
86	self->lnR_min = GSL_NEGINF;	23✔
87	self->lnR_max = GSL_POSINF;	23✔
88	}	23✔
89
90	static void
91	_nc_cluster_mass_ascaso_set_property (GObject object, guint prop_id, const GValue value, GParamSpec *pspec)	92✔
92	{
93	NcClusterMassAscaso *ascaso = NC_CLUSTER_MASS_ASCASO (object);	92✔
94	NcClusterMassAscasoPrivate * const self = ascaso->priv;	92✔
95
96	g_return_if_fail (NC_IS_CLUSTER_MASS_ASCASO (object));	92✔
97
98	switch (prop_id)	92!
99	{
100	case PROP_M0:	23✔
101	self->M0 = g_value_get_double (value);	23✔
102	self->lnM0 = log (self->M0);	23✔
103	break;	23✔
104	case PROP_Z0:	23✔
105	self->z0 = g_value_get_double (value);	23✔
106	self->ln1pz0 = log1p (self->z0);	23✔
107	break;	23✔
108	case PROP_LNRICHNESS_MIN:	23✔
109	self->lnR_min = g_value_get_double (value);	23✔
110	g_assert (self->lnR_min < self->lnR_max);	23✔
111	break;	23✔
112	case PROP_LNRICHNESS_MAX:	23✔
113	self->lnR_max = g_value_get_double (value);	23✔
114	g_assert (self->lnR_min < self->lnR_max);	23✔
115	break;	23✔
116	default:	×
117	G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
118	break;	×
119	}
120	}
121
122	static void
123	_nc_cluster_mass_ascaso_get_property (GObject object, guint prop_id, GValue value, GParamSpec *pspec)	×
124	{
125	NcClusterMassAscaso *ascaso = NC_CLUSTER_MASS_ASCASO (object);	×
126	NcClusterMassAscasoPrivate * const self = ascaso->priv;	×
127
128	g_return_if_fail (NC_IS_CLUSTER_MASS_ASCASO (object));	×
129
130	switch (prop_id)	×
131	{
132	case PROP_M0:	×
133	g_value_set_double (value, self->M0);	×
134	break;	×
135	case PROP_Z0:	×
136	g_value_set_double (value, self->z0);	×
137	break;	×
138	case PROP_LNRICHNESS_MIN:	×
139	g_value_set_double (value, self->lnR_min);	×
140	break;	×
141	case PROP_LNRICHNESS_MAX:	×
142	g_value_set_double (value, self->lnR_max);	×
143	break;	×
144	default:	×
145	G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
146	break;	×
147	}
148	}
149
150	static void
151	_nc_cluster_mass_ascaso_finalize (GObject *object)	23✔
152	{
153	/* Chain up : end */
154	G_OBJECT_CLASS (nc_cluster_mass_ascaso_parent_class)->finalize (object);	23✔
155	}	23✔
156
157	static gdouble _nc_cluster_mass_ascaso_p (NcClusterMass clusterm, NcHICosmo cosmo, gdouble lnM, gdouble z, const gdouble lnM_obs, const gdouble lnM_obs_params);
158	static gdouble _nc_cluster_mass_ascaso_intp (NcClusterMass clusterm, NcHICosmo cosmo, gdouble lnM, gdouble z);
159	static gdouble _nc_cluster_mass_ascaso_intp_bin (NcClusterMass clusterm, NcHICosmo cosmo, gdouble lnM, gdouble z, const gdouble lnM_obs_lower, const gdouble lnM_obs_upper, const gdouble *lnM_obs_params);
160	static gboolean _nc_cluster_mass_ascaso_resample (NcClusterMass clusterm, NcHICosmo cosmo, gdouble lnM, gdouble z, gdouble lnM_obs, const gdouble lnM_obs_params, NcmRNG *rng);
161	static void _nc_cluster_mass_ascaso_p_limits (NcClusterMass clusterm, NcHICosmo cosmo, const gdouble lnM_obs, const gdouble lnM_obs_params, gdouble lnM_lower, gdouble lnM_upper);
162	static void _nc_cluster_mass_ascaso_p_bin_limits (NcClusterMass clusterm, NcHICosmo cosmo, const gdouble lnM_obs_lower, const gdouble lnM_obs_upper, const gdouble lnM_obs_params, gdouble lnM_lower, gdouble *lnM_upper);
163	static void _nc_cluster_mass_ascaso_n_limits (NcClusterMass clusterm, NcHICosmo cosmo, gdouble lnM_lower, gdouble lnM_upper);
164	static gdouble _nc_cluster_mass_ascaso_volume (NcClusterMass *clusterm);
165
166	static void _nc_cluster_mass_ascaso_p_vec_z_lnMobs (NcClusterMass clusterm, NcHICosmo cosmo, const gdouble lnM, const NcmVector z, const NcmMatrix lnM_obs, const NcmMatrix lnM_obs_params, GArray res);
167
168	static void
169	nc_cluster_mass_ascaso_class_init (NcClusterMassAscasoClass *klass)	110✔
170	{
171	GObjectClass *object_class = G_OBJECT_CLASS (klass);	110✔
172	NcClusterMassClass *parent_class = NC_CLUSTER_MASS_CLASS (klass);	110✔
173	NcmModelClass *model_class = NCM_MODEL_CLASS (klass);	110✔
174
175	model_class->set_property = &_nc_cluster_mass_ascaso_set_property;	110✔
176	model_class->get_property = &_nc_cluster_mass_ascaso_get_property;	110✔
177	object_class->finalize = &_nc_cluster_mass_ascaso_finalize;	110✔
178
179	ncm_model_class_set_name_nick (model_class, "Ascaso Ln-normal richness distribution", "Ascaso");	110✔
180	ncm_model_class_add_params (model_class, 6, 0, PROP_SIZE);	110✔
181
182
183	/**
184	* NcClusterMassAscaso:M0:
185	*
186	* Pivot mass FIXME Set correct values (limits)
187	*/
188	g_object_class_install_property (object_class,	110✔
189	PROP_M0,
190	g_param_spec_double ("M0",
191	NULL,
192	"Pivot mass",
193	11.0 * M_LN10, G_MAXDOUBLE, 3.0e14 / 0.71,
194	G_PARAM_READWRITE \| G_PARAM_CONSTRUCT_ONLY \| G_PARAM_STATIC_NAME \| G_PARAM_STATIC_BLURB));
195
196	/**
197	* NcClusterMassAscaso:Z0:
198	*
199	* Pivot redshift FIXME Set correct values (limits)
200	*/
201	g_object_class_install_property (object_class,	110✔
202	PROP_Z0,
203	g_param_spec_double ("z0",
204	NULL,
205	"Pivot redshift",
206	0.0, G_MAXDOUBLE, 0.6,
207	G_PARAM_READWRITE \| G_PARAM_CONSTRUCT_ONLY \| G_PARAM_STATIC_NAME \| G_PARAM_STATIC_BLURB));
208
209
210
211	/**
212	* NcClusterMassAscaso:lnRichness_min:
213	*
214	* FIXME Set correct values (limits)
215	*/
216	g_object_class_install_property (object_class,	110✔
217	PROP_LNRICHNESS_MIN,
218	g_param_spec_double ("lnRichness-min",
219	NULL,
220	"Minimum LnRichness",
221	0.0, G_MAXDOUBLE, M_LN10 * 1.0,
222	G_PARAM_READWRITE \| G_PARAM_CONSTRUCT_ONLY \| G_PARAM_STATIC_NAME \| G_PARAM_STATIC_BLURB));
223
224	/**
225	* NcClusterMassAscaso:lnRichness_max:
226	*
227	* FIXME Set correct values (limits)
228	*/
229	g_object_class_install_property (object_class,	110✔
230	PROP_LNRICHNESS_MAX,
231	g_param_spec_double ("lnRichness-max",
232	NULL,
233	"Maximum LnRichness",
234	0.0, G_MAXDOUBLE, M_LN10 * 2.0,
235	G_PARAM_READWRITE \| G_PARAM_CONSTRUCT_ONLY \| G_PARAM_STATIC_NAME \| G_PARAM_STATIC_BLURB));
236
237	/**
238	* NcClusterMassAscaso:MU_P0:
239	*
240	* Distribution's bias in the mean.
241	* FIXME Set correct values (limits)
242	*/
243	ncm_model_class_set_sparam (model_class, NC_CLUSTER_MASS_ASCASO_MU_P0, "mu_p0", "mup0",	110✔
244	0.0, 6.0, 1.0e-1,
245	NC_CLUSTER_MASS_ASCASO_DEFAULT_PARAMS_ABSTOL, NC_CLUSTER_MASS_ASCASO_DEFAULT_MU_P0,
246	NCM_PARAM_TYPE_FIXED);
247
248	/**
249	* NcClusterMassAscaso:MU_P1:
250	*
251	* Distribution's slope with respect to the mass in the mean.
252	* FIXME Set correct values (limits)
253	*/
254	ncm_model_class_set_sparam (model_class, NC_CLUSTER_MASS_ASCASO_MU_P1, "mu_p1", "mup1",	110✔
255	-10.0, 10.0, 1.0e-2,
256	NC_CLUSTER_MASS_ASCASO_DEFAULT_PARAMS_ABSTOL, NC_CLUSTER_MASS_ASCASO_DEFAULT_MU_P1,
257	NCM_PARAM_TYPE_FIXED);
258
259	/**
260	* NcClusterMassAscaso:MU_P2:
261	*
262	* Distribution's slope with respect to the redshift in the mean.
263	* FIXME Set correct values (limits)
264	*/
265	ncm_model_class_set_sparam (model_class, NC_CLUSTER_MASS_ASCASO_MU_P2, "mu_p2", "mup2",	110✔
266	-10.0, 10.0, 1.0e-2,
267	NC_CLUSTER_MASS_ASCASO_DEFAULT_PARAMS_ABSTOL, NC_CLUSTER_MASS_ASCASO_DEFAULT_MU_P2,
268	NCM_PARAM_TYPE_FIXED);
269
270	/**
271	* NcClusterMassAscaso:sigma_P0:
272	*
273	* Distribution's bias in the standard deviation, $\sigma \in [10^{-4}, 10]$.
274	*
275	*/
276	ncm_model_class_set_sparam (model_class, NC_CLUSTER_MASS_ASCASO_SIGMA_P0, "\\sigma_p0", "sigmap0",	110✔
277	1.0e-4, 10.0, 1.0e-2,
278	NC_CLUSTER_MASS_ASCASO_DEFAULT_PARAMS_ABSTOL, NC_CLUSTER_MASS_ASCASO_DEFAULT_SIGMA_P0,
279	NCM_PARAM_TYPE_FIXED);
280
281	/**
282	* NcClusterMassAscaso:sigma_P1:
283	*
284	* Distribution's slope with respect to the mass in the standard deviation.
285	*
286	*/
287	ncm_model_class_set_sparam (model_class, NC_CLUSTER_MASS_ASCASO_SIGMA_P1, "\\sigma_p1", "sigmap1",	110✔
288	-10.0, 10.0, 1.0e-2,
289	NC_CLUSTER_MASS_ASCASO_DEFAULT_PARAMS_ABSTOL, NC_CLUSTER_MASS_ASCASO_DEFAULT_SIGMA_P1,
290	NCM_PARAM_TYPE_FIXED);
291
292
293	/**
294	* NcClusterMassAscaso:sigma_P2:
295	*
296	* Distribution's slope with respect to the redshift in the standard deviation.
297	*
298	*/
299	ncm_model_class_set_sparam (model_class, NC_CLUSTER_MASS_ASCASO_SIGMA_P2, "\\sigma_p2", "sigmap2",	110✔
300	-10.0, 10.0, 1.0e-2,
301	NC_CLUSTER_MASS_ASCASO_DEFAULT_PARAMS_ABSTOL, NC_CLUSTER_MASS_ASCASO_DEFAULT_SIGMA_P2,
302	NCM_PARAM_TYPE_FIXED);
303
304	/* Check for errors in parameters initialization */
305	ncm_model_class_check_params_info (model_class);	110✔
306
307	parent_class->P = &_nc_cluster_mass_ascaso_p;	110✔
308	parent_class->intP = &_nc_cluster_mass_ascaso_intp;	110✔
309	parent_class->intP_bin = &_nc_cluster_mass_ascaso_intp_bin;	110✔
310	parent_class->resample = &_nc_cluster_mass_ascaso_resample;	110✔
311	parent_class->P_limits = &_nc_cluster_mass_ascaso_p_limits;	110✔
312	parent_class->P_bin_limits = &_nc_cluster_mass_ascaso_p_bin_limits;	110✔
313	parent_class->N_limits = &_nc_cluster_mass_ascaso_n_limits;	110✔
314	parent_class->volume = &_nc_cluster_mass_ascaso_volume;	110✔
315	parent_class->P_vec_z_lnMobs = &_nc_cluster_mass_ascaso_p_vec_z_lnMobs;	110✔
316	parent_class->obs_len = 1;	110✔
317	parent_class->obs_params_len = 0;	110✔
318
319	ncm_model_class_add_impl_flag (model_class, NC_CLUSTER_MASS_IMPL_ALL);	110✔
320	}	110✔
321
322	static void
323	_nc_cluster_mass_ascaso_lnR_sigma (NcClusterMass clusterm, const gdouble lnM, const gdouble z, gdouble lnR, gdouble *sigma)	1,781,067✔
324	{
325	NcClusterMassAscaso *ascaso = NC_CLUSTER_MASS_ASCASO (clusterm);	1,781,067✔
326	NcClusterMassAscasoPrivate * const self = ascaso->priv;	1,781,067✔
327	const gdouble DlnM = lnM - self->lnM0;	1,781,067✔
328	const gdouble Dln1pz = log1p (z) - self->ln1pz0;	1,781,067✔
329
330	lnR[0] = MU_P0 + MU_P1 * DlnM + MU_P2 * Dln1pz;	1,781,067✔
331	sigma[0] = SIGMA_P0 + SIGMA_P1 * DlnM + SIGMA_P2 * Dln1pz;	1,781,067✔
332	}	1,781,067✔
333
334	static gdouble
335	_nc_cluster_mass_ascaso_p (NcClusterMass clusterm, NcHICosmo cosmo, gdouble lnM, gdouble z, const gdouble lnM_obs, const gdouble lnM_obs_params)	13,390✔
336	{
337	/NcClusterMassAscaso ascaso = NC_CLUSTER_MASS_ASCASO (clusterm);*/
338	gdouble lnR_true, sigma;
339
340	_nc_cluster_mass_ascaso_lnR_sigma (clusterm, lnM, z, &lnR_true, &sigma);	13,390✔
341
342	{
343	const gdouble x = (lnM_obs[0] - lnR_true) / sigma;	13,390✔
344
345	return 1.0 / (ncm_c_sqrt_2pi () * sigma) * exp (-0.5 * x * x);	13,390✔
346	}
347	}
348
349	static gdouble
350	_nc_cluster_mass_ascaso_intp (NcClusterMass clusterm, NcHICosmo cosmo, gdouble lnM, gdouble z)	375,627✔
351	{
352	NcClusterMassAscaso *ascaso = NC_CLUSTER_MASS_ASCASO (clusterm);	375,627✔
353	NcClusterMassAscasoPrivate * const self = ascaso->priv;	375,627✔
354	gdouble lnR_true, sigma;
355
356	_nc_cluster_mass_ascaso_lnR_sigma (clusterm, lnM, z, &lnR_true, &sigma);	375,627✔
357
358	{
359	const gdouble x_min = (lnR_true - self->lnR_min) / (M_SQRT2 * sigma);	375,627✔
360	const gdouble x_max = (lnR_true - self->lnR_max) / (M_SQRT2 * sigma);	375,627✔
361
362	if (x_max > 4.0)	375,627!
363	return -(erfc (x_min) - erfc (x_max)) / 2.0;	×
364	else
365	return (erf (x_min) - erf (x_max)) / 2.0;	375,627✔
366	}
367	}
368
369	static gdouble
370	_nc_cluster_mass_ascaso_intp_bin (NcClusterMass clusterm, NcHICosmo cosmo, gdouble lnM, gdouble z, const gdouble lnM_obs_lower, const gdouble lnM_obs_upper, const gdouble *lnM_obs_params)	553,150✔
371	{
372	/NcClusterMassAscaso ascaso = NC_CLUSTER_MASS_ASCASO (clusterm);*/
373	gdouble lnR_true, sigma;
374
375	_nc_cluster_mass_ascaso_lnR_sigma (clusterm, lnM, z, &lnR_true, &sigma);	553,150✔
376
377	{
378	const gdouble x_min = (lnR_true - lnM_obs_lower[0]) / (M_SQRT2 * sigma);	553,150✔
379	const gdouble x_max = (lnR_true - lnM_obs_upper[0]) / (M_SQRT2 * sigma);	553,150✔
380
381	if (x_max > 4.0)	553,150✔
382	return -(erfc (x_min) - erfc (x_max)) / 2.0;	13,415✔
383	else
384	return (erf (x_min) - erf (x_max)) / 2.0;	539,735✔
385	}
386	}
387
388	static gboolean
389	_nc_cluster_mass_ascaso_resample (NcClusterMass clusterm, NcHICosmo cosmo, gdouble lnM, gdouble z, gdouble lnM_obs, const gdouble lnM_obs_params, NcmRNG *rng)	838,900✔
390	{
391	NcClusterMassAscaso *ascaso = NC_CLUSTER_MASS_ASCASO (clusterm);	838,900✔
392	NcClusterMassAscasoPrivate * const self = ascaso->priv;	838,900✔
393	gdouble lnR_true, sigma;
394
395	_nc_cluster_mass_ascaso_lnR_sigma (clusterm, lnM, z, &lnR_true, &sigma);	838,900✔
396
397	ncm_rng_lock (rng);	838,900✔
398	lnM_obs[0] = lnR_true + gsl_ran_gaussian (rng->r, sigma);	838,900✔
399	ncm_rng_unlock (rng);	838,900✔
400
401	return (lnM_obs[0] <= self->lnR_max) && (lnM_obs[0] >= self->lnR_min);	838,900✔
402	}
403
404	static void
405	_nc_cluster_mass_ascaso_p_limits (NcClusterMass clusterm, NcHICosmo cosmo, const gdouble lnM_obs, const gdouble lnM_obs_params, gdouble lnM_lower, gdouble lnM_upper)	2✔
406	{
407	NcClusterMassAscaso *ascaso = NC_CLUSTER_MASS_ASCASO (clusterm);	2✔
408	const gdouble mean = lnM_obs[0] - MU_P0; /* - P2 * log10(1.0 + z); FIX This!!!! What is the mean richeness? */	2✔
409	const gdouble logRichnessl = M_LN10 * log10 (1e13);	2✔
410	const gdouble logRichnessu = M_LN10 * log10 (1e15);	2✔
411
412	*lnM_lower = logRichnessl;	2✔
413	*lnM_upper = logRichnessu;	2✔
414
415	return;	2✔
416	}
417
418	static void
419	_nc_cluster_mass_ascaso_p_bin_limits (NcClusterMass clusterm, NcHICosmo cosmo, const gdouble lnM_obs_lower, const gdouble lnM_obs_upper, const gdouble lnM_obs_params, gdouble lnM_lower, gdouble *lnM_upper)	104✔
420	{
421	NcClusterMassAscaso *ascaso = NC_CLUSTER_MASS_ASCASO (clusterm);	104✔
422	const gdouble lnMl = M_LN10 * log10 (1e13);	104✔
423	const gdouble lnMu = M_LN10 * log10 (1e15);	104✔
424
425	*lnM_lower = lnMl;	104✔
426	*lnM_upper = lnMu;	104✔
427	}	104✔
428
429	static void
430	_nc_cluster_mass_ascaso_n_limits (NcClusterMass clusterm, NcHICosmo cosmo, gdouble lnM_lower, gdouble lnM_upper)	44✔
431	{
432	NcClusterMassAscaso *ascaso = NC_CLUSTER_MASS_ASCASO (clusterm);	44✔
433	const gdouble lnMl = M_LN10 * log10 (1e13);	44✔
434	const gdouble lnMu = M_LN10 * log10 (1e15);	44✔
435
436	*lnM_lower = lnMl;	44✔
437	*lnM_upper = lnMu;	44✔
438
439	return;	44✔
440	}
441
442	static gdouble
443	_nc_cluster_mass_ascaso_volume (NcClusterMass *clusterm)	15✔
444	{
445	NcClusterMassAscaso *ascaso = NC_CLUSTER_MASS_ASCASO (clusterm);	15✔
446	NcClusterMassAscasoPrivate * const self = ascaso->priv;	15✔
447
448	return (self->lnR_max - self->lnR_min);	15✔
449	}
450
451	static void
452	_nc_cluster_mass_ascaso_p_vec_z_lnMobs (NcClusterMass clusterm, NcHICosmo cosmo, const gdouble lnM, const NcmVector z, const NcmMatrix lnM_obs, const NcmMatrix lnM_obs_params, GArray res)	325✔
453	{
454	NcClusterMassAscaso *ascaso = NC_CLUSTER_MASS_ASCASO (clusterm);	325✔
455	NcClusterMassAscasoPrivate * const self = ascaso->priv;	325✔
456
457	const gdouble *lnM_obs_ptr = ncm_matrix_const_data (lnM_obs);	325✔
458	const gdouble *z_ptr = ncm_vector_const_data (z);	325✔
459	const guint tda = ncm_matrix_tda (lnM_obs);	325✔
460	const guint sz = ncm_vector_stride (z);	325✔
461	const guint len = ncm_vector_len (z);	325✔
462	const gdouble DlnM = lnM - self->lnM0;	325✔
463	const gdouble sqrt_2pi = ncm_c_sqrt_2pi ();	325✔
464	const gdouble lnR_pre = MU_P0 + MU_P1 * DlnM;	325✔
465	const gdouble sigma_pre = SIGMA_P0 + SIGMA_P1 * DlnM;	325✔
466	const gdouble mu_p2 = MU_P2;	325✔
467	const gdouble sigma_p2 = SIGMA_P2;	325✔
468	gdouble *res_ptr = &g_array_index (res, gdouble, 0);	325✔
469	gint i;
470
471	if ((tda == 1) && (sz == 1))	325!
472	{
473	for (i = 0; i < len; i++)	223,600✔
474	{
475	const gdouble Dln1pz = log1p (z_ptr[i]) - self->ln1pz0;	223,275✔
476	const gdouble lnR = lnR_pre + mu_p2 * Dln1pz;	223,275✔
477	const gdouble sigma = sigma_pre + sigma_p2 * Dln1pz;	223,275✔
478
479	const gdouble x = (lnM_obs_ptr[i] - lnR) / sigma;	223,275✔
480
481	res_ptr[i] = exp (-0.5 * x * x) / (sqrt_2pi * sigma);	223,275✔
482	}
483	}
484	else
485	{
486	for (i = 0; i < len; i++)	×
487	{
488	const gdouble Dln1pz = log1p (z_ptr[i * sz]) - self->ln1pz0;	×
489	const gdouble lnR = lnR_pre + mu_p2 * Dln1pz;	×
490	const gdouble sigma = sigma_pre + sigma_p2 * Dln1pz;	×
491
492	const gdouble x = (lnM_obs_ptr[i * tda] - lnR) / sigma;	×
493
494	res_ptr[i] = exp (-0.5 * x * x) / (sqrt_2pi * sigma);	×
495	}
496	}
497	}	325✔
498
499	/**
500	* nc_cluster_mass_ascaso_get_mean_richness:
501	* @ascaso: a #NcClusterMassAscaso
502	* @lnM: ln of the mass
503	* @z: redshift
504	*
505	* Computes the mean of the richness distribution.
506	*
507	*/
508	gdouble
NEW 509	nc_cluster_mass_ascaso_get_mean_richness (NcClusterMassAscaso *ascaso, gdouble lnM, gdouble z)	×
510	{
NEW 511	NcClusterMassAscasoPrivate * const self = ascaso->priv;	×
NEW 512	const gdouble DlnM = lnM - self->lnM0;	×
NEW 513	const gdouble Dln1pz = log1p (z) - self->ln1pz0;	×
514
NEW 515	return MU_P0 + MU_P1 * DlnM + MU_P2 * Dln1pz;	×
516	}
517
518	/**
519	* nc_cluster_mass_ascaso_get_std_richness:
520	* @ascaso: a #NcClusterMassAscaso
521	* @lnM: ln of the mass
522	* @z: redshift
523	*
524	* Computes the standard deviation of the richness distribution.
525	*
526	*/
527	gdouble
NEW 528	nc_cluster_mass_ascaso_get_std_richness (NcClusterMassAscaso *ascaso, gdouble lnM, gdouble z)	×
529	{
NEW 530	NcClusterMassAscasoPrivate * const self = ascaso->priv;	×
NEW 531	const gdouble DlnM = lnM - self->lnM0;	×
NEW 532	const gdouble Dln1pz = log1p (z) - self->ln1pz0;	×
533
NEW 534	return SIGMA_P0 + SIGMA_P1 * DlnM + SIGMA_P2 * Dln1pz;	×
535	}
536

NumCosmo / NumCosmo / 6817476894

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