Skip to content

Commit

Permalink
Added support for namespace search in ncm_mset_func_list. Added plot_…
Browse files Browse the repository at this point in the history
…corner helper script.
  • Loading branch information
vitenti committed Nov 7, 2022
1 parent 392bd34 commit 1d6c568
Show file tree
Hide file tree
Showing 3 changed files with 156 additions and 10 deletions.
50 changes: 43 additions & 7 deletions numcosmo/math/ncm_mset_func_list.c
Original file line number Diff line number Diff line change
Expand Up @@ -295,7 +295,7 @@ ncm_mset_func_list_register (const gchar *name, const gchar *symbol, const gchar
* gets from all namespaces, @nvar and/or @dim equals to -1 selects
* any value. The contained strings must not be freed.
*
* Returns: (transfer full) (element-type NcmMSetFuncListStruct): NcmMSetFuncListStruct array.
* Returns: (transfer container) (element-type NcmMSetFuncListStruct): NcmMSetFuncListStruct array.
*/
GArray *
ncm_mset_func_list_select (const gchar *ns, gint nvar, gint dim)
Expand All @@ -309,7 +309,7 @@ ncm_mset_func_list_select (const gchar *ns, gint nvar, gint dim)
for (i = 0; i < flist_class->func_array->len; i++)
{
NcmMSetFuncListStruct *fdata = &g_array_index (flist_class->func_array, NcmMSetFuncListStruct, i);
gboolean in_ns = (ns == NULL) || g_strcmp0 (fdata->ns, ns) == 0;
gboolean in_ns = (ns == NULL) || g_str_has_prefix (fdata->ns, ns);
gboolean in_nvar = (nvar == -1) || nvar == fdata->nvar;
gboolean in_dim = (dim == -1) || dim == fdata->dim;

Expand All @@ -319,7 +319,7 @@ ncm_mset_func_list_select (const gchar *ns, gint nvar, gint dim)
}
}
}

G_UNLOCK (insert_lock);

g_type_class_unref (flist_class);
Expand Down Expand Up @@ -359,13 +359,41 @@ ncm_mset_func_list_new (const gchar *full_name, GObject *obj)
NcmMSetFuncList *
ncm_mset_func_list_new_ns_name (const gchar *ns, const gchar *name, GObject *obj)
{
gchar *full_name = g_strdup_printf ("%s:%s", ns, name);
NcmMSetFuncListClass *flist_class = g_type_class_ref (NCM_TYPE_MSET_FUNC_LIST);
const gchar *full_ns = NULL;

G_LOCK (insert_lock);
{
gint i;
for (i = 0; i < flist_class->func_array->len; i++)
{
NcmMSetFuncListStruct *fdata = &g_array_index (flist_class->func_array, NcmMSetFuncListStruct, i);
if (g_str_has_prefix (fdata->ns, ns))
{
GHashTable *func_hash = g_hash_table_lookup (flist_class->ns_hash, fdata->ns);
if ((func_hash != NULL) && g_hash_table_lookup_extended (func_hash, name, NULL, NULL))
{
full_ns = fdata->ns;
break;
}
}
}
}
G_UNLOCK (insert_lock);
g_type_class_unref (flist_class);

if (full_ns == NULL)
g_error ("ncm_mset_func_list_new_ns_name: namespace `%s' not found.", ns);

{
gchar *full_name = g_strdup_printf ("%s:%s", full_ns, name);
NcmMSetFuncList *flist = g_object_new (NCM_TYPE_MSET_FUNC_LIST,
"full-name", full_name,
"object", obj,
NULL);
g_free (full_name);
return flist;
}
}

