Commit ea8d0693 authored by Theo Steininger's avatar Theo Steininger
Browse files

Fixed Pipeline -> pymultinest now runs with NIFTy-MPI

Added Healpix ring profile function.
parent e9f715ea
...@@ -6,7 +6,7 @@ from magnetic_fields import * ...@@ -6,7 +6,7 @@ from magnetic_fields import *
from observers import * from observers import *
from priors import * from priors import *
from pymultinest import pymultinest from pymultinest_importer import pymultinest
from pipeline import Pipeline from pipeline import Pipeline
......
# -*- coding: utf-8 -*-
import numpy as np
import healpy as hp
def create_ring_profile(input_map):
input_map = np.abs(input_map)
npix = input_map.shape[0]
nside = hp.npix2nside(npix)
rings = hp.pix2ring(nside, np.arange(npix))
rho = np.bincount(rings)[1:]
averages = np.bincount(rings, weights=input_map)[1:]/rho
result = averages[rings]
return result
...@@ -2,28 +2,33 @@ ...@@ -2,28 +2,33 @@
import numpy as np import numpy as np
from imagine.likelihoods.likelihood import Likelihood
from imagine.create_ring_profile import create_ring_profile
class EnsembleLikelihood(object):
class EnsembleLikelihood(Likelihood):
def __init__(self, observable_name, measured_data, def __init__(self, observable_name, measured_data,
data_covariance_operator): data_covariance_operator):
self.observable_name = observable_name self.observable_name = observable_name
self.measured_data = measured_data self.measured_data = measured_data
self.data_covariance_operator = data_covariance_operator self.data_covariance_operator = data_covariance_operator
self.profile = create_ring_profile(
self.measured_data.val.get_full_data())
def __call__(self, observable): def __call__(self, observable):
field = observable[self.observable_name] field = observable[self.observable_name]
return self._process_simple_field(field, return self._process_simple_field(field,
self.measured_data, self.measured_data,
self.data_covariance_operator) self.data_covariance_operator,
self.profile)
def _process_simple_field(self, field, measured_data, def _process_simple_field(self, field, measured_data,
data_covariance_operator): data_covariance_operator, profile):
# https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula#Generalization # https://en.wikipedia.org/wiki/Sherman%E2%80%93Morrison_formula#Generalization
# B = A^{-1} + U U^dagger # B = A^{-1} + U U^dagger
# A = data_covariance # A = data_covariance
# B^{-1} c = (A_inv - # B^{-1} c = (A_inv -
# A_inv U (I_k + U^dagger A_inv U)^{-1} U^dagger A_inv) c # A_inv U (I_k + U^dagger A_inv U)^{-1} U^dagger A_inv) c
observable = field observable = field
k = observable.shape[0] k = observable.shape[0]
...@@ -46,6 +51,7 @@ class EnsembleLikelihood(object): ...@@ -46,6 +51,7 @@ class EnsembleLikelihood(object):
result_array = np.zeros(k) result_array = np.zeros(k)
for i in xrange(k): for i in xrange(k):
c = measured_data - obs_val[i] c = measured_data - obs_val[i]
c /= profile
# assuming that A == A^dagger, this can be shortend # assuming that A == A^dagger, this can be shortend
# a_c = A.inverse_times(c) # a_c = A.inverse_times(c)
......
...@@ -121,7 +121,8 @@ class MagneticFieldFactory(Loggable, object): ...@@ -121,7 +121,8 @@ class MagneticFieldFactory(Loggable, object):
domain = (ensemble, grid_space, vector) domain = (ensemble, grid_space, vector)
result_magnetic_field = self.magnetic_field_class( result_magnetic_field = self.magnetic_field_class(
domain=domain, domain=domain,
parameters=work_parameters) parameters=work_parameters,
distribution_strategy='equal')
return result_magnetic_field return result_magnetic_field
...@@ -136,7 +136,8 @@ class Hammurapy(Observer): ...@@ -136,7 +136,8 @@ class Hammurapy(Observer):
# get the local shape by creating a dummy d2o # get the local shape by creating a dummy d2o
ensemble_number = magnetic_field.shape[0] ensemble_number = magnetic_field.shape[0]
dummy = distributed_data_object(global_shape=(ensemble_number,), dummy = distributed_data_object(global_shape=(ensemble_number,),
distribution_strategy='equal') distribution_strategy='equal',
dtype=np.float)
local_length = dummy.distributor.local_length local_length = dummy.distributor.local_length
for local_ensemble_index in xrange(local_length): for local_ensemble_index in xrange(local_length):
......
...@@ -11,12 +11,15 @@ from likelihoods import Likelihood ...@@ -11,12 +11,15 @@ from likelihoods import Likelihood
from magnetic_fields import MagneticFieldFactory from magnetic_fields import MagneticFieldFactory
from observers import Observer from observers import Observer
from priors import Prior from priors import Prior
from imagine.pymultinest import pymultinest from imagine import pymultinest
comm = MPI.COMM_WORLD comm = MPI.COMM_WORLD
size = comm.size size = comm.size
rank = comm.rank rank = comm.rank
WORK_TAG = 0
DIE_TAG = 1
class Pipeline(Loggable, object): class Pipeline(Loggable, object):
""" """
...@@ -59,6 +62,9 @@ class Pipeline(Loggable, object): ...@@ -59,6 +62,9 @@ class Pipeline(Loggable, object):
def likelihood(self, likelihood): def likelihood(self, likelihood):
self.logger.debug("Setting likelihood.") self.logger.debug("Setting likelihood.")
self._likelihood = () self._likelihood = ()
if not (isinstance(likelihood, list) and
isinstance(likelihood, tuple)):
likelihood = [likelihood]
for l in likelihood: for l in likelihood:
if not isinstance(l, Likelihood): if not isinstance(l, Likelihood):
raise TypeError( raise TypeError(
...@@ -72,12 +78,10 @@ class Pipeline(Loggable, object): ...@@ -72,12 +78,10 @@ class Pipeline(Loggable, object):
@prior.setter @prior.setter
def prior(self, prior): def prior(self, prior):
self.logger.debug("Setting prior.") self.logger.debug("Setting prior.")
self._prior = () if not isinstance(prior, Prior):
for p in prior: raise TypeError(
if not isinstance(p, Prior): "prior must be an instance of Prior-class.")
raise TypeError( self._prior = prior
"prior must be an instance of Prior-class.")
self._prior += (p,)
@property @property
def magnetic_field_factory(self): def magnetic_field_factory(self):
...@@ -129,43 +133,63 @@ class Pipeline(Loggable, object): ...@@ -129,43 +133,63 @@ class Pipeline(Loggable, object):
raise RuntimeError("_multinest_likelihood must only be called on " raise RuntimeError("_multinest_likelihood must only be called on "
"rank==0.") "rank==0.")
for i in xrange(1, size): for i in xrange(1, size):
comm.send(cube_content, dest=i) comm.send(cube_content, dest=i, tag=WORK_TAG)
self.logger.debug("Sent multinest-cube to rank %i" % i)
return self._core_likelihood(cube_content) return self._core_likelihood(cube_content)
def _listen_for_likelihood_calls(self): def _listen_for_likelihood_calls(self):
cube = comm.recv(obj=None, source=0) status = MPI.Status()
self._core_likelihood(cube) while True:
cube = comm.recv(source=0, tag=MPI.ANY_TAG, status=status)
if status == DIE_TAG:
self.logger.debug("Received DIE_TAG from rank 0.")
break
self.logger.debug("Received cube from rank 0.")
self._core_likelihood(cube)
def _core_likelihood(self, cube): def _core_likelihood(self, cube):
self.logger.debug("Beginning Likelihood-calculation for %s." %
str(cube))
# translate cube to variables # translate cube to variables
variables = {} variables = {}
for i, av in enumerate(self.active_variables): for i, av in enumerate(self.active_variables):
variables[av] = cube[i] variables[av] = cube[i]
# create magnetic field # create magnetic field
b_field = self.magnetic_field_factory(variables=variables, self.logger.debug("Creating magnetic field.")
ensemble_size=self.ensemle_size) b_field = self.magnetic_field_factory.generate(
variables=variables,
ensemble_size=self.ensemble_size)
# create observables # create observables
self.logger.debug("Creating observables.")
observables = self.observer(b_field) observables = self.observer(b_field)
# add up individual log-likelihood terms # add up individual log-likelihood terms
self.logger.debug("Evaluating likelihood(s).")
likelihood = 0 likelihood = 0
for like in self.likelihood: for like in self.likelihood:
likelihood += like(observables) likelihood += like(observables)
self.logger.debug("Evaluated likelihood: %f" % likelihood)
return likelihood return likelihood
def __call__(self, variables): def __call__(self, variables):
if rank == 0: if rank == 0:
# kickstart pymultinest # kickstart pymultinest
self.logger.info("starting pymultinest.")
if not os.path.exists("chains"): if not os.path.exists("chains"):
os.mkdir("chains") os.mkdir("chains")
pymultinest.run(self._multinest_likelihood, pymultinest.run(self._multinest_likelihood,
self.prior, self.prior,
len(self.active_variables), len(self.active_variables),
verbose=True) verbose=True)
self.logger.info("pymultinest finished.")
for i in xrange(1, size):
self.logger.debug("Sending DIE_TAG to rank %i." % i)
comm.send(None, dest=i, tag=DIE_TAG)
else: else:
# let all other nodes listen for likelihood evaluations # let all other nodes listen for likelihood evaluations
self._listen_for_likelihood_calls() self._listen_for_likelihood_calls()
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
from pymultinest import pymultinest from pymultinest_importer import pymultinest
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment