Commit e9f715ea by Theo Steininger

### Implemented pymultinest MPI magic in Pipeline; still untested.

parent 374f0b38
 # -*- coding: utf-8 -*- # -*- coding: utf-8 -*- from likelihood import Likelihood from likelihood import Likelihood from ensemble_likelihood import EnsembleLikelihood from ensemble_likelihood import *
 ... @@ -2,24 +2,33 @@ ... @@ -2,24 +2,33 @@ import numpy as np import numpy as np from imagine.likelihoods.likelihood import Likelihood class EnsembleLikelihood(object): class EnsembleLikelihood(Likelihood): def __init__(self, observable_name, measured_data, def __init__(self, measured_data, data_covariance_operator): data_covariance_operator): 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 def __call__(self, observable): def __call__(self, observable): field = observable[self.observable_name] return self._process_simple_field(field, self.measured_data, self.data_covariance_operator) def _process_simple_field(self, field, measured_data, data_covariance_operator): # 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 k = observable.shape[0] k = observable.shape[0] A = self.data_covariance_operator A = data_covariance_operator obs_val = observable.val.get_full_data() obs_val = observable.val.get_full_data() obs_mean = observable.mean(spaces=0).val.get_full_data() obs_mean = observable.mean(spaces=0).val.get_full_data() ... @@ -36,7 +45,7 @@ class EnsembleLikelihood(Likelihood): ... @@ -36,7 +45,7 @@ class EnsembleLikelihood(Likelihood): middle = np.linalg.inv(middle) middle = np.linalg.inv(middle) result_array = np.zeros(k) result_array = np.zeros(k) for i in xrange(k): for i in xrange(k): c = self.measured_data - obs_val[i] c = measured_data - obs_val[i] # 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) ... ...
 # -*- coding: utf-8 -*- # -*- coding: utf-8 -*- import os import numpy as np from mpi4py import MPI from keepers import Loggable from keepers import Loggable from likelihoods import Likelihood 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 comm = MPI.COMM_WORLD size = comm.size rank = comm.rank class Pipeline(Loggable, object): class Pipeline(Loggable, object): ... @@ -111,5 +121,51 @@ class Pipeline(Loggable, object): ... @@ -111,5 +121,51 @@ class Pipeline(Loggable, object): self.logger.debug("Setting ensemble size to %i." % ensemble_size) self.logger.debug("Setting ensemble size to %i." % ensemble_size) self._ensemble_size = ensemble_size self._ensemble_size = ensemble_size def _multinest_likelihood(self, cube, ndim, nparams): cube_content = np.empty(ndim) for i in xrange(ndim): cube_content[i] = cube[i] if rank != 0: raise RuntimeError("_multinest_likelihood must only be called on " "rank==0.") for i in xrange(1, size): comm.send(cube_content, dest=i) return self._core_likelihood(cube_content) def _listen_for_likelihood_calls(self): cube = comm.recv(obj=None, source=0) self._core_likelihood(cube) def _core_likelihood(self, cube): # translate cube to variables variables = {} for i, av in enumerate(self.active_variables): variables[av] = cube[i] # create magnetic field b_field = self.magnetic_field_factory(variables=variables, ensemble_size=self.ensemle_size) # create observables observables = self.observer(b_field) # add up individual log-likelihood terms likelihood = 0 for like in self.likelihood: likelihood += like(observables) return likelihood def __call__(self, variables): def __call__(self, variables): pass if rank == 0: # kickstart pymultinest if not os.path.exists("chains"): os.mkdir("chains") pymultinest.run(self._multinest_likelihood, self.prior, len(self.active_variables), verbose=True) else: # let all other nodes listen for likelihood evaluations self._listen_for_likelihood_calls()
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!