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