5. 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.
5.1. If you are using this notebook file in Google Colab, this block of code can help you install sciunit from PyPI in Colab environment.¶
!pip install -q sciunit
Like what we do before, let’s create some capabilities, tests, models, and test suites.
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
5.2. Capabilities¶
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()
5.3. Models¶
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)
5.4. Tests¶
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}
5.5. Observations¶
5.5.1. The data collected from observations and experiments.¶
# 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
5.6. TestSuites¶
5.6.1. Let’s put the test instances into a poython list, and create the TestSuite instances with them.¶
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')
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]
5.6.2. And then, we can let each suite instance judge each model in the model list.¶
ptolemy = Ptolemy()
copernicus = Copernicus()
kepler = Kepler()
newton = Newton()
einstein = Einstein()
models = [ptolemy,copernicus,kepler,newton,einstein]
for suite in suites:
suite.judge(models)
The model class claimed to implement all methods required by the Test class, but at least one was left unimplemented, so this model will be skipped.
Score: N/A for Ptolemy on Solar Year
The model class claimed to implement all methods required by the Test class, but at least one was left unimplemented, so this model will be skipped.
Score: N/A for Copernicus on Solar Year
The model class claimed to implement all methods required by the Test class, but at least one was left unimplemented, so this model will be skipped.
Score: N/A for Kepler on Solar Year
The model class claimed to implement all methods required by the Test class, but at least one was left unimplemented, so this model will be skipped.
Score: N/A for Newton on Solar Year
The model class claimed to implement all methods required by the Test class, but at least one was left unimplemented, so this model will be skipped.
Score: N/A for Einstein on Solar Year
Score: Fail for Ptolemy on Ecc. Mars
Score: Fail for Ptolemy on Ecc. Saturn
Score: Fail for Copernicus on Ecc. Mars
Score: Fail for Copernicus on Ecc. Saturn
Score: Pass for Kepler on Ecc. Mars
Score: Pass for Kepler on Ecc. Saturn
Score: Pass for Newton on Ecc. Mars
Score: Pass for Newton on Ecc. Saturn
Score: Pass for Einstein on Ecc. Mars
Score: Pass for Einstein on Ecc. Saturn
Score: N/A for Ptolemy on Prlx. Promixa Centauri
Score: N/A for Ptolemy on Prlx. 61 Cygni
Score: Fail for Copernicus on Prlx. Promixa Centauri
Score: Fail for Copernicus on Prlx. 61 Cygni
Score: Fail for Kepler on Prlx. Promixa Centauri
Score: Fail for Kepler on Prlx. 61 Cygni
Score: Pass for Newton on Prlx. Promixa Centauri
Score: Pass for Newton on Prlx. 61 Cygni
Score: Pass for Einstein on Prlx. Promixa Centauri
Score: Pass for Einstein on Prlx. 61 Cygni
The model class claimed to implement all methods required by the Test class, but at least one was left unimplemented, so this model will be skipped.
Score: N/A for Ptolemy on Phln. Mercury
The model class claimed to implement all methods required by the Test class, but at least one was left unimplemented, so this model will be skipped.
Score: N/A for Copernicus on Phln. Mercury
The model class claimed to implement all methods required by the Test class, but at least one was left unimplemented, so this model will be skipped.
Score: N/A for Kepler on Phln. Mercury
The model class claimed to implement all methods required by the Test class, but at least one was left unimplemented, so this model will be skipped.
Score: N/A for Newton on Phln. Mercury
The model class claimed to implement all methods required by the Test class, but at least one was left unimplemented, so this model will be skipped.
Score: N/A for Einstein on Phln. Mercury