Features in potentials (#44)

* add `FeatureFunction` to `potentials.py`, an abstract base class for different types of feature functions such as steps or dips

* add `FeaturePotential` to `potentials.py`, an abstract base class for an inflationary potential with a feature that inherits from `InflationaryPotential` and from `FeatureFunction`

* add `GaussianDip` to `potentials.py`, implementing a Gaussian dip feature function

* add `TanhStep` to `potentials.py`, implementing a step feature function using a hyperbolic tangent

* add `StarobinskyGaussianDipPotential` and `StarobinskyTanhStepPotential` to `potentials.py`, implementing the respective features for a Starobinsky potential

* version bump to 2.9.0; introducing features in potentials

* change amplitude usage for `GaussianDip` such that positive input corresponds to a dip as the name suggests

* rename `GaussDip` to `GaussianDip`

* add test for feature potentials

README.rst

@@ -3,7 +3,7 @@ primpy: calculations for the primordial Universe
 ================================================
 :primpy: calculations for the primordial Universe
 :Author: Lukas Hergt
-:Version: 2.8.0
+:Version: 2.9.0
 :Homepage: https://github.com/lukashergt/primpy
 :Documentation: https://primpy.readthedocs.io
 

primpy/__version__.py

@@ -1,4 +1,4 @@
 #!/usr/bin/env python
 """:mod:`primpy.__version__`: version file for primpy."""
 
-__version__ = '2.8.0'
+__version__ = '2.9.0'

primpy/potentials.py

@@ -1109,3 +1109,155 @@ def sr_As2Lambda(cls, A_s, phi_star, N_star, **pot_kwargs):
 #     def inv_V(self, V):
 #         """`phi(V) = phi0 * (1 - sqrt(V) / Lambda**2)**(1/p)`."""
 #         return self.phi0 * (1 - np.sqrt(V) / self.Lambda**2)**(1/self.p)
+
+
+class FeatureFunction(ABC):
+    """Feature in the inflationary potential."""
+
+    @staticmethod
+    @abstractmethod
+    def F(x, x0, a, b):
+        """Feature function."""
+
+    @staticmethod
+    @abstractmethod
+    def dF(x, x0, a, b):
+        """Feature function derivative."""
+
+    @staticmethod
+    @abstractmethod
+    def d2F(x, x0, a, b):
+        """Feature function 2nd derivative."""
+
+    @staticmethod
+    @abstractmethod
+    def d3F(x, x0, a, b):
+        """Feature function 3rd derivative."""
+
+
+class GaussianDip(FeatureFunction):
+    """Gaussian: `F(x) = -a * exp(-(x-x0)**2 / (2*b**2))`."""
+
+    @staticmethod
+    def F(x, x0, a, b):
+        """`F(x) = -a * exp(-(x-x0)**2 / (2*b**2))`."""
+        return -a * np.exp(-(x - x0)**2 / (2 * b**2))
+
+    @staticmethod
+    def dF(x, x0, a, b):
+        """`F'(x) = a/b**2 * (x-x0) * exp(-(x-x0)**2 / (2*b**2))`."""
+        return a / b**2 * (x - x0) * np.exp(-(x - x0)**2 / (2 * b**2))
+
+    @staticmethod
+    def d2F(x, x0, a, b):
+        """`F''(x) = a/b**4 * (b**2 - (x-x0)**2) * exp(-(x-x0)**2 / (2*b**2))`."""
+        return a / b**4 * (b**2 - (x - x0)**2) * np.exp(-(x - x0)**2 / (2 * b**2))
+
+    @staticmethod
+    def d3F(x, x0, a, b):
+        """`F'''(x) = a/b**6 * (x-x0) * ((x-x0)**2 - 3*b**2) * exp(-(x-x0)**2 / (2*b**2))`."""
+        return a / b**6 * (x - x0) * ((x - x0)**2 - 3 * b**2) * np.exp(-(x - x0)**2 / (2 * b**2))
+
+
+class TanhStep(FeatureFunction):
+    """Tanh step function: `F(x) = a * tanh((x - x0) / b)`."""
+
+    @staticmethod
+    def F(x, x0, a, b):
+        """`F(x) = a * tanh((x-x0)/b)`."""
+        return a * np.tanh((x - x0) / b)
+
+    @staticmethod
+    def dF(x, x0, a, b):
+        """`F'(x) = a/b * (1 - tanh((x-x0)/b)**2)`."""
+        tanh = np.tanh((x - x0) / b)
+        return a / b * (1 - tanh**2)
+
+    @staticmethod
+    def d2F(x, x0, a, b):
+        """`F''(x) = -2*a/b**2 * tanh((x-x0)/b) * (1 - tanh((x-x0)/b)**2)`."""
+        tanh = np.tanh((x - x0) / b)
+        return -2 * a / b**2 * tanh * (1 - tanh**2)
+
+    @staticmethod
+    def d3F(x, x0, a, b):
+        """`F'''(x) = -2*a/b**3 * (1 - 4*tanh((x-x0)/b)**2 + 3*tanh((x-x0)/b)**4)`."""
+        tanh = np.tanh((x - x0) / b)
+        return -2 * a / b**3 * (1 - 4 * tanh**2 + 3 * tanh**4)
+
+
+class FeaturePotential(InflationaryPotential, FeatureFunction):
+    """Inflationary potential with a feature: `V(phi) = V0(phi) * (1+F(phi))`."""
+
+    def __init__(self, **pot_kwargs):
+        self.phi_feature = pot_kwargs.pop('phi_feature')  # position of feature
+        self.a_feature = pot_kwargs.pop('a_feature')      # e.g. height or amplitude of feature
+        self.b_feature = pot_kwargs.pop('b_feature')      # e.g. width or gradient of feature
+        super().__init__(**pot_kwargs)
+
+    def V(self, phi):
+        """Inflationary potential V0(phi) with a feature function F(phi).
+
+        `V(phi) = V0(phi) * (1 + F(phi))`
+        """
+        V0 = super().V(phi)
+        F = super().F(phi, self.phi_feature, self.a_feature, self.b_feature)
+        return V0 * (1 + F)
+
+    def dV(self, phi):
+        """First derivative of the inflationary potential with a feature.
+
+        `V'(phi) = V0'(phi) * (1 + F(phi)) + V0(phi) * F'(phi)`
+        """
+        V0 = super().V(phi)
+        dV0 = super().dV(phi)
+        F = super().F(phi, self.phi_feature, self.a_feature, self.b_feature)
+        dF = super().dF(phi, self.phi_feature, self.a_feature, self.b_feature)
+        return dV0 * (1 + F) + V0 * dF
+
+    def d2V(self, phi):
+        """Second derivative of the inflationary potential with a feature.
+
+        `V''(phi) = V0''(phi) * (1 + F(phi)) + 2*V0'(phi)*F'(phi) + V0(phi)*F''(phi)`
+        """
+        V0 = super().V(phi)
+        dV0 = super().dV(phi)
+        d2V0 = super().d2V(phi)
+        F = super().F(phi, self.phi_feature, self.a_feature, self.b_feature)
+        dF = super().dF(phi, self.phi_feature, self.a_feature, self.b_feature)
+        d2F = super().d2F(phi, self.phi_feature, self.a_feature, self.b_feature)
+        return d2V0 * (1 + F) + 2 * dV0 * dF + V0 * d2F
+
+    def d3V(self, phi):
+        """Third derivative of the inflationary potential with a feature.
+
+        `V'''(phi) = V0'''(phi) * (1 + F(phi))
+                     + 3 * V0''(phi) * F'(phi)
+                     + 3 * V0'(phi) * F''(phi)
+                     + V0(phi) * F'''(phi)`
+        """
+        V0 = super().V(phi)
+        dV0 = super().dV(phi)
+        d2V0 = super().d2V(phi)
+        d3V0 = super().d3V(phi)
+        F = super().F(phi, self.phi_feature, self.a_feature, self.b_feature)
+        dF = super().dF(phi, self.phi_feature, self.a_feature, self.b_feature)
+        d2F = super().d2F(phi, self.phi_feature, self.a_feature, self.b_feature)
+        d3F = super().d3F(phi, self.phi_feature, self.a_feature, self.b_feature)
+        return d3V0 * (1 + F) + 3 * d2V0 * dF + 3 * dV0 * d2F + V0 * d3F
+
+
+class StarobinskyGaussianDipPotential(FeaturePotential, StarobinskyPotential, GaussianDip):
+    """Starobinsky potential with a Gaussian dip."""
+
+    tag = 'sgd'
+    name = 'StarobinskyGaussianDipPotential'
+    tex = r'Starobinsky with a Gaussian dip'
+
+
+class StarobinskyTanhStepPotential(FeaturePotential, StarobinskyPotential, TanhStep):
+    """Starobinsky potential with a hyperbolic tangent step."""
+
+    tag = 'sts'
+    name = 'StarobinskyTanhStepPotential'
+    tex = r'Starobinsky with a tanh step'

tests/test_potentials.py

@@ -1,6 +1,7 @@
 #!/usr/bin/env python
 """Tests for `primpy.potential` module."""
 import pytest
+from pytest import approx
 import numpy as np
 from numpy.testing import assert_array_equal, assert_allclose
 import primpy.potentials as pp
@@ -51,11 +52,11 @@ def test_inflationary_potentials(Pot, pot_kwargs, Lambda, phi):
 def test_quadratic_inflation_V(Lambda, phi):
     """Tests for `QuadraticPotential`."""
     pot1 = pp.QuadraticPotential(Lambda=Lambda)
-    assert pot1.V(phi=phi) == pytest.approx(Lambda**4 * phi**2)
-    assert pot1.dV(phi=phi) == pytest.approx(2 * Lambda**4 * phi)
-    assert pot1.d2V(phi=phi) == pytest.approx(2 * Lambda**4)
-    assert pot1.d3V(phi=phi) == pytest.approx(0)
-    assert pot1.inv_V(V=Lambda**4) == pytest.approx(np.sqrt(1))
+    assert pot1.V(phi=phi) == approx(Lambda**4 * phi**2)
+    assert pot1.dV(phi=phi) == approx(2 * Lambda**4 * phi)
+    assert pot1.d2V(phi=phi) == approx(2 * Lambda**4)
+    assert pot1.d3V(phi=phi) == approx(0)
+    assert pot1.inv_V(V=Lambda**4) == approx(np.sqrt(1))
     with pytest.raises(Exception):
         pp.QuadraticPotential(mass=Lambda**2)
 
@@ -128,12 +129,12 @@ def test_monomial_slow_roll(p, N_star):
     n_s = Pot.sr_Nstar2ns(N_star=N_star, p=p)
     assert 0.8 < n_s < 1
     assert n_s == 1 - p / (2 * N_star) - 1 / N_star
-    assert Pot.sr_ns2Nstar(n_s=n_s, p=p) == pytest.approx(N_star)
+    assert Pot.sr_ns2Nstar(n_s=n_s, p=p) == approx(N_star)
 
     r = Pot.sr_Nstar2r(N_star=N_star, p=p)
     assert 1e-2 < Pot.sr_Nstar2r(N_star=N_star, p=p) < 1
     assert r == 16 * p / (4 * N_star + p)
-    assert Pot.sr_r2Nstar(r=r, p=p) == pytest.approx(N_star)
+    assert Pot.sr_r2Nstar(r=r, p=p) == approx(N_star)
 
 
 @pytest.mark.parametrize('Pot, p', [(pp.LinearPotential, 1),
@@ -145,29 +146,29 @@ def test_specific_monomial_slow_roll(Pot, p, N_star):
     n_s = Pot.sr_Nstar2ns(N_star=N_star)
     assert 0.8 < n_s < 1
     assert n_s == 1 - p / (2 * N_star) - 1 / N_star
-    assert Pot.sr_ns2Nstar(n_s=n_s) == pytest.approx(N_star)
+    assert Pot.sr_ns2Nstar(n_s=n_s) == approx(N_star)
 
     r = Pot.sr_Nstar2r(N_star=N_star)
     assert 1e-2 < Pot.sr_Nstar2r(N_star=N_star) < 1
     assert r == 16 * p / (4 * N_star + p)
-    assert Pot.sr_r2Nstar(r=r) == pytest.approx(N_star)
+    assert Pot.sr_r2Nstar(r=r) == approx(N_star)
 
 
 @pytest.mark.parametrize('N_star', [20, 60, 90])
 def test_starobinsky_slow_roll(N_star):
     Pot = pp.StarobinskyPotential
 
     n_s = Pot.sr_Nstar2ns(N_star=N_star)
-    approx = 1 - 2 / N_star + (-3 + np.sqrt(3)) / N_star**2 + (-3 + 3*np.sqrt(3)) / N_star**3
+    aprx = 1 - 2 / N_star + (-3 + np.sqrt(3)) / N_star**2 + (-3 + 3*np.sqrt(3)) / N_star**3
     assert 0.8 < n_s < 1
-    assert n_s == pytest.approx(approx, rel=1e-3)
-    assert Pot.sr_ns2Nstar(n_s=n_s) == pytest.approx(N_star)
+    assert n_s == approx(aprx, rel=1e-3)
+    assert Pot.sr_ns2Nstar(n_s=n_s) == approx(N_star)
 
     r = Pot.sr_Nstar2r(N_star=N_star)
-    approx = 12 / N_star**2 - 12 * np.sqrt(3) / N_star**3 + 27 / N_star**4
+    aprx = 12 / N_star**2 - 12 * np.sqrt(3) / N_star**3 + 27 / N_star**4
     assert 1e-3 < r < 1
-    assert r == pytest.approx(approx, rel=1e-3)
-    assert Pot.sr_r2Nstar(r=r) == pytest.approx(N_star)
+    assert r == approx(aprx, rel=1e-3)
+    assert Pot.sr_r2Nstar(r=r) == approx(N_star)
 
 
 @pytest.mark.parametrize('pot_kwargs', [dict(phi0=10), dict(phi0=100), dict(phi0=1000)])
@@ -176,8 +177,56 @@ def test_natural_slow_roll(pot_kwargs, N_star):
     Pot = pp.NaturalPotential
     n_s = Pot.sr_Nstar2ns(N_star=N_star, **pot_kwargs)
     assert 0.8 < n_s < 1
-    assert Pot.sr_ns2Nstar(n_s=n_s, **pot_kwargs) == pytest.approx(N_star)
+    assert Pot.sr_ns2Nstar(n_s=n_s, **pot_kwargs) == approx(N_star)
 
     r = Pot.sr_Nstar2r(N_star=N_star, **pot_kwargs)
     assert 1e-4 < r < 1
-    assert Pot.sr_r2Nstar(r=r, **pot_kwargs) == pytest.approx(N_star)
+    assert Pot.sr_r2Nstar(r=r, **pot_kwargs) == approx(N_star)
+
+
+@pytest.mark.parametrize('a_feature', [0.1, 0.01])
+@pytest.mark.parametrize('b_feature', [0.5, 0.05])
+def test_feature_potential(a_feature, b_feature):
+    eps = np.finfo(float).eps
+    Lambda = 1e-3
+    phi_feature = 5
+    phi = np.linspace(0, 10, 200 + 1)[1:]
+    phi_smaller = phi[phi < phi_feature]
+    phi_greater = phi[phi > phi_feature]
+    phi_outside = phi[np.argwhere((phi < phi_feature - 3 * b_feature) |
+                                  (phi > phi_feature + 4 * b_feature)).ravel()]
+    pot = pp.StarobinskyPotential(Lambda=Lambda)
+
+    # Gaussian dip
+    fpot = pp.StarobinskyGaussianDipPotential(Lambda=Lambda, phi_feature=phi_feature,
+                                              a_feature=a_feature, b_feature=b_feature)
+    max_dV = np.max(np.abs(fpot.dV(phi=phi) / pot.dV(phi=phi) - 1))
+    max_d2V = np.max(np.abs(fpot.d2V(phi=phi) / pot.d2V(phi=phi) - 1))
+    max_d3V = np.max(np.abs(fpot.d3V(phi=phi) / pot.d3V(phi=phi) - 1))
+    assert np.all(fpot.V(phi=phi) >= 0)
+    assert np.all(fpot.V(phi=phi) <= pot.V(phi=phi))
+    assert_allclose(fpot.V(phi=phi) / pot.V(phi=phi), 1, rtol=a_feature + eps)
+    assert_allclose(fpot.dV(phi_outside) / pot.dV(phi_outside), 1, rtol=max_dV/10)
+    assert_allclose(fpot.d2V(phi_outside) / pot.d2V(phi_outside), 1, rtol=max_d2V/10)
+    assert_allclose(fpot.d3V(phi_outside) / pot.d3V(phi_outside), 1, rtol=max_d3V/10)
+    assert np.mean(fpot.dV(phi=phi_smaller) / pot.dV(phi=phi_smaller) - 1) < 0
+    assert np.mean(fpot.dV(phi=phi_greater) / pot.dV(phi=phi_greater) - 1) > 0
+    assert np.mean(fpot.d3V(phi=phi_smaller) / pot.d3V(phi=phi_smaller) - 1) > 0
+    assert np.mean(fpot.d3V(phi=phi_greater) / pot.d3V(phi=phi_greater) - 1) < 0
+
+    # Tanh step
+    fpot = pp.StarobinskyTanhStepPotential(Lambda=Lambda, phi_feature=phi_feature,
+                                           a_feature=a_feature, b_feature=b_feature)
+    max_dV = np.max(np.abs(fpot.dV(phi=phi) / pot.dV(phi=phi) - 1))
+    max_d2V = np.max(np.abs(fpot.d2V(phi=phi) / pot.d2V(phi=phi) - 1))
+    max_d3V = np.max(np.abs(fpot.d3V(phi=phi) / pot.d3V(phi=phi) - 1))
+    assert np.all(fpot.V(phi=phi) >= 0)
+    assert np.all(fpot.V(phi=phi_smaller) <= pot.V(phi=phi_smaller))
+    assert np.all(fpot.V(phi=phi_greater) >= pot.V(phi=phi_greater))
+    assert_allclose(fpot.V(phi=phi) / pot.V(phi=phi), 1, rtol=a_feature + eps)
+    assert_allclose(fpot.dV(phi_outside) / pot.dV(phi_outside), 1, rtol=max_dV/10)
+    assert_allclose(fpot.d2V(phi_outside) / pot.d2V(phi_outside), 1, rtol=max_d2V/10)
+    assert_allclose(fpot.d3V(phi_outside) / pot.d3V(phi_outside), 1, rtol=max_d3V/10)
+    assert np.mean(fpot.dV(phi=phi) / pot.dV(phi=phi) - 1) > 0
+    assert np.mean(fpot.d2V(phi=phi_smaller) / pot.d2V(phi=phi_smaller) - 1) < 0
+    assert np.mean(fpot.d2V(phi=phi_greater) / pot.d2V(phi=phi_greater) - 1) > 0