7062014968

Committed 01 Dec 2023 03:33PM UTC coverage: 31.828% (-21.6%) from 53.447%

Build # 7062014968

Build Type

Pull #120

github

Committed by

web-flow

Commit Message

Merge 854e08790 into 24e5f6a9b

Pull Request Pull Request #120: Moving to meson

Run Details

14310 of 71092 branches covered (0.0%)

Branch coverage included in aggregate %.

1956 of 3488 new or added lines in 312 files covered. (56.08%)

3937 existing lines in 210 files now uncovered.

57509 of 154555 relevant lines covered (37.21%)

2445460.13 hits per line

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

12.04

/numcosmo/math/ncm_fit_esmcmc_walker_walk.c

/***************************************************************************
 *            ncm_fit_esmcmc_walker_walk.c
 *
 *  Tue March 29 10:41:49 2016
 *  Copyright  2016  Sandro Dias Pinto Vitenti
 *  <vitenti@uel.br>
 ****************************************************************************/
/*
 * ncm_fit_esmcmc_walker_walk.c
 * Copyright (C) 2016 Sandro Dias Pinto Vitenti <vitenti@uel.br>
 *
 * 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:ncm_fit_esmcmc_walker_walk
 * @title: NcmFitESMCMCWalkerWalk
 * @short_description: Ensemble sampler Markov Chain Monte Carlo walker - walk move.
 * @stability: Unstable
 *
 * Implementing walk move walker for #NcmFitESMCMC (affine invariant).
 *
 */

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

#include "math/ncm_fit_esmcmc_walker.h"
#include "math/ncm_fit_esmcmc_walker_walk.h"

#include "math/ncm_fit_esmcmc.h"

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

enum
{
  PROP_0,
  PROP_SCALE
};

struct _NcmFitESMCMCWalkerWalk
{
  /*< private >*/
  NcmFitESMCMCWalker parent_instance;
  guint size;
  guint size_2;
  guint nparams;
  gdouble a;
  gdouble sqrt_nparams;
  NcmMatrix *z;
  GPtrArray *thetabar;
  GArray *indices;
  GArray *numbers;
};

G_DEFINE_TYPE (NcmFitESMCMCWalkerWalk, ncm_fit_esmcmc_walker_walk, NCM_TYPE_FIT_ESMCMC_WALKER)

static void
ncm_fit_esmcmc_walker_walk_init (NcmFitESMCMCWalkerWalk *walk)
{
  walk->size         = 0;
  walk->size_2       = 0;
  walk->nparams      = 0;
  walk->a            = 0.0;
  walk->sqrt_nparams = 0.0;
  walk->z            = NULL;
  walk->thetabar     = g_ptr_array_new ();
  walk->indices      = g_array_new (TRUE, TRUE, sizeof (guint));
  walk->numbers      = g_array_new (TRUE, TRUE, sizeof (guint));

  g_ptr_array_set_free_func (walk->thetabar, (GDestroyNotify) ncm_vector_free);
}

static void
ncm_fit_esmcmc_walker_walk_dispose (GObject *object)
{
  NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (object);

  ncm_matrix_clear (&walk->z);

  g_clear_pointer (&walk->thetabar, g_ptr_array_unref);
  g_clear_pointer (&walk->indices, g_array_unref);
  g_clear_pointer (&walk->numbers, g_array_unref);

  /* Chain up : end */
  G_OBJECT_CLASS (ncm_fit_esmcmc_walker_walk_parent_class)->dispose (object);
}

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

static void
ncm_fit_esmcmc_walker_walk_set_property (GObject *object, guint prop_id, const GValue *value, GParamSpec *pspec)
{
  NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (object);

  g_return_if_fail (NCM_IS_FIT_ESMCMC_WALKER_WALK (object));

  switch (prop_id)
  {
    case PROP_SCALE:
      ncm_fit_esmcmc_walker_walk_set_scale (walk, g_value_get_double (value));
      break;
    default:
      G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
      break;
  }
}

static void
ncm_fit_esmcmc_walker_walk_get_property (GObject *object, guint prop_id, GValue *value, GParamSpec *pspec)
{
  NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (object);

  g_return_if_fail (NCM_IS_FIT_ESMCMC_WALKER_WALK (object));

  switch (prop_id)
  {
    case PROP_SCALE:
      g_value_set_double (value, ncm_fit_esmcmc_walker_walk_get_scale (walk));
      break;
    default:
      G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
      break;
  }
}

