Minor photonics reformatting (#1002)

facebookresearch · Jan 15, 2021 · 294aed2 · 294aed2
1 parent d58b5c5
commit 294aed2
Show file tree

Hide file tree

Showing 3 changed files with 111 additions and 66 deletions.
diff --git a/nevergrad/functions/photonics/photonics.py b/nevergrad/functions/photonics/photonics.py
@@ -19,27 +19,28 @@
 from pathlib import Path
 import numpy as np
 from scipy.linalg import toeplitz
+
 # pylint: disable=blacklisted-name,too-many-locals,too-many-arguments
 
 
 def bragg(X: np.ndarray) -> float:
     """
-        Cost function for the Bragg mirror problem: maximizing the reflection
-        when the refractive index are given for all the layers.
-        Input: a vector whose components represent each the thickness of each
-        layer.
-        https://hal.archives-ouvertes.fr/hal-02613161
+    Cost function for the Bragg mirror problem: maximizing the reflection
+    when the refractive index are given for all the layers.
+    Input: a vector whose components represent each the thickness of each
+    layer.
+    https://hal.archives-ouvertes.fr/hal-02613161
     """
     lam = 600
     bar = int(np.size(X) / 2)
     n = np.concatenate(([1], np.sqrt(X[0:bar]), [1.7320508075688772]))
-    Type = np.arange(0, bar + 2)
+    type_ = np.arange(0, bar + 2)
     hauteur = np.concatenate(([0], X[bar : 2 * bar], [0]))
-    tmp = np.tan(2 * np.pi * n[Type] * hauteur / lam)
+    tmp = np.tan(2 * np.pi * n[type_] * hauteur / lam)
     # Specific to this substrate.
     Z = n[-1]
-    for k in range(np.size(Type) - 1, 0, -1):
-        Z = (Z - 1j * n[Type[k]] * tmp[k]) / (1 - 1j * tmp[k] * Z / n[Type[k]])
+    for k in range(np.size(type_) - 1, 0, -1):
+        Z = (Z - 1j * n[type_[k]] * tmp[k]) / (1 - 1j * tmp[k] * Z / n[type_[k]])
     # Specific to air.
     r = (1 - Z) / (1 + Z)
     c = np.real(1 - r * np.conj(r))
@@ -49,15 +50,15 @@ def bragg(X: np.ndarray) -> float:
 def chirped(X: np.ndarray) -> float:
     lam = np.linspace(500, 800, 50)
     n = np.array([1, 1.4142135623730951, 1.7320508075688772])
-    Type = np.concatenate(([0], np.tile([2, 1], int(np.size(X) / 2)), [2]))
+    type_ = np.concatenate(([0], np.tile([2, 1], int(np.size(X) / 2)), [2]))
     hauteur = np.concatenate(([0], X, [0]))
     r = np.zeros(np.size(lam)) + 0j
     for m in range(0, np.size(lam)):
         # Specific to this substrate.
-        tmp = np.tan(2 * np.pi * n[Type] * hauteur / lam[m])
+        tmp = np.tan(2 * np.pi * n[type_] * hauteur / lam[m])
         Z = 1.7320508075688772
-        for k in range(np.size(Type) - 1, 0, -1):
-            Z = (Z - 1j * n[Type[k]] * tmp[k]) / (1 - 1j * tmp[k] * Z / n[Type[k]])
+        for k in range(np.size(type_) - 1, 0, -1):
+            Z = (Z - 1j * n[type_[k]] * tmp[k]) / (1 - 1j * tmp[k] * Z / n[type_[k]])
         # Specific to air.
         r[m] = (1 - Z) / (1 + Z)
     # c=1-np.mean(abs(r)**2)
@@ -124,7 +125,9 @@ def marche(a: float, b: float, p: float, n: int, x: float) -> np.ndarray:
     return T  # type: ignore
 
 
-def creneau(k0: float, a0: float, pol: float, e1: float, e2: float, a: float, n: int, x0: float) -> tp.Tuple[np.ndarray, np.ndarray]:
+def creneau(
+    k0: float, a0: float, pol: float, e1: float, e2: float, a: float, n: int, x0: float
+) -> tp.Tuple[np.ndarray, np.ndarray]:
     nmod = int(n / 2)
     alpha = np.diag(a0 + 2 * np.pi * np.arange(-nmod, nmod + 1))
     if pol == 0:
@@ -152,9 +155,7 @@ def creneau(k0: float, a0: float, pol: float, e1: float, e2: float, a: float, n:
 
 def homogene(k0: float, a0: float, pol: float, epsilon: float, n: int) -> tp.Tuple[np.ndarray, np.ndarray]:
     nmod = int(n / 2)
-    valp = np.sqrt(
-        epsilon * k0 * k0 - (a0 + 2 * np.pi * np.arange(-nmod, nmod + 1)) ** 2 + 0j
-    )
+    valp = np.sqrt(epsilon * k0 * k0 - (a0 + 2 * np.pi * np.arange(-nmod, nmod + 1)) ** 2 + 0j)
     valp = valp * (1 - 2 * (valp < 0)) * (pol / epsilon + (1 - pol))
     P = np.block([[np.eye(n)], [np.diag(valp)]])
     return P, valp
@@ -163,11 +164,7 @@ def homogene(k0: float, a0: float, pol: float, epsilon: float, n: int) -> tp.Tup
 def interface(P: np.ndarray, Q: np.ndarray) -> np.ndarray:
     n = int(P.shape[1])
     S = np.matmul(
-        np.linalg.inv(
-            np.block(
-                [[P[0:n, 0:n], -Q[0:n, 0:n]], [P[n : 2 * n, 0:n], Q[n : 2 * n, 0:n]]]
-            )
-        ),
+        np.linalg.inv(np.block([[P[0:n, 0:n], -Q[0:n, 0:n]], [P[n : 2 * n, 0:n], Q[n : 2 * n, 0:n]]])),
         np.block([[-P[0:n, 0:n], Q[0:n, 0:n]], [P[n : 2 * n, 0:n], Q[n : 2 * n, 0:n]]]),
     )
     return S  # type: ignore
@@ -190,9 +187,7 @@ def morpho(X: np.ndarray) -> float:
     l = lam / d  # noqa
     k0 = 2 * np.pi / l
     P, V = homogene(k0, 0, pol, 1, n)
-    S = np.block(
-        [[np.zeros([n, n]), np.eye(n, dtype=np.complex)], [np.eye(n), np.zeros([n, n])]]
-    )
+    S = np.block([[np.zeros([n, n]), np.eye(n, dtype=np.complex)], [np.eye(n), np.zeros([n, n])]])
     for j in range(0, n_motifs):
         Pc, Vc = creneau(k0, 0, pol, e2, 1, a[j], n, x0[j])
         S = cascade(S, interface(P, Pc))
@@ -229,12 +224,15 @@ def morpho(X: np.ndarray) -> float:
     cost += bar / lams.size
     return cost
 
+
 i = complex(0, 1)
 
 
 def epscSi(lam: np.ndarray) -> np.ndarray:
     a = np.arange(250, 1500, 5)
-    e = np.load(Path(__file__).with_name("epsilon_epscSi.npy"))  # saved with np.save(filename, e) and dumped in this folder
+    e = np.load(
+        Path(__file__).with_name("epsilon_epscSi.npy")
+    )  # saved with np.save(filename, e) and dumped in this folder
     y = np.argmin(np.sign(lam - a))
     y = y - 1
     epsilon = (e[y + 1] - e[y]) / (a[y + 1] - a[y]) * (lam - a[y]) + e[y]
@@ -257,42 +255,51 @@ def cascade2(A: np.ndarray, B: np.ndarray) -> np.ndarray:
 
 
 def solar(lam: np.ndarray) -> np.ndarray:
+    # saved with np.save(filename, e) and dumped in this folder
     a = np.load(Path(__file__).with_name("wavelength_solar.npy"))
-    e = np.load(Path(__file__).with_name("epsilon_solar.npy"))  # saved with np.save(filename, e) and dumped in this folder
+    e = np.load(Path(__file__).with_name("epsilon_solar.npy"))
     jsc = np.interp(lam, a, e)
     return jsc  # type: ignore
 
 
-def absorption(lam, Epsilon, Mu, Type, hauteur, pol, theta):
-    if pol == 0:
-        f = Mu
-    else:
-        f = Epsilon
+def absorption(
+    lam: float,
+    epsilon: np.ndarray,
+    mu: np.ndarray,
+    type_: np.ndarray,
+    hauteur: np.ndarray,
+    pol: int,
+    theta: float,
+) -> np.ndarray:
+    f = mu if not pol else epsilon
     k0 = 2 * np.pi / lam
-    g = Type.size
-    alpha = np.sqrt(Epsilon[Type[0]] * Mu[Type[0]]) * k0 * np.sin(theta)
-    gamma = np.sqrt(Epsilon[Type] * Mu[Type] * k0**2 - np.ones(g) * alpha**2)
-    if np.real(Epsilon[Type[0]]) < 0 and np.real(Mu[Type[0]]) < 0:
+    g = type_.size
+    alpha = np.sqrt(epsilon[type_[0]] * mu[type_[0]]) * k0 * np.sin(theta)
+    gamma = np.sqrt(epsilon[type_] * mu[type_] * k0 ** 2 - np.ones(g) * alpha ** 2)
+    if np.real(epsilon[type_[0]]) < 0 and np.real(mu[type_[0]]) < 0:
         gamma[0] = -gamma[0]
     if g > 2:
-        gamma[1:g - 2] = gamma[1:g - 2] * (1 - 2 * (np.imag(gamma[1:g - 2]) < 0))
+        gamma[1 : g - 2] = gamma[1 : g - 2] * (1 - 2 * (np.imag(gamma[1 : g - 2]) < 0))
     if (
-        np.real(Epsilon[Type[g - 1]]) < 0
-        and np.real(Mu[Type[g - 1]]) < 0
-        and np.real(np.sqrt(Epsilon[Type[g - 1]] * Mu[Type[g - 1]] * k0**2 - alpha**2)) != 0
+        np.real(epsilon[type_[g - 1]]) < 0
+        and np.real(mu[type_[g - 1]]) < 0
+        and np.real(np.sqrt(epsilon[type_[g - 1]] * mu[type_[g - 1]] * k0 ** 2 - alpha ** 2)) != 0
     ):
-        gamma[g - 1] = -np.sqrt(Epsilon[Type[g - 1]] * Mu[Type[g - 1]] * k0**2 - alpha**2)
+        gamma[g - 1] = -np.sqrt(epsilon[type_[g - 1]] * mu[type_[g - 1]] * k0 ** 2 - alpha ** 2)
     else:
-        gamma[g - 1] = np.sqrt(Epsilon[Type[g - 1]] * Mu[Type[g - 1]] * k0**2 - alpha**2)
+        gamma[g - 1] = np.sqrt(epsilon[type_[g - 1]] * mu[type_[g - 1]] * k0 ** 2 - alpha ** 2)
     T = np.zeros(((2 * g, 2, 2)), dtype=complex)
     T[0] = [[0, 1], [1, 0]]
     for k2 in range(g - 1):
         t = np.exp(i * gamma[k2] * hauteur[k2])
         T[2 * k2 + 1] = [[0, t], [t, 0]]
-    # Interface scattering matrix
-        b1 = gamma[k2] / f[Type[k2]]
-        b2 = gamma[k2 + 1] / f[Type[k2 + 1]]
-        T[2 * k2 + 2] = [[(b1 - b2) / (b1 + b2), 2 * b2 / (b1 + b2)], [2 * b1 / (b1 + b2), (b2 - b1) / (b1 + b2)]]
+        # Interface scattering matrix
+        b1 = gamma[k2] / f[type_[k2]]
+        b2 = gamma[k2 + 1] / f[type_[k2 + 1]]
+        T[2 * k2 + 2] = [
+            [(b1 - b2) / (b1 + b2), 2 * b2 / (b1 + b2)],
+            [2 * b1 / (b1 + b2), (b2 - b1) / (b1 + b2)],
+        ]
         t = np.exp(i * gamma[g - 1] * hauteur[g - 1])
         T[2 * g - 1] = [[0, t], [t, 0]]
     H = np.zeros(((2 * g - 1, 2, 2)), dtype=complex)
@@ -318,20 +325,31 @@ def absorption(lam, Epsilon, Mu, Type, hauteur, pol, theta):
     poynting = np.zeros(2 * g, dtype=complex)
     if pol == 0:  # TE
         for j in range(2 * g):
-            poynting[j] = np.real((I[j][0, 0] + I[j][1, 0]) * np.conj((I[j][0, 0] - I[j][1, 0])
-                                                                      * gamma[w] / Mu[Type[w]]) * Mu[Type[0]] / (gamma[0]))
+            poynting[j] = np.real(
+                (I[j][0, 0] + I[j][1, 0])
+                * np.conj((I[j][0, 0] - I[j][1, 0]) * gamma[w] / mu[type_[w]])
+                * mu[type_[0]]
+                / (gamma[0])
+            )
             w = w + 1 - np.mod(j + 1, 2)
     else:  # TM
         for j in range(2 * g):
-            poynting[j] = np.real((I[j][0, 0] - I[j][1, 0]) * np.conj((I[j][0, 0] + I[j][1, 0])
-                                                                      * gamma[w] / Epsilon[Type[w]]) * Epsilon[Type[0]] / (gamma[0]))
+            poynting[j] = np.real(
+                (I[j][0, 0] - I[j][1, 0])
+                * np.conj((I[j][0, 0] + I[j][1, 0]) * gamma[w] / epsilon[type_[w]])
+                * epsilon[type_[0]]
+                / (gamma[0])
+            )
             w = w + 1 - np.mod(j + 1, 2)
     tmp = abs(-np.diff(poynting))
     absorb = tmp[np.arange(0, 2 * g, 2)]
-    return absorb
+    return absorb  # type: ignore
 
 
 def cf_photosic_reference(X: np.ndarray) -> float:
+    """vector X is only the thicknesses of each layers, because the materials (so the epislon)
+    are imposed by the function. This is similar in the chirped function.
+    """
     lam_min = 375
     lam_max = 750
     n_lam = 100
@@ -341,12 +359,12 @@ def cf_photosic_reference(X: np.ndarray) -> float:
     Ab = np.zeros(n_lam)
     for k in range(n_lam):
         lam = vlam[k]
-        Epsilon = np.array([1, 2, 3, epscSi(lam)], dtype=complex)
-        Mu = np.ones(Epsilon.size)
-        Type = np.append(0, np.append(np.tile(np.array([1, 2]), int(X.size / 2)), 3))
+        epsilon = np.array([1, 2, 3, epscSi(lam)], dtype=complex)
+        mu = np.ones(epsilon.size)
+        type_ = np.append(0, np.append(np.tile(np.array([1, 2]), int(X.size / 2)), 3))
         hauteur = np.append(0, np.append(X, 30000))
         pol = 0
-        absorb = absorption(lam, Epsilon, Mu, Type, hauteur, pol, theta)
+        absorb = absorption(lam, epsilon, mu, type_, hauteur, pol, theta)
         scc[k] = solar(lam)
         Ab[k] = absorb[len(absorb) - 1]
     max_scc = np.trapz(scc, vlam)
@@ -356,10 +374,17 @@ def cf_photosic_reference(X: np.ndarray) -> float:
     return cost  # type: ignore
 
 
-def cf_photosic_realist(eps_and_d: np.ndarray) -> float:
+def cf_photosic_realistic(eps_and_d: np.ndarray) -> float:
+    """eps_and_d is a vector composed in a first part with the epsilon values
+    (the material used in each one of the layers), and in a second part with the
+    thicknesses of each one of the layers, like in Bragg.
+    Any number of layers can work. Basically I used between 4 and 50 layers,
+    and the best results are generally obtained when the structure has between 10 and 20 layers.
+    The epsilon values are generally comprised between 1.00 and 9.00.
+    """
     dimension = int(eps_and_d.size / 2)
     eps = eps_and_d[0:dimension]
-    d = eps_and_d[dimension:dimension * 2]
+    d = eps_and_d[dimension : dimension * 2]
     epsd = np.array([eps, d])
     lam_min = 375
     lam_max = 750
@@ -372,12 +397,12 @@ def cf_photosic_realist(eps_and_d: np.ndarray) -> float:
     for k in range(n_lam):
         # absorb=absorption(epsd,theta*pi/180,lam[k])
         lam = vlam[k]
-        Epsilon = np.append(1, np.append(epsd[0], epscSi(lam)))
-        Mu = np.ones(Epsilon.size)
-        Type = np.arange(0, epsd[0].size + 2)
+        epsilon = np.append(1, np.append(epsd[0], epscSi(lam)))
+        mu = np.ones(epsilon.size)
+        type_ = np.arange(0, epsd[0].size + 2)
         hauteur = np.append(0, np.append(epsd[1], 30000))
         pol = 0
-        absorb = absorption(lam, Epsilon, Mu, Type, hauteur, pol, theta)
+        absorb = absorption(lam, epsilon, mu, type_, hauteur, pol, theta)
         scc[k] = solar(lam)
         Ab[k] = absorb[len(absorb) - 1]
     max_scc = np.trapz(scc, vlam)

diff --git a/nevergrad/functions/photonics/test_core.py b/nevergrad/functions/photonics/test_core.py
@@ -171,8 +171,29 @@ def test_photosic_reference() -> None:
 
 
 def test_photosic_realist() -> None:
-    eps_and_d = np.array([2.0000, 3.0000, 2.1076, 2.0000, 3.0000, 2.5783, 2.0000, 3.0000,
-                          2.0000, 3.0000, 90.0231, 78.9789, 72.8369, 99.9577, 82.7487,
-                          62.7583, 104.1682, 139.9002, 93.3356, 75.6039])
-    cf_test = photonics.cf_photosic_realist(eps_and_d)
+    eps_and_d = np.array(
+        [
+            2.0000,
+            3.0000,
+            2.1076,
+            2.0000,
+            3.0000,
+            2.5783,
+            2.0000,
+            3.0000,
+            2.0000,
+            3.0000,
+            90.0231,
+            78.9789,
+            72.8369,
+            99.9577,
+            82.7487,
+            62.7583,
+            104.1682,
+            139.9002,
+            93.3356,
+            75.6039,
+        ]
+    )
+    cf_test = photonics.cf_photosic_realistic(eps_and_d)
     np.testing.assert_almost_equal(cf_test, 0.08602574254532869)
diff --git a/pyproject.toml b/pyproject.toml
@@ -9,7 +9,6 @@ exclude = '''
     | dist
     | docs
     | nevergrad/functions/control
-    | nevergrad/functions/photonics
   )
 )
 '''