SciUnit is a framework for validating scientific models by creating experimental-data-driven unit tests.¶
Chapter 5. The Real Example of Using SciUnit To Test Cosmology Models¶
(or back to Chapter 4)
This is a real and executable example of testing 5 different cosmology models with SciUnit.
If you are using this notebook file in Google Colab, this block of code can help you install sciunit from PyPI in Colab environment.¶
In [ ]:
try:
import google.colab
IN_COLAB = True
except:
IN_COLAB = False
if IN_COLAB:
!pip install -q sciunit
Like what we do before, let's create some capabilities, tests, models, and test suites.
In [ ]:
import sciunit
from sciunit import Test, Model, Capability, TestSuite
from sciunit.errors import PredictionError
from sciunit.scores import RatioScore, BooleanScore
from sciunit.converters import RangeToBoolean
import quantities as pq
Capabilities¶
In [ ]:
class HasSun(Capability):
def solar_days(self):
return self.unimplemented()
class HasStars(Capability):
def stellar_parallax(self,star):
return self.unimplemented()
class HasPlanets(Capability):
def orbital_eccentricity(self,planet):
return self.unimplemented()
def num_moons(self,planet):
return self.unimplemented()
def perihelion_precession_rate(self,planet):
return self.unimplemented()
Models¶
In [ ]:
class _CosmoModel(Model):
def solar_year_duration(self):
raise PredictionError(self,self.curr_method(back=1))
def orbital_eccentricity(self, planet):
raise PredictionError(self,self.curr_method(back=1),planet=planet)
def num_moons(self, planet):
raise PredictionError(self,self.curr_method(back=1),planet=planet)
def perihelion_precession_rate(self, planet):
raise PredictionError(self,self.curr_method(back=1),planet=planet)
def stellar_parallax(self, star):
raise PredictionError(self,self.curr_method(back=1),star=star)
class Ptolemy(_CosmoModel, HasSun, HasPlanets):
"""Cladius Ptolemy, "The Almagest", 50 A.D."""
def solar_year_duration(self):
return 365 * pq.day
def orbital_eccentricity(self, planet):
return 0.0
def num_moons(self, planet):
if planet == 'Earth':
return 1
else:
return _CosmoModel.num_moons(self,planet)
def perihelion_precession_rate(self, planet):
return 0.0 * pq.Hz
class Copernicus(Ptolemy, HasStars):
"""Nicholas Copernicus, "De revolutionibus orbium coelestium", 1543"""
def solar_year_duration(self):
return 365.25 * pq.day
def stellar_parallax(self, star):
return 0.0 * pq.arcsecond
class Kepler(Copernicus):
"""Johannes Kepler, "Astronomia nova", 1609"""
def solar_year_duration(self):
return 365.25 * pq.day
def orbital_eccentricity(self, planet):
if planet == 'Mars':
return 0.0934
elif planet == 'Saturn':
return 0.0541
else:
return _CosmoModel.orbital_eccentricity(self,planet)
def num_moons(self, planet):
if planet == 'Jupiter':
return 4
elif planet == 'Earth':
return 1
else:
return _CosmoModel.num_moons(self,planet)
def perihelion_precession_rate(self, planet):
return 0.0 * pq.Hz
class Newton(Kepler):
"""Isaac Newton, "Philosophiae Naturalis Principia Mathematica", 1687"""
def perihelion_precession_rate(self, planet):
if planet == 'Mercury':
return (531.63 * pq.arcsecond)/(100.0 * pq.year)
else:
return _CosmoModel.perihelion_precession_rate(self,planet)
def stellar_parallax(self, star):
if star == '61 Cygni':
return 0.314 * pq.arcsecond
elif star == 'Promixa Centauri':
return 0.769 * pq.arcsecond
else:
raise _CosmoModel.stellar_parallax(self,star)
class Einstein(Newton):
"""Albert Einstein, "The Foundation of the General Theory of Relativity"
Annalen der Physik, 49(7):769-822, 1916."""
def perihelion_precession_rate(self, planet):
if planet == 'Mercury':
return (574.10 * pq.arcsecond)/(100.0 * pq.year)
else:
return _CosmoModel.perihelion_precession_rate(self,planet)
Tests¶
In [ ]:
class _CosmoTest(sciunit.Test):
score_type = BooleanScore
primary_key = None
units = pq.dimensionless
def validate_observation(self, observation):
"""Observation should be a dictionary of containing the length of a
a solar year in units with the dimension of time."""
key = self.primary_key
assert key in observation, "%s not found in %s test observation" % \
(key,self.__class__.__name__)
value = observation[key]
if type(observation[key]) is tuple:
value = value[1]
if self.units is not pq.dimensionless:
assert isinstance(value,pq.Quantity), \
("Key '%s' of observation for '%s' test is not an instance of "
"quantities.Quantity" ) % (key,self.__class__.__name__)
assert value.simplified.units == \
self.units.simplified.units, \
("Key '%s' of observation for '%s' test does not have units of "
"%s" % (key,self.__class__.__name__,self.units))
def compute_score(self, observation, prediction, verbose=True):
key = self.primary_key
obs,pred = observation[key],prediction[key]
if isinstance(self,_CosmoEntityTest):
obs = obs[1]
error = RatioScore.compute(obs,pred)
score = RangeToBoolean(0.97,1.03).convert(error) # +/- 3% of observed
return score
class _CosmoEntityTest(_CosmoTest):
entity_type = None
def validate_observation(self, observation):
super(_CosmoEntityTest,self).validate_observation(observation)
assert type(observation[self.primary_key]) is tuple, \
"Observation for key %s must be a (%s,value) tuple" % \
(self.entity_type,self.primary_key)
class SolarYear(_CosmoTest):
required_capabilities = [HasSun]
primary_key = 'duration'
units = pq.s
def generate_prediction(self, model, verbose=True):
days = model.solar_year_duration()
return {self.primary_key:days}
class OrbitalEccentricity(_CosmoEntityTest):
required_capabilities = [HasPlanets]
primary_key = 'eccentricity'
entity_type = 'planet'
units = pq.dimensionless
def generate_prediction(self, model, verbose=True):
planet,value = self.observation[self.primary_key]
eccentricity = model.orbital_eccentricity(planet)
return {self.primary_key:eccentricity}
class StellarParallax(_CosmoEntityTest):
required_capabilities = [HasStars]
primary_key = 'parallax'
units = pq.arcsecond
entity_type = 'star'
def generate_prediction(self, model, verbose=True):
star,value = self.observation[self.primary_key]
parallax = model.stellar_parallax(star)
return {self.primary_key:parallax}
class PerihelionPrecession(_CosmoEntityTest):
required_capabilities = [HasSun, HasPlanets]
primary_key = 'precession'
entity_type = 'planet'
units = pq.Hz
def generate_prediction(self, model, verbose=True):
planet,value = self.observation[self.primary_key]
precession = model.perihelion_precession_rate(planet)
return {self.primary_key:precession}
In [ ]:
# Orbital eccentricities
eccentricity = {'Mars':0.093,
'Saturn':0.0541506,
}
# Perihelion precessions
precession = {'Mercury':(574.10 * pq.arcsecond)/(100.0 * pq.year),
}
# Stellar parallaxes
parallax = {'61 Cygni':0.3136 * pq.arcsecond,
# Friedrich Bessel in 1838 using a heliometer.
# Bessel, Friedrich
# "Bestimmung der Entfernung des 61sten Sterns des Schwans"
# Astronomische Nachrichten, 16, 65-96 (1838)
'Promixa Centauri':0.7687 * pq.arcsecond,
}
solar_year_duration = 365.25 * pq.day
In [ ]:
planets = ['Mars', 'Saturn']
stars = ['Promixa Centauri', '61 Cygni']
babylon = SolarYear({'duration' : solar_year_duration}, name='Solar Year')
brahe = [OrbitalEccentricity({'eccentricity' : (planet, eccentricity[planet])},
name='Ecc. %s' % planet) \
for planet in planets]
bessel = [StellarParallax({'parallax' : (star, parallax[star])}, name='Prlx. %s' % star) \
for star in stars]
leverrier = PerihelionPrecession({'precession' : ('Mercury', precession['Mercury'])},
name='Phln. Mercury')
In [ ]:
babylon = TestSuite(tests=babylon, name='Babylon')
brahe = TestSuite(tests=brahe, name='Brahe')
bessel = TestSuite(tests=bessel, name='Bessel')
leverrier = TestSuite(tests=leverrier, name='Leverrier')
# Set these test suites to be applied to all models
suites = [babylon, brahe, bessel, leverrier]
And then, we can let each suite instance judge each model in the model list.¶
In [ ]:
ptolemy = Ptolemy()
copernicus = Copernicus()
kepler = Kepler()
newton = Newton()
einstein = Einstein()
models = [ptolemy,copernicus,kepler,newton,einstein]
In [ ]:
for suite in suites:
suite.judge(models)