static void _ncm_fit_esmcmc_walker_walk_set_size (NcmFitESMCMCWalker *walker, guint size);
static guint _ncm_fit_esmcmc_walker_walk_get_size (NcmFitESMCMCWalker *walker);
static void _ncm_fit_esmcmc_walker_walk_set_nparams (NcmFitESMCMCWalker *walker, guint nparams);
static guint _ncm_fit_esmcmc_walker_walk_get_nparams (NcmFitESMCMCWalker *walker);
static void _ncm_fit_esmcmc_walker_walk_setup (NcmFitESMCMCWalker *walker, NcmMSet *mset, GPtrArray *theta, GPtrArray *m2lnL, guint ki, guint kf, NcmRNG *rng);
static void _ncm_fit_esmcmc_walker_walk_step (NcmFitESMCMCWalker *walker, GPtrArray *theta, GPtrArray *m2lnL, NcmVector *thetastar, guint k);
static gdouble _ncm_fit_esmcmc_walker_walk_prob (NcmFitESMCMCWalker *walker, GPtrArray *theta, GPtrArray *m2lnL, NcmVector *thetastar, guint k, const gdouble m2lnL_cur, const gdouble m2lnL_star);
static gdouble _ncm_fit_esmcmc_walker_walk_prob_norm (NcmFitESMCMCWalker *walker, GPtrArray *theta, GPtrArray *m2lnL, NcmVector *thetastar, guint k);
static void _ncm_fit_esmcmc_walker_walk_clean (NcmFitESMCMCWalker *walker, guint ki, guint kf);
static const gchar *_ncm_fit_esmcmc_walker_walk_desc (NcmFitESMCMCWalker *walker);

static void
ncm_fit_esmcmc_walker_walk_class_init (NcmFitESMCMCWalkerWalkClass *klass)
{
  GObjectClass *object_class            = G_OBJECT_CLASS (klass);
  NcmFitESMCMCWalkerClass *walker_class = NCM_FIT_ESMCMC_WALKER_CLASS (klass);

  object_class->set_property = ncm_fit_esmcmc_walker_walk_set_property;
  object_class->get_property = ncm_fit_esmcmc_walker_walk_get_property;
  object_class->dispose      = ncm_fit_esmcmc_walker_walk_dispose;
  object_class->finalize     = ncm_fit_esmcmc_walker_walk_finalize;

  g_object_class_install_property (object_class,
                                   PROP_SCALE,
                                   g_param_spec_double ("scale",
                                                        NULL,
                                                        "Walk scale a",
                                                        1.0e-2, G_MAXDOUBLE, 0.2,
                                                        G_PARAM_READWRITE | G_PARAM_CONSTRUCT | G_PARAM_STATIC_NAME | G_PARAM_STATIC_BLURB));

  walker_class->set_size    = &_ncm_fit_esmcmc_walker_walk_set_size;
  walker_class->get_size    = &_ncm_fit_esmcmc_walker_walk_get_size;
  walker_class->set_nparams = &_ncm_fit_esmcmc_walker_walk_set_nparams;
  walker_class->get_nparams = &_ncm_fit_esmcmc_walker_walk_get_nparams;
  walker_class->setup       = &_ncm_fit_esmcmc_walker_walk_setup;
  walker_class->step        = &_ncm_fit_esmcmc_walker_walk_step;
  walker_class->prob        = &_ncm_fit_esmcmc_walker_walk_prob;
  walker_class->prob_norm   = &_ncm_fit_esmcmc_walker_walk_prob_norm;
  walker_class->clean       = &_ncm_fit_esmcmc_walker_walk_clean;
  walker_class->desc        = &_ncm_fit_esmcmc_walker_walk_desc;
}