/**
Expand All @@ -385,9 +413,17 @@ ncm_mset_func_list_has_ns_name (const gchar *ns, const gchar *name)
G_LOCK (insert_lock);

{
GHashTable *func_hash = g_hash_table_lookup (flist_class->ns_hash, ns);
if ((func_hash != NULL) && g_hash_table_lookup_extended (func_hash, name, NULL, NULL))
has_func = TRUE;
gint i;
for (i = 0; i < flist_class->func_array->len; i++)
{
NcmMSetFuncListStruct *fdata = &g_array_index (flist_class->func_array, NcmMSetFuncListStruct, i);
if (g_str_has_prefix (fdata->ns, ns))
{
GHashTable *func_hash = g_hash_table_lookup (flist_class->ns_hash, fdata->ns);
if ((func_hash != NULL) && g_hash_table_lookup_extended (func_hash, name, NULL, NULL))
has_func = TRUE;
}
}
}

G_UNLOCK (insert_lock);
Expand Down
16 changes: 13 additions & 3 deletions numcosmo/model/nc_hicosmo_de.c
Original file line number Diff line number Diff line change
Expand Up @@ -1154,10 +1154,20 @@ _nc_hicosmo_de_flist_BBN (NcmMSetFuncList *flist, NcmMSet *mset, const gdouble *
f[0] = 1.0 / sqrt (1.0 - Omega_de);
}

static void
_nc_hicosmo_de_flist_w (NcmMSetFuncList *flist, NcmMSet *mset, const gdouble *x, gdouble *f)
{
NcHICosmoDE *cosmo_de = NC_HICOSMO_DE (ncm_mset_peek (mset, nc_hicosmo_id ()));
g_assert (NC_IS_HICOSMO_DE (cosmo_de));
f[0] = nc_hicosmo_de_w_de (cosmo_de, x[0]);
}

void
_nc_hicosmo_de_register_functions (void)
{
ncm_mset_func_list_register ("Omega_x0", "\\Omega_{x0}", "NcHICosmoDE", "Darkenergy density today", G_TYPE_NONE, _nc_hicosmo_de_flist_Omega_x0, 0, 1);
ncm_mset_func_list_register ("wDE", "\\omega_\\mathrm{de}", "NcHICosmoDE", "Darkenergy equation of state today", G_TYPE_NONE, _nc_hicosmo_de_flist_w0, 0, 1);
ncm_mset_func_list_register ("BBN", "BBN", "NcHICosmoDE", "BBN", G_TYPE_NONE, _nc_hicosmo_de_flist_BBN, 0, 1);
ncm_mset_func_list_register ("Omega_x0", "\\Omega_{x0}", "NcHICosmoDE", "Darkenergy density today", G_TYPE_NONE, _nc_hicosmo_de_flist_Omega_x0, 0, 1);
ncm_mset_func_list_register ("wDE", "\\omega_\\mathrm{de}", "NcHICosmoDE", "Darkenergy equation of state today", G_TYPE_NONE, _nc_hicosmo_de_flist_w0, 0, 1);
ncm_mset_func_list_register ("BBN", "BBN", "NcHICosmoDE", "BBN", G_TYPE_NONE, _nc_hicosmo_de_flist_BBN, 0, 1);

ncm_mset_func_list_register ("wDE_z", "\\omega_\\mathrm{de}(z)", "NcHICosmoDE", "Darkenergy equation of state", G_TYPE_NONE, _nc_hicosmo_de_flist_w, 1, 1);
}
100 changes: 100 additions & 0 deletions scripts/plot_corner.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,100 @@
#!/usr/bin/env python

try:
import gi
gi.require_version('NumCosmo', '1.0')
gi.require_version('NumCosmoMath', '1.0')
except:
pass

import math
import argparse
import numpy as np

from gi.repository import NumCosmo as Nc
from gi.repository import NumCosmoMath as Ncm

from scipy.stats import chi2
from chainconsumer import ChainConsumer

parser = argparse.ArgumentParser(description='Process mset catalogs')

parser.add_argument('-C', '--catalog', metavar='file.fits',
help='catalog fits file',
action='append', required = True)

parser.add_argument('-B', '--burnin', metavar='N',
help='catalog burnin', type=int,
action='append')

parser.add_argument('--kde',
help='whether to use kde interpolation',
action='store_true')

parser.add_argument('--col', type=int, nargs="*",
help='Columns to include')

parser.add_argument('--truth', type=float, nargs="*",
help='Columns to include')

parser.add_argument('--sigma', type=int, nargs="+", default = [1, 2],
help='Sigmas to compute the confidence regions')

parser.add_argument('--out', default = "corner.pdf",
help='Output filename')

parser.add_argument('--mode', choices=['corner', 'walks'], default='corner')

args = parser.parse_args()

Ncm.cfg_init ()

c = ChainConsumer()

for cat in args.catalog:

bin = 0
if args.burnin and (len (args.burnin) > 0):
bin = args.burnin.pop (0)

print (f"# Adding {cat} with burnin {bin}")

mcat = Ncm.MSetCatalog.new_from_file_ro (cat, bin)
nwalkers = mcat.nchains ()

m2lnL = mcat.get_m2lnp_var ()

rows = np.array ([mcat.peek_row (i).dup_array () for i in range (mcat.len ())])
params = ["$" + mcat.col_symb (i) + "$" for i in range (mcat.ncols ())]

posterior = -0.5 * rows[:,m2lnL]

rows = np.delete (rows, m2lnL, 1)
params = np.delete (params, m2lnL, 0)

if args.col:
assert max (args.col) < mcat.ncols ()
indices = np.array (args.col)

rows = rows[:,indices]
params = params[indices]

#c.add_chain(rows, posterior = posterior, parameters=list(params), plot_point = True, name = cat.replace ("_", "-"))
c.add_chain(rows, posterior = posterior, parameters=list(params), name = cat.replace ("_", "-"))

c.configure (kde = args.kde, label_font_size=8, sigma2d=False, sigmas = args.sigma, spacing = 0.0, tick_font_size=8)

plot_args = {}

if args.truth is not None:
plot_args['truth'] = args.truth

if args.mode == "corner":
fig = c.plotter.plot(**plot_args)
elif args.mode == "walks":
fig = c.plotter.plot_walks(**plot_args, convolve=100)
else:
assert False

fig.savefig(args.out, bbox_inches="tight")

0 comments on commit 1d6c568

Please sign in to comment.