<jupyter_start><jupyter_text>![SciUnit Logo](https://raw.githubusercontent.com/scidash/assets/master/logos/SciUnit/sci-unit-tag.png) Chapter 5. The Real Example of Using SciUnit To Test Cosmology Models(or [back to Chapter 4](chapter4.ipynb))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.<jupyter_code>!pip install -q sciunit<jupyter_output><empty_output><jupyter_text>Like what we do before, let's create some capabilities, tests, models, and test suites.<jupyter_code>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<jupyter_output><empty_output><jupyter_text>Capabilities<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>Models<jupyter_code>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):
    """<NAME>, "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):
    """<NAME>, "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):
    """<NAME>, "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):
    """<NAME>, "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):
    """<NAME>, "The Foundation of the General Theory of Relativity"
    <NAME> 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)<jupyter_output><empty_output><jupyter_text>Tests<jupyter_code>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}<jupyter_output><empty_output><jupyter_text>Observations The data collected from observations and experiments.<jupyter_code># 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, 
            # <NAME> 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<jupyter_output><empty_output><jupyter_text>TestSuites Let's put the test instances into a poython list, and create the TestSuite instances with them.<jupyter_code>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]<jupyter_output><empty_output><jupyter_text>And then, we can let each suite instance judge each model in the model list.<jupyter_code>ptolemy = Ptolemy()
copernicus = Copernicus()
kepler = Kepler()
newton = Newton()
einstein = Einstein()
models = [ptolemy,copernicus,kepler,newton,einstein]
for suite in suites:
    suite.judge(models)<jupyter_output>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.
StarcoderdataJupyterStructuredCleanDedup