static void
_ncm_fit_esmcmc_walker_walk_set_sys (NcmFitESMCMCWalker *walker, guint size, guint nparams)
{
  NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (walker);

  g_assert_cmpuint (size, >, 0);
  g_assert_cmpuint (nparams, >, 0);

  if ((walk->size != size) || (walk->nparams != nparams))
  {
    guint i;

    ncm_matrix_clear (&walk->z);

    g_assert (size % 2 == 0);

    walk->z = ncm_matrix_new (size, nparams);

    g_array_set_size (walk->indices, size * nparams);

    g_array_set_size (walk->numbers, size);

    for (i = 0; i < size; i++)
      g_array_index (walk->numbers, guint, i) = i;

    g_ptr_array_set_size (walk->thetabar, 0);

    for (i = 0; i < size; i++)
    {
      NcmVector *thetabar_i = ncm_vector_new (nparams);

      g_ptr_array_add (walk->thetabar, thetabar_i);
    }

    walk->size   = size;
    walk->size_2 = size / 2;

    walk->size         = size;
    walk->nparams      = nparams;
    walk->sqrt_nparams = sqrt (1.0 * walk->nparams);
  }
}

static void
_ncm_fit_esmcmc_walker_walk_set_size (NcmFitESMCMCWalker *walker, guint size)
{
  NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (walker);

  g_assert_cmpuint (size, >, 0);

  if (walk->nparams != 0)
    _ncm_fit_esmcmc_walker_walk_set_sys (walker, size, walk->nparams);
  else
    walk->size = size;
}

static guint
_ncm_fit_esmcmc_walker_walk_get_size (NcmFitESMCMCWalker *walker)
{
  NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (walker);

  return walk->size;
}

static void
_ncm_fit_esmcmc_walker_walk_set_nparams (NcmFitESMCMCWalker *walker, guint nparams)
{
  NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (walker);

  g_assert_cmpuint (nparams, >, 0);

  if (walk->size != 0)
    _ncm_fit_esmcmc_walker_walk_set_sys (walker, walk->size, nparams);
  else
    walk->nparams = nparams;
}

static guint
_ncm_fit_esmcmc_walker_walk_get_nparams (NcmFitESMCMCWalker *walker)
{
  NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (walker);

  return walk->nparams;
}

static void
_ncm_fit_esmcmc_walker_walk_setup (NcmFitESMCMCWalker *walker, NcmMSet *mset, GPtrArray *theta, GPtrArray *m2lnL, guint ki, guint kf, NcmRNG *rng)
{
  NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (walker);
  guint k;

  for (k = ki; k < kf; k++)
  {
    const guint subensemble = (k < walk->size_2) ? walk->size_2 : 0;
    guint pi;

    for (pi = 0; pi < walk->nparams; pi++)
    {
      const gdouble zi = gsl_ran_ugaussian (rng->r);

      ncm_matrix_set (walk->z, k, pi, zi);
    }

    gsl_ran_choose (rng->r,
                    &g_array_index (walk->indices, guint, k * walk->nparams),
                    walk->nparams,
                    &g_array_index (walk->numbers, guint, subensemble),
                    walk->size_2,
                    g_array_get_element_size (walk->numbers));

/*
 *   for (pi = 0; pi < walk->nparams; pi++)
 *   {
 *     printf ("Choosing for walker %u, walker %u | nparams %u size_2 %u subensemble %u | numbers[subensemble + %u] = %u\n",
 *             k, g_array_index (walk->indices, guint, k * walk->nparams + pi),
 *             walk->nparams, walk->size_2, subensemble,
 *             pi, g_array_index (walk->numbers, guint, subensemble + pi));
 *   }
 */
  }
}

static void
_ncm_fit_esmcmc_walker_walk_step (NcmFitESMCMCWalker *walker, GPtrArray *theta, GPtrArray *m2lnL, NcmVector *thetastar, guint k)
{
  NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (walker);
  NcmVector *theta_k           = g_ptr_array_index (theta, k);
  NcmVector *thetabar_k        = g_ptr_array_index (walk->thetabar, k);
  guint i;

  ncm_vector_set_zero (thetabar_k);

  for (i = 0; i < walk->nparams; i++)
  {
    const guint j      = g_array_index (walk->indices, guint, k * walk->nparams + i);
    NcmVector *theta_j = g_ptr_array_index (theta, j);

    /*printf ("walker %u using walker %u to calculate mean\n", k, j);*/

    ncm_vector_add (thetabar_k, theta_j);
  }

  ncm_vector_scale (thetabar_k, 1.0 / (1.0 * walk->nparams));

  ncm_vector_memcpy (thetastar, theta_k);

  for (i = 0; i < walk->nparams; i++)
  {
    const gdouble z    = ncm_matrix_get (walk->z, k, i);
    const guint j      = g_array_index (walk->indices, guint, k * walk->nparams + i);
    NcmVector *theta_j = g_ptr_array_index (theta, j);
    guint m;

    for (m = 0; m < walk->nparams; m++)
    {
      const gdouble thetabar_k_i = ncm_vector_get (thetabar_k, i);
      const gdouble theta_j_i    = ncm_vector_get (theta_j, i);

      ncm_vector_addto (thetastar, i, walk->a * z * (theta_j_i - thetabar_k_i) / walk->sqrt_nparams);
    }
  }
}

static gdouble
_ncm_fit_esmcmc_walker_walk_prob (NcmFitESMCMCWalker *walker, GPtrArray *theta, GPtrArray *m2lnL, NcmVector *thetastar, guint k, const gdouble m2lnL_cur, const gdouble m2lnL_star)
{
  return exp ((m2lnL_cur - m2lnL_star) * 0.5);
}

static gdouble
_ncm_fit_esmcmc_walker_walk_prob_norm (NcmFitESMCMCWalker *walker, GPtrArray *theta, GPtrArray *m2lnL, NcmVector *thetastar, guint k)
{
  return 0.0;
}

static void
_ncm_fit_esmcmc_walker_walk_clean (NcmFitESMCMCWalker *walker, guint ki, guint kf)
{
  /* Nothing to do. */
}

const gchar *
_ncm_fit_esmcmc_walker_walk_desc (NcmFitESMCMCWalker *walker)
{
  return "Walk-Move";
}

/**
 * ncm_fit_esmcmc_walker_walk_new:
 * @nwalkers: number of walkers
 *
 * Creates a new #NcmFitESMCMCWalkerWalk to be used
 * with @nwalkers.
 *
 * Returns: (transfer full): a new #NcmFitESMCMCWalkerWalk.
 */
NcmFitESMCMCWalkerWalk *
ncm_fit_esmcmc_walker_walk_new (guint nwalkers)
{
  NcmFitESMCMCWalkerWalk *walk = g_object_new (NCM_TYPE_FIT_ESMCMC_WALKER_WALK,
                                               "size", nwalkers,
                                               NULL);

  return walk;
}

/**
 * ncm_fit_esmcmc_walker_walk_set_scale:
 * @walk: a #NcmFitESMCMCWalkerWalk
 * @a: new scale $a > 1$
 *
 * Sets the value of the scale $a > 1$.
 *
 */
void
ncm_fit_esmcmc_walker_walk_set_scale (NcmFitESMCMCWalkerWalk *walk, const gdouble a)
{
  g_assert_cmpfloat (a, >, 1.0e-2);
  walk->a = a;
}

/**
 * ncm_fit_esmcmc_walker_walk_get_scale:
 * @walk: a #NcmFitESMCMCWalkerWalk
 *
 * Gets the value of the scale $a > 1$.
 *
 * Returns: current value of $a$.
 */
gdouble
ncm_fit_esmcmc_walker_walk_get_scale (NcmFitESMCMCWalkerWalk *walk)
{
  return walk->a;
}


1	/***************************************************************************
2	* ncm_fit_esmcmc_walker_walk.c
3	*
4	* Tue March 29 10:41:49 2016
5	* Copyright 2016 Sandro Dias Pinto Vitenti
6	* <vitenti@uel.br>
7	****************************************************************************/
8	/*
9	* ncm_fit_esmcmc_walker_walk.c
10	* Copyright (C) 2016 Sandro Dias Pinto Vitenti <vitenti@uel.br>
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:ncm_fit_esmcmc_walker_walk
28	* @title: NcmFitESMCMCWalkerWalk
29	* @short_description: Ensemble sampler Markov Chain Monte Carlo walker - walk move.
30	* @stability: Unstable
31	*
32	* Implementing walk move walker for #NcmFitESMCMC (affine invariant).
33	*
34	*/
35
36	#ifdef HAVE_CONFIG_H
37	# include "config.h"
38	#endif /* HAVE_CONFIG_H */
39	#include "build_cfg.h"
40
41	#include "math/ncm_fit_esmcmc_walker.h"
42	#include "math/ncm_fit_esmcmc_walker_walk.h"
43
44	#include "math/ncm_fit_esmcmc.h"
45
46	#ifndef NUMCOSMO_GIR_SCAN
47	#include <gsl/gsl_math.h>
48	#endif /* NUMCOSMO_GIR_SCAN */
49
50	enum
51	{
52	PROP_0,
53	PROP_SCALE
54	};
55
56	struct _NcmFitESMCMCWalkerWalk
57	{
58	/< private >/
59	NcmFitESMCMCWalker parent_instance;
60	guint size;
61	guint size_2;
62	guint nparams;
63	gdouble a;
64	gdouble sqrt_nparams;
65	NcmMatrix *z;
66	GPtrArray *thetabar;
67	GArray *indices;
68	GArray *numbers;
69	};
70
71	G_DEFINE_TYPE (NcmFitESMCMCWalkerWalk, ncm_fit_esmcmc_walker_walk, NCM_TYPE_FIT_ESMCMC_WALKER)	9!
72
73	static void
74	ncm_fit_esmcmc_walker_walk_init (NcmFitESMCMCWalkerWalk *walk)	×
75	{
76	walk->size = 0;	×
77	walk->size_2 = 0;	×
78	walk->nparams = 0;	×
79	walk->a = 0.0;	×
80	walk->sqrt_nparams = 0.0;	×
81	walk->z = NULL;	×
82	walk->thetabar = g_ptr_array_new ();	×
83	walk->indices = g_array_new (TRUE, TRUE, sizeof (guint));	×
84	walk->numbers = g_array_new (TRUE, TRUE, sizeof (guint));	×
85
86	g_ptr_array_set_free_func (walk->thetabar, (GDestroyNotify) ncm_vector_free);	×
87	}	×
88
89	static void
90	ncm_fit_esmcmc_walker_walk_dispose (GObject *object)	×
91	{
92	NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (object);	×
93
94	ncm_matrix_clear (&walk->z);	×
95
96	g_clear_pointer (&walk->thetabar, g_ptr_array_unref);	×
97	g_clear_pointer (&walk->indices, g_array_unref);	×
98	g_clear_pointer (&walk->numbers, g_array_unref);	×
99
100	/* Chain up : end */
101	G_OBJECT_CLASS (ncm_fit_esmcmc_walker_walk_parent_class)->dispose (object);	×
102	}	×
103
104	static void
105	ncm_fit_esmcmc_walker_walk_finalize (GObject *object)	×
106	{
107	/* Chain up : end */
108	G_OBJECT_CLASS (ncm_fit_esmcmc_walker_walk_parent_class)->finalize (object);	×
109	}	×
110
111	static void
112	ncm_fit_esmcmc_walker_walk_set_property (GObject object, guint prop_id, const GValue value, GParamSpec *pspec)	×
113	{
114	NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (object);	×
115
116	g_return_if_fail (NCM_IS_FIT_ESMCMC_WALKER_WALK (object));	×
117
118	switch (prop_id)	×
119	{
120	case PROP_SCALE:	×
121	ncm_fit_esmcmc_walker_walk_set_scale (walk, g_value_get_double (value));	×
122	break;	×
123	default:	×
124	G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
125	break;	×
126	}
127	}
128
129	static void
130	ncm_fit_esmcmc_walker_walk_get_property (GObject object, guint prop_id, GValue value, GParamSpec *pspec)	×
131	{
132	NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (object);	×
133
134	g_return_if_fail (NCM_IS_FIT_ESMCMC_WALKER_WALK (object));	×
135
136	switch (prop_id)	×
137	{
138	case PROP_SCALE:	×
139	g_value_set_double (value, ncm_fit_esmcmc_walker_walk_get_scale (walk));	×
140	break;	×
141	default:	×
142	G_OBJECT_WARN_INVALID_PROPERTY_ID (object, prop_id, pspec);
143	break;	×
144	}
145	}
146
147	static void _ncm_fit_esmcmc_walker_walk_set_size (NcmFitESMCMCWalker *walker, guint size);
148	static guint _ncm_fit_esmcmc_walker_walk_get_size (NcmFitESMCMCWalker *walker);
149	static void _ncm_fit_esmcmc_walker_walk_set_nparams (NcmFitESMCMCWalker *walker, guint nparams);
150	static guint _ncm_fit_esmcmc_walker_walk_get_nparams (NcmFitESMCMCWalker *walker);
151	static void _ncm_fit_esmcmc_walker_walk_setup (NcmFitESMCMCWalker walker, NcmMSet mset, GPtrArray theta, GPtrArray m2lnL, guint ki, guint kf, NcmRNG *rng);
152	static void _ncm_fit_esmcmc_walker_walk_step (NcmFitESMCMCWalker walker, GPtrArray theta, GPtrArray m2lnL, NcmVector thetastar, guint k);
153	static gdouble _ncm_fit_esmcmc_walker_walk_prob (NcmFitESMCMCWalker walker, GPtrArray theta, GPtrArray m2lnL, NcmVector thetastar, guint k, const gdouble m2lnL_cur, const gdouble m2lnL_star);
154	static gdouble _ncm_fit_esmcmc_walker_walk_prob_norm (NcmFitESMCMCWalker walker, GPtrArray theta, GPtrArray m2lnL, NcmVector thetastar, guint k);
155	static void _ncm_fit_esmcmc_walker_walk_clean (NcmFitESMCMCWalker *walker, guint ki, guint kf);
156	static const gchar _ncm_fit_esmcmc_walker_walk_desc (NcmFitESMCMCWalker walker);
157
158	static void
159	ncm_fit_esmcmc_walker_walk_class_init (NcmFitESMCMCWalkerWalkClass *klass)	3✔
160	{
161	GObjectClass *object_class = G_OBJECT_CLASS (klass);	3✔
162	NcmFitESMCMCWalkerClass *walker_class = NCM_FIT_ESMCMC_WALKER_CLASS (klass);	3✔
163
164	object_class->set_property = ncm_fit_esmcmc_walker_walk_set_property;	3✔
165	object_class->get_property = ncm_fit_esmcmc_walker_walk_get_property;	3✔
166	object_class->dispose = ncm_fit_esmcmc_walker_walk_dispose;	3✔
167	object_class->finalize = ncm_fit_esmcmc_walker_walk_finalize;	3✔
168
169	g_object_class_install_property (object_class,	3✔
170	PROP_SCALE,
171	g_param_spec_double ("scale",
172	NULL,
173	"Walk scale a",
174	1.0e-2, G_MAXDOUBLE, 0.2,
175	G_PARAM_READWRITE \| G_PARAM_CONSTRUCT \| G_PARAM_STATIC_NAME \| G_PARAM_STATIC_BLURB));
176
177	walker_class->set_size = &_ncm_fit_esmcmc_walker_walk_set_size;	3✔
178	walker_class->get_size = &_ncm_fit_esmcmc_walker_walk_get_size;	3✔
179	walker_class->set_nparams = &_ncm_fit_esmcmc_walker_walk_set_nparams;	3✔
180	walker_class->get_nparams = &_ncm_fit_esmcmc_walker_walk_get_nparams;	3✔
181	walker_class->setup = &_ncm_fit_esmcmc_walker_walk_setup;	3✔
182	walker_class->step = &_ncm_fit_esmcmc_walker_walk_step;	3✔
183	walker_class->prob = &_ncm_fit_esmcmc_walker_walk_prob;	3✔
184	walker_class->prob_norm = &_ncm_fit_esmcmc_walker_walk_prob_norm;	3✔
185	walker_class->clean = &_ncm_fit_esmcmc_walker_walk_clean;	3✔
186	walker_class->desc = &_ncm_fit_esmcmc_walker_walk_desc;	3✔
187	}	3✔
188
189	static void
190	_ncm_fit_esmcmc_walker_walk_set_sys (NcmFitESMCMCWalker *walker, guint size, guint nparams)	×
191	{
192	NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (walker);	×
193
194	g_assert_cmpuint (size, >, 0);	×
195	g_assert_cmpuint (nparams, >, 0);	×
196
197	if ((walk->size != size) \|\| (walk->nparams != nparams))	×
198	{
UNCOV 199	guint i;	×
200
201	ncm_matrix_clear (&walk->z);	×
202
203	g_assert (size % 2 == 0);	×
204
205	walk->z = ncm_matrix_new (size, nparams);	×
206
207	g_array_set_size (walk->indices, size * nparams);	×
208
209	g_array_set_size (walk->numbers, size);	×
210
211	for (i = 0; i < size; i++)	×
212	g_array_index (walk->numbers, guint, i) = i;	×
213
214	g_ptr_array_set_size (walk->thetabar, 0);	×
215
216	for (i = 0; i < size; i++)	×
217	{
218	NcmVector *thetabar_i = ncm_vector_new (nparams);	×
219
220	g_ptr_array_add (walk->thetabar, thetabar_i);	×
221	}
222
223	walk->size = size;	×
224	walk->size_2 = size / 2;	×
225
226	walk->size = size;	×
227	walk->nparams = nparams;	×
228	walk->sqrt_nparams = sqrt (1.0 * walk->nparams);	×
229	}
230	}	×
231
232	static void
233	_ncm_fit_esmcmc_walker_walk_set_size (NcmFitESMCMCWalker *walker, guint size)	×
234	{
235	NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (walker);	×
236
237	g_assert_cmpuint (size, >, 0);	×
238
239	if (walk->nparams != 0)	×
240	_ncm_fit_esmcmc_walker_walk_set_sys (walker, size, walk->nparams);	×
241	else
242	walk->size = size;	×
243	}	×
244
245	static guint
246	_ncm_fit_esmcmc_walker_walk_get_size (NcmFitESMCMCWalker *walker)	×
247	{
248	NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (walker);	×
249
250	return walk->size;	×
251	}
252
253	static void
254	_ncm_fit_esmcmc_walker_walk_set_nparams (NcmFitESMCMCWalker *walker, guint nparams)	×
255	{
256	NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (walker);	×
257
258	g_assert_cmpuint (nparams, >, 0);	×
259
260	if (walk->size != 0)	×
261	_ncm_fit_esmcmc_walker_walk_set_sys (walker, walk->size, nparams);	×
262	else
263	walk->nparams = nparams;	×
264	}	×
265
266	static guint
267	_ncm_fit_esmcmc_walker_walk_get_nparams (NcmFitESMCMCWalker *walker)	×
268	{
269	NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (walker);	×
270
271	return walk->nparams;	×
272	}
273
274	static void
275	_ncm_fit_esmcmc_walker_walk_setup (NcmFitESMCMCWalker walker, NcmMSet mset, GPtrArray theta, GPtrArray m2lnL, guint ki, guint kf, NcmRNG *rng)	×
276	{
277	NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (walker);	×
UNCOV 278	guint k;	×
279
280	for (k = ki; k < kf; k++)	×
281	{
282	const guint subensemble = (k < walk->size_2) ? walk->size_2 : 0;	×
UNCOV 283	guint pi;	×
284
285	for (pi = 0; pi < walk->nparams; pi++)	×
286	{
287	const gdouble zi = gsl_ran_ugaussian (rng->r);	×
288
289	ncm_matrix_set (walk->z, k, pi, zi);	×
290	}
291
292	gsl_ran_choose (rng->r,	×
293	&g_array_index (walk->indices, guint, k * walk->nparams),	×
294	walk->nparams,	×
295	&g_array_index (walk->numbers, guint, subensemble),	×
296	walk->size_2,	×
297	g_array_get_element_size (walk->numbers));	×
298
299	/*
300	* for (pi = 0; pi < walk->nparams; pi++)
301	* {
302	* printf ("Choosing for walker %u, walker %u \| nparams %u size_2 %u subensemble %u \| numbers[subensemble + %u] = %u\n",
303	* k, g_array_index (walk->indices, guint, k * walk->nparams + pi),
304	* walk->nparams, walk->size_2, subensemble,
305	* pi, g_array_index (walk->numbers, guint, subensemble + pi));
306	* }
307	*/
308	}
309	}	×
310
311	static void
312	_ncm_fit_esmcmc_walker_walk_step (NcmFitESMCMCWalker walker, GPtrArray theta, GPtrArray m2lnL, NcmVector thetastar, guint k)	×
313	{
314	NcmFitESMCMCWalkerWalk *walk = NCM_FIT_ESMCMC_WALKER_WALK (walker);	×
315	NcmVector *theta_k = g_ptr_array_index (theta, k);	×
316	NcmVector *thetabar_k = g_ptr_array_index (walk->thetabar, k);	×
UNCOV 317	guint i;	×
318
319	ncm_vector_set_zero (thetabar_k);	×
320
321	for (i = 0; i < walk->nparams; i++)	×
322	{
323	const guint j = g_array_index (walk->indices, guint, k * walk->nparams + i);	×
324	NcmVector *theta_j = g_ptr_array_index (theta, j);	×
325
326	/printf ("walker %u using walker %u to calculate mean\n", k, j);/
327
328	ncm_vector_add (thetabar_k, theta_j);	×
329	}
330
331	ncm_vector_scale (thetabar_k, 1.0 / (1.0 * walk->nparams));	×
332
333	ncm_vector_memcpy (thetastar, theta_k);	×
334
335	for (i = 0; i < walk->nparams; i++)	×
336	{
337	const gdouble z = ncm_matrix_get (walk->z, k, i);	×
338	const guint j = g_array_index (walk->indices, guint, k * walk->nparams + i);	×
339	NcmVector *theta_j = g_ptr_array_index (theta, j);	×
UNCOV 340	guint m;	×
341
342	for (m = 0; m < walk->nparams; m++)	×
343	{
344	const gdouble thetabar_k_i = ncm_vector_get (thetabar_k, i);	×
345	const gdouble theta_j_i = ncm_vector_get (theta_j, i);	×
346
347	ncm_vector_addto (thetastar, i, walk->a * z * (theta_j_i - thetabar_k_i) / walk->sqrt_nparams);	×
348	}
349	}
350	}	×
351
352	static gdouble
353	_ncm_fit_esmcmc_walker_walk_prob (NcmFitESMCMCWalker walker, GPtrArray theta, GPtrArray m2lnL, NcmVector thetastar, guint k, const gdouble m2lnL_cur, const gdouble m2lnL_star)	×
354	{
355	return exp ((m2lnL_cur - m2lnL_star) * 0.5);	×
356	}
357
358	static gdouble
359	_ncm_fit_esmcmc_walker_walk_prob_norm (NcmFitESMCMCWalker walker, GPtrArray theta, GPtrArray m2lnL, NcmVector thetastar, guint k)	×
360	{
361	return 0.0;	×
362	}
363
364	static void
365	_ncm_fit_esmcmc_walker_walk_clean (NcmFitESMCMCWalker *walker, guint ki, guint kf)	×
366	{
367	/* Nothing to do. */
368	}	×
369
370	const gchar *
371	_ncm_fit_esmcmc_walker_walk_desc (NcmFitESMCMCWalker *walker)	×
372	{
373	return "Walk-Move";	×
374	}
375
376	/**
377	* ncm_fit_esmcmc_walker_walk_new:
378	* @nwalkers: number of walkers
379	*
380	* Creates a new #NcmFitESMCMCWalkerWalk to be used
381	* with @nwalkers.
382	*
383	* Returns: (transfer full): a new #NcmFitESMCMCWalkerWalk.
384	*/
385	NcmFitESMCMCWalkerWalk *
386	ncm_fit_esmcmc_walker_walk_new (guint nwalkers)	×
387	{
388	NcmFitESMCMCWalkerWalk *walk = g_object_new (NCM_TYPE_FIT_ESMCMC_WALKER_WALK,	×
389	"size", nwalkers,
390	NULL);
391
392	return walk;	×
393	}
394
395	/**
396	* ncm_fit_esmcmc_walker_walk_set_scale:
397	* @walk: a #NcmFitESMCMCWalkerWalk
398	* @a: new scale $a > 1$
399	*
400	* Sets the value of the scale $a > 1$.
401	*
402	*/
403	void
404	ncm_fit_esmcmc_walker_walk_set_scale (NcmFitESMCMCWalkerWalk *walk, const gdouble a)	×
405	{
406	g_assert_cmpfloat (a, >, 1.0e-2);	×
407	walk->a = a;	×
408	}	×
409
410	/**
411	* ncm_fit_esmcmc_walker_walk_get_scale:
412	* @walk: a #NcmFitESMCMCWalkerWalk
413	*
414	* Gets the value of the scale $a > 1$.
415	*
416	* Returns: current value of $a$.
417	*/
418	gdouble
419	ncm_fit_esmcmc_walker_walk_get_scale (NcmFitESMCMCWalkerWalk *walk)	×
420	{
421	return walk->a;	×
422	}
423

NumCosmo / NumCosmo / 7062014968

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